Line data Source code
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 : ! **************************************************************************************************
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, 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 :
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 50246 : SUBROUTINE operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
109 100492 : 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 50246 : CALL timeset(routineN, handle)
130 :
131 268959408 : v = 0.0_dp
132 :
133 50246 : IF (PRESENT(maxder)) THEN
134 24837440 : vac_plus = 0.0_dp
135 : END IF
136 :
137 50246 : nmax = la_max + lc_max + 1
138 :
139 50246 : na = 0
140 50246 : nap = 0
141 :
142 101996 : DO ipgf = 1, npgfa
143 :
144 51750 : nc = 0
145 51750 : ncp = 0
146 :
147 115984 : DO jpgf = 1, npgfc
148 :
149 64234 : zetp = 1.0_dp/zeta(ipgf)
150 64234 : zetq = 1.0_dp/zetc(jpgf)
151 64234 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
152 :
153 64234 : rho = zeta(ipgf)*zetc(jpgf)*zetw
154 :
155 64234 : 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 156218 : maxder=maxder, vac_plus=vac_plus)
159 : END DO
160 :
161 51750 : IF (PRESENT(maxder)) THEN
162 24000 : na = na + ncoset(la_max - maxder)
163 : ELSE
164 27750 : na = na + ncoset(la_max)
165 : END IF
166 101996 : nap = nap + ncoset(la_max)
167 : END DO
168 :
169 50246 : CALL timestop(handle)
170 :
171 50246 : 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 68038 : SUBROUTINE operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta_a, zeta_c, zetp, zetq, &
197 68038 : 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 68038 : maxder_local = 0
217 68038 : IF (PRESENT(maxder)) maxder_local = maxder
218 68038 : nmax = la_max + lc_max + 1
219 : ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
220 :
221 68038 : IF (lc_max > 0) THEN
222 :
223 58933 : f1 = 0.5_dp*zetq
224 58933 : f2 = -rho*zetq
225 :
226 235732 : rcw(:) = -zeta_a*zetw*rac(:)
227 :
228 : ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
229 :
230 362727 : DO n = 1, nmax - 1
231 303794 : v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
232 303794 : v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
233 362727 : 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 165573 : DO lc = 2, lc_max
241 :
242 674101 : 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 508528 : f2*v(1, coset(0, 0, lc - 2), n + 1))
248 :
249 508528 : cz = lc - 1
250 508528 : v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
251 :
252 1516952 : DO cy = 2, lc
253 1008424 : 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 1516952 : f2*v(1, coset(0, cy - 2, cz), n + 1))
258 : END DO
259 :
260 2025480 : DO cy = 0, lc - 1
261 1516952 : cz = lc - 1 - cy
262 2025480 : v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
263 : END DO
264 :
265 1623592 : DO cx = 2, lc
266 1008424 : f6 = f1*REAL(cx - 1, dp)
267 3294994 : DO cy = 0, lc - cx
268 1778042 : 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 2786466 : 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 68038 : IF (la_max > 0) THEN
285 :
286 58148 : f3 = 0.5_dp*zetp
287 58148 : f4 = -rho*zetp
288 58148 : f5 = 0.5_dp*zetw
289 :
290 232592 : raw(:) = zeta_c*zetw*rac(:)
291 :
292 : ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
293 :
294 359492 : DO n = 1, nmax - 1
295 301344 : v(2, 1, n) = raw(1)*v(1, 1, n + 1)
296 301344 : v(3, 1, n) = raw(2)*v(1, 1, n + 1)
297 359492 : 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 149354 : DO la = 2, la_max
305 :
306 585161 : 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 435807 : f4*v(coset(0, 0, la - 2), 1, n + 1))
312 :
313 435807 : az = la - 1
314 435807 : v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
315 :
316 1175002 : DO ay = 2, la
317 739195 : 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 1175002 : f4*v(coset(0, ay - 2, az), 1, n + 1))
322 : END DO
323 :
324 1610809 : DO ay = 0, la - 1
325 1175002 : az = la - 1 - ay
326 1610809 : v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
327 : END DO
328 :
329 1266208 : DO ax = 2, la
330 739195 : f6 = f3*REAL(ax - 1, dp)
331 2305703 : DO ay = 0, la - ax
332 1130701 : 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 1869896 : 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 210138 : DO lc = 1, lc_max
345 :
346 717304 : DO cx = 0, lc
347 1888912 : DO cy = 0, lc - cx
348 1229756 : cz = lc - cx - cy
349 :
350 1229756 : coc = coset(cx, cy, cz)
351 1229756 : cocx = coset(MAX(0, cx - 1), cy, cz)
352 1229756 : cocy = coset(cx, MAX(0, cy - 1), cz)
353 1229756 : cocz = coset(cx, cy, MAX(0, cz - 1))
354 :
355 1229756 : fcx = f5*REAL(cx, dp)
356 1229756 : fcy = f5*REAL(cy, dp)
357 1229756 : 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 6086801 : DO n = 1, nmax - 1 - lc
363 4857045 : v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
364 4857045 : v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
365 6086801 : 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 4188229 : DO la = 2, la_max
374 :
375 10635712 : 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 6954649 : fcz*v(coset(0, 0, la - 1), cocz, n + 1)
382 :
383 6954649 : az = la - 1
384 : v(coset(0, 1, az), coc, n) = &
385 : raw(2)*v(coset(0, 0, az), coc, n + 1) + &
386 6954649 : fcy*v(coset(0, 0, az), cocy, n + 1)
387 :
388 18658636 : DO ay = 2, la
389 11703987 : 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 18658636 : fcy*v(coset(0, ay - 1, az), cocy, n + 1)
395 : END DO
396 :
397 25613285 : DO ay = 0, la - 1
398 18658636 : az = la - 1 - ay
399 : v(coset(1, ay, az), coc, n) = &
400 : raw(1)*v(coset(0, ay, az), coc, n + 1) + &
401 25613285 : fcx*v(coset(0, ay, az), cocx, n + 1)
402 : END DO
403 :
404 21109943 : DO ax = 2, la
405 11703987 : f6 = f3*REAL(ax - 1, dp)
406 36429407 : DO ay = 0, la - ax
407 17770771 : 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 29474758 : 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 908084 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
428 11073870 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
429 11005832 : vac(na + i, nc + j) = v(i, j, 1)
430 : END DO
431 : END DO
432 :
433 68038 : IF (PRESENT(maxder)) THEN
434 929600 : DO j = 1, ncoset(lc_max)
435 24837440 : DO i = 1, ncoset(la_max)
436 24813440 : vac_plus(nap + i, ncp + j) = v(i, j, 1)
437 : END DO
438 : END DO
439 : END IF
440 :
441 68038 : nc = nc + ncoset(lc_max - maxder_local)
442 68038 : ncp = ncp + ncoset(lc_max)
443 68038 : 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 30602 : 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 30602 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
468 :
469 : MARK_USED(omega)
470 : MARK_USED(r_cutoff)
471 :
472 91806 : ALLOCATE (f(0:nmax))
473 30602 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
474 :
475 : ! *** Calculate the incomplete Gamma/Boys function ***
476 :
477 30602 : t = rho*rac2
478 30602 : CALL fgamma(nmax - 1, t, f)
479 :
480 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
481 :
482 187422 : DO n = 1, nmax
483 187422 : v(1, 1, n) = f0*f(n - 1)
484 : END DO
485 :
486 30602 : DEALLOCATE (f)
487 30602 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
511 :
512 : MARK_USED(r_cutoff)
513 :
514 14400 : ALLOCATE (f(0:nmax))
515 4800 : comega = omega**2/(omega**2 + rho)
516 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw*comega)*zetp*zetq
517 :
518 : ! *** Calculate the incomplete Gamma/Boys function ***
519 :
520 4800 : t = rho*rac2
521 4800 : arg = comega*t
522 4800 : CALL fgamma(nmax - 1, arg, f)
523 :
524 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
525 :
526 43680 : DO n = 1, nmax
527 43680 : v(1, 1, n) = f0*f(n - 1)*comega**(n - 1)
528 : END DO
529 :
530 4800 : DEALLOCATE (f)
531 :
532 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
556 :
557 : MARK_USED(r_cutoff)
558 :
559 19200 : ALLOCATE (fv(0:nmax), fverf(0:nmax))
560 4800 : comega = omega**2/(omega**2 + rho)
561 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
562 :
563 : ! *** Calculate the incomplete Gamma/Boys function ***
564 :
565 4800 : t = rho*rac2
566 4800 : argerf = comega*t
567 :
568 4800 : CALL fgamma(nmax - 1, t, fv)
569 4800 : CALL fgamma(nmax - 1, argerf, fverf)
570 :
571 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
572 :
573 43680 : DO n = 1, nmax
574 43680 : v(1, 1, n) = f0*(fv(n - 1) - SQRT(comega)*comega**(n - 1)*fverf(n - 1))
575 : END DO
576 :
577 4800 : DEALLOCATE (fv, fverf)
578 :
579 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
603 :
604 : MARK_USED(r_cutoff)
605 :
606 14400 : ALLOCATE (f(0:nmax))
607 :
608 : dummy = zetp
609 4800 : dummy = zetq
610 4800 : eta = rho/(rho + omega)
611 4800 : tau = omega/(rho + omega)
612 :
613 : ! *** Calculate the incomplete Gamma/Boys function ***
614 :
615 4800 : t = rho*rac2
616 4800 : arg = eta*t
617 :
618 4800 : CALL fgamma(nmax - 1, arg, f)
619 :
620 4800 : expT = EXP(-omega/(omega + rho)*t)
621 4800 : f0 = 2.0_dp*SQRT(pi**5*zetw**3)/(rho + omega)*expT
622 :
623 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
624 43680 : v(1, 1, 1:nmax) = 0.0_dp
625 43680 : DO n = 1, nmax
626 38880 : fsign = (-1.0_dp)**(n - 1)
627 228512 : DO j = 0, n - 1
628 : v(1, 1, n) = v(1, 1, n) + f0*fsign* &
629 223712 : 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 4800 : DEALLOCATE (f)
634 :
635 4800 : 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 4800 : 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 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
659 :
660 : MARK_USED(r_cutoff)
661 :
662 14400 : ALLOCATE (f(0:nmax))
663 :
664 : dummy = zetp
665 4800 : dummy = zetq
666 4800 : tau = omega/(rho + omega)
667 4800 : t = rho*rac2
668 4800 : expT = EXP(-tau*t)
669 4800 : f0 = pi**3*SQRT(zetw**3/(rho + omega)**3)*expT
670 :
671 : ! *** Calculate the basic two-center integrals [s||s]{n} ***
672 :
673 43680 : DO n = 1, nmax
674 43680 : v(1, 1, n) = f0*tau**(n - 1)
675 : END DO
676 :
677 4800 : DEALLOCATE (f)
678 :
679 4800 : 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 14432 : 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 14432 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
705 :
706 : MARK_USED(omega)
707 :
708 43296 : ALLOCATE (f(nmax + 1)) !t_c_g0 needs to start at index 1
709 :
710 14432 : r = r_cutoff*SQRT(rho)
711 14432 : t = rho*rac2
712 14432 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
713 :
714 : !check that the operator has been init from file
715 14432 : CPASSERT(get_lmax_init() >= nmax)
716 :
717 14432 : CALL t_c_g0_n(f, use_gamma, r, t, nmax)
718 14432 : IF (use_gamma) CALL fgamma(nmax, t, f)
719 :
720 75494 : DO n = 1, nmax
721 75494 : v(1, 1, n) = f0*f(n)
722 : END DO
723 :
724 14432 : DEALLOCATE (f)
725 :
726 14432 : END SUBROUTINE cps_truncated2
727 :
728 : END MODULE ai_operators_r12
|
