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