(git:374b731)
Loading...
Searching...
No Matches
ai_onecenter.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8!-----------------------------------------------------------------------------!
9! 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!-----------------------------------------------------------------------------!
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
49
50CONTAINS
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 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 n = SIZE(pa)
74 m = SIZE(pb)
75
76 cpassert(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
77
78 spi = rootpi/2.0_dp**(l + 2)*dfac(2*l + 1)
79 el = real(l, dp) + 1.5_dp
80
81 DO iq = 1, m
82 DO ip = 1, n
83 smat(ip, iq) = spi/(pa(ip) + pb(iq))**el
84 END DO
85 END DO
86
87 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 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 n = SIZE(pa)
111 m = SIZE(pb)
112
113 cpassert(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
114
115 spi = dfac(2*l + 3)*rootpi/2.0_dp**(l + 2)
116 DO iq = 1, m
117 DO ip = 1, n
118 kmat(ip, iq) = spi*pa(ip)*pb(iq)/(pa(ip) + pb(iq))**(l + 2.5_dp)
119 END DO
120 END DO
121
122 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 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 n = SIZE(pa)
146 m = SIZE(pb)
147
148 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
149
150 tld = 0.5_dp*fac(l)
151 DO iq = 1, m
152 DO ip = 1, n
153 umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1)
154 END DO
155 END DO
156
157 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 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 n = SIZE(pa)
181 m = SIZE(pb)
182
183 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
184
185 IF (l > 0) THEN
186 tld = 0.5_dp*fac(l)
187 DO iq = 1, m
188 DO ip = 1, n
189 ppq = pa(ip) + pb(iq)
190 pq = pa(ip)*pb(iq)
191 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 DO iq = 1, m
196 DO ip = 1, n
197 ppq = pa(ip) + pb(iq)
198 pq = pa(ip)*pb(iq)
199 umat(ip, iq) = 2.0_dp*pq/ppq**2
200 END DO
201 END DO
202 END IF
203
204 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 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 CALL timeset("sg_erf", handle)
241
242 n = SIZE(pa)
243 m = SIZE(pb)
244
245 cpassert(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
246
247 a2 = a*a
248 tld = a/2._dp**(l + 1)
249 DO iq = 1, m
250 DO ip = 1, n
251 pq = pa(ip) + pb(iq)
252 z = a2/pq
253 upmat(ip, iq) = tld/(1._dp + z)**(l + 0.5_dp)/pq**(l + 1.5_dp)
254 END DO
255 END DO
256
257 DO iq = 1, m
258 SELECT CASE (l)
259 CASE DEFAULT
260 cpabort("")
261 CASE (0)
262 ! nothing left to do
263 CASE (1)
264 DO ip = 1, n
265 pq = pa(ip) + pb(iq)
266 z = a2/pq
267 fpol = 2.0_dp*z + 3.0_dp
268 upmat(ip, iq) = upmat(ip, iq)*fpol
269 END DO
270 CASE (2)
271 DO ip = 1, n
272 pq = pa(ip) + pb(iq)
273 z = a2/pq
274 fpol = 8.0_dp*z*z + 20.0_dp*z + 15.0_dp
275 upmat(ip, iq) = upmat(ip, iq)*fpol
276 END DO
277 CASE (3)
278 DO ip = 1, n
279 pq = pa(ip) + pb(iq)
280 z = a2/pq
281 fpol = 16.0_dp*z*z*z + 56.0_dp*z*z + 70.0_dp*z + 35.0_dp
282 fpol = 3._dp*fpol
283 upmat(ip, iq) = upmat(ip, iq)*fpol
284 END DO
285 CASE (4)
286 DO ip = 1, n
287 pq = pa(ip) + pb(iq)
288 z = a2/pq
289 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 fpol = 3._dp*fpol
291 upmat(ip, iq) = upmat(ip, iq)*fpol
292 END DO
293 CASE (5)
294 DO ip = 1, n
295 pq = pa(ip) + pb(iq)
296 z = a2/pq
297 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 fpol = 15._dp*fpol
299 upmat(ip, iq) = upmat(ip, iq)*fpol
300 END DO
301 CASE (6)
302 DO ip = 1, n
303 pq = pa(ip) + pb(iq)
304 z = a2/pq
305 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 24024.0_dp*z2 + 12012.0_dp*z + 3003.0_dp
308 fpol = 45._dp*fpol
309 upmat(ip, iq) = upmat(ip, iq)*fpol
310 END DO
311 END SELECT
312 END DO
313
314 CALL timestop(handle)
315
316 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 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 cpassert(SIZE(spmat) >= SIZE(p))
347
348 pf = 2._dp**(l + k + 1)*gamma1(l + k + 1)/rc**(l + 2*k + 1.5_dp)/sqrt(gamma1(l + 2*k + 1))
349 orc = 1._dp/(rc*rc)
350 spmat(:) = pf/(2._dp*p(:) + orc)**(l + k + 1.5_dp)
351
352 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 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 n = SIZE(pa)
384 m = SIZE(pb)
385
386 cpassert(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
387
388 tld = gamma1(l + k + 1)*rc**(2*l + 3)*2._dp**(l + k + 0.5)
389
390 DO iq = 1, m
391 DO ip = 1, n
392 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 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 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 n = SIZE(pa)
424 m = SIZE(pb)
425
426 cpassert(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
427
428 tld = 0.5_dp/rc**(2*k)*gamma1(l + k + 1)
429 DO iq = 1, m
430 DO ip = 1, n
431 gmat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + k + 1.5_dp)
432 END DO
433 END DO
434
435 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 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 CALL timeset("sg_coulomb", handle)
465
466 na = SIZE(pa)
467 nc = SIZE(pc)
468 ee = sqrt(2.0_dp*pi)/2.0_dp**(2*(lab + lcd) + 6)
469 jab = 0
470 DO ia = 1, na
471 p = pa(ia)
472 DO ib = ia, na
473 jab = jab + 1
474 q = pa(ib)
475 xab = 0.5_dp*(p + q)
476 vab1 = vgau(2*lab - nu + 1, xab)
477 vab3 = vgau(2*lab + nu + 2, xab)
478 jcd = 0
479 DO ic = 1, nc
480 r = pc(ic)
481 DO id = ic, nc
482 jcd = jcd + 1
483 s = pc(id)
484 xcd = 0.5_dp*(r + s)
485 vcd1 = vgau(2*lcd + nu + 2, xcd)
486 vcd3 = vgau(2*lcd - nu + 1, xcd)
487 cc1 = cgau(2*lab - nu + 1, 2*lcd + nu + 2, xab/xcd)
488 cc2 = cgau(2*lcd - nu + 1, 2*lab + nu + 2, xcd/xab)
489
490 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 CALL timestop(handle)
498
499 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 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 CALL timeset("sg_exchange", handle)
530
531 eri = 0.0_dp
532 na = SIZE(pa)
533 nb = SIZE(pb)
534 ee = sqrt(2.0_dp*pi)/2.0_dp**(2*(lac + lbd) + 7)
535 jac = 0
536 DO ia = 1, na
537 p = pa(ia)
538 DO ic = ia, na
539 jac = jac + 1
540 q = pa(ic)
541 jbd = 0
542 DO ib = 1, nb
543 r = pb(ib)
544 xab = 0.5_dp*(p + r)
545 xbc = 0.5_dp*(q + r)
546 DO id = ib, nb
547 jbd = jbd + 1
548 s = pb(id)
549 xcd = 0.5_dp*(q + s)
550 xad = 0.5_dp*(p + s)
551 v1pr = vgau(lac + lbd - nu + 1, xab)
552 v1qs = vgau(lac + lbd - nu + 1, xcd)
553 v1ps = vgau(lac + lbd - nu + 1, xad)
554 v1qr = vgau(lac + lbd - nu + 1, xbc)
555 v2qs = vgau(lac + lbd + nu + 2, xcd)
556 v2pr = vgau(lac + lbd + nu + 2, xab)
557 v2qr = vgau(lac + lbd + nu + 2, xbc)
558 v2ps = vgau(lac + lbd + nu + 2, xad)
559 cc1 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xab/xcd)
560 cc2 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xcd/xab)
561 cc3 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xad/xbc)
562 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 v1ps*v2qr*cc3 + v1qr*v2ps*cc4)
566
567 END DO
568 END DO
569 END DO
570 END DO
571
572 CALL timestop(handle)
573
574 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 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 n = SIZE(pa)
600 m = SIZE(pb)
601
602 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
603
604 CALL sg_erf(umat, l, a, pa, pb)
605
606 tld = 0.5_dp*fac(l)
607 DO iq = 1, m
608 DO ip = 1, n
609 umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1) - umat(ip, iq)
610 END DO
611 END DO
612
613 END SUBROUTINE sg_erfc
614
615! ***************************************************************************************************
616
617! **************************************************************************************************
618!> \brief ...
619!> \param n ...
620!> \param x ...
621!> \return ...
622! **************************************************************************************************
623 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 v = dfac(n - 1)/x**(0.5_dp*(n + 1))
629
630 END FUNCTION vgau
631
632! **************************************************************************************************
633!> \brief ...
634!> \param a ...
635!> \param b ...
636!> \param t ...
637!> \return ...
638! **************************************************************************************************
639 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 c = 0.0_dp
647 DO l = 0, (a - 1)/2
648 c = c + (t/(1.0_dp + t))**l*dfac(2*l + b - 1)/dfac(2*l)
649 END DO
650 c = c*(1.0_dp + t)**(-0.5_dp*(b + 1))/dfac(b - 1)
651
652 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 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 n = SIZE(pa)
679 m = SIZE(pb)
680
681 cpassert(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
682
683 DO iq = 1, m
684 vq = vsto(2*nb(iq), pb(iq))
685 DO ip = 1, n
686 vp = vsto(2*na(ip), pa(ip))
687 vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
688 smat(ip, iq) = vpq/sqrt(vp*vq)
689 END DO
690 END DO
691
692 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 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 n = SIZE(pa)
723 m = SIZE(pb)
724
725 cpassert(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
726
727 DO iq = 1, m
728 vq = vsto(2*nb(iq), pb(iq))
729 wq = wsto(l, nb(iq), pb(iq))
730 DO ip = 1, n
731 vp = vsto(2*na(ip), pa(ip))
732 vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
733 vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
734 vpq2 = vsto(na(ip) + nb(iq) - 2, 0.5_dp*(pa(ip) + pb(iq)))
735 wp = wsto(l, na(ip), pa(ip))
736 kmat(ip, iq) = 0.5_dp*pa(ip)*pb(iq)/sqrt(vp*vq)* &
737 (vpq - (wp + wq)*vpq1 + wp*wq*vpq2)
738 END DO
739 END DO
740
741 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 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 n = SIZE(pa)
768 m = SIZE(pb)
769
770 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
771
772 DO iq = 1, m
773 vq = vsto(2*nb(iq), pb(iq))
774 DO ip = 1, n
775 vp = vsto(2*na(ip), pa(ip))
776 vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
777 umat(ip, iq) = 2._dp/sqrt(vp*vq)*vpq1
778 END DO
779 END DO
780
781 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 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 n = SIZE(pa)
812 m = SIZE(pb)
813
814 cpassert(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
815
816 DO iq = 1, m
817 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 *gamma0(na(ip) + nb(iq) + 2*k + 1)
822 END DO
823 END DO
824
825 END SUBROUTINE sto_conf
826
827! **************************************************************************************************
828
829! **************************************************************************************************
830!> \brief ...
831!> \param n ...
832!> \param x ...
833!> \return ...
834! **************************************************************************************************
835 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 v = fac(n)/x**(n + 1)
841
842 END FUNCTION vsto
843
844! **************************************************************************************************
845!> \brief ...
846!> \param n ...
847!> \param m ...
848!> \param x ...
849!> \return ...
850! **************************************************************************************************
851 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 w = 2._dp*real(m - n - 1, dp)/x
857
858 END FUNCTION wsto
859!------------------------------------------------------------------------------
860!
861! S(l,pq) = INT(u^2 Ra(u) Rb(u))du
862!
863!------------------------------------------------------------------------------
864! **************************************************************************************************
865!> \brief ...
866!> \param smat ...
867!> \param ra ...
868!> \param rb ...
869!> \param wr ...
870! **************************************************************************************************
871 SUBROUTINE num_overlap(smat, ra, rb, wr)
872
873 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: smat
874 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
875 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wr
876
877 INTEGER :: ip, iq, m, n
878
879 n = SIZE(ra, 2)
880 m = SIZE(rb, 2)
881
882 cpassert(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
883
884 DO iq = 1, m
885 DO ip = 1, n
886 smat(ip, iq) = sum(ra(:, ip)*rb(:, iq)*wr(:))
887 END DO
888 END DO
889
890 END SUBROUTINE num_overlap
891
892!------------------------------------------------------------------------------
893!
894! T(l,pq) = 0.5 INT( u^2 dRa(u) dRb(u) + l(l+1) Ra(u) Rb(u))du
895!
896!------------------------------------------------------------------------------
897! **************************************************************************************************
898!> \brief ...
899!> \param kmat ...
900!> \param l ...
901!> \param ra ...
902!> \param dra ...
903!> \param rb ...
904!> \param drb ...
905!> \param r ...
906!> \param wr ...
907! **************************************************************************************************
908 SUBROUTINE num_kinetic(kmat, l, ra, dra, rb, drb, r, wr)
909
910 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
911 INTEGER, INTENT(IN) :: l
912 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, dra, rb, drb
913 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
914
915 INTEGER :: ip, iq, m, n
916
917 n = SIZE(ra, 2)
918 m = SIZE(rb, 2)
919
920 cpassert(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
921
922 DO iq = 1, m
923 DO ip = 1, n
924 kmat(ip, iq) = 0.5_dp*sum(wr(:)*dra(:, ip)*drb(:, iq) &
925 + wr(:)*real(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**2)
926 END DO
927 END DO
928
929 END SUBROUTINE num_kinetic
930
931!------------------------------------------------------------------------------
932!
933! U(l,pq) = INT(u Ra(u) Rb(u))du
934!
935!------------------------------------------------------------------------------
936! **************************************************************************************************
937!> \brief ...
938!> \param umat ...
939!> \param ra ...
940!> \param rb ...
941!> \param r ...
942!> \param wr ...
943! **************************************************************************************************
944 SUBROUTINE num_nuclear(umat, ra, rb, r, wr)
945
946 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: umat
947 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
948 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
949
950 INTEGER :: ip, iq, m, n
951
952 n = SIZE(ra, 2)
953 m = SIZE(rb, 2)
954
955 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
956
957 DO iq = 1, m
958 DO ip = 1, n
959 umat(ip, iq) = sum(wr(:)*ra(:, ip)*rb(:, iq)/r(:))
960 END DO
961 END DO
962
963 END SUBROUTINE num_nuclear
964
965!------------------------------------------------------------------------------
966!
967! U(l,pq) = INT(u dRa(u) dRb(u))du + l(l+1) * INT(Ra(u) Rb(u) / u)du
968!
969!------------------------------------------------------------------------------
970! **************************************************************************************************
971!> \brief ...
972!> \param umat ...
973!> \param l ...
974!> \param ra ...
975!> \param dra ...
976!> \param rb ...
977!> \param drb ...
978!> \param r ...
979!> \param wr ...
980! **************************************************************************************************
981 SUBROUTINE num_kinnuc(umat, l, ra, dra, rb, drb, r, wr)
982
983 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: umat
984 INTEGER, INTENT(IN) :: l
985 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, dra, rb, drb
986 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
987
988 INTEGER :: ip, iq, m, n
989
990 n = SIZE(ra, 2)
991 m = SIZE(rb, 2)
992
993 cpassert(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
994
995 DO iq = 1, m
996 DO ip = 1, n
997 umat(ip, iq) = sum(wr(:)*dra(:, ip)*drb(:, iq)/r(:) &
998 + wr(:)*real(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**3)
999 END DO
1000 END DO
1001
1002 END SUBROUTINE num_kinnuc
1003
1004!------------------------------------------------------------------------------
1005!
1006! U(l,pq) = INT(u erf(a*u) Ra(u) Rb(u))du
1007!
1008!------------------------------------------------------------------------------
1009! **************************************************************************************************
1010!> \brief ...
1011!> \param upmat ...
1012!> \param a ...
1013!> \param ra ...
1014!> \param rb ...
1015!> \param r ...
1016!> \param wr ...
1017! **************************************************************************************************
1018 SUBROUTINE num_erf(upmat, a, ra, rb, r, wr)
1019
1020 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: upmat
1021 REAL(kind=dp), INTENT(IN) :: a
1022 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1023 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
1024
1025 INTEGER :: ip, iq, k, m, n
1026
1027 n = SIZE(ra, 2)
1028 m = SIZE(rb, 2)
1029
1030 cpassert(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
1031
1032 DO iq = 1, m
1033 DO ip = 1, n
1034 upmat(ip, iq) = 0._dp
1035 DO k = 1, SIZE(r)
1036 upmat(ip, iq) = upmat(ip, iq) + &
1037 (wr(k)*ra(k, ip)*rb(k, iq)*erf(a*r(k))/r(k))
1038 END DO
1039 END DO
1040 END DO
1041
1042 END SUBROUTINE num_erf
1043
1044!------------------------------------------------------------------------------
1045!
1046! Overlap with Projectors P(l,k,rc) for k=0,1,..
1047!
1048! 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]
1049!
1050! SP(l,k,p,rc) = INT(u^2 R(u) P(u))du
1051!
1052!------------------------------------------------------------------------------
1053! **************************************************************************************************
1054!> \brief ...
1055!> \param spmat ...
1056!> \param l ...
1057!> \param ra ...
1058!> \param k ...
1059!> \param rc ...
1060!> \param r ...
1061!> \param wr ...
1062! **************************************************************************************************
1063 SUBROUTINE num_proj_ol(spmat, l, ra, k, rc, r, wr)
1064
1065 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: spmat
1066 INTEGER, INTENT(IN) :: l
1067 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra
1068 INTEGER, INTENT(IN) :: k
1069 REAL(kind=dp), INTENT(IN) :: rc
1070 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
1071
1072 INTEGER :: ip, n
1073 REAL(kind=dp) :: pf
1074 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pro
1075
1076 n = SIZE(ra, 2)
1077 cpassert(SIZE(spmat) >= n)
1078
1079 ALLOCATE (pro(n))
1080
1081 pf = sqrt(2._dp)/sqrt(gamma1(l + 2*k + 1))/rc**(l + 2*k + 1.5_dp)
1082 pro(:) = pf*r(:)**(l + 2*k)*exp(-0.5_dp*(r(:)/rc)**2)
1083
1084 DO ip = 1, n
1085 spmat(ip) = sum(wr(:)*pro(:)*ra(:, ip))
1086 END DO
1087
1088 DEALLOCATE (pro)
1089
1090 END SUBROUTINE num_proj_ol
1091
1092!------------------------------------------------------------------------------
1093!
1094! Matrix elements for Gaussian potentials
1095!
1096! V(k,rc) = (r/rc)^2k exp[-1/2(r/rc)^2]
1097!
1098! VP(l,k,p+q,rc) = INT(u^2 V(u) Ra(u) Rb(u))du
1099!
1100!------------------------------------------------------------------------------
1101! **************************************************************************************************
1102!> \brief ...
1103!> \param vpmat ...
1104!> \param k ...
1105!> \param rc ...
1106!> \param ra ...
1107!> \param rb ...
1108!> \param r ...
1109!> \param wr ...
1110! **************************************************************************************************
1111 SUBROUTINE num_gpot(vpmat, k, rc, ra, rb, r, wr)
1112
1113 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: vpmat
1114 INTEGER, INTENT(IN) :: k
1115 REAL(kind=dp), INTENT(IN) :: rc
1116 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1117 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
1118
1119 INTEGER :: ip, iq, m, n
1120 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: op
1121
1122 n = SIZE(ra, 2)
1123 m = SIZE(rb, 2)
1124
1125 cpassert(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
1126
1127 ALLOCATE (op(n))
1128
1129 op(:) = (r(:)/rc)**(2*k)*exp(-0.5_dp*(r(:)/rc)**2)
1130
1131 DO iq = 1, m
1132 DO ip = 1, n
1133 vpmat(ip, iq) = sum(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1134 END DO
1135 END DO
1136
1137 DEALLOCATE (op)
1138
1139 END SUBROUTINE num_gpot
1140
1141!------------------------------------------------------------------------------
1142!
1143! G(l,k,pq) = <a|[r/rc]^2k|b>
1144! = INT(u^2 [u/rc]^2k Ra(u) Rb(u))du
1145!
1146!------------------------------------------------------------------------------
1147! **************************************************************************************************
1148!> \brief ...
1149!> \param gmat ...
1150!> \param rc ...
1151!> \param k ...
1152!> \param ra ...
1153!> \param rb ...
1154!> \param r ...
1155!> \param wr ...
1156! **************************************************************************************************
1157 SUBROUTINE num_conf(gmat, rc, k, ra, rb, r, wr)
1158
1159 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
1160 REAL(kind=dp), INTENT(IN) :: rc
1161 INTEGER, INTENT(IN) :: k
1162 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1163 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, wr
1164
1165 INTEGER :: ip, iq, m, n
1166 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: op
1167
1168 n = SIZE(ra, 2)
1169 m = SIZE(rb, 2)
1170
1171 cpassert(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
1172
1173 ALLOCATE (op(n))
1174
1175 op(:) = (r(:)/rc)**(2*k)
1176
1177 DO iq = 1, m
1178 DO ip = 1, n
1179 gmat(ip, iq) = sum(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1180 END DO
1181 END DO
1182
1183 DEALLOCATE (op)
1184
1185 END SUBROUTINE num_conf
1186
1187END MODULE ai_onecenter
subroutine, public sg_kinnuc(umat, l, pa, pb)
...
subroutine, public sg_nuclear(umat, l, pa, pb)
...
subroutine, public sg_coulomb(eri, nu, pa, lab, pc, lcd)
...
subroutine, public sg_exchange(eri, nu, pa, lac, pb, lbd)
...
subroutine, public sg_erfc(umat, l, a, pa, pb)
...
subroutine, public sg_kinetic(kmat, l, pa, pb)
...
subroutine, public sto_nuclear(umat, na, pa, nb, pb)
...
subroutine, public sto_overlap(smat, na, pa, nb, pb)
...
subroutine, public sto_kinetic(kmat, l, na, pa, nb, pb)
...
subroutine, public sg_gpot(vpmat, k, rc, l, pa, pb)
...
subroutine, public sg_proj_ol(spmat, l, p, k, rc)
...
subroutine, public sg_erf(upmat, l, a, pa, pb)
...
subroutine, public sg_conf(gmat, rc, k, l, pa, pb)
...
subroutine, public sg_overlap(smat, l, pa, pb)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), dimension(0:maxfac), parameter, public gamma1
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public gamma0
real(kind=dp), dimension(0:maxfac), parameter, public fac