(git:e7e05ae)
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 ! **************************************************************************************************
15 MODULE 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 
44 CONTAINS
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 
738 END 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