(git:3f0d10b)
Loading...
Searching...
No Matches
ai_operators_r12.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of integrals over Cartesian Gaussian-type functions for different r12
10!> operators: 1/r12, erf(omega*r12/r12), erfc(omega*r12/r12), exp(-omega*r12^2)/r12 and
11!> exp(-omega*r12^2)
12!> \par Literature
13!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
14!> R. Ahlrichs, PCCP, 8, 3072 (2006)
15!> \par History
16!> 05.2019 Added the truncated Coulomb operator (A. Bussy)
17!> \par Parameters
18!> - ax,ay,az : Angular momentum index numbers of orbital a.
19!> - cx,cy,cz : Angular momentum index numbers of orbital c.
20!> - coset : Cartesian orbital set pointer.
21!> - dac : Distance between the atomic centers a and c.
22!> - l{a,c} : Angular momentum quantum number of shell a or c.
23!> - l{a,c}_max : Maximum angular momentum quantum number of shell a or c.
24!> - l{a,c}_min : Minimum angular momentum quantum number of shell a or c.
25!> - ncoset : Number of orbitals in a Cartesian orbital set.
26!> - npgf{a,c} : Degree of contraction of shell a or c.
27!> - rac : Distance vector between the atomic centers a and c.
28!> - rac2 : Square of the distance between the atomic centers a and c.
29!> - zet{a,c} : Exponents of the Gaussian-type functions a or c.
30!> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
31!> - zetw : Reciprocal of the sum of the exponents of orbital a and c.
32!> - omega : Parameter in the operator
33!> - r_cutoff : The cutoff radius for the truncated Coulomb operator
34!> \author Dorothea Golze (05.2016)
35! **************************************************************************************************
36MODULE ai_operators_r12
37
38 USE gamma, ONLY: fgamma => fgamma_0
39 USE kinds, ONLY: dp
40 USE mathconstants, ONLY: fac,&
41 pi
42 USE orbital_pointers, ONLY: coset,&
43 ncoset
44 USE t_c_g0, ONLY: get_lmax_init,&
45 t_c_g0_n
46#include "../base/base_uses.f90"
47
48 IMPLICIT NONE
49 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_operators_r12'
50 PRIVATE
51
52 ! *** Public subroutines ***
53
54 PUBLIC :: operator2, operator2_recurrence, cps_coulomb2, cps_verf2, cps_verfc2, cps_vgauss2, &
55 cps_gauss2, ab_sint_os, cps_truncated2
56
57 abstract INTERFACE
58! **************************************************************************************************
59!> \brief Interface for the calculation of integrals over s-functions and the s-type auxiliary
60!> integrals using the Obara-Saika (OS) scheme
61!> \param v ...
62!> \param nmax ...
63!> \param zetp ...
64!> \param zetq ...
65!> \param zetw ...
66!> \param rho ...
67!> \param rac2 ...
68!> \param omega ...
69!> \param r_cutoff ...
70! **************************************************************************************************
71 SUBROUTINE ab_sint_os(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
72 USE kinds, ONLY: dp
73 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
74 INTEGER, INTENT(IN) :: nmax
75 REAL(KIND=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
76 r_cutoff
77
78 END SUBROUTINE ab_sint_os
79 END INTERFACE
80
81CONTAINS
82
83! **************************************************************************************************
84!> \brief Calculation of the primitive two-center integrals over Cartesian Gaussian-type
85!> functions for different r12 operators.
86!> \param cps_operator2 procedure pointer for the respective operator. The integrals evaluation
87!> differs only in the evaluation of the cartesian primitive s (cps) integrals [s|O(r12)|s]
88!> and auxiliary integrals [s|O(r12)|s]^n. This pointer selects the correct routine.
89!> \param la_max ...
90!> \param npgfa ...
91!> \param zeta ...
92!> \param la_min ...
93!> \param lc_max ...
94!> \param npgfc ...
95!> \param zetc ...
96!> \param lc_min ...
97!> \param omega ...
98!> \param r_cutoff ...
99!> \param rac ...
100!> \param rac2 ...
101!> \param vac matrix storing the integrals
102!> \param v temporary work array
103!> \param maxder maximal derivative
104!> \param vac_plus matrix storing the integrals for highler l-quantum numbers; used to
105!> construct the derivatives
106! **************************************************************************************************
107
108 SUBROUTINE operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
109 omega, r_cutoff, rac, rac2, vac, v, maxder, vac_plus)
110 PROCEDURE(ab_sint_os), POINTER :: cps_operator2
111 INTEGER, INTENT(IN) :: la_max, npgfa
112 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
113 INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
114 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc
115 INTEGER, INTENT(IN) :: lc_min
116 REAL(kind=dp), INTENT(IN) :: omega, r_cutoff
117 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
118 REAL(kind=dp), INTENT(IN) :: rac2
119 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
120 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
121 INTEGER, INTENT(IN), OPTIONAL :: maxder
122 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
123
124 CHARACTER(len=*), PARAMETER :: routinen = 'operator2'
125
126 INTEGER :: ipgf, jpgf, na, nap, nc, ncp, nmax, handle
127 REAL(kind=dp) :: rho, zetp, zetq, zetw
128
129 CALL timeset(routinen, handle)
130
131 v = 0.0_dp
132
133 IF (PRESENT(maxder)) THEN
134 vac_plus = 0.0_dp
135 END IF
136
137 nmax = la_max + lc_max + 1
138
139 na = 0
140 nap = 0
141
142 DO ipgf = 1, npgfa
143
144 nc = 0
145 ncp = 0
146
147 DO jpgf = 1, npgfc
148
149 zetp = 1.0_dp/zeta(ipgf)
150 zetq = 1.0_dp/zetc(jpgf)
151 zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
152
153 rho = zeta(ipgf)*zetc(jpgf)*zetw
154
155 CALL cps_operator2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
156 CALL operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta(ipgf), zetc(jpgf), &
157 zetp, zetq, zetw, rho, rac, vac, v, na, nc, nap, ncp, &
158 maxder=maxder, vac_plus=vac_plus)
159 END DO
160
161 IF (PRESENT(maxder)) THEN
162 na = na + ncoset(la_max - maxder)
163 ELSE
164 na = na + ncoset(la_max)
165 END IF
166 nap = nap + ncoset(la_max)
167 END DO
168
169 CALL timestop(handle)
170
171 END SUBROUTINE operator2
172
173! **************************************************************************************************
174!> \brief Apply the common two-center OS recurrence to one primitive pair.
175!> The caller initializes v(1,1,:) with the operator-specific [s||s] values.
176!> \param la_max ...
177!> \param la_min ...
178!> \param lc_max ...
179!> \param lc_min ...
180!> \param zeta_a ...
181!> \param zeta_c ...
182!> \param zetp ...
183!> \param zetq ...
184!> \param zetw ...
185!> \param rho ...
186!> \param rac ...
187!> \param vac ...
188!> \param v ...
189!> \param na ...
190!> \param nc ...
191!> \param nap ...
192!> \param ncp ...
193!> \param maxder ...
194!> \param vac_plus ...
195! **************************************************************************************************
196 SUBROUTINE operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta_a, zeta_c, zetp, zetq, &
197 zetw, rho, rac, vac, v, na, nc, nap, ncp, maxder, vac_plus)
198 INTEGER, INTENT(IN) :: la_max, la_min, lc_max, lc_min
199 REAL(kind=dp), INTENT(IN) :: zeta_a, zeta_c, zetp, zetq, zetw, rho
200 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
201 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
202 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
203 INTEGER, INTENT(IN) :: na
204 INTEGER, INTENT(INOUT) :: nc
205 INTEGER, INTENT(IN) :: nap
206 INTEGER, INTENT(INOUT) :: ncp
207 INTEGER, INTENT(IN), OPTIONAL :: maxder
208 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
209
210 INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
211 cy, cz, i, j, la, lc, maxder_local, n, &
212 nmax
213 REAL(kind=dp) :: f1, f2, f3, f4, f5, f6, fcx, fcy, fcz
214 REAL(kind=dp), DIMENSION(3) :: raw, rcw
215
216 maxder_local = 0
217 IF (PRESENT(maxder)) maxder_local = maxder
218 nmax = la_max + lc_max + 1
219 ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
220
221 IF (lc_max > 0) THEN
222
223 f1 = 0.5_dp*zetq
224 f2 = -rho*zetq
225
226 rcw(:) = -zeta_a*zetw*rac(:)
227
228 ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
229
230 DO n = 1, nmax - 1
231 v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
232 v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
233 v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
234 END DO
235
236 ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
237 ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
238 ! ** f2*[s||c-2i]{n+1} ***
239
240 DO lc = 2, lc_max
241
242 DO n = 1, nmax - lc
243
244 v(1, coset(0, 0, lc), n) = &
245 rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
246 f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
247 f2*v(1, coset(0, 0, lc - 2), n + 1))
248
249 cz = lc - 1
250 v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
251
252 DO cy = 2, lc
253 cz = lc - cy
254 v(1, coset(0, cy, cz), n) = &
255 rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
256 f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
257 f2*v(1, coset(0, cy - 2, cz), n + 1))
258 END DO
259
260 DO cy = 0, lc - 1
261 cz = lc - 1 - cy
262 v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
263 END DO
264
265 DO cx = 2, lc
266 f6 = f1*real(cx - 1, dp)
267 DO cy = 0, lc - cx
268 cz = lc - cx - cy
269 v(1, coset(cx, cy, cz), n) = &
270 rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
271 f6*(v(1, coset(cx - 2, cy, cz), n) + &
272 f2*v(1, coset(cx - 2, cy, cz), n + 1))
273 END DO
274 END DO
275
276 END DO
277
278 END DO
279
280 END IF
281
282 ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
283
284 IF (la_max > 0) THEN
285
286 f3 = 0.5_dp*zetp
287 f4 = -rho*zetp
288 f5 = 0.5_dp*zetw
289
290 raw(:) = zeta_c*zetw*rac(:)
291
292 ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
293
294 DO n = 1, nmax - 1
295 v(2, 1, n) = raw(1)*v(1, 1, n + 1)
296 v(3, 1, n) = raw(2)*v(1, 1, n + 1)
297 v(4, 1, n) = raw(3)*v(1, 1, n + 1)
298 END DO
299
300 ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
301 ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
302 ! *** f4*[a-2i||s]{n+1}) ***
303
304 DO la = 2, la_max
305
306 DO n = 1, nmax - la
307
308 v(coset(0, 0, la), 1, n) = &
309 raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
310 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
311 f4*v(coset(0, 0, la - 2), 1, n + 1))
312
313 az = la - 1
314 v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
315
316 DO ay = 2, la
317 az = la - ay
318 v(coset(0, ay, az), 1, n) = &
319 raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
320 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
321 f4*v(coset(0, ay - 2, az), 1, n + 1))
322 END DO
323
324 DO ay = 0, la - 1
325 az = la - 1 - ay
326 v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
327 END DO
328
329 DO ax = 2, la
330 f6 = f3*real(ax - 1, dp)
331 DO ay = 0, la - ax
332 az = la - ax - ay
333 v(coset(ax, ay, az), 1, n) = &
334 raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
335 f6*(v(coset(ax - 2, ay, az), 1, n) + &
336 f4*v(coset(ax - 2, ay, az), 1, n + 1))
337 END DO
338 END DO
339
340 END DO
341
342 END DO
343
344 DO lc = 1, lc_max
345
346 DO cx = 0, lc
347 DO cy = 0, lc - cx
348 cz = lc - cx - cy
349
350 coc = coset(cx, cy, cz)
351 cocx = coset(max(0, cx - 1), cy, cz)
352 cocy = coset(cx, max(0, cy - 1), cz)
353 cocz = coset(cx, cy, max(0, cz - 1))
354
355 fcx = f5*real(cx, dp)
356 fcy = f5*real(cy, dp)
357 fcz = f5*real(cz, dp)
358
359 ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
360 ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
361
362 DO n = 1, nmax - 1 - lc
363 v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
364 v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
365 v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
366 END DO
367
368 ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
369 ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
370 ! *** f4*[a-2i||c]{n+1}) + ***
371 ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
372
373 DO la = 2, la_max
374
375 DO n = 1, nmax - la - lc
376
377 v(coset(0, 0, la), coc, n) = &
378 raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
379 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
380 f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
381 fcz*v(coset(0, 0, la - 1), cocz, n + 1)
382
383 az = la - 1
384 v(coset(0, 1, az), coc, n) = &
385 raw(2)*v(coset(0, 0, az), coc, n + 1) + &
386 fcy*v(coset(0, 0, az), cocy, n + 1)
387
388 DO ay = 2, la
389 az = la - ay
390 v(coset(0, ay, az), coc, n) = &
391 raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
392 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
393 f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
394 fcy*v(coset(0, ay - 1, az), cocy, n + 1)
395 END DO
396
397 DO ay = 0, la - 1
398 az = la - 1 - ay
399 v(coset(1, ay, az), coc, n) = &
400 raw(1)*v(coset(0, ay, az), coc, n + 1) + &
401 fcx*v(coset(0, ay, az), cocx, n + 1)
402 END DO
403
404 DO ax = 2, la
405 f6 = f3*real(ax - 1, dp)
406 DO ay = 0, la - ax
407 az = la - ax - ay
408 v(coset(ax, ay, az), coc, n) = &
409 raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
410 f6*(v(coset(ax - 2, ay, az), coc, n) + &
411 f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
412 fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
413 END DO
414 END DO
415
416 END DO
417
418 END DO
419
420 END DO
421 END DO
422
423 END DO
424
425 END IF
426
427 DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
428 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
429 vac(na + i, nc + j) = v(i, j, 1)
430 END DO
431 END DO
432
433 IF (PRESENT(maxder)) THEN
434 DO j = 1, ncoset(lc_max)
435 DO i = 1, ncoset(la_max)
436 vac_plus(nap + i, ncp + j) = v(i, j, 1)
437 END DO
438 END DO
439 END IF
440
441 nc = nc + ncoset(lc_max - maxder_local)
442 ncp = ncp + ncoset(lc_max)
443 END SUBROUTINE operator2_recurrence
444
445! **************************************************************************************************
446!> \brief Calculation of Coulomb integrals for s-function, i.e, [s|1/r12|s], and the auxiliary
447!> integrals [s|1/r12|s]^n
448!> \param v matrix storing the integrals
449!> \param nmax maximal n in the auxiliary integrals [s|1/r12|s]^n
450!> \param zetp = 1/zeta
451!> \param zetq = 1/zetc
452!> \param zetw = 1/(zeta+zetc)
453!> \param rho = zeta*zetc*zetw
454!> \param rac2 square distance between center A and C, |Ra-Rc|^2
455!> \param omega this parameter is actually not used, but included for the sake of the abstract
456!> interface
457!> \param r_cutoff same as above
458! **************************************************************************************************
459 SUBROUTINE cps_coulomb2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
460 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
461 INTEGER, INTENT(IN) :: nmax
462 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
463 r_cutoff
464
465 INTEGER :: n
466 REAL(kind=dp) :: f0, t
467 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
468
469 mark_used(omega)
470 mark_used(r_cutoff)
471
472 ALLOCATE (f(0:nmax))
473 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
474
475 ! *** Calculate the incomplete Gamma/Boys function ***
476
477 t = rho*rac2
478 CALL fgamma(nmax - 1, t, f)
479
480 ! *** Calculate the basic two-center integrals [s||s]{n} ***
481
482 DO n = 1, nmax
483 v(1, 1, n) = f0*f(n - 1)
484 END DO
485
486 DEALLOCATE (f)
487 END SUBROUTINE cps_coulomb2
488
489! **************************************************************************************************
490!> \brief Calculation of verf integrals for s-function, i.e, [s|erf(omega*r12)/r12|s], and the
491!> auxiliary integrals [s|erf(omega*r12)/r12|s]^n
492!> \param v matrix storing the integrals
493!> \param nmax maximal n in the auxiliary integrals [s|erf(omega*r12)/r12|s]^n
494!> \param zetp = 1/zeta
495!> \param zetq = 1/zetc
496!> \param zetw = 1/(zeta+zetc)
497!> \param rho = zeta*zetc*zetw
498!> \param rac2 square distance between center A and C, |Ra-Rc|^2
499!> \param omega parameter in the operator
500!> \param r_cutoff dummy argument for the sake of generality
501! **************************************************************************************************
502 SUBROUTINE cps_verf2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
503 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
504 INTEGER, INTENT(IN) :: nmax
505 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
506 r_cutoff
507
508 INTEGER :: n
509 REAL(kind=dp) :: arg, comega, f0, t
510 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
511
512 mark_used(r_cutoff)
513
514 ALLOCATE (f(0:nmax))
515 comega = omega**2/(omega**2 + rho)
516 f0 = 2.0_dp*sqrt(pi**5*zetw*comega)*zetp*zetq
517
518 ! *** Calculate the incomplete Gamma/Boys function ***
519
520 t = rho*rac2
521 arg = comega*t
522 CALL fgamma(nmax - 1, arg, f)
523
524 ! *** Calculate the basic two-center integrals [s||s]{n} ***
525
526 DO n = 1, nmax
527 v(1, 1, n) = f0*f(n - 1)*comega**(n - 1)
528 END DO
529
530 DEALLOCATE (f)
531
532 END SUBROUTINE cps_verf2
533
534! **************************************************************************************************
535!> \brief Calculation of verfc integrals for s-function, i.e, [s|erfc(omega*r12)/r12|s], and
536!> the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
537!> \param v matrix storing the integrals
538!> \param nmax maximal n in the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
539!> \param zetp = 1/zeta
540!> \param zetq = 1/zetc
541!> \param zetw = 1/(zeta+zetc)
542!> \param rho = zeta*zetc*zetw
543!> \param rac2 square distance between center A and C, |Ra-Rc|^2
544!> \param omega parameter in the operator
545!> \param r_cutoff dummy argument for the sake of generality
546! **************************************************************************************************
547 SUBROUTINE cps_verfc2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
548 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
549 INTEGER, INTENT(IN) :: nmax
550 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
551 r_cutoff
552
553 INTEGER :: n
554 REAL(kind=dp) :: argerf, comega, f0, t
555 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
556
557 mark_used(r_cutoff)
558
559 ALLOCATE (fv(0:nmax), fverf(0:nmax))
560 comega = omega**2/(omega**2 + rho)
561 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
562
563 ! *** Calculate the incomplete Gamma/Boys function ***
564
565 t = rho*rac2
566 argerf = comega*t
567
568 CALL fgamma(nmax - 1, t, fv)
569 CALL fgamma(nmax - 1, argerf, fverf)
570
571 ! *** Calculate the basic two-center integrals [s||s]{n} ***
572
573 DO n = 1, nmax
574 v(1, 1, n) = f0*(fv(n - 1) - sqrt(comega)*comega**(n - 1)*fverf(n - 1))
575 END DO
576
577 DEALLOCATE (fv, fverf)
578
579 END SUBROUTINE cps_verfc2
580
581! **************************************************************************************************
582!> \brief Calculation of vgauss integrals for s-function, i.e, [s|exp(-omega*r12^2)/r12|s], and
583!> the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
584!> \param v matrix storing the integrals
585!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
586!> \param zetp = 1/zeta
587!> \param zetq = 1/zetc
588!> \param zetw = 1/(zeta+zetc)
589!> \param rho = zeta*zetc*zetw
590!> \param rac2 square distance between center A and C, |Ra-Rc|^2
591!> \param omega parameter in the operator
592!> \param r_cutoff dummy argument for the sake of generality
593! **************************************************************************************************
594 SUBROUTINE cps_vgauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
595 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
596 INTEGER, INTENT(IN) :: nmax
597 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
598 r_cutoff
599
600 INTEGER :: j, n
601 REAL(kind=dp) :: arg, dummy, eta, expt, f0, fsign, t, tau
602 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
603
604 mark_used(r_cutoff)
605
606 ALLOCATE (f(0:nmax))
607
608 dummy = zetp
609 dummy = zetq
610 eta = rho/(rho + omega)
611 tau = omega/(rho + omega)
612
613 ! *** Calculate the incomplete Gamma/Boys function ***
614
615 t = rho*rac2
616 arg = eta*t
617
618 CALL fgamma(nmax - 1, arg, f)
619
620 expt = exp(-omega/(omega + rho)*t)
621 f0 = 2.0_dp*sqrt(pi**5*zetw**3)/(rho + omega)*expt
622
623 ! *** Calculate the basic two-center integrals [s||s]{n} ***
624 v(1, 1, 1:nmax) = 0.0_dp
625 DO n = 1, nmax
626 fsign = (-1.0_dp)**(n - 1)
627 DO j = 0, n - 1
628 v(1, 1, n) = v(1, 1, n) + f0*fsign* &
629 fac(n - 1)/fac(n - j - 1)/fac(j)*(-tau)**(n - j - 1)*(-eta)**j*f(j)
630 END DO
631 END DO
632
633 DEALLOCATE (f)
634
635 END SUBROUTINE cps_vgauss2
636
637! **************************************************************************************************
638!> \brief Calculation of gauss integrals for s-function, i.e, [s|exp(-omega*r12^2)|s], and
639!> the auxiliary integrals [s|exp(-omega*r12^2)|s]
640!> \param v matrix storing the integrals
641!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)|s]
642!> \param zetp = 1/zeta
643!> \param zetq = 1/zetc
644!> \param zetw = 1/(zeta+zetc)
645!> \param rho = zeta*zetc*zetw
646!> \param rac2 square distance between center A and C, |Ra-Rc|^2
647!> \param omega parameter in the operator
648!> \param r_cutoff dummy argument for the sake of generality
649! **************************************************************************************************
650 SUBROUTINE cps_gauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
651 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
652 INTEGER, INTENT(IN) :: nmax
653 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
654 r_cutoff
655
656 INTEGER :: n
657 REAL(kind=dp) :: dummy, expt, f0, t, tau
658 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
659
660 mark_used(r_cutoff)
661
662 ALLOCATE (f(0:nmax))
663
664 dummy = zetp
665 dummy = zetq
666 tau = omega/(rho + omega)
667 t = rho*rac2
668 expt = exp(-tau*t)
669 f0 = pi**3*sqrt(zetw**3/(rho + omega)**3)*expt
670
671 ! *** Calculate the basic two-center integrals [s||s]{n} ***
672
673 DO n = 1, nmax
674 v(1, 1, n) = f0*tau**(n - 1)
675 END DO
676
677 DEALLOCATE (f)
678
679 END SUBROUTINE cps_gauss2
680
681! **************************************************************************************************
682!> \brief Calculation of truncated Coulomb integrals for s-function, i.e, [s|TC|s] where TC = 1/r12
683!> if r12 <= r_cutoff and 0 otherwise
684!> \param v matrix storing the integrals
685!> \param nmax maximal n in the auxiliary integrals [s|TC|s]
686!> \param zetp = 1/zeta
687!> \param zetq = 1/zetc
688!> \param zetw = 1/(zeta+zetc)
689!> \param rho = zeta*zetc*zetw
690!> \param rac2 square distance between center A and C, |Ra-Rc|^2
691!> \param omega dummy argument for the sake of generality
692!> \param r_cutoff the radius at which the operator is cut
693!> \note The truncated operator must have been initialized from the data file prior to this call
694! **************************************************************************************************
695 SUBROUTINE cps_truncated2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
696 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
697 INTEGER, INTENT(IN) :: nmax
698 REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
699 r_cutoff
700
701 INTEGER :: n
702 LOGICAL :: use_gamma
703 REAL(kind=dp) :: f0, r, t
704 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
705
706 mark_used(omega)
707
708 ALLOCATE (f(nmax + 1)) !t_c_g0 needs to start at index 1
709
710 r = r_cutoff*sqrt(rho)
711 t = rho*rac2
712 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
713
714 !check that the operator has been init from file
715 cpassert(get_lmax_init() >= nmax)
716
717 CALL t_c_g0_n(f, use_gamma, r, t, nmax)
718 IF (use_gamma) CALL fgamma(nmax, t, f)
719
720 DO n = 1, nmax
721 v(1, 1, n) = f0*f(n)
722 END DO
723
724 DEALLOCATE (f)
725
726 END SUBROUTINE cps_truncated2
727
728END MODULE ai_operators_r12
ai_operators_r12::ab_sint_os
Interface for the calculation of integrals over s-functions and the s-type auxiliary integrals using ...
Definition ai_operators_r12.F:71
ai_operators_r12
Calculation of integrals over Cartesian Gaussian-type functions for different r12 operators: 1/r12,...
Definition ai_operators_r12.F:36
ai_operators_r12::cps_verfc2
subroutine, public cps_verfc2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of verfc integrals for s-function, i.e, [s|erfc(omega*r12)/r12|s], and the auxiliary inte...
Definition ai_operators_r12.F:548
ai_operators_r12::cps_vgauss2
subroutine, public cps_vgauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of vgauss integrals for s-function, i.e, [s|exp(-omega*r12^2)/r12|s], and the auxiliary i...
Definition ai_operators_r12.F:595
ai_operators_r12::cps_gauss2
subroutine, public cps_gauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of gauss integrals for s-function, i.e, [s|exp(-omega*r12^2)|s], and the auxiliary integr...
Definition ai_operators_r12.F:651
ai_operators_r12::cps_coulomb2
subroutine, public cps_coulomb2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of Coulomb integrals for s-function, i.e, [s|1/r12|s], and the auxiliary integrals [s|1/r...
Definition ai_operators_r12.F:460
ai_operators_r12::operator2_recurrence
subroutine, public operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta_a, zeta_c, zetp, zetq, zetw, rho, rac, vac, v, na, nc, nap, ncp, maxder, vac_plus)
Apply the common two-center OS recurrence to one primitive pair. The caller initializes v(1,...
Definition ai_operators_r12.F:198
ai_operators_r12::cps_truncated2
subroutine, public cps_truncated2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of truncated Coulomb integrals for s-function, i.e, [s|TC|s] where TC = 1/r12 if r12 <= r...
Definition ai_operators_r12.F:696
ai_operators_r12::cps_verf2
subroutine, public cps_verf2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
Calculation of verf integrals for s-function, i.e, [s|erf(omega*r12)/r12|s], and the auxiliary integr...
Definition ai_operators_r12.F:503
ai_operators_r12::operator2
subroutine, public operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, omega, r_cutoff, rac, rac2, vac, v, maxder, vac_plus)
Calculation of the primitive two-center integrals over Cartesian Gaussian-type functions for differen...
Definition ai_operators_r12.F:110
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:58
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(res, use_gamma, r, t, nderiv)
...
Definition t_c_g0.F:88
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1468