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 : ! Calculates atomic integrals over unnormalized spherical Gaussian functions
10 : !-----------------------------------------------------------------------------!
11 : !
12 : ! phi(r) = r^l * exp[-p*r^2] Ylm
13 : !
14 : !-----------------------------------------------------------------------------!
15 : ! Calculates atomic integrals over normalized Slater type functions
16 : !-----------------------------------------------------------------------------!
17 : !
18 : ! phi(r) = N(nlm) r^(n-1) * exp[-p*r] Ylm
19 : ! N(nlm) = [(2n)!]^(-1/2) (2p)^(n+1/2)
20 : !
21 : !-----------------------------------------------------------------------------!
22 : ! Calculates atomic integrals over spherical numerical functions
23 : !-----------------------------------------------------------------------------!
24 : !
25 : ! phi(r) = R(r) Ylm
26 : !
27 : !-----------------------------------------------------------------------------!
28 : MODULE ai_onecenter
29 :
30 : USE kinds, ONLY: dp
31 : USE mathconstants, ONLY: dfac,&
32 : fac,&
33 : gamma0,&
34 : gamma1,&
35 : pi,&
36 : rootpi
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_onecenter'
44 :
45 : PUBLIC :: sg_overlap, sg_kinetic, sg_nuclear, sg_erf, sg_gpot, &
46 : sg_proj_ol, sg_coulomb, sg_exchange, sg_kinnuc, &
47 : sg_erfc, sg_conf
48 : PUBLIC :: sto_overlap, sto_kinetic, sto_nuclear
49 :
50 : CONTAINS
51 :
52 : !------------------------------------------------------------------------------
53 : !
54 : ! S(l,pq) = pi^(1/2)*(2*l+1)!! / 2^(l+2) / (p+q)^(l+1.5)
55 : !
56 : !------------------------------------------------------------------------------
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param smat ...
60 : !> \param l ...
61 : !> \param pa ...
62 : !> \param pb ...
63 : ! **************************************************************************************************
64 33324 : SUBROUTINE sg_overlap(smat, l, pa, pb)
65 :
66 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: smat
67 : INTEGER, INTENT(IN) :: l
68 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
69 :
70 : INTEGER :: ip, iq, m, n
71 : REAL(KIND=dp) :: el, spi
72 :
73 33324 : n = SIZE(pa)
74 33324 : m = SIZE(pb)
75 :
76 33324 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
77 :
78 33324 : spi = rootpi/2.0_dp**(l + 2)*dfac(2*l + 1)
79 33324 : el = REAL(l, dp) + 1.5_dp
80 :
81 195466 : DO iq = 1, m
82 2442604 : DO ip = 1, n
83 2409280 : smat(ip, iq) = spi/(pa(ip) + pb(iq))**el
84 : END DO
85 : END DO
86 :
87 33324 : END SUBROUTINE sg_overlap
88 :
89 : !------------------------------------------------------------------------------
90 : !
91 : ! T(l,pq) = (2l+3)!! pi^(1/2)/2^(l+2) [pq/(p+q)^(l+2.5)]
92 : !
93 : !------------------------------------------------------------------------------
94 : ! **************************************************************************************************
95 : !> \brief ...
96 : !> \param kmat ...
97 : !> \param l ...
98 : !> \param pa ...
99 : !> \param pb ...
100 : ! **************************************************************************************************
101 33077 : SUBROUTINE sg_kinetic(kmat, l, pa, pb)
102 :
103 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
104 : INTEGER, INTENT(IN) :: l
105 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
106 :
107 : INTEGER :: ip, iq, m, n
108 : REAL(KIND=dp) :: spi
109 :
110 33077 : n = SIZE(pa)
111 33077 : m = SIZE(pb)
112 :
113 33077 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
114 :
115 33077 : spi = dfac(2*l + 3)*rootpi/2.0_dp**(l + 2)
116 193960 : DO iq = 1, m
117 2430611 : DO ip = 1, n
118 2397534 : kmat(ip, iq) = spi*pa(ip)*pb(iq)/(pa(ip) + pb(iq))**(l + 2.5_dp)
119 : END DO
120 : END DO
121 :
122 33077 : END SUBROUTINE sg_kinetic
123 :
124 : !------------------------------------------------------------------------------
125 : !
126 : ! U(l,pq) = l!/2 / (p+q)^(l+1)
127 : !
128 : !------------------------------------------------------------------------------
129 : ! **************************************************************************************************
130 : !> \brief ...
131 : !> \param umat ...
132 : !> \param l ...
133 : !> \param pa ...
134 : !> \param pb ...
135 : ! **************************************************************************************************
136 11420 : SUBROUTINE sg_nuclear(umat, l, pa, pb)
137 :
138 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
139 : INTEGER, INTENT(IN) :: l
140 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
141 :
142 : INTEGER :: ip, iq, m, n
143 : REAL(KIND=dp) :: tld
144 :
145 11420 : n = SIZE(pa)
146 11420 : m = SIZE(pb)
147 :
148 11420 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
149 :
150 11420 : tld = 0.5_dp*fac(l)
151 69120 : DO iq = 1, m
152 1650876 : DO ip = 1, n
153 1639456 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1)
154 : END DO
155 : END DO
156 :
157 11420 : END SUBROUTINE sg_nuclear
158 :
159 : !------------------------------------------------------------------------------
160 : !
161 : ! U(l,pq) = l!/2 / (p+q)^l * [4/(p+q)^2 *pq*(l+1) + 1]
162 : !
163 : !------------------------------------------------------------------------------
164 : ! **************************************************************************************************
165 : !> \brief ...
166 : !> \param umat ...
167 : !> \param l ...
168 : !> \param pa ...
169 : !> \param pb ...
170 : ! **************************************************************************************************
171 278 : SUBROUTINE sg_kinnuc(umat, l, pa, pb)
172 :
173 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
174 : INTEGER, INTENT(IN) :: l
175 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
176 :
177 : INTEGER :: ip, iq, m, n
178 : REAL(KIND=dp) :: ppq, pq, tld
179 :
180 278 : n = SIZE(pa)
181 278 : m = SIZE(pb)
182 :
183 278 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
184 :
185 278 : IF (l > 0) THEN
186 204 : tld = 0.5_dp*fac(l)
187 5136 : DO iq = 1, m
188 168052 : DO ip = 1, n
189 162916 : ppq = pa(ip) + pb(iq)
190 162916 : pq = pa(ip)*pb(iq)
191 167848 : umat(ip, iq) = tld/ppq**l*(4.0_dp/ppq**2*pq*REAL(l + 1, dp) + 1.0_dp)
192 : END DO
193 : END DO
194 : ELSE
195 1862 : DO iq = 1, m
196 56070 : DO ip = 1, n
197 54208 : ppq = pa(ip) + pb(iq)
198 54208 : pq = pa(ip)*pb(iq)
199 55996 : umat(ip, iq) = 2.0_dp*pq/ppq**2
200 : END DO
201 : END DO
202 : END IF
203 :
204 278 : END SUBROUTINE sg_kinnuc
205 :
206 : !------------------------------------------------------------------------------
207 : !
208 : ! z = a/(p+q)
209 : !
210 : ! UP(l,pq,a) = Gamma(l+3/2)*a/SQRT(Pi)/(p+q)^(l+3/2)*
211 : ! Hypergeom([1/2, 3/2 + l], [3/2], -z)
212 : !
213 : ! UP(l,pq,a) = a/2^(l+1)/(p+q)^(l+3/2)/(1+z)^(l+1/2) * F(z,l)
214 : !
215 : ! F(z,0) = 1
216 : ! F(z,1) = 3 + 2*z
217 : ! F(z,2) = 15 + 20*z + 8*z^2
218 : ! F(z,3) = 35 + 70*z + 56*z^2 + 16*z^3
219 : ! F(z,4) = 315 + 840*z + 1008*z^2 + 576*z^3 + 128*z^4
220 : !
221 : !------------------------------------------------------------------------------
222 : ! **************************************************************************************************
223 : !> \brief ...
224 : !> \param upmat ...
225 : !> \param l ...
226 : !> \param a ...
227 : !> \param pa ...
228 : !> \param pb ...
229 : ! **************************************************************************************************
230 28881 : SUBROUTINE sg_erf(upmat, l, a, pa, pb)
231 :
232 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: upmat
233 : INTEGER, INTENT(IN) :: l
234 : REAL(KIND=dp), INTENT(IN) :: a
235 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
236 :
237 : INTEGER :: handle, ip, iq, m, n
238 : REAL(KIND=dp) :: a2, fpol, pq, tld, z, z2
239 :
240 28881 : CALL timeset("sg_erf", handle)
241 :
242 28881 : n = SIZE(pa)
243 28881 : m = SIZE(pb)
244 :
245 28881 : CPASSERT(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
246 :
247 28881 : a2 = a*a
248 28881 : tld = a/2._dp**(l + 1)
249 158797 : DO iq = 1, m
250 1571529 : DO ip = 1, n
251 1412732 : pq = pa(ip) + pb(iq)
252 1412732 : z = a2/pq
253 1542648 : upmat(ip, iq) = tld/(1._dp + z)**(l + 0.5_dp)/pq**(l + 1.5_dp)
254 : END DO
255 : END DO
256 :
257 158797 : DO iq = 1, m
258 28881 : SELECT CASE (l)
259 : CASE DEFAULT
260 0 : CPABORT("Only 0, 1, 2, 3, 4, 5, 6 are supported as the value of l in sg_erf")
261 : CASE (0)
262 : ! nothing left to do
263 : CASE (1)
264 364034 : DO ip = 1, n
265 324065 : pq = pa(ip) + pb(iq)
266 324065 : z = a2/pq
267 324065 : fpol = 2.0_dp*z + 3.0_dp
268 364034 : upmat(ip, iq) = upmat(ip, iq)*fpol
269 : END DO
270 : CASE (2)
271 232230 : DO ip = 1, n
272 213605 : pq = pa(ip) + pb(iq)
273 213605 : z = a2/pq
274 213605 : fpol = 8.0_dp*z*z + 20.0_dp*z + 15.0_dp
275 232230 : upmat(ip, iq) = upmat(ip, iq)*fpol
276 : END DO
277 : CASE (3)
278 159016 : DO ip = 1, n
279 153106 : pq = pa(ip) + pb(iq)
280 153106 : z = a2/pq
281 153106 : fpol = 16.0_dp*z*z*z + 56.0_dp*z*z + 70.0_dp*z + 35.0_dp
282 153106 : fpol = 3._dp*fpol
283 159016 : upmat(ip, iq) = upmat(ip, iq)*fpol
284 : END DO
285 : CASE (4)
286 156316 : DO ip = 1, n
287 150788 : pq = pa(ip) + pb(iq)
288 150788 : z = a2/pq
289 150788 : fpol = 128.0_dp*z*z*z*z + 576.0_dp*z*z*z + 1008.0_dp*z*z + 840.0_dp*z + 315.0_dp
290 150788 : fpol = 3._dp*fpol
291 156316 : upmat(ip, iq) = upmat(ip, iq)*fpol
292 : END DO
293 : CASE (5)
294 156200 : DO ip = 1, n
295 150688 : pq = pa(ip) + pb(iq)
296 150688 : z = a2/pq
297 150688 : fpol = 256.0_dp*z*z*z*z*z + 1408.0_dp*z*z*z*z + 3168.0_dp*z*z*z + 3696.0_dp*z*z + 2310.0_dp*z + 693.0_dp
298 150688 : fpol = 15._dp*fpol
299 156200 : upmat(ip, iq) = upmat(ip, iq)*fpol
300 : END DO
301 : CASE (6)
302 129916 : DO ip = 1, n
303 0 : pq = pa(ip) + pb(iq)
304 0 : z = a2/pq
305 0 : z2 = z*z
306 : fpol = 1024.0_dp*z2*z2*z2 + 6656.0_dp*z*z2*z2 + 18304.0_dp*z2*z2 + 27456.0_dp*z2*z + &
307 0 : 24024.0_dp*z2 + 12012.0_dp*z + 3003.0_dp
308 0 : fpol = 45._dp*fpol
309 0 : upmat(ip, iq) = upmat(ip, iq)*fpol
310 : END DO
311 : END SELECT
312 : END DO
313 :
314 28881 : CALL timestop(handle)
315 :
316 28881 : END SUBROUTINE sg_erf
317 :
318 : !------------------------------------------------------------------------------
319 : !
320 : ! Overlap with Projectors P(l,k,rc) for k=0,1,..
321 : !
322 : ! P(l,k,rc) = SQRT(2)/SQRT(Gamma[l+2k+1.5])/rc^(l+2k+1.5) r^(l+2k) exp[-0.5(r/rc)^2]
323 : !
324 : ! SP(l,k,p,rc) = 2^(l+k+1) / SQRT(gamma[l+2k+1.5]) / rc^(l+2k+1.5)
325 : ! * Gamma(l+k+1.5) / (2p+1/rc^2)^(l+k+1.5)
326 : !
327 : !------------------------------------------------------------------------------
328 : ! **************************************************************************************************
329 : !> \brief ...
330 : !> \param spmat ...
331 : !> \param l ...
332 : !> \param p ...
333 : !> \param k ...
334 : !> \param rc ...
335 : ! **************************************************************************************************
336 12279 : SUBROUTINE sg_proj_ol(spmat, l, p, k, rc)
337 :
338 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: spmat
339 : INTEGER, INTENT(IN) :: l
340 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p
341 : INTEGER, INTENT(IN) :: k
342 : REAL(KIND=dp), INTENT(IN) :: rc
343 :
344 : REAL(KIND=dp) :: orc, pf
345 :
346 12279 : CPASSERT(SIZE(spmat) >= SIZE(p))
347 :
348 12279 : pf = 2._dp**(l + k + 1)*gamma1(l + k + 1)/rc**(l + 2*k + 1.5_dp)/SQRT(gamma1(l + 2*k + 1))
349 12279 : orc = 1._dp/(rc*rc)
350 83659 : spmat(:) = pf/(2._dp*p(:) + orc)**(l + k + 1.5_dp)
351 :
352 12279 : END SUBROUTINE sg_proj_ol
353 :
354 : !------------------------------------------------------------------------------
355 : !
356 : ! Matrix elements for Gaussian potentials
357 : !
358 : ! V(k,rc) = (r/rc)^2k exp[-1/2(r/rc)^2]
359 : !
360 : ! VP(l,k,p+q,rc) = 2^(l+k+0.5) * rc^(2l+3) * Gamma(l+k+1.5) / (1+2rc^2(p+q))^(l+k+1.5)
361 : !
362 : !------------------------------------------------------------------------------
363 : ! **************************************************************************************************
364 : !> \brief ...
365 : !> \param vpmat ...
366 : !> \param k ...
367 : !> \param rc ...
368 : !> \param l ...
369 : !> \param pa ...
370 : !> \param pb ...
371 : ! **************************************************************************************************
372 50836 : SUBROUTINE sg_gpot(vpmat, k, rc, l, pa, pb)
373 :
374 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: vpmat
375 : INTEGER, INTENT(IN) :: k
376 : REAL(KIND=dp), INTENT(IN) :: rc
377 : INTEGER, INTENT(IN) :: l
378 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
379 :
380 : INTEGER :: ip, iq, m, n
381 : REAL(KIND=dp) :: tld
382 :
383 50836 : n = SIZE(pa)
384 50836 : m = SIZE(pb)
385 :
386 50836 : CPASSERT(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
387 :
388 50836 : tld = gamma1(l + k + 1)*rc**(2*l + 3)*2._dp**(l + k + 0.5)
389 :
390 289734 : DO iq = 1, m
391 3009000 : DO ip = 1, n
392 2958164 : vpmat(ip, iq) = tld/(1._dp + 2._dp*rc*rc*(pa(ip) + pb(iq)))**(l + k + 1.5_dp)
393 : END DO
394 : END DO
395 :
396 50836 : END SUBROUTINE sg_gpot
397 :
398 : !------------------------------------------------------------------------------
399 : !
400 : ! G(l,k,pq) = <a|[r/rc]^2k|b>
401 : ! = 0.5*Gamma(l+k+1.5)/rc^(2k)/(p+q)^(l+k+1.5)
402 : !
403 : !------------------------------------------------------------------------------
404 : ! **************************************************************************************************
405 : !> \brief ...
406 : !> \param gmat ...
407 : !> \param rc ...
408 : !> \param k ...
409 : !> \param l ...
410 : !> \param pa ...
411 : !> \param pb ...
412 : ! **************************************************************************************************
413 84 : SUBROUTINE sg_conf(gmat, rc, k, l, pa, pb)
414 :
415 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
416 : REAL(KIND=dp), INTENT(IN) :: rc
417 : INTEGER, INTENT(IN) :: k, l
418 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
419 :
420 : INTEGER :: ip, iq, m, n
421 : REAL(KIND=dp) :: tld
422 :
423 84 : n = SIZE(pa)
424 84 : m = SIZE(pb)
425 :
426 84 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
427 :
428 84 : tld = 0.5_dp/rc**(2*k)*gamma1(l + k + 1)
429 420 : DO iq = 1, m
430 1764 : DO ip = 1, n
431 1680 : gmat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + k + 1.5_dp)
432 : END DO
433 : END DO
434 :
435 84 : END SUBROUTINE sg_conf
436 :
437 : !------------------------------------------------------------------------------
438 : !
439 : ! (plql,rl'sl')
440 : !
441 : !------------------------------------------------------------------------------
442 : ! **************************************************************************************************
443 : !> \brief ...
444 : !> \param eri ...
445 : !> \param nu ...
446 : !> \param pa ...
447 : !> \param lab ...
448 : !> \param pc ...
449 : !> \param lcd ...
450 : ! **************************************************************************************************
451 1218 : SUBROUTINE sg_coulomb(eri, nu, pa, lab, pc, lcd)
452 :
453 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eri
454 : INTEGER, INTENT(IN) :: nu
455 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
456 : INTEGER, INTENT(IN) :: lab
457 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pc
458 : INTEGER, INTENT(IN) :: lcd
459 :
460 : INTEGER :: handle, ia, ib, ic, id, jab, jcd, na, nc
461 : REAL(KIND=dp) :: cc1, cc2, ee, p, q, r, s, vab1, vab3, &
462 : vcd1, vcd3, xab, xcd
463 :
464 1218 : CALL timeset("sg_coulomb", handle)
465 :
466 1218 : na = SIZE(pa)
467 1218 : nc = SIZE(pc)
468 1218 : ee = SQRT(2.0_dp*pi)/2.0_dp**(2*(lab + lcd) + 6)
469 1218 : jab = 0
470 12782 : DO ia = 1, na
471 11564 : p = pa(ia)
472 141744 : DO ib = ia, na
473 128962 : jab = jab + 1
474 128962 : q = pa(ib)
475 128962 : xab = 0.5_dp*(p + q)
476 128962 : vab1 = vgau(2*lab - nu + 1, xab)
477 128962 : vab3 = vgau(2*lab + nu + 2, xab)
478 128962 : jcd = 0
479 3272892 : DO ic = 1, nc
480 3132366 : r = pc(ic)
481 45753440 : DO id = ic, nc
482 42492112 : jcd = jcd + 1
483 42492112 : s = pc(id)
484 42492112 : xcd = 0.5_dp*(r + s)
485 42492112 : vcd1 = vgau(2*lcd + nu + 2, xcd)
486 42492112 : vcd3 = vgau(2*lcd - nu + 1, xcd)
487 42492112 : cc1 = cgau(2*lab - nu + 1, 2*lcd + nu + 2, xab/xcd)
488 42492112 : cc2 = cgau(2*lcd - nu + 1, 2*lab + nu + 2, xcd/xab)
489 :
490 45624478 : eri(jab, jcd) = ee*(cc1*vab1*vcd1 + cc2*vab3*vcd3)
491 :
492 : END DO
493 : END DO
494 : END DO
495 : END DO
496 :
497 1218 : CALL timestop(handle)
498 :
499 1218 : END SUBROUTINE sg_coulomb
500 :
501 : !------------------------------------------------------------------------------
502 : !
503 : ! (plql',rlsl')
504 : !
505 : !------------------------------------------------------------------------------
506 : ! **************************************************************************************************
507 : !> \brief ...
508 : !> \param eri ...
509 : !> \param nu ...
510 : !> \param pa ...
511 : !> \param lac ...
512 : !> \param pb ...
513 : !> \param lbd ...
514 : ! **************************************************************************************************
515 2128 : SUBROUTINE sg_exchange(eri, nu, pa, lac, pb, lbd)
516 :
517 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eri
518 : INTEGER, INTENT(IN) :: nu
519 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
520 : INTEGER, INTENT(IN) :: lac
521 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
522 : INTEGER, INTENT(IN) :: lbd
523 :
524 : INTEGER :: handle, ia, ib, ic, id, jac, jbd, na, nb
525 : REAL(KIND=dp) :: cc1, cc2, cc3, cc4, ee, p, q, r, s, &
526 : v1pr, v1ps, v1qr, v1qs, v2pr, v2ps, &
527 : v2qr, v2qs, xab, xad, xbc, xcd
528 :
529 2128 : CALL timeset("sg_exchange", handle)
530 :
531 40295595 : eri = 0.0_dp
532 2128 : na = SIZE(pa)
533 2128 : nb = SIZE(pb)
534 2128 : ee = SQRT(2.0_dp*pi)/2.0_dp**(2*(lac + lbd) + 7)
535 2128 : jac = 0
536 15842 : DO ia = 1, na
537 13714 : p = pa(ia)
538 162070 : DO ic = ia, na
539 146228 : jac = jac + 1
540 146228 : q = pa(ic)
541 146228 : jbd = 0
542 3424067 : DO ib = 1, nb
543 3264125 : r = pb(ib)
544 3264125 : xab = 0.5_dp*(p + r)
545 3264125 : xbc = 0.5_dp*(q + r)
546 43482700 : DO id = ib, nb
547 40072347 : jbd = jbd + 1
548 40072347 : s = pb(id)
549 40072347 : xcd = 0.5_dp*(q + s)
550 40072347 : xad = 0.5_dp*(p + s)
551 40072347 : v1pr = vgau(lac + lbd - nu + 1, xab)
552 40072347 : v1qs = vgau(lac + lbd - nu + 1, xcd)
553 40072347 : v1ps = vgau(lac + lbd - nu + 1, xad)
554 40072347 : v1qr = vgau(lac + lbd - nu + 1, xbc)
555 40072347 : v2qs = vgau(lac + lbd + nu + 2, xcd)
556 40072347 : v2pr = vgau(lac + lbd + nu + 2, xab)
557 40072347 : v2qr = vgau(lac + lbd + nu + 2, xbc)
558 40072347 : v2ps = vgau(lac + lbd + nu + 2, xad)
559 40072347 : cc1 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xab/xcd)
560 40072347 : cc2 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xcd/xab)
561 40072347 : cc3 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xad/xbc)
562 40072347 : cc4 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xbc/xad)
563 :
564 : eri(jac, jbd) = ee*(v1pr*v2qs*cc1 + v1qs*v2pr*cc2 + &
565 43336472 : v1ps*v2qr*cc3 + v1qr*v2ps*cc4)
566 :
567 : END DO
568 : END DO
569 : END DO
570 : END DO
571 :
572 2128 : CALL timestop(handle)
573 :
574 2128 : END SUBROUTINE sg_exchange
575 :
576 : !------------------------------------------------------------------------------
577 : !
578 : ! erfc(a*x)/x = 1/x - erf(a*x)/x
579 : !
580 : !------------------------------------------------------------------------------
581 : ! **************************************************************************************************
582 : !> \brief ...
583 : !> \param umat ...
584 : !> \param l ...
585 : !> \param a ...
586 : !> \param pa ...
587 : !> \param pb ...
588 : ! **************************************************************************************************
589 210 : SUBROUTINE sg_erfc(umat, l, a, pa, pb)
590 :
591 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
592 : INTEGER, INTENT(IN) :: l
593 : REAL(KIND=dp), INTENT(IN) :: a
594 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
595 :
596 : INTEGER :: ip, iq, m, n
597 : REAL(KIND=dp) :: tld
598 :
599 210 : n = SIZE(pa)
600 210 : m = SIZE(pb)
601 :
602 210 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
603 :
604 210 : CALL sg_erf(umat, l, a, pa, pb)
605 :
606 210 : tld = 0.5_dp*fac(l)
607 1210 : DO iq = 1, m
608 13042 : DO ip = 1, n
609 12832 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1) - umat(ip, iq)
610 : END DO
611 : END DO
612 :
613 210 : END SUBROUTINE sg_erfc
614 :
615 : ! ***************************************************************************************************
616 :
617 : ! **************************************************************************************************
618 : !> \brief ...
619 : !> \param n ...
620 : !> \param x ...
621 : !> \return ...
622 : ! **************************************************************************************************
623 405820924 : ELEMENTAL FUNCTION vgau(n, x) RESULT(v)
624 : INTEGER, INTENT(IN) :: n
625 : REAL(KIND=dp), INTENT(IN) :: x
626 : REAL(KIND=dp) :: v
627 :
628 405820924 : v = dfac(n - 1)/x**(0.5_dp*(n + 1))
629 :
630 405820924 : END FUNCTION vgau
631 :
632 : ! **************************************************************************************************
633 : !> \brief ...
634 : !> \param a ...
635 : !> \param b ...
636 : !> \param t ...
637 : !> \return ...
638 : ! **************************************************************************************************
639 245273612 : ELEMENTAL FUNCTION cgau(a, b, t) RESULT(c)
640 : INTEGER, INTENT(IN) :: a, b
641 : REAL(KIND=dp), INTENT(IN) :: t
642 : REAL(KIND=dp) :: c
643 :
644 : INTEGER :: l
645 :
646 245273612 : c = 0.0_dp
647 761069194 : DO l = 0, (a - 1)/2
648 761069194 : c = c + (t/(1.0_dp + t))**l*dfac(2*l + b - 1)/dfac(2*l)
649 : END DO
650 245273612 : c = c*(1.0_dp + t)**(-0.5_dp*(b + 1))/dfac(b - 1)
651 :
652 245273612 : END FUNCTION cgau
653 :
654 : !------------------------------------------------------------------------------
655 : !
656 : ! S(l,pn,qm) = ( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m,(p+q)/2]
657 : !
658 : !------------------------------------------------------------------------------
659 : ! **************************************************************************************************
660 : !> \brief ...
661 : !> \param smat ...
662 : !> \param na ...
663 : !> \param pa ...
664 : !> \param nb ...
665 : !> \param pb ...
666 : ! **************************************************************************************************
667 6720 : SUBROUTINE sto_overlap(smat, na, pa, nb, pb)
668 :
669 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: smat
670 : INTEGER, DIMENSION(:), INTENT(IN) :: na
671 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
672 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
673 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
674 :
675 : INTEGER :: ip, iq, m, n
676 : REAL(KIND=dp) :: vp, vpq, vq
677 :
678 6720 : n = SIZE(pa)
679 6720 : m = SIZE(pb)
680 :
681 6720 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
682 :
683 10686 : DO iq = 1, m
684 3966 : vq = vsto(2*nb(iq), pb(iq))
685 21712 : DO ip = 1, n
686 11026 : vp = vsto(2*na(ip), pa(ip))
687 11026 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
688 14992 : smat(ip, iq) = vpq/SQRT(vp*vq)
689 : END DO
690 : END DO
691 :
692 6720 : END SUBROUTINE sto_overlap
693 :
694 : !------------------------------------------------------------------------------
695 : !
696 : ! T(l,pn,qm) = 0.5*p*q*( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m,(p+q)/2]
697 : ! -(W[l,n,p]+W[l,m,q]) * V[n+m-1,(p+q)/2]
698 : ! +W[l,n,p]*W[l,m,q] * V[n+m-2,(p+q)/2]
699 : !
700 : !------------------------------------------------------------------------------
701 : ! **************************************************************************************************
702 : !> \brief ...
703 : !> \param kmat ...
704 : !> \param l ...
705 : !> \param na ...
706 : !> \param pa ...
707 : !> \param nb ...
708 : !> \param pb ...
709 : ! **************************************************************************************************
710 6720 : SUBROUTINE sto_kinetic(kmat, l, na, pa, nb, pb)
711 :
712 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
713 : INTEGER, INTENT(IN) :: l
714 : INTEGER, DIMENSION(:), INTENT(IN) :: na
715 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
716 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
717 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
718 :
719 : INTEGER :: ip, iq, m, n
720 : REAL(KIND=dp) :: vp, vpq, vpq1, vpq2, vq, wp, wq
721 :
722 6720 : n = SIZE(pa)
723 6720 : m = SIZE(pb)
724 :
725 6720 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
726 :
727 10686 : DO iq = 1, m
728 3966 : vq = vsto(2*nb(iq), pb(iq))
729 3966 : wq = wsto(l, nb(iq), pb(iq))
730 21712 : DO ip = 1, n
731 11026 : vp = vsto(2*na(ip), pa(ip))
732 11026 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
733 11026 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
734 11026 : vpq2 = vsto(na(ip) + nb(iq) - 2, 0.5_dp*(pa(ip) + pb(iq)))
735 11026 : wp = wsto(l, na(ip), pa(ip))
736 : kmat(ip, iq) = 0.5_dp*pa(ip)*pb(iq)/SQRT(vp*vq)* &
737 14992 : (vpq - (wp + wq)*vpq1 + wp*wq*vpq2)
738 : END DO
739 : END DO
740 :
741 6720 : END SUBROUTINE sto_kinetic
742 :
743 : !------------------------------------------------------------------------------
744 : !
745 : ! U(l,pq) = 2( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m-1,(p+q)/2]
746 : !
747 : !------------------------------------------------------------------------------
748 : ! **************************************************************************************************
749 : !> \brief ...
750 : !> \param umat ...
751 : !> \param na ...
752 : !> \param pa ...
753 : !> \param nb ...
754 : !> \param pb ...
755 : ! **************************************************************************************************
756 7668 : SUBROUTINE sto_nuclear(umat, na, pa, nb, pb)
757 :
758 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
759 : INTEGER, DIMENSION(:), INTENT(IN) :: na
760 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
761 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
762 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
763 :
764 : INTEGER :: ip, iq, m, n
765 : REAL(KIND=dp) :: vp, vpq1, vq
766 :
767 7668 : n = SIZE(pa)
768 7668 : m = SIZE(pb)
769 :
770 7668 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
771 :
772 11886 : DO iq = 1, m
773 4218 : vq = vsto(2*nb(iq), pb(iq))
774 23164 : DO ip = 1, n
775 11278 : vp = vsto(2*na(ip), pa(ip))
776 11278 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
777 15496 : umat(ip, iq) = 2._dp/SQRT(vp*vq)*vpq1
778 : END DO
779 : END DO
780 :
781 7668 : END SUBROUTINE sto_nuclear
782 :
783 : !------------------------------------------------------------------------------
784 : !
785 : ! G(l,k,pq) = <aln|[r/rc]^2k|blm>
786 : ! = N(aln)*N(blm) (a+b)^(-(n+m+2k+1))/rc^2k * GAMMA(n+m+2k+1)
787 : !
788 : !------------------------------------------------------------------------------
789 : ! **************************************************************************************************
790 : !> \brief ...
791 : !> \param gmat ...
792 : !> \param rc ...
793 : !> \param k ...
794 : !> \param na ...
795 : !> \param pa ...
796 : !> \param nb ...
797 : !> \param pb ...
798 : ! **************************************************************************************************
799 0 : SUBROUTINE sto_conf(gmat, rc, k, na, pa, nb, pb)
800 :
801 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
802 : REAL(KIND=dp), INTENT(IN) :: rc
803 : INTEGER, INTENT(IN) :: k
804 : INTEGER, DIMENSION(:), INTENT(IN) :: na
805 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
806 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
807 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
808 :
809 : INTEGER :: ip, iq, m, n
810 :
811 0 : n = SIZE(pa)
812 0 : m = SIZE(pb)
813 :
814 0 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
815 :
816 0 : DO iq = 1, m
817 0 : DO ip = 1, n
818 : gmat(ip, iq) = (2._dp*pa(ip))**(na(ip) + 0.5_dp)/SQRT(fac(2*na(ip))) &
819 : *(2._dp*pb(iq))**(nb(iq) + 0.5_dp)/SQRT(fac(2*nb(iq))) &
820 : /rc**(2*k)/(pa(ip) + pb(iq))**(na(ip) + nb(iq) + 2*k + 1) &
821 0 : *gamma0(na(ip) + nb(iq) + 2*k + 1)
822 : END DO
823 : END DO
824 :
825 0 : END SUBROUTINE sto_conf
826 :
827 : ! **************************************************************************************************
828 :
829 : ! **************************************************************************************************
830 : !> \brief ...
831 : !> \param n ...
832 : !> \param x ...
833 : !> \return ...
834 : ! **************************************************************************************************
835 100862 : ELEMENTAL FUNCTION vsto(n, x) RESULT(v)
836 : INTEGER, INTENT(IN) :: n
837 : REAL(KIND=dp), INTENT(IN) :: x
838 : REAL(KIND=dp) :: v
839 :
840 100862 : v = fac(n)/x**(n + 1)
841 :
842 100862 : END FUNCTION vsto
843 :
844 : ! **************************************************************************************************
845 : !> \brief ...
846 : !> \param n ...
847 : !> \param m ...
848 : !> \param x ...
849 : !> \return ...
850 : ! **************************************************************************************************
851 14992 : ELEMENTAL FUNCTION wsto(n, m, x) RESULT(w)
852 : INTEGER, INTENT(IN) :: n, m
853 : REAL(KIND=dp), INTENT(IN) :: x
854 : REAL(KIND=dp) :: w
855 :
856 14992 : w = 2._dp*REAL(m - n - 1, dp)/x
857 :
858 14992 : END FUNCTION wsto
859 : END MODULE ai_onecenter
|
