Line data Source code
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 : ! **************************************************************************************************
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 292 : 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 292 : INTEGER, DIMENSION(:), POINTER :: econf
71 : REAL(KIND=dp) :: xv, zval
72 292 : 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 292 : 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 292 : 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 292 : 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 292 : CPASSERT(.NOT. ASSOCIATED(ri_aux_basis_set))
87 292 : NULLIFY (orb_basis_set)
88 292 : IF (.NOT. PRESENT(basis_type)) THEN
89 234 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
90 : ELSE
91 58 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type=basis_type)
92 : END IF
93 292 : 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 1310 : 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 0 : "you should state BASIS_SET RI_AUX NONE explicitly in the input.")
102 : END IF
103 292 : 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 292 : CALL init_orbital_pointers(2*lmax)
107 292 : CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
108 292 : CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
109 292 : CALL get_ptable_info(element_symbol, ielement=z)
110 292 : lval = 0
111 1436 : DO l = 0, MAXVAL(UBOUND(econf))
112 1144 : IF (econf(l) > 0) lval = l
113 : END DO
114 1144 : IF (SUM(econf) /= NINT(zval)) THEN
115 0 : CPWARN("Valence charge and electron configuration not consistent")
116 : END IF
117 292 : pend = 0.0_dp
118 292 : linc = 1
119 292 : IF (z > 18) linc = 2
120 440 : SELECT CASE (basis_cntrl)
121 : CASE (0)
122 148 : laux = MAX(2*lval, lmax + linc)
123 : CASE (1)
124 140 : laux = MAX(2*lval, lmax + linc)
125 : CASE (2)
126 0 : laux = MAX(2*lval, lmax + linc + 1)
127 : CASE (3)
128 4 : laux = MAX(2*lmax, lmax + linc + 2)
129 : CASE DEFAULT
130 292 : CPABORT("Invalid value of control variable")
131 : END SELECT
132 : !
133 346 : DO l = 2*lmax + 1, laux
134 54 : xv = peff(2*lmax)
135 54 : pmin(l) = xv
136 54 : pmax(l) = xv
137 54 : peff(l) = xv
138 346 : pend(l) = xv
139 : END DO
140 : !
141 1256 : DO l = 0, laux
142 964 : IF (l <= 2*lval) THEN
143 660 : pend(l) = MIN(fv(l)*peff(l), pmax(l))
144 660 : bval(l) = 1.8_dp
145 : ELSE
146 304 : pend(l) = peff(l)
147 304 : bval(l) = bv(l)
148 : END IF
149 964 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
150 1256 : nval(l) = MAX(CEILING(xv), 0)
151 : END DO
152 : ! first set include valence only
153 292 : nsets = 1
154 292 : ls1(1) = 0
155 292 : ls2(1) = lval
156 442 : DO l = lval + 1, laux
157 416 : IF (nval(l) < nval(lval) - 1) EXIT
158 442 : ls2(1) = l
159 : END DO
160 : ! second set up to 2*lval
161 292 : IF (laux > ls2(1)) THEN
162 266 : IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN
163 266 : nsets = 2
164 266 : ls1(2) = ls2(1) + 1
165 266 : ls2(2) = laux
166 : ELSE
167 0 : nsets = 2
168 0 : ls1(2) = ls2(1) + 1
169 0 : ls2(2) = MIN(2*lval, laux)
170 0 : lx = ls2(2)
171 0 : DO l = lx + 1, laux
172 0 : IF (nval(l) < nval(lx) - 1) EXIT
173 0 : ls2(2) = l
174 : END DO
175 0 : IF (laux > ls2(2)) THEN
176 0 : nsets = 3
177 0 : ls1(3) = ls2(2) + 1
178 0 : ls2(3) = laux
179 : END IF
180 : END IF
181 : END IF
182 : !
183 292 : amax = 0.0
184 1168 : amin = HUGE(0.0_dp)
185 1168 : bmin = HUGE(0.0_dp)
186 850 : DO i = 1, nsets
187 1522 : DO j = ls1(i), ls2(i)
188 964 : amax(i) = MAX(amax(i), pend(j))
189 964 : amin(i) = MIN(amin(i), pmin(j))
190 1522 : bmin(i) = MIN(bmin(i), bval(j))
191 : END DO
192 558 : xv = LOG(amax(i)/amin(i))/LOG(bmin(i)) + 1.e-10_dp
193 850 : npgf(i) = MAX(CEILING(xv), 0)
194 : END DO
195 850 : nx = MAXVAL(npgf(1:nsets))
196 1168 : ALLOCATE (zet(nx, nsets))
197 6018 : zet = 0.0_dp
198 292 : nl = 0
199 850 : DO i = 1, nsets
200 3712 : DO j = 1, npgf(i)
201 3154 : jj = npgf(i) - j + 1
202 3712 : zet(jj, i) = amin(i)*bmin(i)**(j - 1)
203 : END DO
204 1814 : DO l = ls1(i), ls2(i)
205 1522 : nl(l, i) = nval(l)
206 : END DO
207 : END DO
208 292 : bsname = TRIM(element_symbol)//"-RI-AUX-"//TRIM(orb_basis_set%name)
209 : !
210 292 : CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
211 :
212 292 : DEALLOCATE (zet)
213 :
214 876 : IF (PRESENT(basis_sort)) THEN
215 186 : CALL sort_gto_basis_set(ri_aux_basis_set, basis_sort)
216 : END IF
217 :
218 : END IF
219 :
220 584 : 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 30 : 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 30 : INTEGER, DIMENSION(:), POINTER :: econf
246 : LOGICAL :: e1terms, kernel_basis
247 : REAL(KIND=dp) :: xv, zval
248 30 : 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 30 : IF (PRESENT(exact_1c_terms)) THEN
256 30 : e1terms = exact_1c_terms
257 : ELSE
258 : e1terms = .FALSE.
259 : END IF
260 30 : IF (PRESENT(tda_kernel)) THEN
261 12 : kernel_basis = tda_kernel
262 : ELSE
263 : kernel_basis = .FALSE.
264 : END IF
265 30 : IF (kernel_basis .AND. e1terms) THEN
266 0 : CALL cp_warn(__LOCATION__, "LRI Kernel basis generation will ignore exact 1C term option.")
267 : END IF
268 : !
269 30 : CPASSERT(.NOT. ASSOCIATED(lri_aux_basis_set))
270 30 : NULLIFY (orb_basis_set)
271 30 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
272 30 : IF (ASSOCIATED(orb_basis_set)) THEN
273 30 : CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
274 30 : CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
275 30 : CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
276 30 : CALL get_ptable_info(element_symbol, ielement=z)
277 30 : lval = 0
278 106 : DO l = 0, MAXVAL(UBOUND(econf))
279 76 : IF (econf(l) > 0) lval = l
280 : END DO
281 76 : IF (SUM(econf) /= NINT(zval)) THEN
282 0 : CPWARN("Valence charge and electron configuration not consistent")
283 : END IF
284 : !
285 30 : linc = 1
286 30 : IF (z > 18) linc = 2
287 30 : pend = 0.0_dp
288 30 : IF (kernel_basis) THEN
289 12 : bv(1:4) = (/3.20_dp, 2.80_dp, 2.40_dp, 2.00_dp/)
290 12 : bx(1:4) = (/4.00_dp, 3.50_dp, 3.00_dp, 2.50_dp/)
291 : !
292 12 : SELECT CASE (basis_cntrl)
293 : CASE (0)
294 0 : laux = lval + 1
295 : CASE (1)
296 12 : laux = MAX(lval + 1, lmax)
297 : CASE (2)
298 0 : laux = MAX(lval + 2, lmax + 1)
299 : CASE (3)
300 0 : laux = MAX(lval + 3, lmax + 2)
301 0 : laux = MIN(laux, 2 + linc)
302 : CASE DEFAULT
303 12 : CPABORT("Invalid value of control variable")
304 : END SELECT
305 : ELSE
306 18 : bv(1:4) = (/2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp/)
307 18 : bx(1:4) = (/2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp/)
308 : !
309 18 : SELECT CASE (basis_cntrl)
310 : CASE (0)
311 0 : laux = MAX(2*lval, lmax + linc)
312 0 : laux = MIN(laux, 2 + linc)
313 : CASE (1)
314 18 : laux = MAX(2*lval, lmax + linc)
315 18 : laux = MIN(laux, 3 + linc)
316 : CASE (2)
317 0 : laux = MAX(2*lval, lmax + linc + 1)
318 0 : laux = MIN(laux, 4 + linc)
319 : CASE (3)
320 0 : laux = MAX(2*lval, lmax + linc + 1)
321 0 : laux = MIN(laux, 4 + linc)
322 : CASE DEFAULT
323 18 : CPABORT("Invalid value of control variable")
324 : END SELECT
325 : END IF
326 : !
327 30 : DO l = 2*lmax + 1, laux
328 0 : pmin(l) = pmin(2*lmax)
329 0 : pmax(l) = pmax(2*lmax)
330 30 : peff(l) = peff(2*lmax)
331 : END DO
332 : !
333 30 : nval = 0
334 30 : IF (exact_1c_terms) THEN
335 0 : DO l = 0, laux
336 0 : IF (l <= lval + 1) THEN
337 0 : pend(l) = zmax(l) + 1.0_dp
338 0 : bval(l) = bv(basis_cntrl + 1)
339 : ELSE
340 0 : pend(l) = 2.0_dp*peff(l)
341 0 : bval(l) = bx(basis_cntrl + 1)
342 : END IF
343 0 : pmin(l) = zmin(l)
344 0 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
345 0 : nval(l) = MAX(CEILING(xv), 0)
346 0 : bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
347 : END DO
348 : ELSE
349 124 : DO l = 0, laux
350 94 : IF (l <= lval + 1) THEN
351 76 : pend(l) = pmax(l)
352 76 : bval(l) = bv(basis_cntrl + 1)
353 76 : pmin(l) = zmin(l)
354 : ELSE
355 18 : pend(l) = 4.0_dp*peff(l)
356 18 : bval(l) = bx(basis_cntrl + 1)
357 : END IF
358 94 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
359 94 : nval(l) = MAX(CEILING(xv), 0)
360 124 : bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
361 : END DO
362 : END IF
363 : !
364 30 : lm = MIN(2*lval, 3)
365 92 : n1 = MAXVAL(nval(0:lm))
366 30 : IF (laux < lm + 1) THEN
367 : n2 = 0
368 : ELSE
369 56 : n2 = MAXVAL(nval(lm + 1:laux))
370 : END IF
371 : !
372 30 : nsets = n1 + n2
373 90 : ALLOCATE (zet(1, nsets))
374 830 : zet = 0.0_dp
375 30 : nl = 0
376 152 : j = MAXVAL(MAXLOC(nval(0:lm)))
377 280 : DO i = 1, n1
378 250 : ls1(i) = 0
379 250 : ls2(i) = lm
380 250 : npgf(i) = 1
381 250 : zet(1, i) = pmin(j)*bval(j)**(i - 1)
382 786 : DO l = 0, lm
383 756 : nl(l, i) = 1
384 : END DO
385 : END DO
386 30 : j = lm + 1
387 180 : DO i = n1 + 1, nsets
388 150 : ls1(i) = lm + 1
389 150 : ls2(i) = laux
390 150 : npgf(i) = 1
391 150 : zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1)
392 418 : DO l = lm + 1, laux
393 388 : nl(l, i) = 1
394 : END DO
395 : END DO
396 : !
397 30 : bsname = TRIM(element_symbol)//"-LRI-AUX-"//TRIM(orb_basis_set%name)
398 : !
399 30 : CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
400 : !
401 90 : DEALLOCATE (zet)
402 : END IF
403 :
404 60 : 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 12 : 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 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lmin, lset, npgf
420 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl
421 12 : INTEGER, DIMENSION(:), POINTER :: npgf_orb
422 : REAL(KIND=dp) :: cval, x, z0, z1
423 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
424 : REAL(KIND=dp), DIMENSION(0:9) :: zeff, zmax, zmin
425 :
426 12 : CALL get_basis_keyfigures(orb_basis, lmax, zmin, zmax, zeff)
427 12 : IF (nbas_oce < 1) THEN
428 12 : CALL get_gto_basis_set(gto_basis_set=orb_basis, nset=nset, npgf=npgf_orb)
429 54 : nx = SUM(npgf_orb(1:nset))
430 : ELSE
431 : nx = 0
432 : END IF
433 12 : nset = MAX(nbas_oce, nx)
434 12 : lx = MAX(lmax_oce, lmax)
435 : !
436 12 : bsname = "OCE-"//TRIM(orb_basis%name)
437 108 : ALLOCATE (lmin(nset), lset(nset), nl(0:9, nset), npgf(nset), zet(1, nset))
438 114 : lmin = 0
439 114 : lset = 0
440 1134 : nl = 1
441 114 : npgf = 1
442 216 : zet = 0.0_dp
443 : !
444 56 : z0 = MINVAL(zmin(0:lmax))
445 56 : z1 = MAXVAL(zmax(0:lmax))
446 12 : x = 1.0_dp/REAL(nset - 1, KIND=dp)
447 12 : cval = (z1/z0)**x
448 12 : zet(1, nset) = z0
449 102 : DO i = nset - 1, 1, -1
450 102 : zet(1, i) = zet(1, i + 1)*cval
451 : END DO
452 114 : DO i = 1, nset
453 102 : x = zet(1, i)
454 284 : DO l = 1, lmax
455 182 : z1 = 1.05_dp*zmax(l)
456 284 : IF (x < z1) lset(i) = l
457 : END DO
458 114 : IF (lset(i) == lmax) lset(i) = lx
459 : END DO
460 : !
461 12 : CALL create_aux_basis(oce_basis, bsname, nset, lmin, lset, nl, npgf, zet)
462 : !
463 12 : DEALLOCATE (lmin, lset, nl, npgf, zet)
464 :
465 12 : 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 334 : 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 334 : INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
481 334 : INTEGER, DIMENSION(:, :), POINTER :: lshell
482 : REAL(KIND=dp) :: aeff, gcca, gccb, kval, rexp, rint, rno, &
483 : zeta
484 334 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
485 334 : 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 334 : gcc=gcc)
495 :
496 1432 : lmax = MAXVAL(lm)
497 334 : CPASSERT(lmax <= 9)
498 :
499 334 : zmax = 0.0_dp
500 3674 : zmin = HUGE(0.0_dp)
501 334 : zeff = 0.0_dp
502 :
503 1432 : DO iset = 1, nset
504 : ! zmin zmax
505 3564 : DO ipgf = 1, npgf(iset)
506 7602 : DO ishell = 1, nshell(iset)
507 4038 : l = lshell(ishell, iset)
508 4038 : zeta = zet(ipgf, iset)
509 4038 : zmax(l) = MAX(zmax(l), zeta)
510 6504 : zmin(l) = MIN(zmin(l), zeta)
511 : END DO
512 : END DO
513 : ! zeff
514 2932 : DO ishell = 1, nshell(iset)
515 1500 : l = lshell(ishell, iset)
516 1500 : kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2)
517 1500 : rexp = 0.0_dp
518 1500 : rno = 0.0_dp
519 5538 : DO i = 1, npgf(iset)
520 4038 : gcca = gcc(i, ishell, iset)
521 22036 : DO j = 1, npgf(iset)
522 16498 : zeta = zet(i, iset) + zet(j, iset)
523 16498 : gccb = gcc(j, ishell, iset)
524 16498 : rint = 0.5_dp*fac(l + 1)/zeta**(l + 2)
525 16498 : rexp = rexp + gcca*gccb*rint
526 16498 : rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp)
527 20536 : rno = rno + gcca*gccb*rint
528 : END DO
529 : END DO
530 1500 : rexp = rexp/rno
531 1500 : aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2)
532 2598 : zeff(l) = MAX(zeff(l), aeff)
533 : END DO
534 : END DO
535 :
536 334 : 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 322 : 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 6440 : pmin = HUGE(0.0_dp)
556 322 : pmax = 0.0_dp
557 322 : peff = 0.0_dp
558 :
559 1062 : DO l1 = 0, lmax
560 2368 : DO l2 = l1, lmax
561 4788 : DO la = l2 - l1, l2 + l1
562 2742 : pmax(la) = MAX(pmax(la), zmax(l1) + zmax(l2))
563 2742 : pmin(la) = MIN(pmin(la), zmin(l1) + zmin(l2))
564 4048 : peff(la) = MAX(peff(la), zeff(l1) + zeff(l2))
565 : END DO
566 : END DO
567 : END DO
568 :
569 322 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fx, tx, x2, xx
594 :
595 : ! SUM_i(fi M fi)
596 0 : fij = 0.0_dp
597 0 : DO ia = 1, npgf
598 0 : DO ib = 1, npgf
599 0 : p = zet(ia) + zet(ib) + amet
600 0 : intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
601 0 : DO i = 1, nfun
602 0 : 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 0 : ALLOCATE (fx(nfit, nfun), tx(nfit, nfun))
609 0 : fx = 0.0_dp
610 0 : DO ia = 1, npgf
611 0 : DO ib = 1, nfit
612 0 : p = zet(ia) + afit(ib) + amet
613 0 : intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
614 0 : DO i = 1, nfun
615 0 : 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 0 : ALLOCATE (xx(nfit, nfit), x2(nfit, nfit))
622 0 : DO ia = 1, nfit
623 0 : DO ib = 1, nfit
624 0 : p = afit(ia) + afit(ib) + amet
625 0 : 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 0 : tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun)
631 0 : x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit)
632 0 : CALL dposv("U", nfit, nfun, x2, nfit, tx, nfit, info)
633 0 : IF (info == 0) THEN
634 : ! value t*xx*t
635 : xij = 0.0_dp
636 0 : DO i = 1, nfun
637 0 : xij = xij + DOT_PRODUCT(tx(:, i), MATMUL(xx, tx(:, i)))
638 : END DO
639 : ! value t*fx
640 : fxij = 0.0_dp
641 0 : DO i = 1, nfun
642 0 : fxij = fxij + DOT_PRODUCT(tx(:, i), fx(:, i))
643 : END DO
644 : !
645 0 : eval = fij - 2.0_dp*fxij + xij
646 : ELSE
647 : ! error in solving for max overlap
648 0 : eval = 1.0e10_dp
649 : END IF
650 :
651 0 : DEALLOCATE (fx, xx, x2, tx)
652 :
653 0 : END SUBROUTINE overlap_maximum
654 : ! **************************************************************************************************
655 : !> \brief ...
656 : !> \param x ...
657 : !> \param n ...
658 : !> \param eval ...
659 : ! **************************************************************************************************
660 0 : 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 0 : DO i = 2, n
668 0 : IF (x(i) < 1.5_dp) THEN
669 0 : eval = eval + 10.0_dp*(1.5_dp - x(i))**2
670 : END IF
671 : END DO
672 :
673 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
692 0 : INTEGER, DIMENSION(:, :), POINTER :: lshell
693 : LOGICAL :: toadd
694 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
695 0 : 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 0 : gcc=gcc)
705 :
706 0 : np = 0
707 0 : nf = 0
708 0 : DO iset = 1, nset
709 0 : toadd = .TRUE.
710 0 : DO ishell = 1, nshell(iset)
711 0 : l = lshell(ishell, iset)
712 0 : IF (l == lin) THEN
713 0 : nf = nf + 1
714 0 : IF (toadd) THEN
715 0 : np = np + npgf(iset)
716 0 : toadd = .FALSE.
717 : END IF
718 : END IF
719 : END DO
720 : END DO
721 0 : ALLOCATE (zval(np), gcval(np, nf))
722 0 : zval = 0.0_dp
723 0 : gcval = 0.0_dp
724 : !
725 : jp = 0
726 : jf = 0
727 0 : DO iset = 1, nset
728 0 : toadd = .TRUE.
729 0 : DO ishell = 1, nshell(iset)
730 0 : l = lshell(ishell, iset)
731 0 : IF (l == lin) THEN
732 0 : jf = jf + 1
733 0 : IF (toadd) THEN
734 0 : j1 = jp + 1
735 0 : j2 = jp + npgf(iset)
736 0 : zval(j1:j2) = zet(1:npgf(iset), iset)
737 0 : jp = jp + npgf(iset)
738 0 : toadd = .FALSE.
739 : END IF
740 0 : gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset)
741 : END IF
742 : END DO
743 : END DO
744 :
745 0 : END SUBROUTINE get_basis_functions
746 :
747 0 : END MODULE auto_basis
|