(git:e5b1968)
Loading...
Searching...
No Matches
auto_basis.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Automatic generation of auxiliary basis sets of different kind
10!> \author JGH
11!>
12!> <b>Modification history:</b>
13!> - 11.2017 creation [JGH]
14! **************************************************************************************************
15MODULE auto_basis
16 USE aux_basis_set, ONLY: create_aux_basis
17 USE basis_set_types, ONLY: get_gto_basis_set,&
18 gto_basis_set_type,&
19 sort_gto_basis_set
20 USE bibliography, ONLY: stoychev2016,&
21 cite_reference
22 USE kinds, ONLY: default_string_length,&
23 dp
24 USE mathconstants, ONLY: dfac,&
25 fac,&
26 gamma1,&
27 pi,&
28 rootpi
29 USE orbital_pointers, ONLY: init_orbital_pointers
30 USE periodic_table, ONLY: get_ptable_info
31 USE qs_kind_types, ONLY: get_qs_kind,&
32 qs_kind_type
33#include "./base/base_uses.f90"
34
35 IMPLICIT NONE
36
37 PRIVATE
38
39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'auto_basis'
40
41 PUBLIC :: create_ri_aux_basis_set, create_lri_aux_basis_set, &
42 create_oce_basis
43
44CONTAINS
45
46! **************************************************************************************************
47!> \brief Create a RI_AUX basis set using some heuristics
48!> \param ri_aux_basis_set ...
49!> \param qs_kind ...
50!> \param basis_cntrl ...
51!> \param basis_type ...
52!> \param basis_sort ...
53!> \date 01.11.2017
54!> \author JGH
55! **************************************************************************************************
56 SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_type, basis_sort)
57 TYPE(gto_basis_set_type), POINTER :: ri_aux_basis_set
58 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
59 INTEGER, INTENT(IN) :: basis_cntrl
60 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
61 INTEGER, INTENT(IN), OPTIONAL :: basis_sort
62
63 CHARACTER(LEN=2) :: element_symbol
64 CHARACTER(LEN=default_string_length) :: bsname
65 INTEGER :: i, j, jj, l, laux, linc, lmax, lval, lx, &
66 nsets, nx, z
67 INTEGER, DIMENSION(0:18) :: nval
68 INTEGER, DIMENSION(0:9, 1:20) :: nl
69 INTEGER, DIMENSION(1:3) :: ls1, ls2, npgf
70 INTEGER, DIMENSION(:), POINTER :: econf
71 REAL(kind=dp) :: xv, zval
72 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
73 REAL(kind=dp), DIMENSION(0:18) :: bv, bval, fv, peff, pend, pmax, pmin
74 REAL(kind=dp), DIMENSION(0:9) :: zeff, zmax, zmin
75 REAL(kind=dp), DIMENSION(3) :: amax, amin, bmin
76 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
77
78 !
79 CALL cite_reference(stoychev2016)
80 !
81 bv(0:18) = (/1.8_dp, 2.0_dp, 2.2_dp, 2.2_dp, 2.3_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, &
82 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp/)
83 fv(0:18) = (/20.0_dp, 4.0_dp, 4.0_dp, 3.5_dp, 2.5_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, &
84 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp/)
85 !
86 cpassert(.NOT. ASSOCIATED(ri_aux_basis_set))
87 NULLIFY (orb_basis_set)
88 IF (.NOT. PRESENT(basis_type)) THEN
89 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
90 ELSE
91 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type=basis_type)
92 END IF
93 IF (ASSOCIATED(orb_basis_set)) THEN
94 ! BASIS_SET ORB NONE associates the pointer orb_basis_set, but does not contain
95 ! any actual basis functions. Therefore, we catch it here to avoid spurious autogenerated
96 ! RI_AUX basis sets.
97 IF (sum(orb_basis_set%nsgf_set) == 0) THEN
98 CALL cp_abort(__location__, &
99 "Cannot autocreate RI_AUX basis set for at least one of the given "// &
100 "primary basis sets due to missing exponents. If you have invoked BASIS_SET NONE, "// &
101 "you should state BASIS_SET RI_AUX NONE explicitly in the input.")
102 END IF
103 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
104 !Note: RI basis coud require lmax up to 2*orb_lmax. This ensures that all orbital pointers
105 ! are properly initialized before building the basis
106 CALL init_orbital_pointers(2*lmax)
107 CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
108 CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
109 CALL get_ptable_info(element_symbol, ielement=z)
110 lval = 0
111 DO l = 0, maxval(ubound(econf))
112 IF (econf(l) > 0) lval = l
113 END DO
114 IF (sum(econf) /= nint(zval)) THEN
115 cpwarn("Valence charge and electron configuration not consistent")
116 END IF
117 pend = 0.0_dp
118 linc = 1
119 IF (z > 18) linc = 2
120 SELECT CASE (basis_cntrl)
121 CASE (0)
122 laux = max(2*lval, lmax + linc)
123 CASE (1)
124 laux = max(2*lval, lmax + linc)
125 CASE (2)
126 laux = max(2*lval, lmax + linc + 1)
127 CASE (3)
128 laux = max(2*lmax, lmax + linc + 2)
129 CASE DEFAULT
130 cpabort("Invalid value of control variable")
131 END SELECT
132 !
133 DO l = 2*lmax + 1, laux
134 xv = peff(2*lmax)
135 pmin(l) = xv
136 pmax(l) = xv
137 peff(l) = xv
138 pend(l) = xv
139 END DO
140 !
141 DO l = 0, laux
142 IF (l <= 2*lval) THEN
143 pend(l) = min(fv(l)*peff(l), pmax(l))
144 bval(l) = 1.8_dp
145 ELSE
146 pend(l) = peff(l)
147 bval(l) = bv(l)
148 END IF
149 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
150 nval(l) = max(ceiling(xv), 0)
151 END DO
152 ! first set include valence only
153 nsets = 1
154 ls1(1) = 0
155 ls2(1) = lval
156 DO l = lval + 1, laux
157 IF (nval(l) < nval(lval) - 1) EXIT
158 ls2(1) = l
159 END DO
160 ! second set up to 2*lval
161 IF (laux > ls2(1)) THEN
162 IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN
163 nsets = 2
164 ls1(2) = ls2(1) + 1
165 ls2(2) = laux
166 ELSE
167 nsets = 2
168 ls1(2) = ls2(1) + 1
169 ls2(2) = min(2*lval, laux)
170 lx = ls2(2)
171 DO l = lx + 1, laux
172 IF (nval(l) < nval(lx) - 1) EXIT
173 ls2(2) = l
174 END DO
175 IF (laux > ls2(2)) THEN
176 nsets = 3
177 ls1(3) = ls2(2) + 1
178 ls2(3) = laux
179 END IF
180 END IF
181 END IF
182 !
183 amax = 0.0
184 amin = huge(0.0_dp)
185 bmin = huge(0.0_dp)
186 DO i = 1, nsets
187 DO j = ls1(i), ls2(i)
188 amax(i) = max(amax(i), pend(j))
189 amin(i) = min(amin(i), pmin(j))
190 bmin(i) = min(bmin(i), bval(j))
191 END DO
192 xv = log(amax(i)/amin(i))/log(bmin(i)) + 1.e-10_dp
193 npgf(i) = max(ceiling(xv), 0)
194 END DO
195 nx = maxval(npgf(1:nsets))
196 ALLOCATE (zet(nx, nsets))
197 zet = 0.0_dp
198 nl = 0
199 DO i = 1, nsets
200 DO j = 1, npgf(i)
201 jj = npgf(i) - j + 1
202 zet(jj, i) = amin(i)*bmin(i)**(j - 1)
203 END DO
204 DO l = ls1(i), ls2(i)
205 nl(l, i) = nval(l)
206 END DO
207 END DO
208 bsname = trim(element_symbol)//"-RI-AUX-"//trim(orb_basis_set%name)
209 !
210 CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
211
212 DEALLOCATE (zet)
213
214 IF (PRESENT(basis_sort)) THEN
215 CALL sort_gto_basis_set(ri_aux_basis_set, basis_sort)
216 END IF
217
218 END IF
219
220 END SUBROUTINE create_ri_aux_basis_set
221! **************************************************************************************************
222!> \brief Create a LRI_AUX basis set using some heuristics
223!> \param lri_aux_basis_set ...
224!> \param qs_kind ...
225!> \param basis_cntrl ...
226!> \param exact_1c_terms ...
227!> \param tda_kernel ...
228!> \date 01.11.2017
229!> \author JGH
230! **************************************************************************************************
231 SUBROUTINE create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, &
232 exact_1c_terms, tda_kernel)
233 TYPE(gto_basis_set_type), POINTER :: lri_aux_basis_set
234 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
235 INTEGER, INTENT(IN) :: basis_cntrl
236 LOGICAL, INTENT(IN), OPTIONAL :: exact_1c_terms, tda_kernel
237
238 CHARACTER(LEN=2) :: element_symbol
239 CHARACTER(LEN=default_string_length) :: bsname
240 INTEGER :: i, j, l, laux, linc, lm, lmax, lval, n1, &
241 n2, nsets, z
242 INTEGER, DIMENSION(0:18) :: nval
243 INTEGER, DIMENSION(0:9, 1:50) :: nl
244 INTEGER, DIMENSION(1:50) :: ls1, ls2, npgf
245 INTEGER, DIMENSION(:), POINTER :: econf
246 LOGICAL :: e1terms, kernel_basis
247 REAL(kind=dp) :: xv, zval
248 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
249 REAL(kind=dp), DIMENSION(0:18) :: bval, peff, pend, pmax, pmin
250 REAL(kind=dp), DIMENSION(0:9) :: zeff, zmax, zmin
251 REAL(kind=dp), DIMENSION(4) :: bv, bx
252 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
253
254 !
255 IF (PRESENT(exact_1c_terms)) THEN
256 e1terms = exact_1c_terms
257 ELSE
258 e1terms = .false.
259 END IF
260 IF (PRESENT(tda_kernel)) THEN
261 kernel_basis = tda_kernel
262 ELSE
263 kernel_basis = .false.
264 END IF
265 IF (kernel_basis .AND. e1terms) THEN
266 CALL cp_warn(__location__, "LRI Kernel basis generation will ignore exact 1C term option.")
267 END IF
268 !
269 cpassert(.NOT. ASSOCIATED(lri_aux_basis_set))
270 NULLIFY (orb_basis_set)
271 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
272 IF (ASSOCIATED(orb_basis_set)) THEN
273 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
274 CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
275 CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
276 CALL get_ptable_info(element_symbol, ielement=z)
277 lval = 0
278 DO l = 0, maxval(ubound(econf))
279 IF (econf(l) > 0) lval = l
280 END DO
281 IF (sum(econf) /= nint(zval)) THEN
282 cpwarn("Valence charge and electron configuration not consistent")
283 END IF
284 !
285 linc = 1
286 IF (z > 18) linc = 2
287 pend = 0.0_dp
288 IF (kernel_basis) THEN
289 bv(1:4) = (/3.20_dp, 2.80_dp, 2.40_dp, 2.00_dp/)
290 bx(1:4) = (/4.00_dp, 3.50_dp, 3.00_dp, 2.50_dp/)
291 !
292 SELECT CASE (basis_cntrl)
293 CASE (0)
294 laux = lval + 1
295 CASE (1)
296 laux = max(lval + 1, lmax)
297 CASE (2)
298 laux = max(lval + 2, lmax + 1)
299 CASE (3)
300 laux = max(lval + 3, lmax + 2)
301 laux = min(laux, 2 + linc)
302 CASE DEFAULT
303 cpabort("Invalid value of control variable")
304 END SELECT
305 ELSE
306 bv(1:4) = (/2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp/)
307 bx(1:4) = (/2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp/)
308 !
309 SELECT CASE (basis_cntrl)
310 CASE (0)
311 laux = max(2*lval, lmax + linc)
312 laux = min(laux, 2 + linc)
313 CASE (1)
314 laux = max(2*lval, lmax + linc)
315 laux = min(laux, 3 + linc)
316 CASE (2)
317 laux = max(2*lval, lmax + linc + 1)
318 laux = min(laux, 4 + linc)
319 CASE (3)
320 laux = max(2*lval, lmax + linc + 1)
321 laux = min(laux, 4 + linc)
322 CASE DEFAULT
323 cpabort("Invalid value of control variable")
324 END SELECT
325 END IF
326 !
327 DO l = 2*lmax + 1, laux
328 pmin(l) = pmin(2*lmax)
329 pmax(l) = pmax(2*lmax)
330 peff(l) = peff(2*lmax)
331 END DO
332 !
333 nval = 0
334 IF (exact_1c_terms) THEN
335 DO l = 0, laux
336 IF (l <= lval + 1) THEN
337 pend(l) = zmax(l) + 1.0_dp
338 bval(l) = bv(basis_cntrl + 1)
339 ELSE
340 pend(l) = 2.0_dp*peff(l)
341 bval(l) = bx(basis_cntrl + 1)
342 END IF
343 pmin(l) = zmin(l)
344 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
345 nval(l) = max(ceiling(xv), 0)
346 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
347 END DO
348 ELSE
349 DO l = 0, laux
350 IF (l <= lval + 1) THEN
351 pend(l) = pmax(l)
352 bval(l) = bv(basis_cntrl + 1)
353 pmin(l) = zmin(l)
354 ELSE
355 pend(l) = 4.0_dp*peff(l)
356 bval(l) = bx(basis_cntrl + 1)
357 END IF
358 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
359 nval(l) = max(ceiling(xv), 0)
360 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
361 END DO
362 END IF
363 !
364 lm = min(2*lval, 3)
365 n1 = maxval(nval(0:lm))
366 IF (laux < lm + 1) THEN
367 n2 = 0
368 ELSE
369 n2 = maxval(nval(lm + 1:laux))
370 END IF
371 !
372 nsets = n1 + n2
373 ALLOCATE (zet(1, nsets))
374 zet = 0.0_dp
375 nl = 0
376 j = maxval(maxloc(nval(0:lm)))
377 DO i = 1, n1
378 ls1(i) = 0
379 ls2(i) = lm
380 npgf(i) = 1
381 zet(1, i) = pmin(j)*bval(j)**(i - 1)
382 DO l = 0, lm
383 nl(l, i) = 1
384 END DO
385 END DO
386 j = lm + 1
387 DO i = n1 + 1, nsets
388 ls1(i) = lm + 1
389 ls2(i) = laux
390 npgf(i) = 1
391 zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1)
392 DO l = lm + 1, laux
393 nl(l, i) = 1
394 END DO
395 END DO
396 !
397 bsname = trim(element_symbol)//"-LRI-AUX-"//trim(orb_basis_set%name)
398 !
399 CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
400 !
401 DEALLOCATE (zet)
402 END IF
403
404 END SUBROUTINE create_lri_aux_basis_set
405
406! **************************************************************************************************
407!> \brief ...
408!> \param oce_basis ...
409!> \param orb_basis ...
410!> \param lmax_oce ...
411!> \param nbas_oce ...
412! **************************************************************************************************
413 SUBROUTINE create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
414 TYPE(gto_basis_set_type), POINTER :: oce_basis, orb_basis
415 INTEGER, INTENT(IN) :: lmax_oce, nbas_oce
416
417 CHARACTER(LEN=default_string_length) :: bsname
418 INTEGER :: i, l, lmax, lx, nset, nx
419 INTEGER, ALLOCATABLE, DIMENSION(:) :: lmin, lset, npgf
420 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl
421 INTEGER, DIMENSION(:), POINTER :: npgf_orb
422 REAL(kind=dp) :: cval, x, z0, z1
423 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
424 REAL(kind=dp), DIMENSION(0:9) :: zeff, zmax, zmin
425
426 CALL get_basis_keyfigures(orb_basis, lmax, zmin, zmax, zeff)
427 IF (nbas_oce < 1) THEN
428 CALL get_gto_basis_set(gto_basis_set=orb_basis, nset=nset, npgf=npgf_orb)
429 nx = sum(npgf_orb(1:nset))
430 ELSE
431 nx = 0
432 END IF
433 nset = max(nbas_oce, nx)
434 lx = max(lmax_oce, lmax)
435 !
436 bsname = "OCE-"//trim(orb_basis%name)
437 ALLOCATE (lmin(nset), lset(nset), nl(0:9, nset), npgf(nset), zet(1, nset))
438 lmin = 0
439 lset = 0
440 nl = 1
441 npgf = 1
442 zet = 0.0_dp
443 !
444 z0 = minval(zmin(0:lmax))
445 z1 = maxval(zmax(0:lmax))
446 x = 1.0_dp/real(nset - 1, kind=dp)
447 cval = (z1/z0)**x
448 zet(1, nset) = z0
449 DO i = nset - 1, 1, -1
450 zet(1, i) = zet(1, i + 1)*cval
451 END DO
452 DO i = 1, nset
453 x = zet(1, i)
454 DO l = 1, lmax
455 z1 = 1.05_dp*zmax(l)
456 IF (x < z1) lset(i) = l
457 END DO
458 IF (lset(i) == lmax) lset(i) = lx
459 END DO
460 !
461 CALL create_aux_basis(oce_basis, bsname, nset, lmin, lset, nl, npgf, zet)
462 !
463 DEALLOCATE (lmin, lset, nl, npgf, zet)
464
465 END SUBROUTINE create_oce_basis
466! **************************************************************************************************
467!> \brief ...
468!> \param basis_set ...
469!> \param lmax ...
470!> \param zmin ...
471!> \param zmax ...
472!> \param zeff ...
473! **************************************************************************************************
474 SUBROUTINE get_basis_keyfigures(basis_set, lmax, zmin, zmax, zeff)
475 TYPE(gto_basis_set_type), POINTER :: basis_set
476 INTEGER, INTENT(OUT) :: lmax
477 REAL(kind=dp), DIMENSION(0:9), INTENT(OUT) :: zmin, zmax, zeff
478
479 INTEGER :: i, ipgf, iset, ishell, j, l, nset
480 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
481 INTEGER, DIMENSION(:, :), POINTER :: lshell
482 REAL(kind=dp) :: aeff, gcca, gccb, kval, rexp, rint, rno, &
483 zeta
484 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
485 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
486
487 CALL get_gto_basis_set(gto_basis_set=basis_set, &
488 nset=nset, &
489 nshell=nshell, &
490 npgf=npgf, &
491 l=lshell, &
492 lmax=lm, &
493 zet=zet, &
494 gcc=gcc)
495
496 lmax = maxval(lm)
497 cpassert(lmax <= 9)
498
499 zmax = 0.0_dp
500 zmin = huge(0.0_dp)
501 zeff = 0.0_dp
502
503 DO iset = 1, nset
504 ! zmin zmax
505 DO ipgf = 1, npgf(iset)
506 DO ishell = 1, nshell(iset)
507 l = lshell(ishell, iset)
508 zeta = zet(ipgf, iset)
509 zmax(l) = max(zmax(l), zeta)
510 zmin(l) = min(zmin(l), zeta)
511 END DO
512 END DO
513 ! zeff
514 DO ishell = 1, nshell(iset)
515 l = lshell(ishell, iset)
516 kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2)
517 rexp = 0.0_dp
518 rno = 0.0_dp
519 DO i = 1, npgf(iset)
520 gcca = gcc(i, ishell, iset)
521 DO j = 1, npgf(iset)
522 zeta = zet(i, iset) + zet(j, iset)
523 gccb = gcc(j, ishell, iset)
524 rint = 0.5_dp*fac(l + 1)/zeta**(l + 2)
525 rexp = rexp + gcca*gccb*rint
526 rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp)
527 rno = rno + gcca*gccb*rint
528 END DO
529 END DO
530 rexp = rexp/rno
531 aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2)
532 zeff(l) = max(zeff(l), aeff)
533 END DO
534 END DO
535
536 END SUBROUTINE get_basis_keyfigures
537
538! **************************************************************************************************
539!> \brief ...
540!> \param lmax ...
541!> \param zmin ...
542!> \param zmax ...
543!> \param zeff ...
544!> \param pmin ...
545!> \param pmax ...
546!> \param peff ...
547! **************************************************************************************************
548 SUBROUTINE get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
549 INTEGER, INTENT(IN) :: lmax
550 REAL(kind=dp), DIMENSION(0:9), INTENT(IN) :: zmin, zmax, zeff
551 REAL(kind=dp), DIMENSION(0:18), INTENT(OUT) :: pmin, pmax, peff
552
553 INTEGER :: l1, l2, la
554
555 pmin = huge(0.0_dp)
556 pmax = 0.0_dp
557 peff = 0.0_dp
558
559 DO l1 = 0, lmax
560 DO l2 = l1, lmax
561 DO la = l2 - l1, l2 + l1
562 pmax(la) = max(pmax(la), zmax(l1) + zmax(l2))
563 pmin(la) = min(pmin(la), zmin(l1) + zmin(l2))
564 peff(la) = max(peff(la), zeff(l1) + zeff(l2))
565 END DO
566 END DO
567 END DO
568
569 END SUBROUTINE get_basis_products
570! **************************************************************************************************
571!> \brief ...
572!> \param lm ...
573!> \param npgf ...
574!> \param nfun ...
575!> \param zet ...
576!> \param gcc ...
577!> \param nfit ...
578!> \param afit ...
579!> \param amet ...
580!> \param eval ...
581! **************************************************************************************************
582 SUBROUTINE overlap_maximum(lm, npgf, nfun, zet, gcc, nfit, afit, amet, eval)
583 INTEGER, INTENT(IN) :: lm, npgf, nfun
584 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zet
585 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gcc
586 INTEGER, INTENT(IN) :: nfit
587 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: afit
588 REAL(kind=dp), INTENT(IN) :: amet
589 REAL(kind=dp), INTENT(OUT) :: eval
590
591 INTEGER :: i, ia, ib, info
592 REAL(kind=dp) :: fij, fxij, intab, p, xij
593 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fx, tx, x2, xx
594
595 ! SUM_i(fi M fi)
596 fij = 0.0_dp
597 DO ia = 1, npgf
598 DO ib = 1, npgf
599 p = zet(ia) + zet(ib) + amet
600 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
601 DO i = 1, nfun
602 fij = fij + gcc(ia, i)*gcc(ib, i)*intab
603 END DO
604 END DO
605 END DO
606
607 !Integrals (fi M xj)
608 ALLOCATE (fx(nfit, nfun), tx(nfit, nfun))
609 fx = 0.0_dp
610 DO ia = 1, npgf
611 DO ib = 1, nfit
612 p = zet(ia) + afit(ib) + amet
613 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
614 DO i = 1, nfun
615 fx(ib, i) = fx(ib, i) + gcc(ia, i)*intab
616 END DO
617 END DO
618 END DO
619
620 !Integrals (xi M xj)
621 ALLOCATE (xx(nfit, nfit), x2(nfit, nfit))
622 DO ia = 1, nfit
623 DO ib = 1, nfit
624 p = afit(ia) + afit(ib) + amet
625 xx(ia, ib) = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
626 END DO
627 END DO
628
629 !Solve for tab
630 tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun)
631 x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit)
632 CALL dposv("U", nfit, nfun, x2, nfit, tx, nfit, info)
633 IF (info == 0) THEN
634 ! value t*xx*t
635 xij = 0.0_dp
636 DO i = 1, nfun
637 xij = xij + dot_product(tx(:, i), matmul(xx, tx(:, i)))
638 END DO
639 ! value t*fx
640 fxij = 0.0_dp
641 DO i = 1, nfun
642 fxij = fxij + dot_product(tx(:, i), fx(:, i))
643 END DO
644 !
645 eval = fij - 2.0_dp*fxij + xij
646 ELSE
647 ! error in solving for max overlap
648 eval = 1.0e10_dp
649 END IF
650
651 DEALLOCATE (fx, xx, x2, tx)
652
653 END SUBROUTINE overlap_maximum
654! **************************************************************************************************
655!> \brief ...
656!> \param x ...
657!> \param n ...
658!> \param eval ...
659! **************************************************************************************************
660 SUBROUTINE neb_potential(x, n, eval)
661 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: x
662 INTEGER, INTENT(IN) :: n
663 REAL(kind=dp), INTENT(INOUT) :: eval
664
665 INTEGER :: i
666
667 DO i = 2, n
668 IF (x(i) < 1.5_dp) THEN
669 eval = eval + 10.0_dp*(1.5_dp - x(i))**2
670 END IF
671 END DO
672
673 END SUBROUTINE neb_potential
674! **************************************************************************************************
675!> \brief ...
676!> \param basis_set ...
677!> \param lin ...
678!> \param np ...
679!> \param nf ...
680!> \param zval ...
681!> \param gcval ...
682! **************************************************************************************************
683 SUBROUTINE get_basis_functions(basis_set, lin, np, nf, zval, gcval)
684 TYPE(gto_basis_set_type), POINTER :: basis_set
685 INTEGER, INTENT(IN) :: lin
686 INTEGER, INTENT(OUT) :: np, nf
687 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: zval
688 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gcval
689
690 INTEGER :: iset, ishell, j1, j2, jf, jp, l, nset
691 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
692 INTEGER, DIMENSION(:, :), POINTER :: lshell
693 LOGICAL :: toadd
694 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
695 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
696
697 CALL get_gto_basis_set(gto_basis_set=basis_set, &
698 nset=nset, &
699 nshell=nshell, &
700 npgf=npgf, &
701 l=lshell, &
702 lmax=lm, &
703 zet=zet, &
704 gcc=gcc)
705
706 np = 0
707 nf = 0
708 DO iset = 1, nset
709 toadd = .true.
710 DO ishell = 1, nshell(iset)
711 l = lshell(ishell, iset)
712 IF (l == lin) THEN
713 nf = nf + 1
714 IF (toadd) THEN
715 np = np + npgf(iset)
716 toadd = .false.
717 END IF
718 END IF
719 END DO
720 END DO
721 ALLOCATE (zval(np), gcval(np, nf))
722 zval = 0.0_dp
723 gcval = 0.0_dp
724 !
725 jp = 0
726 jf = 0
727 DO iset = 1, nset
728 toadd = .true.
729 DO ishell = 1, nshell(iset)
730 l = lshell(ishell, iset)
731 IF (l == lin) THEN
732 jf = jf + 1
733 IF (toadd) THEN
734 j1 = jp + 1
735 j2 = jp + npgf(iset)
736 zval(j1:j2) = zet(1:npgf(iset), iset)
737 jp = jp + npgf(iset)
738 toadd = .false.
739 END IF
740 gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset)
741 END IF
742 END DO
743 END DO
744
745 END SUBROUTINE get_basis_functions
746
747END MODULE auto_basis
auto_basis
Automatic generation of auxiliary basis sets of different kind.
Definition auto_basis.F:15
auto_basis::create_lri_aux_basis_set
subroutine, public create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, exact_1c_terms, tda_kernel)
Create a LRI_AUX basis set using some heuristics.
Definition auto_basis.F:233
auto_basis::create_ri_aux_basis_set
subroutine, public create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl, basis_type, basis_sort)
Create a RI_AUX basis set using some heuristics.
Definition auto_basis.F:57
auto_basis::create_oce_basis
subroutine, public create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
...
Definition auto_basis.F:414
aux_basis_set
Definition aux_basis_set.F:13
aux_basis_set::create_aux_basis
subroutine, public create_aux_basis(aux_basis, bsname, nsets, lmin, lmax, nl, npgf, zet)
create a basis in GTO form
Definition aux_basis_set.F:53
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
basis_set_types::sort_gto_basis_set
subroutine, public sort_gto_basis_set(basis_set, sort_method)
sort basis sets w.r.t. radius
Definition basis_set_types.F:2618
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::stoychev2016
integer, save, public stoychev2016
Definition bibliography.F:36
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::gamma1
real(kind=dp), dimension(0:maxfac), parameter, public gamma1
Definition mathconstants.F:95
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171