(git:b279b6b)
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-2024 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 ! **************************************************************************************************
36 MODULE 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, cps_coulomb2, cps_verf2, cps_verfc2, cps_vgauss2, cps_gauss2, ab_sint_os, &
55  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 
81 CONTAINS
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 :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
127  cy, cz, i, ipgf, j, jpgf, la, lc, &
128  maxder_local, n, na, nap, nc, ncp, &
129  nmax, handle
130  REAL(kind=dp) :: dac, f1, f2, f3, f4, f5, f6, fcx, &
131  fcy, fcz, rho, zetp, zetq, zetw
132  REAL(kind=dp), DIMENSION(3) :: raw, rcw
133 
134  CALL timeset(routinen, handle)
135 
136  v = 0.0_dp
137 
138  maxder_local = 0
139  IF (PRESENT(maxder)) THEN
140  maxder_local = maxder
141  vac_plus = 0.0_dp
142  END IF
143 
144  nmax = la_max + lc_max + 1
145 
146  ! *** Calculate the distance of the centers a and c ***
147 
148  dac = sqrt(rac2)
149 
150  ! *** Loop over all pairs of primitive Gaussian-type functions ***
151 
152  na = 0
153  nap = 0
154 
155  DO ipgf = 1, npgfa
156 
157  nc = 0
158  ncp = 0
159 
160  DO jpgf = 1, npgfc
161 
162  ! *** Calculate some prefactors ***
163 
164  zetp = 1.0_dp/zeta(ipgf)
165  zetq = 1.0_dp/zetc(jpgf)
166  zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
167 
168  rho = zeta(ipgf)*zetc(jpgf)*zetw
169 
170  ! *** Calculate the basic two-center integrals [s||s]{n} ***
171 
172  CALL cps_operator2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
173 
174  ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
175 
176  IF (lc_max > 0) THEN
177 
178  f1 = 0.5_dp*zetq
179  f2 = -rho*zetq
180 
181  rcw(:) = -zeta(ipgf)*zetw*rac(:)
182 
183  ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
184 
185  DO n = 1, nmax - 1
186  v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
187  v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
188  v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
189  END DO
190 
191  ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
192  ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
193  ! ** f2*[s||c-2i]{n+1} ***
194 
195  DO lc = 2, lc_max
196 
197  DO n = 1, nmax - lc
198 
199  ! **** Increase the angular momentum component z of c ***
200 
201  v(1, coset(0, 0, lc), n) = &
202  rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
203  f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
204  f2*v(1, coset(0, 0, lc - 2), n + 1))
205 
206  ! *** Increase the angular momentum component y of c ***
207 
208  cz = lc - 1
209  v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
210 
211  DO cy = 2, lc
212  cz = lc - cy
213  v(1, coset(0, cy, cz), n) = &
214  rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
215  f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
216  f2*v(1, coset(0, cy - 2, cz), n + 1))
217  END DO
218 
219  ! *** Increase the angular momentum component x of c ***
220 
221  DO cy = 0, lc - 1
222  cz = lc - 1 - cy
223  v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
224  END DO
225 
226  DO cx = 2, lc
227  f6 = f1*real(cx - 1, dp)
228  DO cy = 0, lc - cx
229  cz = lc - cx - cy
230  v(1, coset(cx, cy, cz), n) = &
231  rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
232  f6*(v(1, coset(cx - 2, cy, cz), n) + &
233  f2*v(1, coset(cx - 2, cy, cz), n + 1))
234  END DO
235  END DO
236 
237  END DO
238 
239  END DO
240 
241  END IF
242 
243  ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
244 
245  IF (la_max > 0) THEN
246 
247  f3 = 0.5_dp*zetp
248  f4 = -rho*zetp
249  f5 = 0.5_dp*zetw
250 
251  raw(:) = zetc(jpgf)*zetw*rac(:)
252 
253  ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
254 
255  DO n = 1, nmax - 1
256  v(2, 1, n) = raw(1)*v(1, 1, n + 1)
257  v(3, 1, n) = raw(2)*v(1, 1, n + 1)
258  v(4, 1, n) = raw(3)*v(1, 1, n + 1)
259  END DO
260 
261  ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
262  ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
263  ! *** f4*[a-2i||s]{n+1}) ***
264 
265  DO la = 2, la_max
266 
267  DO n = 1, nmax - la
268 
269  ! *** Increase the angular momentum component z of a ***
270 
271  v(coset(0, 0, la), 1, n) = &
272  raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
273  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
274  f4*v(coset(0, 0, la - 2), 1, n + 1))
275 
276  ! *** Increase the angular momentum component y of a ***
277 
278  az = la - 1
279  v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
280 
281  DO ay = 2, la
282  az = la - ay
283  v(coset(0, ay, az), 1, n) = &
284  raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
285  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
286  f4*v(coset(0, ay - 2, az), 1, n + 1))
287  END DO
288 
289  ! *** Increase the angular momentum component x of a ***
290 
291  DO ay = 0, la - 1
292  az = la - 1 - ay
293  v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
294  END DO
295 
296  DO ax = 2, la
297  f6 = f3*real(ax - 1, dp)
298  DO ay = 0, la - ax
299  az = la - ax - ay
300  v(coset(ax, ay, az), 1, n) = &
301  raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
302  f6*(v(coset(ax - 2, ay, az), 1, n) + &
303  f4*v(coset(ax - 2, ay, az), 1, n + 1))
304  END DO
305  END DO
306 
307  END DO
308 
309  END DO
310 
311  DO lc = 1, lc_max
312 
313  DO cx = 0, lc
314  DO cy = 0, lc - cx
315  cz = lc - cx - cy
316 
317  coc = coset(cx, cy, cz)
318  cocx = coset(max(0, cx - 1), cy, cz)
319  cocy = coset(cx, max(0, cy - 1), cz)
320  cocz = coset(cx, cy, max(0, cz - 1))
321 
322  fcx = f5*real(cx, dp)
323  fcy = f5*real(cy, dp)
324  fcz = f5*real(cz, dp)
325 
326  ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
327  ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
328 
329  DO n = 1, nmax - 1 - lc
330  v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
331  v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
332  v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
333  END DO
334 
335  ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
336  ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
337  ! *** f4*[a-2i||c]{n+1}) + ***
338  ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
339 
340  DO la = 2, la_max
341 
342  DO n = 1, nmax - la - lc
343 
344  ! *** Increase the angular momentum component z of a ***
345 
346  v(coset(0, 0, la), coc, n) = &
347  raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
348  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
349  f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
350  fcz*v(coset(0, 0, la - 1), cocz, n + 1)
351 
352  ! *** Increase the angular momentum component y of a ***
353 
354  az = la - 1
355  v(coset(0, 1, az), coc, n) = &
356  raw(2)*v(coset(0, 0, az), coc, n + 1) + &
357  fcy*v(coset(0, 0, az), cocy, n + 1)
358 
359  DO ay = 2, la
360  az = la - ay
361  v(coset(0, ay, az), coc, n) = &
362  raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
363  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
364  f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
365  fcy*v(coset(0, ay - 1, az), cocy, n + 1)
366  END DO
367 
368  ! *** Increase the angular momentum component x of a ***
369 
370  DO ay = 0, la - 1
371  az = la - 1 - ay
372  v(coset(1, ay, az), coc, n) = &
373  raw(1)*v(coset(0, ay, az), coc, n + 1) + &
374  fcx*v(coset(0, ay, az), cocx, n + 1)
375  END DO
376 
377  DO ax = 2, la
378  f6 = f3*real(ax - 1, dp)
379  DO ay = 0, la - ax
380  az = la - ax - ay
381  v(coset(ax, ay, az), coc, n) = &
382  raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
383  f6*(v(coset(ax - 2, ay, az), coc, n) + &
384  f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
385  fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
386  END DO
387  END DO
388 
389  END DO
390 
391  END DO
392 
393  END DO
394  END DO
395 
396  END DO
397 
398  END IF
399 
400  DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
401  DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
402  vac(na + i, nc + j) = v(i, j, 1)
403  END DO
404  END DO
405 
406  IF (PRESENT(maxder)) THEN
407  DO j = 1, ncoset(lc_max)
408  DO i = 1, ncoset(la_max)
409  vac_plus(nap + i, ncp + j) = v(i, j, 1)
410  END DO
411  END DO
412  END IF
413 
414  nc = nc + ncoset(lc_max - maxder_local)
415  ncp = ncp + ncoset(lc_max)
416 
417  END DO
418 
419  na = na + ncoset(la_max - maxder_local)
420  nap = nap + ncoset(la_max)
421 
422  END DO
423 
424  CALL timestop(handle)
425 
426  END SUBROUTINE operator2
427 
428 ! **************************************************************************************************
429 !> \brief Calculation of Coulomb integrals for s-function, i.e, [s|1/r12|s], and the auxiliary
430 !> integrals [s|1/r12|s]^n
431 !> \param v matrix storing the integrals
432 !> \param nmax maximal n in the auxiliary integrals [s|1/r12|s]^n
433 !> \param zetp = 1/zeta
434 !> \param zetq = 1/zetc
435 !> \param zetw = 1/(zeta+zetc)
436 !> \param rho = zeta*zetc*zetw
437 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
438 !> \param omega this parameter is actually not used, but included for the sake of the abstract
439 !> interface
440 !> \param r_cutoff same as above
441 ! **************************************************************************************************
442  SUBROUTINE cps_coulomb2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
443  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
444  INTEGER, INTENT(IN) :: nmax
445  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
446  r_cutoff
447 
448  INTEGER :: n
449  REAL(kind=dp) :: f0, t
450  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
451 
452  mark_used(omega)
453  mark_used(r_cutoff)
454 
455  ALLOCATE (f(0:nmax))
456  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
457 
458  ! *** Calculate the incomplete Gamma/Boys function ***
459 
460  t = rho*rac2
461  CALL fgamma(nmax - 1, t, f)
462 
463  ! *** Calculate the basic two-center integrals [s||s]{n} ***
464 
465  DO n = 1, nmax
466  v(1, 1, n) = f0*f(n - 1)
467  END DO
468 
469  DEALLOCATE (f)
470  END SUBROUTINE cps_coulomb2
471 
472 ! **************************************************************************************************
473 !> \brief Calculation of verf integrals for s-function, i.e, [s|erf(omega*r12)/r12|s], and the
474 !> auxiliary integrals [s|erf(omega*r12)/r12|s]^n
475 !> \param v matrix storing the integrals
476 !> \param nmax maximal n in the auxiliary integrals [s|erf(omega*r12)/r12|s]^n
477 !> \param zetp = 1/zeta
478 !> \param zetq = 1/zetc
479 !> \param zetw = 1/(zeta+zetc)
480 !> \param rho = zeta*zetc*zetw
481 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
482 !> \param omega parameter in the operator
483 !> \param r_cutoff dummy argument for the sake of generality
484 ! **************************************************************************************************
485  SUBROUTINE cps_verf2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
486  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
487  INTEGER, INTENT(IN) :: nmax
488  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
489  r_cutoff
490 
491  INTEGER :: n
492  REAL(kind=dp) :: arg, comega, f0, t
493  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
494 
495  mark_used(r_cutoff)
496 
497  ALLOCATE (f(0:nmax))
498  comega = omega**2/(omega**2 + rho)
499  f0 = 2.0_dp*sqrt(pi**5*zetw*comega)*zetp*zetq
500 
501  ! *** Calculate the incomplete Gamma/Boys function ***
502 
503  t = rho*rac2
504  arg = comega*t
505  CALL fgamma(nmax - 1, arg, f)
506 
507  ! *** Calculate the basic two-center integrals [s||s]{n} ***
508 
509  DO n = 1, nmax
510  v(1, 1, n) = f0*f(n - 1)*comega**(n - 1)
511  END DO
512 
513  DEALLOCATE (f)
514 
515  END SUBROUTINE cps_verf2
516 
517 ! **************************************************************************************************
518 !> \brief Calculation of verfc integrals for s-function, i.e, [s|erfc(omega*r12)/r12|s], and
519 !> the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
520 !> \param v matrix storing the integrals
521 !> \param nmax maximal n in the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
522 !> \param zetp = 1/zeta
523 !> \param zetq = 1/zetc
524 !> \param zetw = 1/(zeta+zetc)
525 !> \param rho = zeta*zetc*zetw
526 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
527 !> \param omega parameter in the operator
528 !> \param r_cutoff dummy argument for the sake of generality
529 ! **************************************************************************************************
530  SUBROUTINE cps_verfc2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
531  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
532  INTEGER, INTENT(IN) :: nmax
533  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
534  r_cutoff
535 
536  INTEGER :: n
537  REAL(kind=dp) :: argerf, comega, f0, t
538  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
539 
540  mark_used(r_cutoff)
541 
542  ALLOCATE (fv(0:nmax), fverf(0:nmax))
543  comega = omega**2/(omega**2 + rho)
544  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
545 
546  ! *** Calculate the incomplete Gamma/Boys function ***
547 
548  t = rho*rac2
549  argerf = comega*t
550 
551  CALL fgamma(nmax - 1, t, fv)
552  CALL fgamma(nmax - 1, argerf, fverf)
553 
554  ! *** Calculate the basic two-center integrals [s||s]{n} ***
555 
556  DO n = 1, nmax
557  v(1, 1, n) = f0*(fv(n - 1) - sqrt(comega)*comega**(n - 1)*fverf(n - 1))
558  END DO
559 
560  DEALLOCATE (fv, fverf)
561 
562  END SUBROUTINE cps_verfc2
563 
564 ! **************************************************************************************************
565 !> \brief Calculation of vgauss integrals for s-function, i.e, [s|exp(-omega*r12^2)/r12|s], and
566 !> the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
567 !> \param v matrix storing the integrals
568 !> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
569 !> \param zetp = 1/zeta
570 !> \param zetq = 1/zetc
571 !> \param zetw = 1/(zeta+zetc)
572 !> \param rho = zeta*zetc*zetw
573 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
574 !> \param omega parameter in the operator
575 !> \param r_cutoff dummy argument for the sake of generality
576 ! **************************************************************************************************
577  SUBROUTINE cps_vgauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
578  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
579  INTEGER, INTENT(IN) :: nmax
580  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
581  r_cutoff
582 
583  INTEGER :: j, n
584  REAL(kind=dp) :: arg, dummy, eta, expt, f0, fsign, t, tau
585  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
586 
587  mark_used(r_cutoff)
588 
589  ALLOCATE (f(0:nmax))
590 
591  dummy = zetp
592  dummy = zetq
593  eta = rho/(rho + omega)
594  tau = omega/(rho + omega)
595 
596  ! *** Calculate the incomplete Gamma/Boys function ***
597 
598  t = rho*rac2
599  arg = eta*t
600 
601  CALL fgamma(nmax - 1, arg, f)
602 
603  expt = exp(-omega/(omega + rho)*t)
604  f0 = 2.0_dp*sqrt(pi**5*zetw**3)/(rho + omega)*expt
605 
606  ! *** Calculate the basic two-center integrals [s||s]{n} ***
607  v(1, 1, 1:nmax) = 0.0_dp
608  DO n = 1, nmax
609  fsign = (-1.0_dp)**(n - 1)
610  DO j = 0, n - 1
611  v(1, 1, n) = v(1, 1, n) + f0*fsign* &
612  fac(n - 1)/fac(n - j - 1)/fac(j)*(-tau)**(n - j - 1)*(-eta)**j*f(j)
613  END DO
614  END DO
615 
616  DEALLOCATE (f)
617 
618  END SUBROUTINE cps_vgauss2
619 
620 ! **************************************************************************************************
621 !> \brief Calculation of gauss integrals for s-function, i.e, [s|exp(-omega*r12^2)|s], and
622 !> the auxiliary integrals [s|exp(-omega*r12^2)|s]
623 !> \param v matrix storing the integrals
624 !> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)|s]
625 !> \param zetp = 1/zeta
626 !> \param zetq = 1/zetc
627 !> \param zetw = 1/(zeta+zetc)
628 !> \param rho = zeta*zetc*zetw
629 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
630 !> \param omega parameter in the operator
631 !> \param r_cutoff dummy argument for the sake of generality
632 ! **************************************************************************************************
633  SUBROUTINE cps_gauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
634  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
635  INTEGER, INTENT(IN) :: nmax
636  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
637  r_cutoff
638 
639  INTEGER :: n
640  REAL(kind=dp) :: dummy, expt, f0, t, tau
641  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
642 
643  mark_used(r_cutoff)
644 
645  ALLOCATE (f(0:nmax))
646 
647  dummy = zetp
648  dummy = zetq
649  tau = omega/(rho + omega)
650  t = rho*rac2
651  expt = exp(-tau*t)
652  f0 = pi**3*sqrt(zetw**3/(rho + omega)**3)*expt
653 
654  ! *** Calculate the basic two-center integrals [s||s]{n} ***
655 
656  DO n = 1, nmax
657  v(1, 1, n) = f0*tau**(n - 1)
658  END DO
659 
660  DEALLOCATE (f)
661 
662  END SUBROUTINE cps_gauss2
663 
664 ! **************************************************************************************************
665 !> \brief Calculation of truncated Coulomb integrals for s-function, i.e, [s|TC|s] where TC = 1/r12
666 !> if r12 <= r_cutoff and 0 otherwise
667 !> \param v matrix storing the integrals
668 !> \param nmax maximal n in the auxiliary integrals [s|TC|s]
669 !> \param zetp = 1/zeta
670 !> \param zetq = 1/zetc
671 !> \param zetw = 1/(zeta+zetc)
672 !> \param rho = zeta*zetc*zetw
673 !> \param rac2 square distance between center A and C, |Ra-Rc|^2
674 !> \param omega dummy argument for the sake of generality
675 !> \param r_cutoff the radius at which the operator is cut
676 !> \note The truncated operator must have been initialized from the data file prior to this call
677 ! **************************************************************************************************
678  SUBROUTINE cps_truncated2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
679  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
680  INTEGER, INTENT(IN) :: nmax
681  REAL(kind=dp), INTENT(IN) :: zetp, zetq, zetw, rho, rac2, omega, &
682  r_cutoff
683 
684  INTEGER :: n
685  LOGICAL :: use_gamma
686  REAL(kind=dp) :: f0, r, t
687  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
688 
689  mark_used(omega)
690 
691  ALLOCATE (f(nmax + 1)) !t_c_g0 needs to start at index 1
692 
693  r = r_cutoff*sqrt(rho)
694  t = rho*rac2
695  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
696 
697  !check that the operator has been init from file
698  cpassert(get_lmax_init() .GE. nmax)
699 
700  CALL t_c_g0_n(f, use_gamma, r, t, nmax)
701  IF (use_gamma) CALL fgamma(nmax, t, f)
702 
703  DO n = 1, nmax
704  v(1, 1, n) = f0*f(n)
705  END DO
706 
707  DEALLOCATE (f)
708 
709  END SUBROUTINE cps_truncated2
710 
711 END MODULE ai_operators_r12
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:531
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:578
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:634
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:443
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:679
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:486
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:57
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(RES, use_gamma, R, T, NDERIV)
...
Definition: t_c_g0.F:84
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:1464