Line data Source code
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 266 : 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 266 : INTEGER, DIMENSION(:), POINTER :: econf
71 : REAL(KIND=dp) :: xv, zval
72 266 : 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 266 : 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 266 : 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 266 : 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 266 : CPASSERT(.NOT. ASSOCIATED(ri_aux_basis_set))
87 266 : NULLIFY (orb_basis_set)
88 266 : IF (.NOT. PRESENT(basis_type)) THEN
89 210 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
90 : ELSE
91 56 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type=basis_type)
92 : END IF
93 266 : IF (ASSOCIATED(orb_basis_set)) THEN
94 266 : 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 266 : CALL init_orbital_pointers(2*lmax)
98 266 : CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
99 266 : CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
100 266 : CALL get_ptable_info(element_symbol, ielement=z)
101 266 : lval = 0
102 1316 : DO l = 0, MAXVAL(UBOUND(econf))
103 1050 : IF (econf(l) > 0) lval = l
104 : END DO
105 1050 : IF (SUM(econf) /= NINT(zval)) THEN
106 0 : CPWARN("Valence charge and electron configuration not consistent")
107 : END IF
108 266 : pend = 0.0_dp
109 266 : linc = 1
110 266 : IF (z > 18) linc = 2
111 398 : SELECT CASE (basis_cntrl)
112 : CASE (0)
113 132 : laux = MAX(2*lval, lmax + linc)
114 : CASE (1)
115 130 : laux = MAX(2*lval, lmax + linc)
116 : CASE (2)
117 0 : laux = MAX(2*lval, lmax + linc + 1)
118 : CASE (3)
119 4 : laux = MAX(2*lmax, lmax + linc + 2)
120 : CASE DEFAULT
121 266 : CPABORT("Invalid value of control variable")
122 : END SELECT
123 : !
124 312 : DO l = 2*lmax + 1, laux
125 46 : xv = peff(2*lmax)
126 46 : pmin(l) = xv
127 46 : pmax(l) = xv
128 46 : peff(l) = xv
129 312 : pend(l) = xv
130 : END DO
131 : !
132 1150 : DO l = 0, laux
133 884 : IF (l <= 2*lval) THEN
134 602 : pend(l) = MIN(fv(l)*peff(l), pmax(l))
135 602 : bval(l) = 1.8_dp
136 : ELSE
137 282 : pend(l) = peff(l)
138 282 : bval(l) = bv(l)
139 : END IF
140 884 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
141 1150 : nval(l) = MAX(CEILING(xv), 0)
142 : END DO
143 : ! first set include valence only
144 266 : nsets = 1
145 266 : ls1(1) = 0
146 266 : ls2(1) = lval
147 406 : DO l = lval + 1, laux
148 380 : IF (nval(l) < nval(lval) - 1) EXIT
149 406 : ls2(1) = l
150 : END DO
151 : ! second set up to 2*lval
152 266 : IF (laux > ls2(1)) THEN
153 240 : IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN
154 240 : nsets = 2
155 240 : ls1(2) = ls2(1) + 1
156 240 : ls2(2) = laux
157 : ELSE
158 0 : nsets = 2
159 0 : ls1(2) = ls2(1) + 1
160 0 : ls2(2) = MIN(2*lval, laux)
161 0 : lx = ls2(2)
162 0 : DO l = lx + 1, laux
163 0 : IF (nval(l) < nval(lx) - 1) EXIT
164 0 : ls2(2) = l
165 : END DO
166 0 : IF (laux > ls2(2)) THEN
167 0 : nsets = 3
168 0 : ls1(3) = ls2(2) + 1
169 0 : ls2(3) = laux
170 : END IF
171 : END IF
172 : END IF
173 : !
174 266 : amax = 0.0
175 1064 : amin = HUGE(0.0_dp)
176 1064 : bmin = HUGE(0.0_dp)
177 772 : DO i = 1, nsets
178 1390 : DO j = ls1(i), ls2(i)
179 884 : amax(i) = MAX(amax(i), pend(j))
180 884 : amin(i) = MIN(amin(i), pmin(j))
181 1390 : bmin(i) = MIN(bmin(i), bval(j))
182 : END DO
183 506 : xv = LOG(amax(i)/amin(i))/LOG(bmin(i)) + 1.e-10_dp
184 772 : npgf(i) = MAX(CEILING(xv), 0)
185 : END DO
186 772 : nx = MAXVAL(npgf(1:nsets))
187 1064 : ALLOCATE (zet(nx, nsets))
188 5428 : zet = 0.0_dp
189 266 : nl = 0
190 772 : DO i = 1, nsets
191 3346 : DO j = 1, npgf(i)
192 2840 : jj = npgf(i) - j + 1
193 3346 : zet(jj, i) = amin(i)*bmin(i)**(j - 1)
194 : END DO
195 1656 : DO l = ls1(i), ls2(i)
196 1390 : nl(l, i) = nval(l)
197 : END DO
198 : END DO
199 266 : bsname = TRIM(element_symbol)//"-RI-AUX-"//TRIM(orb_basis_set%name)
200 : !
201 266 : CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
202 :
203 266 : DEALLOCATE (zet)
204 :
205 798 : IF (PRESENT(basis_sort)) THEN
206 174 : CALL sort_gto_basis_set(ri_aux_basis_set, basis_sort)
207 : END IF
208 :
209 : END IF
210 :
211 532 : 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 30 : 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 30 : INTEGER, DIMENSION(:), POINTER :: econf
237 : LOGICAL :: e1terms, kernel_basis
238 : REAL(KIND=dp) :: xv, zval
239 30 : 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 30 : IF (PRESENT(exact_1c_terms)) THEN
247 30 : e1terms = exact_1c_terms
248 : ELSE
249 : e1terms = .FALSE.
250 : END IF
251 30 : IF (PRESENT(tda_kernel)) THEN
252 12 : kernel_basis = tda_kernel
253 : ELSE
254 : kernel_basis = .FALSE.
255 : END IF
256 30 : IF (kernel_basis .AND. e1terms) THEN
257 0 : CALL cp_warn(__LOCATION__, "LRI Kernel basis generation will ignore exact 1C term option.")
258 : END IF
259 : !
260 30 : CPASSERT(.NOT. ASSOCIATED(lri_aux_basis_set))
261 30 : NULLIFY (orb_basis_set)
262 30 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB")
263 30 : IF (ASSOCIATED(orb_basis_set)) THEN
264 30 : CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff)
265 30 : CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff)
266 30 : CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol)
267 30 : CALL get_ptable_info(element_symbol, ielement=z)
268 30 : lval = 0
269 106 : DO l = 0, MAXVAL(UBOUND(econf))
270 76 : IF (econf(l) > 0) lval = l
271 : END DO
272 76 : IF (SUM(econf) /= NINT(zval)) THEN
273 0 : CPWARN("Valence charge and electron configuration not consistent")
274 : END IF
275 : !
276 30 : linc = 1
277 30 : IF (z > 18) linc = 2
278 30 : pend = 0.0_dp
279 30 : IF (kernel_basis) THEN
280 12 : bv(1:4) = (/3.20_dp, 2.80_dp, 2.40_dp, 2.00_dp/)
281 12 : bx(1:4) = (/4.00_dp, 3.50_dp, 3.00_dp, 2.50_dp/)
282 : !
283 12 : SELECT CASE (basis_cntrl)
284 : CASE (0)
285 0 : laux = lval + 1
286 : CASE (1)
287 12 : laux = MAX(lval + 1, lmax)
288 : CASE (2)
289 0 : laux = MAX(lval + 2, lmax + 1)
290 : CASE (3)
291 0 : laux = MAX(lval + 3, lmax + 2)
292 0 : laux = MIN(laux, 2 + linc)
293 : CASE DEFAULT
294 12 : CPABORT("Invalid value of control variable")
295 : END SELECT
296 : ELSE
297 18 : bv(1:4) = (/2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp/)
298 18 : bx(1:4) = (/2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp/)
299 : !
300 18 : SELECT CASE (basis_cntrl)
301 : CASE (0)
302 0 : laux = MAX(2*lval, lmax + linc)
303 0 : laux = MIN(laux, 2 + linc)
304 : CASE (1)
305 18 : laux = MAX(2*lval, lmax + linc)
306 18 : laux = MIN(laux, 3 + linc)
307 : CASE (2)
308 0 : laux = MAX(2*lval, lmax + linc + 1)
309 0 : laux = MIN(laux, 4 + linc)
310 : CASE (3)
311 0 : laux = MAX(2*lval, lmax + linc + 1)
312 0 : laux = MIN(laux, 4 + linc)
313 : CASE DEFAULT
314 18 : CPABORT("Invalid value of control variable")
315 : END SELECT
316 : END IF
317 : !
318 30 : DO l = 2*lmax + 1, laux
319 0 : pmin(l) = pmin(2*lmax)
320 0 : pmax(l) = pmax(2*lmax)
321 30 : peff(l) = peff(2*lmax)
322 : END DO
323 : !
324 30 : nval = 0
325 30 : IF (exact_1c_terms) THEN
326 0 : DO l = 0, laux
327 0 : IF (l <= lval + 1) THEN
328 0 : pend(l) = zmax(l) + 1.0_dp
329 0 : bval(l) = bv(basis_cntrl + 1)
330 : ELSE
331 0 : pend(l) = 2.0_dp*peff(l)
332 0 : bval(l) = bx(basis_cntrl + 1)
333 : END IF
334 0 : pmin(l) = zmin(l)
335 0 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
336 0 : nval(l) = MAX(CEILING(xv), 0)
337 0 : bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
338 : END DO
339 : ELSE
340 124 : DO l = 0, laux
341 94 : IF (l <= lval + 1) THEN
342 76 : pend(l) = pmax(l)
343 76 : bval(l) = bv(basis_cntrl + 1)
344 76 : pmin(l) = zmin(l)
345 : ELSE
346 18 : pend(l) = 4.0_dp*peff(l)
347 18 : bval(l) = bx(basis_cntrl + 1)
348 : END IF
349 94 : xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp
350 94 : nval(l) = MAX(CEILING(xv), 0)
351 124 : bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l))
352 : END DO
353 : END IF
354 : !
355 30 : lm = MIN(2*lval, 3)
356 92 : n1 = MAXVAL(nval(0:lm))
357 30 : IF (laux < lm + 1) THEN
358 : n2 = 0
359 : ELSE
360 56 : n2 = MAXVAL(nval(lm + 1:laux))
361 : END IF
362 : !
363 30 : nsets = n1 + n2
364 90 : ALLOCATE (zet(1, nsets))
365 830 : zet = 0.0_dp
366 30 : nl = 0
367 152 : j = MAXVAL(MAXLOC(nval(0:lm)))
368 280 : DO i = 1, n1
369 250 : ls1(i) = 0
370 250 : ls2(i) = lm
371 250 : npgf(i) = 1
372 250 : zet(1, i) = pmin(j)*bval(j)**(i - 1)
373 786 : DO l = 0, lm
374 756 : nl(l, i) = 1
375 : END DO
376 : END DO
377 30 : j = lm + 1
378 180 : DO i = n1 + 1, nsets
379 150 : ls1(i) = lm + 1
380 150 : ls2(i) = laux
381 150 : npgf(i) = 1
382 150 : zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1)
383 418 : DO l = lm + 1, laux
384 388 : nl(l, i) = 1
385 : END DO
386 : END DO
387 : !
388 30 : bsname = TRIM(element_symbol)//"-LRI-AUX-"//TRIM(orb_basis_set%name)
389 : !
390 30 : CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet)
391 : !
392 90 : DEALLOCATE (zet)
393 : END IF
394 :
395 60 : 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 12 : 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 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lmin, lset, npgf
411 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl
412 12 : INTEGER, DIMENSION(:), POINTER :: npgf_orb
413 : REAL(KIND=dp) :: cval, x, z0, z1
414 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zet
415 : REAL(KIND=dp), DIMENSION(0:9) :: zeff, zmax, zmin
416 :
417 12 : CALL get_basis_keyfigures(orb_basis, lmax, zmin, zmax, zeff)
418 12 : IF (nbas_oce < 1) THEN
419 12 : CALL get_gto_basis_set(gto_basis_set=orb_basis, nset=nset, npgf=npgf_orb)
420 54 : nx = SUM(npgf_orb(1:nset))
421 : ELSE
422 : nx = 0
423 : END IF
424 12 : nset = MAX(nbas_oce, nx)
425 12 : lx = MAX(lmax_oce, lmax)
426 : !
427 12 : bsname = "OCE-"//TRIM(orb_basis%name)
428 108 : ALLOCATE (lmin(nset), lset(nset), nl(0:9, nset), npgf(nset), zet(1, nset))
429 114 : lmin = 0
430 114 : lset = 0
431 1134 : nl = 1
432 114 : npgf = 1
433 216 : zet = 0.0_dp
434 : !
435 56 : z0 = MINVAL(zmin(0:lmax))
436 56 : z1 = MAXVAL(zmax(0:lmax))
437 12 : x = 1.0_dp/REAL(nset - 1, KIND=dp)
438 12 : cval = (z1/z0)**x
439 12 : zet(1, nset) = z0
440 102 : DO i = nset - 1, 1, -1
441 102 : zet(1, i) = zet(1, i + 1)*cval
442 : END DO
443 114 : DO i = 1, nset
444 102 : x = zet(1, i)
445 284 : DO l = 1, lmax
446 182 : z1 = 1.05_dp*zmax(l)
447 284 : IF (x < z1) lset(i) = l
448 : END DO
449 114 : IF (lset(i) == lmax) lset(i) = lx
450 : END DO
451 : !
452 12 : CALL create_aux_basis(oce_basis, bsname, nset, lmin, lset, nl, npgf, zet)
453 : !
454 12 : DEALLOCATE (lmin, lset, nl, npgf, zet)
455 :
456 12 : 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 308 : 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 308 : INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
472 308 : INTEGER, DIMENSION(:, :), POINTER :: lshell
473 : REAL(KIND=dp) :: aeff, gcca, gccb, kval, rexp, rint, rno, &
474 : zeta
475 308 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
476 308 : 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 308 : gcc=gcc)
486 :
487 1310 : lmax = MAXVAL(lm)
488 308 : CPASSERT(lmax <= 9)
489 :
490 308 : zmax = 0.0_dp
491 3388 : zmin = HUGE(0.0_dp)
492 308 : zeff = 0.0_dp
493 :
494 1310 : DO iset = 1, nset
495 : ! zmin zmax
496 3266 : DO ipgf = 1, npgf(iset)
497 7036 : DO ishell = 1, nshell(iset)
498 3770 : l = lshell(ishell, iset)
499 3770 : zeta = zet(ipgf, iset)
500 3770 : zmax(l) = MAX(zmax(l), zeta)
501 6034 : zmin(l) = MIN(zmin(l), zeta)
502 : END DO
503 : END DO
504 : ! zeff
505 2692 : DO ishell = 1, nshell(iset)
506 1382 : l = lshell(ishell, iset)
507 1382 : kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2)
508 1382 : rexp = 0.0_dp
509 1382 : rno = 0.0_dp
510 5152 : DO i = 1, npgf(iset)
511 3770 : gcca = gcc(i, ishell, iset)
512 20694 : DO j = 1, npgf(iset)
513 15542 : zeta = zet(i, iset) + zet(j, iset)
514 15542 : gccb = gcc(j, ishell, iset)
515 15542 : rint = 0.5_dp*fac(l + 1)/zeta**(l + 2)
516 15542 : rexp = rexp + gcca*gccb*rint
517 15542 : rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp)
518 19312 : rno = rno + gcca*gccb*rint
519 : END DO
520 : END DO
521 1382 : rexp = rexp/rno
522 1382 : aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2)
523 2384 : zeff(l) = MAX(zeff(l), aeff)
524 : END DO
525 : END DO
526 :
527 308 : 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 296 : 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 5920 : pmin = HUGE(0.0_dp)
547 296 : pmax = 0.0_dp
548 296 : peff = 0.0_dp
549 :
550 982 : DO l1 = 0, lmax
551 2196 : DO l2 = l1, lmax
552 4454 : DO la = l2 - l1, l2 + l1
553 2554 : pmax(la) = MAX(pmax(la), zmax(l1) + zmax(l2))
554 2554 : pmin(la) = MIN(pmin(la), zmin(l1) + zmin(l2))
555 3768 : peff(la) = MAX(peff(la), zeff(l1) + zeff(l2))
556 : END DO
557 : END DO
558 : END DO
559 :
560 296 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fx, tx, x2, xx
585 :
586 : ! SUM_i(fi M fi)
587 0 : fij = 0.0_dp
588 0 : DO ia = 1, npgf
589 0 : DO ib = 1, npgf
590 0 : p = zet(ia) + zet(ib) + amet
591 0 : intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
592 0 : DO i = 1, nfun
593 0 : 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 0 : ALLOCATE (fx(nfit, nfun), tx(nfit, nfun))
600 0 : fx = 0.0_dp
601 0 : DO ia = 1, npgf
602 0 : DO ib = 1, nfit
603 0 : p = zet(ia) + afit(ib) + amet
604 0 : intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1)
605 0 : DO i = 1, nfun
606 0 : 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 0 : ALLOCATE (xx(nfit, nfit), x2(nfit, nfit))
613 0 : DO ia = 1, nfit
614 0 : DO ib = 1, nfit
615 0 : p = afit(ia) + afit(ib) + amet
616 0 : 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 0 : tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun)
622 0 : x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit)
623 0 : CALL DPOSV("U", nfit, nfun, x2, nfit, tx, nfit, info)
624 0 : IF (info == 0) THEN
625 : ! value t*xx*t
626 : xij = 0.0_dp
627 0 : DO i = 1, nfun
628 0 : xij = xij + DOT_PRODUCT(tx(:, i), MATMUL(xx, tx(:, i)))
629 : END DO
630 : ! value t*fx
631 : fxij = 0.0_dp
632 0 : DO i = 1, nfun
633 0 : fxij = fxij + DOT_PRODUCT(tx(:, i), fx(:, i))
634 : END DO
635 : !
636 0 : eval = fij - 2.0_dp*fxij + xij
637 : ELSE
638 : ! error in solving for max overlap
639 0 : eval = 1.0e10_dp
640 : END IF
641 :
642 0 : DEALLOCATE (fx, xx, x2, tx)
643 :
644 0 : END SUBROUTINE overlap_maximum
645 : ! **************************************************************************************************
646 : !> \brief ...
647 : !> \param x ...
648 : !> \param n ...
649 : !> \param eval ...
650 : ! **************************************************************************************************
651 0 : 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 0 : DO i = 2, n
659 0 : IF (x(i) < 1.5_dp) THEN
660 0 : eval = eval + 10.0_dp*(1.5_dp - x(i))**2
661 : END IF
662 : END DO
663 :
664 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell
683 0 : INTEGER, DIMENSION(:, :), POINTER :: lshell
684 : LOGICAL :: toadd
685 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
686 0 : 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 0 : gcc=gcc)
696 :
697 0 : np = 0
698 0 : nf = 0
699 0 : DO iset = 1, nset
700 0 : toadd = .TRUE.
701 0 : DO ishell = 1, nshell(iset)
702 0 : l = lshell(ishell, iset)
703 0 : IF (l == lin) THEN
704 0 : nf = nf + 1
705 0 : IF (toadd) THEN
706 0 : np = np + npgf(iset)
707 0 : toadd = .FALSE.
708 : END IF
709 : END IF
710 : END DO
711 : END DO
712 0 : ALLOCATE (zval(np), gcval(np, nf))
713 0 : zval = 0.0_dp
714 0 : gcval = 0.0_dp
715 : !
716 : jp = 0
717 : jf = 0
718 0 : DO iset = 1, nset
719 0 : toadd = .TRUE.
720 0 : DO ishell = 1, nshell(iset)
721 0 : l = lshell(ishell, iset)
722 0 : IF (l == lin) THEN
723 0 : jf = jf + 1
724 0 : IF (toadd) THEN
725 0 : j1 = jp + 1
726 0 : j2 = jp + npgf(iset)
727 0 : zval(j1:j2) = zet(1:npgf(iset), iset)
728 0 : jp = jp + npgf(iset)
729 0 : toadd = .FALSE.
730 : END IF
731 0 : gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset)
732 : END IF
733 : END DO
734 : END DO
735 :
736 0 : END SUBROUTINE get_basis_functions
737 :
738 0 : END MODULE auto_basis
|
