Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !-----------------------------------------------------------------------------!
9 : ! 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 : #include "../base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_onecenter'
43 :
44 : PUBLIC :: sg_overlap, sg_kinetic, sg_nuclear, sg_erf, sg_gpot, &
45 : sg_proj_ol, sg_coulomb, sg_exchange, sg_kinnuc, &
46 : sg_erfc, sg_conf
47 : PUBLIC :: sto_overlap, sto_kinetic, sto_nuclear
48 :
49 : CONTAINS
50 :
51 : !------------------------------------------------------------------------------
52 : !
53 : ! S(l,pq) = pi^(1/2)*(2*l+1)!! / 2^(l+2) / (p+q)^(l+1.5)
54 : !
55 : !------------------------------------------------------------------------------
56 : ! **************************************************************************************************
57 : !> \brief ...
58 : !> \param smat ...
59 : !> \param l ...
60 : !> \param pa ...
61 : !> \param pb ...
62 : ! **************************************************************************************************
63 27683 : SUBROUTINE sg_overlap(smat, l, pa, pb)
64 :
65 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: smat
66 : INTEGER, INTENT(IN) :: l
67 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
68 :
69 : INTEGER :: ip, iq, m, n
70 : REAL(KIND=dp) :: el, spi
71 :
72 27683 : n = SIZE(pa)
73 27683 : m = SIZE(pb)
74 :
75 27683 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
76 :
77 27683 : spi = SQRT(pi)/2.0_dp**(l + 2)*dfac(2*l + 1)
78 27683 : el = REAL(l, dp) + 1.5_dp
79 :
80 167659 : DO iq = 1, m
81 2256059 : DO ip = 1, n
82 2228376 : smat(ip, iq) = spi/(pa(ip) + pb(iq))**el
83 : END DO
84 : END DO
85 :
86 27683 : END SUBROUTINE sg_overlap
87 :
88 : !------------------------------------------------------------------------------
89 : !
90 : ! T(l,pq) = (2l+3)!! pi^(1/2)/2^(l+2) [pq/(p+q)^(l+2.5)]
91 : !
92 : !------------------------------------------------------------------------------
93 : ! **************************************************************************************************
94 : !> \brief ...
95 : !> \param kmat ...
96 : !> \param l ...
97 : !> \param pa ...
98 : !> \param pb ...
99 : ! **************************************************************************************************
100 27620 : SUBROUTINE sg_kinetic(kmat, l, pa, pb)
101 :
102 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
103 : INTEGER, INTENT(IN) :: l
104 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
105 :
106 : INTEGER :: ip, iq, m, n
107 : REAL(KIND=dp) :: spi
108 :
109 27620 : n = SIZE(pa)
110 27620 : m = SIZE(pb)
111 :
112 27620 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
113 :
114 27620 : spi = dfac(2*l + 3)*SQRT(pi)/2.0_dp**(l + 2)
115 166917 : DO iq = 1, m
116 2246994 : DO ip = 1, n
117 2219374 : kmat(ip, iq) = spi*pa(ip)*pb(iq)/(pa(ip) + pb(iq))**(l + 2.5_dp)
118 : END DO
119 : END DO
120 :
121 27620 : END SUBROUTINE sg_kinetic
122 :
123 : !------------------------------------------------------------------------------
124 : !
125 : ! U(l,pq) = l!/2 / (p+q)^(l+1)
126 : !
127 : !------------------------------------------------------------------------------
128 : ! **************************************************************************************************
129 : !> \brief ...
130 : !> \param umat ...
131 : !> \param l ...
132 : !> \param pa ...
133 : !> \param pb ...
134 : ! **************************************************************************************************
135 8819 : SUBROUTINE sg_nuclear(umat, l, pa, pb)
136 :
137 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
138 : INTEGER, INTENT(IN) :: l
139 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
140 :
141 : INTEGER :: ip, iq, m, n
142 : REAL(KIND=dp) :: tld
143 :
144 8819 : n = SIZE(pa)
145 8819 : m = SIZE(pb)
146 :
147 8819 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
148 :
149 8819 : tld = 0.5_dp*fac(l)
150 62231 : DO iq = 1, m
151 1603251 : DO ip = 1, n
152 1594432 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1)
153 : END DO
154 : END DO
155 :
156 8819 : END SUBROUTINE sg_nuclear
157 :
158 : !------------------------------------------------------------------------------
159 : !
160 : ! U(l,pq) = l!/2 / (p+q)^l * [4/(p+q)^2 *pq*(l+1) + 1]
161 : !
162 : !------------------------------------------------------------------------------
163 : ! **************************************************************************************************
164 : !> \brief ...
165 : !> \param umat ...
166 : !> \param l ...
167 : !> \param pa ...
168 : !> \param pb ...
169 : ! **************************************************************************************************
170 272 : SUBROUTINE sg_kinnuc(umat, l, pa, pb)
171 :
172 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
173 : INTEGER, INTENT(IN) :: l
174 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
175 :
176 : INTEGER :: ip, iq, m, n
177 : REAL(KIND=dp) :: ppq, pq, tld
178 :
179 272 : n = SIZE(pa)
180 272 : m = SIZE(pb)
181 :
182 272 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
183 :
184 272 : IF (l > 0) THEN
185 200 : tld = 0.5_dp*fac(l)
186 5114 : DO iq = 1, m
187 167900 : DO ip = 1, n
188 162786 : ppq = pa(ip) + pb(iq)
189 162786 : pq = pa(ip)*pb(iq)
190 167700 : umat(ip, iq) = tld/ppq**l*(4.0_dp/ppq**2*pq*REAL(l + 1, dp) + 1.0_dp)
191 : END DO
192 : END DO
193 : ELSE
194 1836 : DO iq = 1, m
195 55756 : DO ip = 1, n
196 53920 : ppq = pa(ip) + pb(iq)
197 53920 : pq = pa(ip)*pb(iq)
198 55684 : umat(ip, iq) = 2.0_dp*pq/ppq**2
199 : END DO
200 : END DO
201 : END IF
202 :
203 272 : END SUBROUTINE sg_kinnuc
204 :
205 : !------------------------------------------------------------------------------
206 : !
207 : ! z = a/(p+q)
208 : !
209 : ! UP(l,pq,a) = Gamma(l+3/2)*a/SQRT(Pi)/(p+q)^(l+3/2)*
210 : ! Hypergeom([1/2, 3/2 + l], [3/2], -z)
211 : !
212 : ! UP(l,pq,a) = a/2^(l+1)/(p+q)^(l+3/2)/(1+z)^(l+1/2) * F(z,l)
213 : !
214 : ! F(z,0) = 1
215 : ! F(z,1) = 3 + 2*z
216 : ! F(z,2) = 15 + 20*z + 8*z^2
217 : ! F(z,3) = 35 + 70*z + 56*z^2 + 16*z^3
218 : ! F(z,4) = 315 + 840*z + 1008*z^2 + 576*z^3 + 128*z^4
219 : !
220 : !------------------------------------------------------------------------------
221 : ! **************************************************************************************************
222 : !> \brief ...
223 : !> \param upmat ...
224 : !> \param l ...
225 : !> \param a ...
226 : !> \param pa ...
227 : !> \param pb ...
228 : ! **************************************************************************************************
229 24353 : SUBROUTINE sg_erf(upmat, l, a, pa, pb)
230 :
231 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: upmat
232 : INTEGER, INTENT(IN) :: l
233 : REAL(KIND=dp), INTENT(IN) :: a
234 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
235 :
236 : INTEGER :: handle, ip, iq, m, n
237 : REAL(KIND=dp) :: a2, fpol, pq, tld, z, z2
238 :
239 24353 : CALL timeset("sg_erf", handle)
240 :
241 24353 : n = SIZE(pa)
242 24353 : m = SIZE(pb)
243 :
244 24353 : CPASSERT(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
245 :
246 24353 : a2 = a*a
247 24353 : tld = a/2._dp**(l + 1)
248 136281 : DO iq = 1, m
249 1404189 : DO ip = 1, n
250 1267908 : pq = pa(ip) + pb(iq)
251 1267908 : z = a2/pq
252 1379836 : upmat(ip, iq) = tld/(1._dp + z)**(l + 0.5_dp)/pq**(l + 1.5_dp)
253 : END DO
254 : END DO
255 :
256 136281 : DO iq = 1, m
257 24353 : SELECT CASE (l)
258 : CASE DEFAULT
259 0 : CPABORT("")
260 : CASE (0)
261 : ! nothing left to do
262 : CASE (1)
263 320146 : DO ip = 1, n
264 286185 : pq = pa(ip) + pb(iq)
265 286185 : z = a2/pq
266 286185 : fpol = 2.0_dp*z + 3.0_dp
267 320146 : upmat(ip, iq) = upmat(ip, iq)*fpol
268 : END DO
269 : CASE (2)
270 213350 : DO ip = 1, n
271 196679 : pq = pa(ip) + pb(iq)
272 196679 : z = a2/pq
273 196679 : fpol = 8.0_dp*z*z + 20.0_dp*z + 15.0_dp
274 213350 : upmat(ip, iq) = upmat(ip, iq)*fpol
275 : END DO
276 : CASE (3)
277 148464 : DO ip = 1, n
278 142876 : pq = pa(ip) + pb(iq)
279 142876 : z = a2/pq
280 142876 : fpol = 16.0_dp*z*z*z + 56.0_dp*z*z + 70.0_dp*z + 35.0_dp
281 142876 : fpol = 3._dp*fpol
282 148464 : upmat(ip, iq) = upmat(ip, iq)*fpol
283 : END DO
284 : CASE (4)
285 145780 : DO ip = 1, n
286 140578 : pq = pa(ip) + pb(iq)
287 140578 : z = a2/pq
288 140578 : 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
289 140578 : fpol = 3._dp*fpol
290 145780 : upmat(ip, iq) = upmat(ip, iq)*fpol
291 : END DO
292 : CASE (5)
293 145776 : DO ip = 1, n
294 140576 : pq = pa(ip) + pb(iq)
295 140576 : z = a2/pq
296 140576 : 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
297 140576 : fpol = 15._dp*fpol
298 145776 : upmat(ip, iq) = upmat(ip, iq)*fpol
299 : END DO
300 : CASE (6)
301 111928 : DO ip = 1, n
302 0 : pq = pa(ip) + pb(iq)
303 0 : z = a2/pq
304 0 : z2 = z*z
305 : fpol = 1024.0_dp*z2*z2*z2 + 6656.0_dp*z*z2*z2 + 18304.0_dp*z2*z2 + 27456.0_dp*z2*z + &
306 0 : 24024.0_dp*z2 + 12012.0_dp*z + 3003.0_dp
307 0 : fpol = 45._dp*fpol
308 0 : upmat(ip, iq) = upmat(ip, iq)*fpol
309 : END DO
310 : END SELECT
311 : END DO
312 :
313 24353 : CALL timestop(handle)
314 :
315 24353 : END SUBROUTINE sg_erf
316 :
317 : !------------------------------------------------------------------------------
318 : !
319 : ! Overlap with Projectors P(l,k,rc) for k=0,1,..
320 : !
321 : ! 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]
322 : !
323 : ! SP(l,k,p,rc) = 2^(l+k+1) / SQRT(gamma[l+2k+1.5]) / rc^(l+2k+1.5)
324 : ! * Gamma(l+k+1.5) / (2p+1/rc^2)^(l+k+1.5)
325 : !
326 : !------------------------------------------------------------------------------
327 : ! **************************************************************************************************
328 : !> \brief ...
329 : !> \param spmat ...
330 : !> \param l ...
331 : !> \param p ...
332 : !> \param k ...
333 : !> \param rc ...
334 : ! **************************************************************************************************
335 10617 : SUBROUTINE sg_proj_ol(spmat, l, p, k, rc)
336 :
337 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: spmat
338 : INTEGER, INTENT(IN) :: l
339 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p
340 : INTEGER, INTENT(IN) :: k
341 : REAL(KIND=dp), INTENT(IN) :: rc
342 :
343 : REAL(KIND=dp) :: orc, pf
344 :
345 10617 : CPASSERT(SIZE(spmat) >= SIZE(p))
346 :
347 10617 : pf = 2._dp**(l + k + 1)*gamma1(l + k + 1)/rc**(l + 2*k + 1.5_dp)/SQRT(gamma1(l + 2*k + 1))
348 10617 : orc = 1._dp/(rc*rc)
349 73137 : spmat(:) = pf/(2._dp*p(:) + orc)**(l + k + 1.5_dp)
350 :
351 10617 : END SUBROUTINE sg_proj_ol
352 :
353 : !------------------------------------------------------------------------------
354 : !
355 : ! Matrix elements for Gaussian potentials
356 : !
357 : ! V(k,rc) = (r/rc)^2k exp[-1/2(r/rc)^2]
358 : !
359 : ! 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)
360 : !
361 : !------------------------------------------------------------------------------
362 : ! **************************************************************************************************
363 : !> \brief ...
364 : !> \param vpmat ...
365 : !> \param k ...
366 : !> \param rc ...
367 : !> \param l ...
368 : !> \param pa ...
369 : !> \param pb ...
370 : ! **************************************************************************************************
371 42410 : SUBROUTINE sg_gpot(vpmat, k, rc, l, pa, pb)
372 :
373 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: vpmat
374 : INTEGER, INTENT(IN) :: k
375 : REAL(KIND=dp), INTENT(IN) :: rc
376 : INTEGER, INTENT(IN) :: l
377 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
378 :
379 : INTEGER :: ip, iq, m, n
380 : REAL(KIND=dp) :: tld
381 :
382 42410 : n = SIZE(pa)
383 42410 : m = SIZE(pb)
384 :
385 42410 : CPASSERT(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
386 :
387 42410 : tld = gamma1(l + k + 1)*rc**(2*l + 3)*2._dp**(l + k + 0.5)
388 :
389 247578 : DO iq = 1, m
390 2692934 : DO ip = 1, n
391 2650524 : vpmat(ip, iq) = tld/(1._dp + 2._dp*rc*rc*(pa(ip) + pb(iq)))**(l + k + 1.5_dp)
392 : END DO
393 : END DO
394 :
395 42410 : END SUBROUTINE sg_gpot
396 :
397 : !------------------------------------------------------------------------------
398 : !
399 : ! G(l,k,pq) = <a|[r/rc]^2k|b>
400 : ! = 0.5*Gamma(l+k+1.5)/rc^(2k)/(p+q)^(l+k+1.5)
401 : !
402 : !------------------------------------------------------------------------------
403 : ! **************************************************************************************************
404 : !> \brief ...
405 : !> \param gmat ...
406 : !> \param rc ...
407 : !> \param k ...
408 : !> \param l ...
409 : !> \param pa ...
410 : !> \param pb ...
411 : ! **************************************************************************************************
412 75 : SUBROUTINE sg_conf(gmat, rc, k, l, pa, pb)
413 :
414 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
415 : REAL(KIND=dp), INTENT(IN) :: rc
416 : INTEGER, INTENT(IN) :: k, l
417 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
418 :
419 : INTEGER :: ip, iq, m, n
420 : REAL(KIND=dp) :: tld
421 :
422 75 : n = SIZE(pa)
423 75 : m = SIZE(pb)
424 :
425 75 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
426 :
427 75 : tld = 0.5_dp/rc**(2*k)*gamma1(l + k + 1)
428 375 : DO iq = 1, m
429 1575 : DO ip = 1, n
430 1500 : gmat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + k + 1.5_dp)
431 : END DO
432 : END DO
433 :
434 75 : END SUBROUTINE sg_conf
435 :
436 : !------------------------------------------------------------------------------
437 : !
438 : ! (plql,rl'sl')
439 : !
440 : !------------------------------------------------------------------------------
441 : ! **************************************************************************************************
442 : !> \brief ...
443 : !> \param eri ...
444 : !> \param nu ...
445 : !> \param pa ...
446 : !> \param lab ...
447 : !> \param pc ...
448 : !> \param lcd ...
449 : ! **************************************************************************************************
450 1218 : SUBROUTINE sg_coulomb(eri, nu, pa, lab, pc, lcd)
451 :
452 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eri
453 : INTEGER, INTENT(IN) :: nu
454 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
455 : INTEGER, INTENT(IN) :: lab
456 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pc
457 : INTEGER, INTENT(IN) :: lcd
458 :
459 : INTEGER :: handle, ia, ib, ic, id, jab, jcd, na, nc
460 : REAL(KIND=dp) :: cc1, cc2, ee, p, q, r, s, vab1, vab3, &
461 : vcd1, vcd3, xab, xcd
462 :
463 1218 : CALL timeset("sg_coulomb", handle)
464 :
465 1218 : na = SIZE(pa)
466 1218 : nc = SIZE(pc)
467 1218 : ee = SQRT(2.0_dp*pi)/2.0_dp**(2*(lab + lcd) + 6)
468 1218 : jab = 0
469 12782 : DO ia = 1, na
470 11564 : p = pa(ia)
471 141744 : DO ib = ia, na
472 128962 : jab = jab + 1
473 128962 : q = pa(ib)
474 128962 : xab = 0.5_dp*(p + q)
475 128962 : vab1 = vgau(2*lab - nu + 1, xab)
476 128962 : vab3 = vgau(2*lab + nu + 2, xab)
477 128962 : jcd = 0
478 3272892 : DO ic = 1, nc
479 3132366 : r = pc(ic)
480 45753440 : DO id = ic, nc
481 42492112 : jcd = jcd + 1
482 42492112 : s = pc(id)
483 42492112 : xcd = 0.5_dp*(r + s)
484 42492112 : vcd1 = vgau(2*lcd + nu + 2, xcd)
485 42492112 : vcd3 = vgau(2*lcd - nu + 1, xcd)
486 42492112 : cc1 = cgau(2*lab - nu + 1, 2*lcd + nu + 2, xab/xcd)
487 42492112 : cc2 = cgau(2*lcd - nu + 1, 2*lab + nu + 2, xcd/xab)
488 :
489 45624478 : eri(jab, jcd) = ee*(cc1*vab1*vcd1 + cc2*vab3*vcd3)
490 :
491 : END DO
492 : END DO
493 : END DO
494 : END DO
495 :
496 1218 : CALL timestop(handle)
497 :
498 1218 : END SUBROUTINE sg_coulomb
499 :
500 : !------------------------------------------------------------------------------
501 : !
502 : ! (plql',rlsl')
503 : !
504 : !------------------------------------------------------------------------------
505 : ! **************************************************************************************************
506 : !> \brief ...
507 : !> \param eri ...
508 : !> \param nu ...
509 : !> \param pa ...
510 : !> \param lac ...
511 : !> \param pb ...
512 : !> \param lbd ...
513 : ! **************************************************************************************************
514 2128 : SUBROUTINE sg_exchange(eri, nu, pa, lac, pb, lbd)
515 :
516 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eri
517 : INTEGER, INTENT(IN) :: nu
518 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
519 : INTEGER, INTENT(IN) :: lac
520 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
521 : INTEGER, INTENT(IN) :: lbd
522 :
523 : INTEGER :: handle, ia, ib, ic, id, jac, jbd, na, nb
524 : REAL(KIND=dp) :: cc1, cc2, cc3, cc4, ee, p, q, r, s, &
525 : v1pr, v1ps, v1qr, v1qs, v2pr, v2ps, &
526 : v2qr, v2qs, xab, xad, xbc, xcd
527 :
528 2128 : CALL timeset("sg_exchange", handle)
529 :
530 40295595 : eri = 0.0_dp
531 2128 : na = SIZE(pa)
532 2128 : nb = SIZE(pb)
533 2128 : ee = SQRT(2.0_dp*pi)/2.0_dp**(2*(lac + lbd) + 7)
534 2128 : jac = 0
535 15842 : DO ia = 1, na
536 13714 : p = pa(ia)
537 162070 : DO ic = ia, na
538 146228 : jac = jac + 1
539 146228 : q = pa(ic)
540 146228 : jbd = 0
541 3424067 : DO ib = 1, nb
542 3264125 : r = pb(ib)
543 3264125 : xab = 0.5_dp*(p + r)
544 3264125 : xbc = 0.5_dp*(q + r)
545 43482700 : DO id = ib, nb
546 40072347 : jbd = jbd + 1
547 40072347 : s = pb(id)
548 40072347 : xcd = 0.5_dp*(q + s)
549 40072347 : xad = 0.5_dp*(p + s)
550 40072347 : v1pr = vgau(lac + lbd - nu + 1, xab)
551 40072347 : v1qs = vgau(lac + lbd - nu + 1, xcd)
552 40072347 : v1ps = vgau(lac + lbd - nu + 1, xad)
553 40072347 : v1qr = vgau(lac + lbd - nu + 1, xbc)
554 40072347 : v2qs = vgau(lac + lbd + nu + 2, xcd)
555 40072347 : v2pr = vgau(lac + lbd + nu + 2, xab)
556 40072347 : v2qr = vgau(lac + lbd + nu + 2, xbc)
557 40072347 : v2ps = vgau(lac + lbd + nu + 2, xad)
558 40072347 : cc1 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xab/xcd)
559 40072347 : cc2 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xcd/xab)
560 40072347 : cc3 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xad/xbc)
561 40072347 : cc4 = cgau(lac + lbd - nu + 1, lac + lbd + nu + 2, xbc/xad)
562 :
563 : eri(jac, jbd) = ee*(v1pr*v2qs*cc1 + v1qs*v2pr*cc2 + &
564 43336472 : v1ps*v2qr*cc3 + v1qr*v2ps*cc4)
565 :
566 : END DO
567 : END DO
568 : END DO
569 : END DO
570 :
571 2128 : CALL timestop(handle)
572 :
573 2128 : END SUBROUTINE sg_exchange
574 :
575 : !------------------------------------------------------------------------------
576 : !
577 : ! erfc(a*x)/x = 1/x - erf(a*x)/x
578 : !
579 : !------------------------------------------------------------------------------
580 : ! **************************************************************************************************
581 : !> \brief ...
582 : !> \param umat ...
583 : !> \param l ...
584 : !> \param a ...
585 : !> \param pa ...
586 : !> \param pb ...
587 : ! **************************************************************************************************
588 192 : SUBROUTINE sg_erfc(umat, l, a, pa, pb)
589 :
590 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
591 : INTEGER, INTENT(IN) :: l
592 : REAL(KIND=dp), INTENT(IN) :: a
593 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa, pb
594 :
595 : INTEGER :: ip, iq, m, n
596 : REAL(KIND=dp) :: tld
597 :
598 192 : n = SIZE(pa)
599 192 : m = SIZE(pb)
600 :
601 192 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
602 :
603 192 : CALL sg_erf(umat, l, a, pa, pb)
604 :
605 192 : tld = 0.5_dp*fac(l)
606 1108 : DO iq = 1, m
607 12104 : DO ip = 1, n
608 11912 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1) - umat(ip, iq)
609 : END DO
610 : END DO
611 :
612 192 : END SUBROUTINE sg_erfc
613 :
614 : ! ***************************************************************************************************
615 :
616 : ! **************************************************************************************************
617 : !> \brief ...
618 : !> \param n ...
619 : !> \param x ...
620 : !> \return ...
621 : ! **************************************************************************************************
622 405820924 : ELEMENTAL FUNCTION vgau(n, x) RESULT(v)
623 : INTEGER, INTENT(IN) :: n
624 : REAL(KIND=dp), INTENT(IN) :: x
625 : REAL(KIND=dp) :: v
626 :
627 405820924 : v = dfac(n - 1)/x**(0.5_dp*(n + 1))
628 :
629 405820924 : END FUNCTION vgau
630 :
631 : ! **************************************************************************************************
632 : !> \brief ...
633 : !> \param a ...
634 : !> \param b ...
635 : !> \param t ...
636 : !> \return ...
637 : ! **************************************************************************************************
638 245273612 : ELEMENTAL FUNCTION cgau(a, b, t) RESULT(c)
639 : INTEGER, INTENT(IN) :: a, b
640 : REAL(KIND=dp), INTENT(IN) :: t
641 : REAL(KIND=dp) :: c
642 :
643 : INTEGER :: l
644 :
645 245273612 : c = 0.0_dp
646 761069194 : DO l = 0, (a - 1)/2
647 761069194 : c = c + (t/(1.0_dp + t))**l*dfac(2*l + b - 1)/dfac(2*l)
648 : END DO
649 245273612 : c = c*(1.0_dp + t)**(-0.5_dp*(b + 1))/dfac(b - 1)
650 :
651 245273612 : END FUNCTION cgau
652 :
653 : !------------------------------------------------------------------------------
654 : !
655 : ! S(l,pn,qm) = ( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m,(p+q)/2]
656 : !
657 : !------------------------------------------------------------------------------
658 : ! **************************************************************************************************
659 : !> \brief ...
660 : !> \param smat ...
661 : !> \param na ...
662 : !> \param pa ...
663 : !> \param nb ...
664 : !> \param pb ...
665 : ! **************************************************************************************************
666 7296 : SUBROUTINE sto_overlap(smat, na, pa, nb, pb)
667 :
668 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: smat
669 : INTEGER, DIMENSION(:), INTENT(IN) :: na
670 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
671 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
672 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
673 :
674 : INTEGER :: ip, iq, m, n
675 : REAL(KIND=dp) :: vp, vpq, vq
676 :
677 7296 : n = SIZE(pa)
678 7296 : m = SIZE(pb)
679 :
680 7296 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
681 :
682 11558 : DO iq = 1, m
683 4262 : vq = vsto(2*nb(iq), pb(iq))
684 23464 : DO ip = 1, n
685 11906 : vp = vsto(2*na(ip), pa(ip))
686 11906 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
687 16168 : smat(ip, iq) = vpq/SQRT(vp*vq)
688 : END DO
689 : END DO
690 :
691 7296 : END SUBROUTINE sto_overlap
692 :
693 : !------------------------------------------------------------------------------
694 : !
695 : ! T(l,pn,qm) = 0.5*p*q*( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m,(p+q)/2]
696 : ! -(W[l,n,p]+W[l,m,q]) * V[n+m-1,(p+q)/2]
697 : ! +W[l,n,p]*W[l,m,q] * V[n+m-2,(p+q)/2]
698 : !
699 : !------------------------------------------------------------------------------
700 : ! **************************************************************************************************
701 : !> \brief ...
702 : !> \param kmat ...
703 : !> \param l ...
704 : !> \param na ...
705 : !> \param pa ...
706 : !> \param nb ...
707 : !> \param pb ...
708 : ! **************************************************************************************************
709 7296 : SUBROUTINE sto_kinetic(kmat, l, na, pa, nb, pb)
710 :
711 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
712 : INTEGER, INTENT(IN) :: l
713 : INTEGER, DIMENSION(:), INTENT(IN) :: na
714 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
715 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
716 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
717 :
718 : INTEGER :: ip, iq, m, n
719 : REAL(KIND=dp) :: vp, vpq, vpq1, vpq2, vq, wp, wq
720 :
721 7296 : n = SIZE(pa)
722 7296 : m = SIZE(pb)
723 :
724 7296 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
725 :
726 11558 : DO iq = 1, m
727 4262 : vq = vsto(2*nb(iq), pb(iq))
728 4262 : wq = wsto(l, nb(iq), pb(iq))
729 23464 : DO ip = 1, n
730 11906 : vp = vsto(2*na(ip), pa(ip))
731 11906 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
732 11906 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
733 11906 : vpq2 = vsto(na(ip) + nb(iq) - 2, 0.5_dp*(pa(ip) + pb(iq)))
734 11906 : wp = wsto(l, na(ip), pa(ip))
735 : kmat(ip, iq) = 0.5_dp*pa(ip)*pb(iq)/SQRT(vp*vq)* &
736 16168 : (vpq - (wp + wq)*vpq1 + wp*wq*vpq2)
737 : END DO
738 : END DO
739 :
740 7296 : END SUBROUTINE sto_kinetic
741 :
742 : !------------------------------------------------------------------------------
743 : !
744 : ! U(l,pq) = 2( V[2n,p]*V[2m,q] )^(-1/2) * V[n+m-1,(p+q)/2]
745 : !
746 : !------------------------------------------------------------------------------
747 : ! **************************************************************************************************
748 : !> \brief ...
749 : !> \param umat ...
750 : !> \param na ...
751 : !> \param pa ...
752 : !> \param nb ...
753 : !> \param pb ...
754 : ! **************************************************************************************************
755 8244 : SUBROUTINE sto_nuclear(umat, na, pa, nb, pb)
756 :
757 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
758 : INTEGER, DIMENSION(:), INTENT(IN) :: na
759 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
760 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
761 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
762 :
763 : INTEGER :: ip, iq, m, n
764 : REAL(KIND=dp) :: vp, vpq1, vq
765 :
766 8244 : n = SIZE(pa)
767 8244 : m = SIZE(pb)
768 :
769 8244 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
770 :
771 12758 : DO iq = 1, m
772 4514 : vq = vsto(2*nb(iq), pb(iq))
773 24916 : DO ip = 1, n
774 12158 : vp = vsto(2*na(ip), pa(ip))
775 12158 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
776 16672 : umat(ip, iq) = 2._dp/SQRT(vp*vq)*vpq1
777 : END DO
778 : END DO
779 :
780 8244 : END SUBROUTINE sto_nuclear
781 :
782 : !------------------------------------------------------------------------------
783 : !
784 : ! G(l,k,pq) = <aln|[r/rc]^2k|blm>
785 : ! = N(aln)*N(blm) (a+b)^(-(n+m+2k+1))/rc^2k * GAMMA(n+m+2k+1)
786 : !
787 : !------------------------------------------------------------------------------
788 : ! **************************************************************************************************
789 : !> \brief ...
790 : !> \param gmat ...
791 : !> \param rc ...
792 : !> \param k ...
793 : !> \param na ...
794 : !> \param pa ...
795 : !> \param nb ...
796 : !> \param pb ...
797 : ! **************************************************************************************************
798 0 : SUBROUTINE sto_conf(gmat, rc, k, na, pa, nb, pb)
799 :
800 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
801 : REAL(KIND=dp), INTENT(IN) :: rc
802 : INTEGER, INTENT(IN) :: k
803 : INTEGER, DIMENSION(:), INTENT(IN) :: na
804 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pa
805 : INTEGER, DIMENSION(:), INTENT(IN) :: nb
806 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pb
807 :
808 : INTEGER :: ip, iq, m, n
809 :
810 0 : n = SIZE(pa)
811 0 : m = SIZE(pb)
812 :
813 0 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
814 :
815 0 : DO iq = 1, m
816 0 : DO ip = 1, n
817 : gmat(ip, iq) = (2._dp*pa(ip))**(na(ip) + 0.5_dp)/SQRT(fac(2*na(ip))) &
818 : *(2._dp*pb(iq))**(nb(iq) + 0.5_dp)/SQRT(fac(2*nb(iq))) &
819 : /rc**(2*k)/(pa(ip) + pb(iq))**(na(ip) + nb(iq) + 2*k + 1) &
820 0 : *gamma0(na(ip) + nb(iq) + 2*k + 1)
821 : END DO
822 : END DO
823 :
824 0 : END SUBROUTINE sto_conf
825 :
826 : ! **************************************************************************************************
827 :
828 : ! **************************************************************************************************
829 : !> \brief ...
830 : !> \param n ...
831 : !> \param x ...
832 : !> \return ...
833 : ! **************************************************************************************************
834 108790 : ELEMENTAL FUNCTION vsto(n, x) RESULT(v)
835 : INTEGER, INTENT(IN) :: n
836 : REAL(KIND=dp), INTENT(IN) :: x
837 : REAL(KIND=dp) :: v
838 :
839 108790 : v = fac(n)/x**(n + 1)
840 :
841 108790 : END FUNCTION vsto
842 :
843 : ! **************************************************************************************************
844 : !> \brief ...
845 : !> \param n ...
846 : !> \param m ...
847 : !> \param x ...
848 : !> \return ...
849 : ! **************************************************************************************************
850 16168 : ELEMENTAL FUNCTION wsto(n, m, x) RESULT(w)
851 : INTEGER, INTENT(IN) :: n, m
852 : REAL(KIND=dp), INTENT(IN) :: x
853 : REAL(KIND=dp) :: w
854 :
855 16168 : w = 2._dp*REAL(m - n - 1, dp)/x
856 :
857 16168 : END FUNCTION wsto
858 : !------------------------------------------------------------------------------
859 : !
860 : ! S(l,pq) = INT(u^2 Ra(u) Rb(u))du
861 : !
862 : !------------------------------------------------------------------------------
863 : ! **************************************************************************************************
864 : !> \brief ...
865 : !> \param smat ...
866 : !> \param ra ...
867 : !> \param rb ...
868 : !> \param wr ...
869 : ! **************************************************************************************************
870 0 : SUBROUTINE num_overlap(smat, ra, rb, wr)
871 :
872 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: smat
873 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
874 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wr
875 :
876 : INTEGER :: ip, iq, m, n
877 :
878 0 : n = SIZE(ra, 2)
879 0 : m = SIZE(rb, 2)
880 :
881 0 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
882 :
883 0 : DO iq = 1, m
884 0 : DO ip = 1, n
885 0 : smat(ip, iq) = SUM(ra(:, ip)*rb(:, iq)*wr(:))
886 : END DO
887 : END DO
888 :
889 0 : END SUBROUTINE num_overlap
890 :
891 : !------------------------------------------------------------------------------
892 : !
893 : ! T(l,pq) = 0.5 INT( u^2 dRa(u) dRb(u) + l(l+1) Ra(u) Rb(u))du
894 : !
895 : !------------------------------------------------------------------------------
896 : ! **************************************************************************************************
897 : !> \brief ...
898 : !> \param kmat ...
899 : !> \param l ...
900 : !> \param ra ...
901 : !> \param dra ...
902 : !> \param rb ...
903 : !> \param drb ...
904 : !> \param r ...
905 : !> \param wr ...
906 : ! **************************************************************************************************
907 0 : SUBROUTINE num_kinetic(kmat, l, ra, dra, rb, drb, r, wr)
908 :
909 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: kmat
910 : INTEGER, INTENT(IN) :: l
911 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, dra, rb, drb
912 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
913 :
914 : INTEGER :: ip, iq, m, n
915 :
916 0 : n = SIZE(ra, 2)
917 0 : m = SIZE(rb, 2)
918 :
919 0 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
920 :
921 0 : DO iq = 1, m
922 0 : DO ip = 1, n
923 : kmat(ip, iq) = 0.5_dp*SUM(wr(:)*dra(:, ip)*drb(:, iq) &
924 0 : + wr(:)*REAL(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**2)
925 : END DO
926 : END DO
927 :
928 0 : END SUBROUTINE num_kinetic
929 :
930 : !------------------------------------------------------------------------------
931 : !
932 : ! U(l,pq) = INT(u Ra(u) Rb(u))du
933 : !
934 : !------------------------------------------------------------------------------
935 : ! **************************************************************************************************
936 : !> \brief ...
937 : !> \param umat ...
938 : !> \param ra ...
939 : !> \param rb ...
940 : !> \param r ...
941 : !> \param wr ...
942 : ! **************************************************************************************************
943 0 : SUBROUTINE num_nuclear(umat, ra, rb, r, wr)
944 :
945 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
946 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
947 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
948 :
949 : INTEGER :: ip, iq, m, n
950 :
951 0 : n = SIZE(ra, 2)
952 0 : m = SIZE(rb, 2)
953 :
954 0 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
955 :
956 0 : DO iq = 1, m
957 0 : DO ip = 1, n
958 0 : umat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)/r(:))
959 : END DO
960 : END DO
961 :
962 0 : END SUBROUTINE num_nuclear
963 :
964 : !------------------------------------------------------------------------------
965 : !
966 : ! U(l,pq) = INT(u dRa(u) dRb(u))du + l(l+1) * INT(Ra(u) Rb(u) / u)du
967 : !
968 : !------------------------------------------------------------------------------
969 : ! **************************************************************************************************
970 : !> \brief ...
971 : !> \param umat ...
972 : !> \param l ...
973 : !> \param ra ...
974 : !> \param dra ...
975 : !> \param rb ...
976 : !> \param drb ...
977 : !> \param r ...
978 : !> \param wr ...
979 : ! **************************************************************************************************
980 0 : SUBROUTINE num_kinnuc(umat, l, ra, dra, rb, drb, r, wr)
981 :
982 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: umat
983 : INTEGER, INTENT(IN) :: l
984 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, dra, rb, drb
985 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
986 :
987 : INTEGER :: ip, iq, m, n
988 :
989 0 : n = SIZE(ra, 2)
990 0 : m = SIZE(rb, 2)
991 :
992 0 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
993 :
994 0 : DO iq = 1, m
995 0 : DO ip = 1, n
996 : umat(ip, iq) = SUM(wr(:)*dra(:, ip)*drb(:, iq)/r(:) &
997 0 : + wr(:)*REAL(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**3)
998 : END DO
999 : END DO
1000 :
1001 0 : END SUBROUTINE num_kinnuc
1002 :
1003 : !------------------------------------------------------------------------------
1004 : !
1005 : ! U(l,pq) = INT(u erf(a*u) Ra(u) Rb(u))du
1006 : !
1007 : !------------------------------------------------------------------------------
1008 : ! **************************************************************************************************
1009 : !> \brief ...
1010 : !> \param upmat ...
1011 : !> \param a ...
1012 : !> \param ra ...
1013 : !> \param rb ...
1014 : !> \param r ...
1015 : !> \param wr ...
1016 : ! **************************************************************************************************
1017 0 : SUBROUTINE num_erf(upmat, a, ra, rb, r, wr)
1018 :
1019 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: upmat
1020 : REAL(KIND=dp), INTENT(IN) :: a
1021 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1022 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
1023 :
1024 : INTEGER :: ip, iq, k, m, n
1025 :
1026 0 : n = SIZE(ra, 2)
1027 0 : m = SIZE(rb, 2)
1028 :
1029 0 : CPASSERT(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
1030 :
1031 0 : DO iq = 1, m
1032 0 : DO ip = 1, n
1033 0 : upmat(ip, iq) = 0._dp
1034 0 : DO k = 1, SIZE(r)
1035 : upmat(ip, iq) = upmat(ip, iq) + &
1036 0 : (wr(k)*ra(k, ip)*rb(k, iq)*erf(a*r(k))/r(k))
1037 : END DO
1038 : END DO
1039 : END DO
1040 :
1041 0 : END SUBROUTINE num_erf
1042 :
1043 : !------------------------------------------------------------------------------
1044 : !
1045 : ! Overlap with Projectors P(l,k,rc) for k=0,1,..
1046 : !
1047 : ! 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]
1048 : !
1049 : ! SP(l,k,p,rc) = INT(u^2 R(u) P(u))du
1050 : !
1051 : !------------------------------------------------------------------------------
1052 : ! **************************************************************************************************
1053 : !> \brief ...
1054 : !> \param spmat ...
1055 : !> \param l ...
1056 : !> \param ra ...
1057 : !> \param k ...
1058 : !> \param rc ...
1059 : !> \param r ...
1060 : !> \param wr ...
1061 : ! **************************************************************************************************
1062 0 : SUBROUTINE num_proj_ol(spmat, l, ra, k, rc, r, wr)
1063 :
1064 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: spmat
1065 : INTEGER, INTENT(IN) :: l
1066 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra
1067 : INTEGER, INTENT(IN) :: k
1068 : REAL(KIND=dp), INTENT(IN) :: rc
1069 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
1070 :
1071 : INTEGER :: ip, n
1072 : REAL(KIND=dp) :: pf
1073 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pro
1074 :
1075 0 : n = SIZE(ra, 2)
1076 0 : CPASSERT(SIZE(spmat) >= n)
1077 :
1078 0 : ALLOCATE (pro(n))
1079 :
1080 0 : pf = SQRT(2._dp)/SQRT(gamma1(l + 2*k + 1))/rc**(l + 2*k + 1.5_dp)
1081 0 : pro(:) = pf*r(:)**(l + 2*k)*EXP(-0.5_dp*(r(:)/rc)**2)
1082 :
1083 0 : DO ip = 1, n
1084 0 : spmat(ip) = SUM(wr(:)*pro(:)*ra(:, ip))
1085 : END DO
1086 :
1087 0 : DEALLOCATE (pro)
1088 :
1089 0 : END SUBROUTINE num_proj_ol
1090 :
1091 : !------------------------------------------------------------------------------
1092 : !
1093 : ! Matrix elements for Gaussian potentials
1094 : !
1095 : ! V(k,rc) = (r/rc)^2k exp[-1/2(r/rc)^2]
1096 : !
1097 : ! VP(l,k,p+q,rc) = INT(u^2 V(u) Ra(u) Rb(u))du
1098 : !
1099 : !------------------------------------------------------------------------------
1100 : ! **************************************************************************************************
1101 : !> \brief ...
1102 : !> \param vpmat ...
1103 : !> \param k ...
1104 : !> \param rc ...
1105 : !> \param ra ...
1106 : !> \param rb ...
1107 : !> \param r ...
1108 : !> \param wr ...
1109 : ! **************************************************************************************************
1110 0 : SUBROUTINE num_gpot(vpmat, k, rc, ra, rb, r, wr)
1111 :
1112 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: vpmat
1113 : INTEGER, INTENT(IN) :: k
1114 : REAL(KIND=dp), INTENT(IN) :: rc
1115 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1116 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
1117 :
1118 : INTEGER :: ip, iq, m, n
1119 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: op
1120 :
1121 0 : n = SIZE(ra, 2)
1122 0 : m = SIZE(rb, 2)
1123 :
1124 0 : CPASSERT(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
1125 :
1126 0 : ALLOCATE (op(n))
1127 :
1128 0 : op(:) = (r(:)/rc)**(2*k)*EXP(-0.5_dp*(r(:)/rc)**2)
1129 :
1130 0 : DO iq = 1, m
1131 0 : DO ip = 1, n
1132 0 : vpmat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1133 : END DO
1134 : END DO
1135 :
1136 0 : DEALLOCATE (op)
1137 :
1138 0 : END SUBROUTINE num_gpot
1139 :
1140 : !------------------------------------------------------------------------------
1141 : !
1142 : ! G(l,k,pq) = <a|[r/rc]^2k|b>
1143 : ! = INT(u^2 [u/rc]^2k Ra(u) Rb(u))du
1144 : !
1145 : !------------------------------------------------------------------------------
1146 : ! **************************************************************************************************
1147 : !> \brief ...
1148 : !> \param gmat ...
1149 : !> \param rc ...
1150 : !> \param k ...
1151 : !> \param ra ...
1152 : !> \param rb ...
1153 : !> \param r ...
1154 : !> \param wr ...
1155 : ! **************************************************************************************************
1156 0 : SUBROUTINE num_conf(gmat, rc, k, ra, rb, r, wr)
1157 :
1158 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gmat
1159 : REAL(KIND=dp), INTENT(IN) :: rc
1160 : INTEGER, INTENT(IN) :: k
1161 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
1162 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r, wr
1163 :
1164 : INTEGER :: ip, iq, m, n
1165 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: op
1166 :
1167 0 : n = SIZE(ra, 2)
1168 0 : m = SIZE(rb, 2)
1169 :
1170 0 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
1171 :
1172 0 : ALLOCATE (op(n))
1173 :
1174 0 : op(:) = (r(:)/rc)**(2*k)
1175 :
1176 0 : DO iq = 1, m
1177 0 : DO ip = 1, n
1178 0 : gmat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1179 : END DO
1180 : END DO
1181 :
1182 0 : DEALLOCATE (op)
1183 :
1184 0 : END SUBROUTINE num_conf
1185 :
1186 : END MODULE ai_onecenter
|
