(git:374b731)
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-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! **************************************************************************************************
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, 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
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 :: 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
711END 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: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