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 : !> \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 25018 : SUBROUTINE operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
109 50036 : 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 25018 : CALL timeset(routineN, handle)
135 :
136 242004620 : v = 0.0_dp
137 :
138 25018 : maxder_local = 0
139 25018 : IF (PRESENT(maxder)) THEN
140 24000 : maxder_local = maxder
141 24837440 : vac_plus = 0.0_dp
142 : END IF
143 :
144 25018 : 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 25018 : na = 0
153 25018 : nap = 0
154 :
155 50992 : DO ipgf = 1, npgfa
156 :
157 25974 : nc = 0
158 25974 : ncp = 0
159 :
160 60006 : DO jpgf = 1, npgfc
161 :
162 : ! *** Calculate some prefactors ***
163 :
164 34032 : zetp = 1.0_dp/zeta(ipgf)
165 34032 : zetq = 1.0_dp/zetc(jpgf)
166 34032 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
167 :
168 34032 : rho = zeta(ipgf)*zetc(jpgf)*zetw
169 :
170 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
171 :
172 34032 : CALL cps_operator2(v, nmax, zetp, zetq, zetw, rho, rac2, omega, r_cutoff)
173 :
174 : ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
175 :
176 34032 : IF (lc_max > 0) THEN
177 :
178 33780 : f1 = 0.5_dp*zetq
179 33780 : f2 = -rho*zetq
180 :
181 135120 : rcw(:) = -zeta(ipgf)*zetw*rac(:)
182 :
183 : ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
184 :
185 243768 : DO n = 1, nmax - 1
186 209988 : v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
187 209988 : v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
188 243768 : 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 111204 : DO lc = 2, lc_max
196 :
197 524020 : 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 412816 : f2*v(1, coset(0, 0, lc - 2), n + 1))
205 :
206 : ! *** Increase the angular momentum component y of c ***
207 :
208 412816 : cz = lc - 1
209 412816 : v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
210 :
211 1270126 : DO cy = 2, lc
212 857310 : 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 1270126 : 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 1682942 : DO cy = 0, lc - 1
222 1270126 : cz = lc - 1 - cy
223 1682942 : v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
224 : END DO
225 :
226 1347550 : DO cx = 2, lc
227 857310 : f6 = f1*REAL(cx - 1, dp)
228 2822024 : DO cy = 0, lc - cx
229 1551898 : 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 2409208 : 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 34032 : IF (la_max > 0) THEN
246 :
247 33780 : f3 = 0.5_dp*zetp
248 33780 : f4 = -rho*zetp
249 33780 : f5 = 0.5_dp*zetw
250 :
251 135120 : raw(:) = zetc(jpgf)*zetw*rac(:)
252 :
253 : ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
254 :
255 243768 : DO n = 1, nmax - 1
256 209988 : v(2, 1, n) = raw(1)*v(1, 1, n + 1)
257 209988 : v(3, 1, n) = raw(2)*v(1, 1, n + 1)
258 243768 : 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 99204 : DO la = 2, la_max
266 :
267 448020 : 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 348816 : f4*v(coset(0, 0, la - 2), 1, n + 1))
275 :
276 : ! *** Increase the angular momentum component y of a ***
277 :
278 348816 : az = la - 1
279 348816 : v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
280 :
281 954126 : DO ay = 2, la
282 605310 : 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 954126 : 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 1302942 : DO ay = 0, la - 1
292 954126 : az = la - 1 - ay
293 1302942 : v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
294 : END DO
295 :
296 1019550 : DO ax = 2, la
297 605310 : f6 = f3*REAL(ax - 1, dp)
298 1888424 : DO ay = 0, la - ax
299 934298 : 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 1539608 : 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 144564 : DO lc = 1, lc_max
312 :
313 534348 : DO cx = 0, lc
314 1483920 : DO cy = 0, lc - cx
315 983352 : cz = lc - cx - cy
316 :
317 983352 : coc = coset(cx, cy, cz)
318 983352 : cocx = coset(MAX(0, cx - 1), cy, cz)
319 983352 : cocy = coset(cx, MAX(0, cy - 1), cz)
320 983352 : cocz = coset(cx, cy, MAX(0, cz - 1))
321 :
322 983352 : fcx = f5*REAL(cx, dp)
323 983352 : fcy = f5*REAL(cy, dp)
324 983352 : 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 5130058 : DO n = 1, nmax - 1 - lc
330 4146706 : v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
331 4146706 : v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
332 5130058 : 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 3524612 : DO la = 2, la_max
341 :
342 9373274 : 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 6238446 : fcz*v(coset(0, 0, la - 1), cocz, n + 1)
351 :
352 : ! *** Increase the angular momentum component y of a ***
353 :
354 6238446 : az = la - 1
355 : v(coset(0, 1, az), coc, n) = &
356 : raw(2)*v(coset(0, 0, az), coc, n + 1) + &
357 6238446 : fcy*v(coset(0, 0, az), cocy, n + 1)
358 :
359 16780592 : DO ay = 2, la
360 10542146 : 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 16780592 : 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 23019038 : DO ay = 0, la - 1
371 16780592 : az = la - 1 - ay
372 : v(coset(1, ay, az), coc, n) = &
373 : raw(1)*v(coset(0, ay, az), coc, n + 1) + &
374 23019038 : fcx*v(coset(0, ay, az), cocx, n + 1)
375 : END DO
376 :
377 18932068 : DO ax = 2, la
378 10542146 : f6 = f3*REAL(ax - 1, dp)
379 32757172 : DO ay = 0, la - ax
380 15976580 : 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 26518726 : 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 672836 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
401 9582232 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
402 9548200 : vac(na + i, nc + j) = v(i, j, 1)
403 : END DO
404 : END DO
405 :
406 34032 : IF (PRESENT(maxder)) THEN
407 929600 : DO j = 1, ncoset(lc_max)
408 24837440 : DO i = 1, ncoset(la_max)
409 24813440 : vac_plus(nap + i, ncp + j) = v(i, j, 1)
410 : END DO
411 : END DO
412 : END IF
413 :
414 34032 : nc = nc + ncoset(lc_max - maxder_local)
415 60006 : ncp = ncp + ncoset(lc_max)
416 :
417 : END DO
418 :
419 25974 : na = na + ncoset(la_max - maxder_local)
420 50992 : nap = nap + ncoset(la_max)
421 :
422 : END DO
423 :
424 25018 : CALL timestop(handle)
425 :
426 25018 : 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 14832 : 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 14832 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
451 :
452 : MARK_USED(omega)
453 : MARK_USED(r_cutoff)
454 :
455 44496 : ALLOCATE (f(0:nmax))
456 14832 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
457 :
458 : ! *** Calculate the incomplete Gamma/Boys function ***
459 :
460 14832 : t = rho*rac2
461 14832 : CALL fgamma(nmax - 1, t, f)
462 :
463 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
464 :
465 103752 : DO n = 1, nmax
466 103752 : v(1, 1, n) = f0*f(n - 1)
467 : END DO
468 :
469 14832 : DEALLOCATE (f)
470 14832 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
494 :
495 : MARK_USED(r_cutoff)
496 :
497 14400 : ALLOCATE (f(0:nmax))
498 4800 : comega = omega**2/(omega**2 + rho)
499 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw*comega)*zetp*zetq
500 :
501 : ! *** Calculate the incomplete Gamma/Boys function ***
502 :
503 4800 : t = rho*rac2
504 4800 : arg = comega*t
505 4800 : CALL fgamma(nmax - 1, arg, f)
506 :
507 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
508 :
509 43680 : DO n = 1, nmax
510 43680 : v(1, 1, n) = f0*f(n - 1)*comega**(n - 1)
511 : END DO
512 :
513 4800 : DEALLOCATE (f)
514 :
515 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
539 :
540 : MARK_USED(r_cutoff)
541 :
542 19200 : ALLOCATE (fv(0:nmax), fverf(0:nmax))
543 4800 : comega = omega**2/(omega**2 + rho)
544 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
545 :
546 : ! *** Calculate the incomplete Gamma/Boys function ***
547 :
548 4800 : t = rho*rac2
549 4800 : argerf = comega*t
550 :
551 4800 : CALL fgamma(nmax - 1, t, fv)
552 4800 : CALL fgamma(nmax - 1, argerf, fverf)
553 :
554 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
555 :
556 43680 : DO n = 1, nmax
557 43680 : v(1, 1, n) = f0*(fv(n - 1) - SQRT(comega)*comega**(n - 1)*fverf(n - 1))
558 : END DO
559 :
560 4800 : DEALLOCATE (fv, fverf)
561 :
562 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
586 :
587 : MARK_USED(r_cutoff)
588 :
589 14400 : ALLOCATE (f(0:nmax))
590 :
591 4800 : dummy = zetp
592 4800 : dummy = zetq
593 4800 : eta = rho/(rho + omega)
594 4800 : tau = omega/(rho + omega)
595 :
596 : ! *** Calculate the incomplete Gamma/Boys function ***
597 :
598 4800 : t = rho*rac2
599 4800 : arg = eta*t
600 :
601 4800 : CALL fgamma(nmax - 1, arg, f)
602 :
603 4800 : expT = EXP(-omega/(omega + rho)*t)
604 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw**3)/(rho + omega)*expT
605 :
606 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
607 43680 : v(1, 1, 1:nmax) = 0.0_dp
608 43680 : DO n = 1, nmax
609 38880 : fsign = (-1.0_dp)**(n - 1)
610 228512 : DO j = 0, n - 1
611 : v(1, 1, n) = v(1, 1, n) + f0*fsign* &
612 223712 : 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 4800 : DEALLOCATE (f)
617 :
618 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
642 :
643 : MARK_USED(r_cutoff)
644 :
645 14400 : ALLOCATE (f(0:nmax))
646 :
647 4800 : dummy = zetp
648 4800 : dummy = zetq
649 4800 : tau = omega/(rho + omega)
650 4800 : t = rho*rac2
651 4800 : expT = EXP(-tau*t)
652 4800 : f0 = pi**3*SQRT(zetw**3/(rho + omega)**3)*expT
653 :
654 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
655 :
656 43680 : DO n = 1, nmax
657 43680 : v(1, 1, n) = f0*tau**(n - 1)
658 : END DO
659 :
660 4800 : DEALLOCATE (f)
661 :
662 4800 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
688 :
689 : MARK_USED(omega)
690 :
691 0 : ALLOCATE (f(nmax + 1)) !t_c_g0 needs to start at index 1
692 :
693 0 : r = r_cutoff*SQRT(rho)
694 0 : t = rho*rac2
695 0 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
696 :
697 : !check that the operator has been init from file
698 0 : CPASSERT(get_lmax_init() .GE. nmax)
699 :
700 0 : CALL t_c_g0_n(f, use_gamma, r, t, nmax)
701 0 : IF (use_gamma) CALL fgamma(nmax, t, f)
702 :
703 0 : DO n = 1, nmax
704 0 : v(1, 1, n) = f0*f(n)
705 : END DO
706 :
707 0 : DEALLOCATE (f)
708 :
709 0 : END SUBROUTINE cps_truncated2
710 :
711 : END MODULE ai_operators_r12
|
