(git:21ef868)
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-2026 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, kname
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, econf)
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 IF (.NOT. ASSOCIATED(econf)) THEN
110 CALL get_qs_kind(qs_kind, name=kname)
111 CALL cp_abort(__location__, &
112 "AUTO_BASIS RI_AUX cannot process atom kind "// &
113 "<"//trim(adjustl(kname))//"> due to missing "// &
114 "definition of potential or electron configuration; "// &
115 "consider setting keyword ELEC_CONF explicitly for "// &
116 "GHOST atom kind that has assigned a basis set")
117 END IF
118 CALL get_ptable_info(element_symbol, ielement=z)
119 lval = 0
120 DO l = 0, maxval(ubound(econf))
121 IF (econf(l) > 0) lval = l
122 END DO
123 IF (sum(econf) /= nint(zval)) THEN
124 cpwarn("Valence charge and electron configuration not consistent")
125 END IF
126 pend = 0.0_dp
127 linc = 1
128 IF (z > 18) linc = 2
129 SELECT CASE (basis_cntrl)
130 CASE (0)
131 laux = max(2*lval, lmax + linc)
132 CASE (1)
133 laux = max(2*lval, lmax + linc)
134 CASE (2)
135 laux = max(2*lval, lmax + linc + 1)
136 CASE (3)
137 laux = max(2*lmax, lmax + linc + 2)
138 CASE DEFAULT
139 cpabort("Invalid value of control variable")
140 END SELECT
141 !
142 DO l = 2*lmax + 1, laux
143 xv = peff(2*lmax)
144 pmin(l) = xv
145 pmax(l) = xv
146 peff(l) = xv
147 pend(l) = xv
148 END DO
149 !
150 DO l = 0, laux
151 IF (l <= 2*lval) THEN
152 pend(l) = min(fv(l)*peff(l), pmax(l))
153 bval(l) = 1.8_dp
154 ELSE
155 pend(l) = peff(l)
156 bval(l) = bv(l)
157 END IF
158 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
159 nval(l) = max(ceiling(xv), 0)
160 END DO
161 ! first set include valence only
162 nsets = 1
163 ls1(1) = 0
164 ls2(1) = lval
165 DO l = lval + 1, laux
166 IF (nval(l) < nval(lval) - 1) EXIT
167 ls2(1) = l
168 END DO
169 ! second set up to 2*lval
170 IF (laux > ls2(1)) THEN
171 IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN
172 nsets = 2
173 ls1(2) = ls2(1) + 1
174 ls2(2) = laux
175 ELSE
176 nsets = 2
177 ls1(2) = ls2(1) + 1
178 ls2(2) = min(2*lval, laux)
179 lx = ls2(2)
180 DO l = lx + 1, laux
181 IF (nval(l) < nval(lx) - 1) EXIT
182 ls2(2) = l
183 END DO
184 IF (laux > ls2(2)) THEN
185 nsets = 3
186 ls1(3) = ls2(2) + 1
187 ls2(3) = laux
188 END IF
189 END IF
190 END IF
191 !
192 amax = 0.0
193 amin = huge(0.0_dp)
194 bmin = huge(0.0_dp)
195 DO i = 1, nsets
196 DO j = ls1(i), ls2(i)
197 amax(i) = max(amax(i), pend(j))
198 amin(i) = min(amin(i), pmin(j))
199 bmin(i) = min(bmin(i), bval(j))
200 END DO
201 xv = log(amax(i)/amin(i))/log(bmin(i)) + 1.e-10_dp
202 npgf(i) = max(ceiling(xv), 0)
203 END DO
204 nx = maxval(npgf(1:nsets))
205 ALLOCATE (zet(nx, nsets))
206 zet = 0.0_dp
207 nl = 0
208 DO i = 1, nsets
209 DO j = 1, npgf(i)
210 jj = npgf(i) - j + 1
211 zet(jj, i) = amin(i)*bmin(i)**(j - 1)
212 END DO
213 DO l = ls1(i), ls2(i)
214 nl(l, i) = nval(l)
215 END DO
216 END DO
217 bsname = trim(element_symbol)//"-RI-AUX-"//trim(orb_basis_set%name)
218 !
219 CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
220
221 DEALLOCATE (zet)
222
223 IF (PRESENT(basis_sort)) THEN
224 CALL sort_gto_basis_set(ri_aux_basis_set, basis_sort)
225 END IF
226
227 END IF
228
229 END SUBROUTINE create_ri_aux_basis_set
230! **************************************************************************************************
231!> \brief Create a LRI_AUX basis set using some heuristics
232!> \param lri_aux_basis_set ...
233!> \param qs_kind ...
234!> \param basis_cntrl ...
235!> \param exact_1c_terms ...
236!> \param tda_kernel ...
237!> \date 01.11.2017
238!> \author JGH
239! **************************************************************************************************
240 SUBROUTINE create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, &
241 exact_1c_terms, tda_kernel)
242 TYPE(gto_basis_set_type), POINTER :: lri_aux_basis_set
243 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
244 INTEGER, INTENT(IN) :: basis_cntrl
245 LOGICAL, INTENT(IN), OPTIONAL :: exact_1c_terms, tda_kernel
246
247 CHARACTER(LEN=2) :: element_symbol
248 CHARACTER(LEN=default_string_length) :: bsname, kname
249 INTEGER :: i, j, l, laux, linc, lm, lmax, lval, n1, &
250 n2, nsets, z
251 INTEGER, DIMENSION(0:18) :: nval
252 INTEGER, DIMENSION(0:9, 1:50) :: nl
253 INTEGER, DIMENSION(1:50) :: ls1, ls2, npgf
254 INTEGER, DIMENSION(:), POINTER :: econf
255 LOGICAL :: e1terms, kernel_basis
256 REAL(kind=dp) :: xv, zval
257 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
258 REAL(kind=dp), DIMENSION(0:18) :: bval, peff, pend, pmax, pmin
259 REAL(kind=dp), DIMENSION(0:9) :: zeff, zmax, zmin
260 REAL(kind=dp), DIMENSION(4) :: bv, bx
261 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
262
263 !
264 IF (PRESENT(exact_1c_terms)) THEN
265 e1terms = exact_1c_terms
266 ELSE
267 e1terms = .false.
268 END IF
269 IF (PRESENT(tda_kernel)) THEN
270 kernel_basis = tda_kernel
271 ELSE
272 kernel_basis = .false.
273 END IF
274 IF (kernel_basis .AND. e1terms) THEN
275 CALL cp_warn(__location__, "LRI Kernel basis generation will ignore exact 1C term option.")
276 END IF
277 !
278 cpassert(.NOT. ASSOCIATED(lri_aux_basis_set))
279 NULLIFY (orb_basis_set, econf)
280 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
281 IF (ASSOCIATED(orb_basis_set)) THEN
282 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
283 CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
284 CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
285 IF (.NOT. ASSOCIATED(econf)) THEN
286 CALL get_qs_kind(qs_kind, name=kname)
287 CALL cp_abort(__location__, &
288 "AUTO_BASIS LRI_AUX cannot process atom kind "// &
289 "<"//trim(adjustl(kname))//"> due to missing "// &
290 "definition of potential or electron configuration; "// &
291 "consider setting keyword ELEC_CONF explicitly for "// &
292 "GHOST atom kind that has assigned a basis set")
293 END IF
294 CALL get_ptable_info(element_symbol, ielement=z)
295 lval = 0
296 DO l = 0, maxval(ubound(econf))
297 IF (econf(l) > 0) lval = l
298 END DO
299 IF (sum(econf) /= nint(zval)) THEN
300 cpwarn("Valence charge and electron configuration not consistent")
301 END IF
302 !
303 linc = 1
304 IF (z > 18) linc = 2
305 pend = 0.0_dp
306 IF (kernel_basis) THEN
307 bv(1:4) = [3.20_dp, 2.80_dp, 2.40_dp, 2.00_dp]
308 bx(1:4) = [4.00_dp, 3.50_dp, 3.00_dp, 2.50_dp]
309 !
310 SELECT CASE (basis_cntrl)
311 CASE (0)
312 laux = lval + 1
313 CASE (1)
314 laux = max(lval + 1, lmax)
315 CASE (2)
316 laux = max(lval + 2, lmax + 1)
317 CASE (3)
318 laux = max(lval + 3, lmax + 2)
319 laux = min(laux, 2 + linc)
320 CASE DEFAULT
321 cpabort("Invalid value of control variable")
322 END SELECT
323 ELSE
324 bv(1:4) = [2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp]
325 bx(1:4) = [2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp]
326 !
327 SELECT CASE (basis_cntrl)
328 CASE (0)
329 laux = max(2*lval, lmax + linc)
330 laux = min(laux, 2 + linc)
331 CASE (1)
332 laux = max(2*lval, lmax + linc)
333 laux = min(laux, 3 + linc)
334 CASE (2)
335 laux = max(2*lval, lmax + linc + 1)
336 laux = min(laux, 4 + linc)
337 CASE (3)
338 laux = max(2*lval, lmax + linc + 1)
339 laux = min(laux, 4 + linc)
340 CASE DEFAULT
341 cpabort("Invalid value of control variable")
342 END SELECT
343 END IF
344 !
345 DO l = 2*lmax + 1, laux
346 pmin(l) = pmin(2*lmax)
347 pmax(l) = pmax(2*lmax)
348 peff(l) = peff(2*lmax)
349 END DO
350 !
351 nval = 0
352 IF (exact_1c_terms) THEN
353 DO l = 0, laux
354 IF (l <= lval + 1) THEN
355 pend(l) = zmax(l) + 1.0_dp
356 bval(l) = bv(basis_cntrl + 1)
357 ELSE
358 pend(l) = 2.0_dp*peff(l)
359 bval(l) = bx(basis_cntrl + 1)
360 END IF
361 pmin(l) = zmin(l)
362 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
363 nval(l) = max(ceiling(xv), 0)
364 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
365 END DO
366 ELSE
367 DO l = 0, laux
368 IF (l <= lval + 1) THEN
369 pend(l) = pmax(l)
370 bval(l) = bv(basis_cntrl + 1)
371 pmin(l) = zmin(l)
372 ELSE
373 pend(l) = 4.0_dp*peff(l)
374 bval(l) = bx(basis_cntrl + 1)
375 END IF
376 xv = log(pend(l)/pmin(l))/log(bval(l)) + 1.e-10_dp
377 nval(l) = max(ceiling(xv), 0)
378 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
379 END DO
380 END IF
381 !
382 lm = min(2*lval, 3)
383 n1 = maxval(nval(0:lm))
384 IF (laux < lm + 1) THEN
385 n2 = 0
386 ELSE
387 n2 = maxval(nval(lm + 1:laux))
388 END IF
389 !
390 nsets = n1 + n2
391 ALLOCATE (zet(1, nsets))
392 zet = 0.0_dp
393 nl = 0
394 j = maxval(maxloc(nval(0:lm)))
395 DO i = 1, n1
396 ls1(i) = 0
397 ls2(i) = lm
398 npgf(i) = 1
399 zet(1, i) = pmin(j)*bval(j)**(i - 1)
400 DO l = 0, lm
401 nl(l, i) = 1
402 END DO
403 END DO
404 j = lm + 1
405 DO i = n1 + 1, nsets
406 ls1(i) = lm + 1
407 ls2(i) = laux
408 npgf(i) = 1
409 zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1)
410 DO l = lm + 1, laux
411 nl(l, i) = 1
412 END DO
413 END DO
414 !
415 bsname = trim(element_symbol)//"-LRI-AUX-"//trim(orb_basis_set%name)
416 !
417 CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
418 !
419 DEALLOCATE (zet)
420 END IF
421
422 END SUBROUTINE create_lri_aux_basis_set
423
424! **************************************************************************************************
425!> \brief ...
426!> \param oce_basis ...
427!> \param orb_basis ...
428!> \param lmax_oce ...
429!> \param nbas_oce ...
430! **************************************************************************************************
431 SUBROUTINE create_oce_basis(oce_basis, orb_basis, lmax_oce, nbas_oce)
432 TYPE(gto_basis_set_type), POINTER :: oce_basis, orb_basis
433 INTEGER, INTENT(IN) :: lmax_oce, nbas_oce
434
435 CHARACTER(LEN=default_string_length) :: bsname
436 INTEGER :: i, l, lmax, lx, nset, nx
437 INTEGER, ALLOCATABLE, DIMENSION(:) :: lmin, lset, npgf
438 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl
439 INTEGER, DIMENSION(:), POINTER :: npgf_orb
440 REAL(kind=dp) :: cval, x, z0, z1
441 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
442 REAL(kind=dp), DIMENSION(0:9) :: zeff, zmax, zmin
443
444 CALL get_basis_keyfigures(orb_basis, lmax, zmin, zmax, zeff)
445 IF (nbas_oce < 1) THEN
446 CALL get_gto_basis_set(gto_basis_set=orb_basis, nset=nset, npgf=npgf_orb)
447 nx = sum(npgf_orb(1:nset))
448 ELSE
449 nx = 0
450 END IF
451 nset = max(nbas_oce, nx)
452 lx = max(lmax_oce, lmax)
453 !
454 bsname = "OCE-"//trim(orb_basis%name)
455 ALLOCATE (lmin(nset), lset(nset), nl(0:9, nset), npgf(nset), zet(1, nset))
456 lmin = 0
457 lset = 0
458 nl = 1
459 npgf = 1
460 zet = 0.0_dp
461 !
462 z0 = minval(zmin(0:lmax))
463 z1 = maxval(zmax(0:lmax))
464 x = 1.0_dp/real(nset - 1, kind=dp)
465 cval = (z1/z0)**x
466 zet(1, nset) = z0
467 DO i = nset - 1, 1, -1
468 zet(1, i) = zet(1, i + 1)*cval
469 END DO
470 DO i = 1, nset
471 x = zet(1, i)
472 DO l = 1, lmax
473 z1 = 1.05_dp*zmax(l)
474 IF (x < z1) lset(i) = l
475 END DO
476 IF (lset(i) == lmax) lset(i) = lx
477 END DO
478 !
479 CALL create_aux_basis(oce_basis, bsname, nset, lmin, lset, nl, npgf, zet)
480 !
481 DEALLOCATE (lmin, lset, nl, npgf, zet)
482
483 END SUBROUTINE create_oce_basis
484! **************************************************************************************************
485!> \brief ...
486!> \param basis_set ...
487!> \param lmax ...
488!> \param zmin ...
489!> \param zmax ...
490!> \param zeff ...
491! **************************************************************************************************
492 SUBROUTINE get_basis_keyfigures(basis_set, lmax, zmin, zmax, zeff)
493 TYPE(gto_basis_set_type), POINTER :: basis_set
494 INTEGER, INTENT(OUT) :: lmax
495 REAL(kind=dp), DIMENSION(0:9), INTENT(OUT) :: zmin, zmax, zeff
496
497 INTEGER :: i, ipgf, iset, ishell, j, l, nset
498 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
499 INTEGER, DIMENSION(:, :), POINTER :: lshell
500 REAL(kind=dp) :: aeff, gcca, gccb, kval, rexp, rint, rno, &
501 zeta
502 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
503 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
504
505 CALL get_gto_basis_set(gto_basis_set=basis_set, &
506 nset=nset, &
507 nshell=nshell, &
508 npgf=npgf, &
509 l=lshell, &
510 lmax=lm, &
511 zet=zet, &
512 gcc=gcc)
513
514 lmax = maxval(lm)
515 cpassert(lmax <= 9)
516
517 zmax = 0.0_dp
518 zmin = huge(0.0_dp)
519 zeff = 0.0_dp
520
521 DO iset = 1, nset
522 ! zmin zmax
523 DO ipgf = 1, npgf(iset)
524 DO ishell = 1, nshell(iset)
525 l = lshell(ishell, iset)
526 zeta = zet(ipgf, iset)
527 zmax(l) = max(zmax(l), zeta)
528 zmin(l) = min(zmin(l), zeta)
529 END DO
530 END DO
531 ! zeff
532 DO ishell = 1, nshell(iset)
533 l = lshell(ishell, iset)
534 kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2)
535 rexp = 0.0_dp
536 rno = 0.0_dp
537 DO i = 1, npgf(iset)
538 gcca = gcc(i, ishell, iset)
539 DO j = 1, npgf(iset)
540 zeta = zet(i, iset) + zet(j, iset)
541 gccb = gcc(j, ishell, iset)
542 rint = 0.5_dp*fac(l + 1)/zeta**(l + 2)
543 rexp = rexp + gcca*gccb*rint
544 rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp)
545 rno = rno + gcca*gccb*rint
546 END DO
547 END DO
548 rexp = rexp/rno
549 aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2)
550 zeff(l) = max(zeff(l), aeff)
551 END DO
552 END DO
553
554 END SUBROUTINE get_basis_keyfigures
555
556! **************************************************************************************************
557!> \brief ...
558!> \param lmax ...
559!> \param zmin ...
560!> \param zmax ...
561!> \param zeff ...
562!> \param pmin ...
563!> \param pmax ...
564!> \param peff ...
565! **************************************************************************************************
566 SUBROUTINE get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
567 INTEGER, INTENT(IN) :: lmax
568 REAL(kind=dp), DIMENSION(0:9), INTENT(IN) :: zmin, zmax, zeff
569 REAL(kind=dp), DIMENSION(0:18), INTENT(OUT) :: pmin, pmax, peff
570
571 INTEGER :: l1, l2, la
572
573 pmin = huge(0.0_dp)
574 pmax = 0.0_dp
575 peff = 0.0_dp
576
577 DO l1 = 0, lmax
578 DO l2 = l1, lmax
579 DO la = l2 - l1, l2 + l1
580 pmax(la) = max(pmax(la), zmax(l1) + zmax(l2))
581 pmin(la) = min(pmin(la), zmin(l1) + zmin(l2))
582 peff(la) = max(peff(la), zeff(l1) + zeff(l2))
583 END DO
584 END DO
585 END DO
586
587 END SUBROUTINE get_basis_products
588! **************************************************************************************************
589!> \brief ...
590!> \param lm ...
591!> \param npgf ...
592!> \param nfun ...
593!> \param zet ...
594!> \param gcc ...
595!> \param nfit ...
596!> \param afit ...
597!> \param amet ...
598!> \param eval ...
599! **************************************************************************************************
600 SUBROUTINE overlap_maximum(lm, npgf, nfun, zet, gcc, nfit, afit, amet, eval)
601 INTEGER, INTENT(IN) :: lm, npgf, nfun
602 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zet
603 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gcc
604 INTEGER, INTENT(IN) :: nfit
605 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: afit
606 REAL(kind=dp), INTENT(IN) :: amet
607 REAL(kind=dp), INTENT(OUT) :: eval
608
609 INTEGER :: i, ia, ib, info
610 REAL(kind=dp) :: fij, fxij, intab, p, xij
611 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fx, tx, x2, xx
612
613 ! SUM_i(fi M fi)
614 fij = 0.0_dp
615 DO ia = 1, npgf
616 DO ib = 1, npgf
617 p = zet(ia) + zet(ib) + amet
618 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
619 DO i = 1, nfun
620 fij = fij + gcc(ia, i)*gcc(ib, i)*intab
621 END DO
622 END DO
623 END DO
624
625 !Integrals (fi M xj)
626 ALLOCATE (fx(nfit, nfun), tx(nfit, nfun))
627 fx = 0.0_dp
628 DO ia = 1, npgf
629 DO ib = 1, nfit
630 p = zet(ia) + afit(ib) + amet
631 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
632 DO i = 1, nfun
633 fx(ib, i) = fx(ib, i) + gcc(ia, i)*intab
634 END DO
635 END DO
636 END DO
637
638 !Integrals (xi M xj)
639 ALLOCATE (xx(nfit, nfit), x2(nfit, nfit))
640 DO ia = 1, nfit
641 DO ib = 1, nfit
642 p = afit(ia) + afit(ib) + amet
643 xx(ia, ib) = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
644 END DO
645 END DO
646
647 !Solve for tab
648 tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun)
649 x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit)
650 CALL dposv("U", nfit, nfun, x2, nfit, tx, nfit, info)
651 IF (info == 0) THEN
652 ! value t*xx*t
653 xij = 0.0_dp
654 DO i = 1, nfun
655 xij = xij + dot_product(tx(:, i), matmul(xx, tx(:, i)))
656 END DO
657 ! value t*fx
658 fxij = 0.0_dp
659 DO i = 1, nfun
660 fxij = fxij + dot_product(tx(:, i), fx(:, i))
661 END DO
662 !
663 eval = fij - 2.0_dp*fxij + xij
664 ELSE
665 ! error in solving for max overlap
666 eval = 1.0e10_dp
667 END IF
668
669 DEALLOCATE (fx, xx, x2, tx)
670
671 END SUBROUTINE overlap_maximum
672! **************************************************************************************************
673!> \brief ...
674!> \param x ...
675!> \param n ...
676!> \param eval ...
677! **************************************************************************************************
678 SUBROUTINE neb_potential(x, n, eval)
679 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: x
680 INTEGER, INTENT(IN) :: n
681 REAL(kind=dp), INTENT(INOUT) :: eval
682
683 INTEGER :: i
684
685 DO i = 2, n
686 IF (x(i) < 1.5_dp) THEN
687 eval = eval + 10.0_dp*(1.5_dp - x(i))**2
688 END IF
689 END DO
690
691 END SUBROUTINE neb_potential
692! **************************************************************************************************
693!> \brief ...
694!> \param basis_set ...
695!> \param lin ...
696!> \param np ...
697!> \param nf ...
698!> \param zval ...
699!> \param gcval ...
700! **************************************************************************************************
701 SUBROUTINE get_basis_functions(basis_set, lin, np, nf, zval, gcval)
702 TYPE(gto_basis_set_type), POINTER :: basis_set
703 INTEGER, INTENT(IN) :: lin
704 INTEGER, INTENT(OUT) :: np, nf
705 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: zval
706 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gcval
707
708 INTEGER :: iset, ishell, j1, j2, jf, jp, l, nset
709 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
710 INTEGER, DIMENSION(:, :), POINTER :: lshell
711 LOGICAL :: toadd
712 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
713 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
714
715 CALL get_gto_basis_set(gto_basis_set=basis_set, &
716 nset=nset, &
717 nshell=nshell, &
718 npgf=npgf, &
719 l=lshell, &
720 lmax=lm, &
721 zet=zet, &
722 gcc=gcc)
723
724 np = 0
725 nf = 0
726 DO iset = 1, nset
727 toadd = .true.
728 DO ishell = 1, nshell(iset)
729 l = lshell(ishell, iset)
730 IF (l == lin) THEN
731 nf = nf + 1
732 IF (toadd) THEN
733 np = np + npgf(iset)
734 toadd = .false.
735 END IF
736 END IF
737 END DO
738 END DO
739 ALLOCATE (zval(np), gcval(np, nf))
740 zval = 0.0_dp
741 gcval = 0.0_dp
742 !
743 jp = 0
744 jf = 0
745 DO iset = 1, nset
746 toadd = .true.
747 DO ishell = 1, nshell(iset)
748 l = lshell(ishell, iset)
749 IF (l == lin) THEN
750 jf = jf + 1
751 IF (toadd) THEN
752 j1 = jp + 1
753 j2 = jp + npgf(iset)
754 zval(j1:j2) = zet(1:npgf(iset), iset)
755 jp = jp + npgf(iset)
756 toadd = .false.
757 END IF
758 gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset)
759 END IF
760 END DO
761 END DO
762
763 END SUBROUTINE get_basis_functions
764
765END 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:242
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:432
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, ccon)
...
Definition basis_set_types.F:643
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:2765
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, cneo_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, 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:491
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:190