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 : 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 28511 : 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 28511 : n = SIZE(pa)
74 28511 : m = SIZE(pb)
75 :
76 28511 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
77 :
78 28511 : spi = rootpi/2.0_dp**(l + 2)*dfac(2*l + 1)
79 28511 : el = REAL(l, dp) + 1.5_dp
80 :
81 171385 : DO iq = 1, m
82 2282027 : DO ip = 1, n
83 2253516 : smat(ip, iq) = spi/(pa(ip) + pb(iq))**el
84 : END DO
85 : END DO
86 :
87 28511 : 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 28448 : 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 28448 : n = SIZE(pa)
111 28448 : m = SIZE(pb)
112 :
113 28448 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
114 :
115 28448 : spi = dfac(2*l + 3)*rootpi/2.0_dp**(l + 2)
116 170643 : DO iq = 1, m
117 2272962 : DO ip = 1, n
118 2244514 : kmat(ip, iq) = spi*pa(ip)*pb(iq)/(pa(ip) + pb(iq))**(l + 2.5_dp)
119 : END DO
120 : END DO
121 :
122 28448 : 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 8775 : 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 8775 : n = SIZE(pa)
146 8775 : m = SIZE(pb)
147 :
148 8775 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
149 :
150 8775 : tld = 0.5_dp*fac(l)
151 62091 : DO iq = 1, m
152 1602367 : DO ip = 1, n
153 1593592 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1)
154 : END DO
155 : END DO
156 :
157 8775 : 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 272 : 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 272 : n = SIZE(pa)
181 272 : m = SIZE(pb)
182 :
183 272 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
184 :
185 272 : IF (l > 0) THEN
186 200 : tld = 0.5_dp*fac(l)
187 5114 : DO iq = 1, m
188 167900 : DO ip = 1, n
189 162786 : ppq = pa(ip) + pb(iq)
190 162786 : pq = pa(ip)*pb(iq)
191 167700 : 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 1836 : DO iq = 1, m
196 55756 : DO ip = 1, n
197 53920 : ppq = pa(ip) + pb(iq)
198 53920 : pq = pa(ip)*pb(iq)
199 55684 : umat(ip, iq) = 2.0_dp*pq/ppq**2
200 : END DO
201 : END DO
202 : END IF
203 :
204 272 : 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 25197 : 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 25197 : CALL timeset("sg_erf", handle)
241 :
242 25197 : n = SIZE(pa)
243 25197 : m = SIZE(pb)
244 :
245 25197 : CPASSERT(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
246 :
247 25197 : a2 = a*a
248 25197 : tld = a/2._dp**(l + 1)
249 140119 : DO iq = 1, m
250 1431013 : DO ip = 1, n
251 1290894 : pq = pa(ip) + pb(iq)
252 1290894 : z = a2/pq
253 1405816 : upmat(ip, iq) = tld/(1._dp + z)**(l + 0.5_dp)/pq**(l + 1.5_dp)
254 : END DO
255 : END DO
256 :
257 140119 : DO iq = 1, m
258 25197 : SELECT CASE (l)
259 : CASE DEFAULT
260 0 : CPABORT("")
261 : CASE (0)
262 : ! nothing left to do
263 : CASE (1)
264 327650 : DO ip = 1, n
265 292725 : pq = pa(ip) + pb(iq)
266 292725 : z = a2/pq
267 292725 : fpol = 2.0_dp*z + 3.0_dp
268 327650 : upmat(ip, iq) = upmat(ip, iq)*fpol
269 : END DO
270 : CASE (2)
271 216622 : DO ip = 1, n
272 199449 : pq = pa(ip) + pb(iq)
273 199449 : z = a2/pq
274 199449 : fpol = 8.0_dp*z*z + 20.0_dp*z + 15.0_dp
275 216622 : upmat(ip, iq) = upmat(ip, iq)*fpol
276 : END DO
277 : CASE (3)
278 148464 : DO ip = 1, n
279 142876 : pq = pa(ip) + pb(iq)
280 142876 : z = a2/pq
281 142876 : fpol = 16.0_dp*z*z*z + 56.0_dp*z*z + 70.0_dp*z + 35.0_dp
282 142876 : fpol = 3._dp*fpol
283 148464 : upmat(ip, iq) = upmat(ip, iq)*fpol
284 : END DO
285 : CASE (4)
286 145780 : DO ip = 1, n
287 140578 : pq = pa(ip) + pb(iq)
288 140578 : z = a2/pq
289 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
290 140578 : fpol = 3._dp*fpol
291 145780 : upmat(ip, iq) = upmat(ip, iq)*fpol
292 : END DO
293 : CASE (5)
294 145776 : DO ip = 1, n
295 140576 : pq = pa(ip) + pb(iq)
296 140576 : z = a2/pq
297 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
298 140576 : fpol = 15._dp*fpol
299 145776 : upmat(ip, iq) = upmat(ip, iq)*fpol
300 : END DO
301 : CASE (6)
302 114922 : 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 25197 : CALL timestop(handle)
315 :
316 25197 : 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 11121 : 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 11121 : CPASSERT(SIZE(spmat) >= SIZE(p))
347 :
348 11121 : pf = 2._dp**(l + k + 1)*gamma1(l + k + 1)/rc**(l + 2*k + 1.5_dp)/SQRT(gamma1(l + 2*k + 1))
349 11121 : orc = 1._dp/(rc*rc)
350 76665 : spmat(:) = pf/(2._dp*p(:) + orc)**(l + k + 1.5_dp)
351 :
352 11121 : 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 43686 : 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 43686 : n = SIZE(pa)
384 43686 : m = SIZE(pb)
385 :
386 43686 : CPASSERT(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
387 :
388 43686 : tld = gamma1(l + k + 1)*rc**(2*l + 3)*2._dp**(l + k + 0.5)
389 :
390 253586 : DO iq = 1, m
391 2737470 : DO ip = 1, n
392 2693784 : 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 43686 : 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 79 : 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 79 : n = SIZE(pa)
424 79 : m = SIZE(pb)
425 :
426 79 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
427 :
428 79 : tld = 0.5_dp/rc**(2*k)*gamma1(l + k + 1)
429 395 : DO iq = 1, m
430 1659 : DO ip = 1, n
431 1580 : gmat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + k + 1.5_dp)
432 : END DO
433 : END DO
434 :
435 79 : 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 192 : 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 192 : n = SIZE(pa)
600 192 : m = SIZE(pb)
601 :
602 192 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
603 :
604 192 : CALL sg_erf(umat, l, a, pa, pb)
605 :
606 192 : tld = 0.5_dp*fac(l)
607 1108 : DO iq = 1, m
608 12104 : DO ip = 1, n
609 11912 : umat(ip, iq) = tld/(pa(ip) + pb(iq))**(l + 1) - umat(ip, iq)
610 : END DO
611 : END DO
612 :
613 192 : 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 6780 : 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 6780 : n = SIZE(pa)
679 6780 : m = SIZE(pb)
680 :
681 6780 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
682 :
683 10776 : DO iq = 1, m
684 3996 : vq = vsto(2*nb(iq), pb(iq))
685 21892 : DO ip = 1, n
686 11116 : vp = vsto(2*na(ip), pa(ip))
687 11116 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
688 15112 : smat(ip, iq) = vpq/SQRT(vp*vq)
689 : END DO
690 : END DO
691 :
692 6780 : 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 6780 : 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 6780 : n = SIZE(pa)
723 6780 : m = SIZE(pb)
724 :
725 6780 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
726 :
727 10776 : DO iq = 1, m
728 3996 : vq = vsto(2*nb(iq), pb(iq))
729 3996 : wq = wsto(l, nb(iq), pb(iq))
730 21892 : DO ip = 1, n
731 11116 : vp = vsto(2*na(ip), pa(ip))
732 11116 : vpq = vsto(na(ip) + nb(iq), 0.5_dp*(pa(ip) + pb(iq)))
733 11116 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
734 11116 : vpq2 = vsto(na(ip) + nb(iq) - 2, 0.5_dp*(pa(ip) + pb(iq)))
735 11116 : wp = wsto(l, na(ip), pa(ip))
736 : kmat(ip, iq) = 0.5_dp*pa(ip)*pb(iq)/SQRT(vp*vq)* &
737 15112 : (vpq - (wp + wq)*vpq1 + wp*wq*vpq2)
738 : END DO
739 : END DO
740 :
741 6780 : 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 7728 : 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 7728 : n = SIZE(pa)
768 7728 : m = SIZE(pb)
769 :
770 7728 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
771 :
772 11976 : DO iq = 1, m
773 4248 : vq = vsto(2*nb(iq), pb(iq))
774 23344 : DO ip = 1, n
775 11368 : vp = vsto(2*na(ip), pa(ip))
776 11368 : vpq1 = vsto(na(ip) + nb(iq) - 1, 0.5_dp*(pa(ip) + pb(iq)))
777 15616 : umat(ip, iq) = 2._dp/SQRT(vp*vq)*vpq1
778 : END DO
779 : END DO
780 :
781 7728 : 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 101672 : 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 101672 : v = fac(n)/x**(n + 1)
841 :
842 101672 : END FUNCTION vsto
843 :
844 : ! **************************************************************************************************
845 : !> \brief ...
846 : !> \param n ...
847 : !> \param m ...
848 : !> \param x ...
849 : !> \return ...
850 : ! **************************************************************************************************
851 15112 : 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 15112 : w = 2._dp*REAL(m - n - 1, dp)/x
857 :
858 15112 : 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 0 : 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 0 : n = SIZE(ra, 2)
880 0 : m = SIZE(rb, 2)
881 :
882 0 : CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
883 :
884 0 : DO iq = 1, m
885 0 : DO ip = 1, n
886 0 : smat(ip, iq) = SUM(ra(:, ip)*rb(:, iq)*wr(:))
887 : END DO
888 : END DO
889 :
890 0 : 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 0 : 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 0 : n = SIZE(ra, 2)
918 0 : m = SIZE(rb, 2)
919 :
920 0 : CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
921 :
922 0 : DO iq = 1, m
923 0 : DO ip = 1, n
924 : kmat(ip, iq) = 0.5_dp*SUM(wr(:)*dra(:, ip)*drb(:, iq) &
925 0 : + wr(:)*REAL(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**2)
926 : END DO
927 : END DO
928 :
929 0 : 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 0 : 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 0 : n = SIZE(ra, 2)
953 0 : m = SIZE(rb, 2)
954 :
955 0 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
956 :
957 0 : DO iq = 1, m
958 0 : DO ip = 1, n
959 0 : umat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)/r(:))
960 : END DO
961 : END DO
962 :
963 0 : 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 0 : 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 0 : n = SIZE(ra, 2)
991 0 : m = SIZE(rb, 2)
992 :
993 0 : CPASSERT(.NOT. (n > SIZE(umat, 1) .OR. m > SIZE(umat, 2)))
994 :
995 0 : DO iq = 1, m
996 0 : DO ip = 1, n
997 : umat(ip, iq) = SUM(wr(:)*dra(:, ip)*drb(:, iq)/r(:) &
998 0 : + wr(:)*REAL(l*(l + 1), dp)*ra(:, ip)*rb(:, iq)/r(:)**3)
999 : END DO
1000 : END DO
1001 :
1002 0 : 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 0 : 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 0 : n = SIZE(ra, 2)
1028 0 : m = SIZE(rb, 2)
1029 :
1030 0 : CPASSERT(.NOT. (n > SIZE(upmat, 1) .OR. m > SIZE(upmat, 2)))
1031 :
1032 0 : DO iq = 1, m
1033 0 : DO ip = 1, n
1034 0 : upmat(ip, iq) = 0._dp
1035 0 : DO k = 1, SIZE(r)
1036 : upmat(ip, iq) = upmat(ip, iq) + &
1037 0 : (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 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pro
1075 :
1076 0 : n = SIZE(ra, 2)
1077 0 : CPASSERT(SIZE(spmat) >= n)
1078 :
1079 0 : ALLOCATE (pro(n))
1080 :
1081 0 : pf = SQRT(2._dp)/SQRT(gamma1(l + 2*k + 1))/rc**(l + 2*k + 1.5_dp)
1082 0 : pro(:) = pf*r(:)**(l + 2*k)*EXP(-0.5_dp*(r(:)/rc)**2)
1083 :
1084 0 : DO ip = 1, n
1085 0 : spmat(ip) = SUM(wr(:)*pro(:)*ra(:, ip))
1086 : END DO
1087 :
1088 0 : DEALLOCATE (pro)
1089 :
1090 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: op
1121 :
1122 0 : n = SIZE(ra, 2)
1123 0 : m = SIZE(rb, 2)
1124 :
1125 0 : CPASSERT(.NOT. (n > SIZE(vpmat, 1) .OR. m > SIZE(vpmat, 2)))
1126 :
1127 0 : ALLOCATE (op(n))
1128 :
1129 0 : op(:) = (r(:)/rc)**(2*k)*EXP(-0.5_dp*(r(:)/rc)**2)
1130 :
1131 0 : DO iq = 1, m
1132 0 : DO ip = 1, n
1133 0 : vpmat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1134 : END DO
1135 : END DO
1136 :
1137 0 : DEALLOCATE (op)
1138 :
1139 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: op
1167 :
1168 0 : n = SIZE(ra, 2)
1169 0 : m = SIZE(rb, 2)
1170 :
1171 0 : CPASSERT(.NOT. (n > SIZE(gmat, 1) .OR. m > SIZE(gmat, 2)))
1172 :
1173 0 : ALLOCATE (op(n))
1174 :
1175 0 : op(:) = (r(:)/rc)**(2*k)
1176 :
1177 0 : DO iq = 1, m
1178 0 : DO ip = 1, n
1179 0 : gmat(ip, iq) = SUM(wr(:)*ra(:, ip)*rb(:, iq)*op(:))
1180 : END DO
1181 : END DO
1182 :
1183 0 : DEALLOCATE (op)
1184 :
1185 0 : END SUBROUTINE num_conf
1186 :
1187 : END MODULE ai_onecenter
|
