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