(git:cdd058b)
Loading...
Searching...
No Matches
ai_coulomb.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of Coulomb integrals over Cartesian Gaussian-type functions
10!> (electron repulsion integrals, ERIs).
11!> \par Literature
12!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13!> \par History
14!> none
15!> \par Parameters
16!> - ax,ay,az : Angular momentum index numbers of orbital a.
17!> - bx,by,bz : Angular momentum index numbers of orbital b.
18!> - cx,cy,cz : Angular momentum index numbers of orbital c.
19!> - coset : Cartesian orbital set pointer.
20!> - dab : Distance between the atomic centers a and b.
21!> - dac : Distance between the atomic centers a and c.
22!> - dbc : Distance between the atomic centers b and c.
23!> - gccc : Prefactor of the primitive Gaussian function c.
24!> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
25!> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
26!> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
27!> - ncoset : Number of orbitals in a Cartesian orbital set.
28!> - npgf{a,b} : Degree of contraction of shell a or b.
29!> - rab : Distance vector between the atomic centers a and b.
30!> - rab2 : Square of the distance between the atomic centers a and b.
31!> - rac : Distance vector between the atomic centers a and c.
32!> - rac2 : Square of the distance between the atomic centers a and c.
33!> - rbc : Distance vector between the atomic centers b and c.
34!> - rbc2 : Square of the distance between the atomic centers b and c.
35!> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
36!> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
37!> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
38!> - zetw : Reciprocal of the sum of the exponents of orbital a, b and c.
39!> \author Matthias Krack (22.08.2000)
40! **************************************************************************************************
41MODULE ai_coulomb
42
43 USE ai_operators_r12, ONLY: operator2_recurrence
44 USE gamma, ONLY: fgamma => fgamma_0
45 USE kinds, ONLY: dp
46 USE mathconstants, ONLY: pi
47 USE orbital_pointers, ONLY: coset,&
48 ncoset
49#include "../base/base_uses.f90"
50
51 IMPLICIT NONE
52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_coulomb'
53 PRIVATE
54
55 ! *** Public subroutines ***
56
57 PUBLIC :: coulomb2, coulomb3
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief Calculation of the primitive two-center Coulomb integrals over
63!> Cartesian Gaussian-type functions.
64!> \param la_max ...
65!> \param npgfa ...
66!> \param zeta ...
67!> \param rpgfa ...
68!> \param la_min ...
69!> \param lc_max ...
70!> \param npgfc ...
71!> \param zetc ...
72!> \param rpgfc ...
73!> \param lc_min ...
74!> \param rac ...
75!> \param rac2 ...
76!> \param vac ...
77!> \param v ...
78!> \param f ...
79!> \param screening optional primitive-pair screening switch
80!> \date 05.12.2000
81!> \author Matthias Krack
82!> \version 1.0
83! **************************************************************************************************
84 SUBROUTINE coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, &
85 rac, rac2, vac, v, f, screening)
86 INTEGER, INTENT(IN) :: la_max, npgfa
87 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
88 INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
89 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
90 INTEGER, INTENT(IN) :: lc_min
91 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
92 REAL(kind=dp), INTENT(IN) :: rac2
93 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
94 REAL(kind=dp), DIMENSION(:, :, :) :: v
95 REAL(kind=dp), DIMENSION(0:) :: f
96 LOGICAL, INTENT(IN), OPTIONAL :: screening
97
98 INTEGER :: i, ipgf, j, jpgf, n, na, nap, nc, ncp, &
99 nmax
100 LOGICAL :: do_screening
101 REAL(kind=dp) :: dac, f0, rho, t, zetp, zetq, zetw
102
103 do_screening = .true.
104 IF (PRESENT(screening)) do_screening = screening
105
106 v = 0.0_dp
107
108 nmax = la_max + lc_max + 1
109
110 dac = sqrt(rac2)
111
112 na = 0
113 nap = 0
114 DO ipgf = 1, npgfa
115
116 nc = 0
117 ncp = 0
118
119 DO jpgf = 1, npgfc
120
121 IF (do_screening .AND. rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
122 DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
123 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
124 vac(i, j) = 0.0_dp
125 END DO
126 END DO
127 nc = nc + ncoset(lc_max)
128 cycle
129 END IF
130
131 zetp = 1.0_dp/zeta(ipgf)
132 zetq = 1.0_dp/zetc(jpgf)
133 zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
134
135 rho = zeta(ipgf)*zetc(jpgf)*zetw
136
137 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
138
139 t = rho*rac2
140
141 CALL fgamma(nmax - 1, t, f)
142
143 DO n = 1, nmax
144 v(1, 1, n) = f0*f(n - 1)
145 END DO
146
147 CALL operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta(ipgf), zetc(jpgf), &
148 zetp, zetq, zetw, rho, rac, vac, v, na, nc, nap, ncp)
149
150 END DO
151
152 na = na + ncoset(la_max)
153 nap = nap + ncoset(la_max)
154
155 END DO
156
157 END SUBROUTINE coulomb2
158! **************************************************************************************************
159!> \brief Calculation of the primitive three-center Coulomb integrals over
160!> Cartesian Gaussian-type functions (electron repulsion integrals,
161!> ERIs).
162!> \param la_max ...
163!> \param npgfa ...
164!> \param zeta ...
165!> \param rpgfa ...
166!> \param la_min ...
167!> \param lb_max ...
168!> \param npgfb ...
169!> \param zetb ...
170!> \param rpgfb ...
171!> \param lb_min ...
172!> \param lc_max ...
173!> \param zetc ...
174!> \param rpgfc ...
175!> \param lc_min ...
176!> \param gccc ...
177!> \param rab ...
178!> \param rab2 ...
179!> \param rac ...
180!> \param rac2 ...
181!> \param rbc2 ...
182!> \param vabc ...
183!> \param int_abc ...
184!> \param v ...
185!> \param f ...
186!> \param maxder ...
187!> \param vabc_plus ...
188!> \date 06.11.2000
189!> \author Matthias Krack
190!> \version 1.0
191! **************************************************************************************************
192 SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
193 lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
194 v, f, maxder, vabc_plus)
195 INTEGER, INTENT(IN) :: la_max, npgfa
196 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
197 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
198 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
199 INTEGER, INTENT(IN) :: lb_min, lc_max
200 REAL(kind=dp), INTENT(IN) :: zetc, rpgfc
201 INTEGER, INTENT(IN) :: lc_min
202 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: gccc
203 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
204 REAL(kind=dp), INTENT(IN) :: rab2
205 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
206 REAL(kind=dp), INTENT(IN) :: rac2, rbc2
207 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
208 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: int_abc
209 REAL(kind=dp), DIMENSION(:, :, :, :) :: v
210 REAL(kind=dp), DIMENSION(0:) :: f
211 INTEGER, INTENT(IN), OPTIONAL :: maxder
212 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus
213
214 INTEGER :: ax, ay, az, bx, by, bz, coc, cocx, cocy, &
215 cocz, cx, cy, cz, i, ipgf, j, jpgf, k, &
216 kk, la, la_start, lb, lc, &
217 maxder_local, n, na, nap, nb, nmax
218 REAL(kind=dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, f5, &
219 f6, f7, fcx, fcy, fcz, fx, fy, fz, t, &
220 zetp, zetq, zetw
221 REAL(kind=dp), DIMENSION(3) :: rap, rbp, rcp, rcw, rpw
222
223 v = 0.0_dp
224
225 maxder_local = 0
226 IF (PRESENT(maxder)) THEN
227 maxder_local = maxder
228 END IF
229
230 nmax = la_max + lb_max + lc_max + 1
231
232 ! *** Calculate the distances of the centers a, b and c ***
233
234 dab = sqrt(rab2)
235 dac = sqrt(rac2)
236 dbc = sqrt(rbc2)
237
238 ! *** Initialize integrals array
239 int_abc = 0.0_dp
240
241 ! *** Loop over all pairs of primitive Gaussian-type functions ***
242
243 na = 0
244 nap = 0
245
246 DO ipgf = 1, npgfa
247
248 ! *** Screening ***
249 IF (rpgfa(ipgf) + rpgfc < dac) THEN
250 na = na + ncoset(la_max - maxder_local)
251 nap = nap + ncoset(la_max)
252 cycle
253 END IF
254
255 nb = 0
256
257 DO jpgf = 1, npgfb
258
259 ! *** Screening ***
260 IF ( &
261 (rpgfb(jpgf) + rpgfc < dbc) .OR. &
262 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
263 nb = nb + ncoset(lb_max)
264 cycle
265 END IF
266
267 ! *** Calculate some prefactors ***
268
269 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
270 zetq = 1.0_dp/zetc
271 zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
272
273 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
274 f1 = zetb(jpgf)*zetp
275 f2 = 0.5_dp*zetp
276 f4 = -zetc*zetw
277
278 f0 = f0*exp(-zeta(ipgf)*f1*rab2)
279
280 rap(:) = f1*rab(:)
281 rcp(:) = rap(:) - rac(:)
282 rpw(:) = f4*rcp(:)
283
284 ! *** Calculate the incomplete Gamma function ***
285
286 t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
287
288 CALL fgamma(nmax - 1, t, f)
289
290 ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
291
292 DO n = 1, nmax
293 v(1, 1, 1, n) = f0*f(n - 1)
294 END DO
295
296 ! *** Recurrence steps: [ss||s] -> [as||s] ***
297
298 IF (la_max > 0) THEN
299
300 ! *** Vertical recurrence steps: [ss||s] -> [as||s] ***
301
302 ! *** [ps||s]{n} = (Pi - Ai)*[ss||s]{n} + ***
303 ! *** (Wi - Pi)*[ss||s]{n+1} (i = x,y,z) ***
304
305 DO n = 1, nmax - 1
306 v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
307 v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
308 v(4, 1, 1, n) = rap(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
309 END DO
310
311 ! *** [as||s]{n} = (Pi - Ai)*[(a-1i)s||s]{n} + ***
312 ! *** (Wi - Pi)*[(a-1i)s||s]{n+1} + ***
313 ! *** f2*Ni(a-1i)*( [(a-2i)s||s]{n} + ***
314 ! *** f4*[(a-2i)s||s]{n+1}) ***
315
316 DO la = 2, la_max
317
318 DO n = 1, nmax - la
319
320 ! *** Increase the angular momentum component z of a ***
321
322 v(coset(0, 0, la), 1, 1, n) = &
323 rap(3)*v(coset(0, 0, la - 1), 1, 1, n) + &
324 rpw(3)*v(coset(0, 0, la - 1), 1, 1, n + 1) + &
325 f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, 1, n) + &
326 f4*v(coset(0, 0, la - 2), 1, 1, n + 1))
327
328 ! *** Increase the angular momentum component y of a ***
329
330 az = la - 1
331 v(coset(0, 1, az), 1, 1, n) = &
332 rap(2)*v(coset(0, 0, az), 1, 1, n) + &
333 rpw(2)*v(coset(0, 0, az), 1, 1, n + 1)
334
335 DO ay = 2, la
336 az = la - ay
337 v(coset(0, ay, az), 1, 1, n) = &
338 rap(2)*v(coset(0, ay - 1, az), 1, 1, n) + &
339 rpw(2)*v(coset(0, ay - 1, az), 1, 1, n + 1) + &
340 f2*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, 1, n) + &
341 f4*v(coset(0, ay - 2, az), 1, 1, n + 1))
342 END DO
343
344 ! *** Increase the angular momentum component x of a ***
345
346 DO ay = 0, la - 1
347 az = la - 1 - ay
348 v(coset(1, ay, az), 1, 1, n) = &
349 rap(1)*v(coset(0, ay, az), 1, 1, n) + &
350 rpw(1)*v(coset(0, ay, az), 1, 1, n + 1)
351 END DO
352
353 DO ax = 2, la
354 f3 = f2*real(ax - 1, dp)
355 DO ay = 0, la - ax
356 az = la - ax - ay
357 v(coset(ax, ay, az), 1, 1, n) = &
358 rap(1)*v(coset(ax - 1, ay, az), 1, 1, n) + &
359 rpw(1)*v(coset(ax - 1, ay, az), 1, 1, n + 1) + &
360 f3*(v(coset(ax - 2, ay, az), 1, 1, n) + &
361 f4*v(coset(ax - 2, ay, az), 1, 1, n + 1))
362 END DO
363 END DO
364
365 END DO
366
367 END DO
368
369 ! *** Recurrence steps: [as||s] -> [ab||s] ***
370
371 IF (lb_max > 0) THEN
372
373 ! *** Horizontal recurrence steps ***
374
375 rbp(:) = rap(:) - rab(:)
376
377 ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
378
379 la_start = max(0, la_min - 1)
380
381 DO la = la_start, la_max - 1
382 DO n = 1, nmax - la - 1
383 DO ax = 0, la
384 DO ay = 0, la - ax
385 az = la - ax - ay
386 v(coset(ax, ay, az), 2, 1, n) = &
387 v(coset(ax + 1, ay, az), 1, 1, n) - &
388 rab(1)*v(coset(ax, ay, az), 1, 1, n)
389 v(coset(ax, ay, az), 3, 1, n) = &
390 v(coset(ax, ay + 1, az), 1, 1, n) - &
391 rab(2)*v(coset(ax, ay, az), 1, 1, n)
392 v(coset(ax, ay, az), 4, 1, n) = &
393 v(coset(ax, ay, az + 1), 1, 1, n) - &
394 rab(3)*v(coset(ax, ay, az), 1, 1, n)
395 END DO
396 END DO
397 END DO
398 END DO
399
400 ! *** Vertical recurrence step ***
401
402 ! *** [ap||s]{n} = (Pi - Bi)*[as||s]{n} + ***
403 ! *** (Wi - Pi)*[as||s]{n+1} + ***
404 ! *** f2*Ni(a)*( [(a-1i)s||s]{n} + ***
405 ! *** f4*[(a-1i)s||s]{n+1}) ***
406
407 DO n = 1, nmax - la_max - 1
408 DO ax = 0, la_max
409 fx = f2*real(ax, dp)
410 DO ay = 0, la_max - ax
411 fy = f2*real(ay, dp)
412 az = la_max - ax - ay
413 fz = f2*real(az, dp)
414
415 IF (ax == 0) THEN
416 v(coset(ax, ay, az), 2, 1, n) = &
417 rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
418 rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1)
419 ELSE
420 v(coset(ax, ay, az), 2, 1, n) = &
421 rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
422 rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1) + &
423 fx*(v(coset(ax - 1, ay, az), 1, 1, n) + &
424 f4*v(coset(ax - 1, ay, az), 1, 1, n + 1))
425 END IF
426
427 IF (ay == 0) THEN
428 v(coset(ax, ay, az), 3, 1, n) = &
429 rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
430 rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1)
431 ELSE
432 v(coset(ax, ay, az), 3, 1, n) = &
433 rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
434 rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1) + &
435 fy*(v(coset(ax, ay - 1, az), 1, 1, n) + &
436 f4*v(coset(ax, ay - 1, az), 1, 1, n + 1))
437 END IF
438
439 IF (az == 0) THEN
440 v(coset(ax, ay, az), 4, 1, n) = &
441 rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
442 rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1)
443 ELSE
444 v(coset(ax, ay, az), 4, 1, n) = &
445 rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
446 rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1) + &
447 fz*(v(coset(ax, ay, az - 1), 1, 1, n) + &
448 f4*v(coset(ax, ay, az - 1), 1, 1, n + 1))
449 END IF
450
451 END DO
452 END DO
453 END DO
454
455 ! *** Recurrence steps: [ap||s] -> [ab||s] ***
456
457 DO lb = 2, lb_max
458
459 ! *** Horizontal recurrence steps ***
460
461 ! *** [ab||s]{n} = [(a+1i)(b-1i)||s]{n} - ***
462 ! *** (Bi - Ai)*[a(b-1i)||s]{n} ***
463
464 la_start = max(0, la_min - 1)
465
466 DO la = la_start, la_max - 1
467 DO n = 1, nmax - la - lb
468 DO ax = 0, la
469 DO ay = 0, la - ax
470 az = la - ax - ay
471
472 ! *** Shift of angular momentum component z from a to b ***
473
474 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
475 v(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1, n) - &
476 rab(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n)
477
478 ! *** Shift of angular momentum component y from a to b ***
479
480 DO by = 1, lb
481 bz = lb - by
482 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
483 v(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1, n) - &
484 rab(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n)
485 END DO
486
487 ! *** Shift of angular momentum component x from a to b ***
488
489 DO bx = 1, lb
490 DO by = 0, lb - bx
491 bz = lb - bx - by
492 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
493 v(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1, n) - &
494 rab(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n)
495 END DO
496 END DO
497
498 END DO
499 END DO
500 END DO
501 END DO
502
503 ! *** Vertical recurrence step ***
504
505 ! *** [ab||s]{n} = (Pi - Bi)*[a(b-1i)||s]{n} + ***
506 ! *** (Wi - Pi)*[a(b-1i)||s]{n+1} + ***
507 ! *** f2*Ni(a)*( [(a-1i)(b-1i)||s]{n} + ***
508 ! *** f4*[(a-1i)(b-1i)||s]{n+1}) ***
509 ! *** f2*Ni(b-1i)*( [a(b-2i)||s]{n} + ***
510 ! *** f4*[a(b-2i)||s]{n+1}) ***
511
512 DO n = 1, nmax - la_max - lb
513 DO ax = 0, la_max
514 fx = f2*real(ax, dp)
515 DO ay = 0, la_max - ax
516 fy = f2*real(ay, dp)
517 az = la_max - ax - ay
518 fz = f2*real(az, dp)
519
520 ! *** Shift of angular momentum component z from a to b ***
521
522 f3 = f2*real(lb - 1, dp)
523
524 IF (az == 0) THEN
525 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
526 rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
527 rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
528 f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
529 f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
530 ELSE
531 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
532 rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
533 rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
534 fz*(v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n) + &
535 f4*v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n + 1)) + &
536 f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
537 f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
538 END IF
539
540 ! *** Shift of angular momentum component y from a to b ***
541
542 IF (ay == 0) THEN
543 bz = lb - 1
544 v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
545 rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
546 rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1)
547 DO by = 2, lb
548 bz = lb - by
549 f3 = f2*real(by - 1, dp)
550 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
551 rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
552 rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
553 f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
554 f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
555 END DO
556 ELSE
557 bz = lb - 1
558 v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
559 rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
560 rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1) + &
561 fy*(v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n) + &
562 f4*v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n + 1))
563 DO by = 2, lb
564 bz = lb - by
565 f3 = f2*real(by - 1, dp)
566 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
567 rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
568 rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
569 fy*(v(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1, n) + &
570 f4*v(coset(ax, ay - 1, az), &
571 coset(0, by - 1, bz), 1, n + 1)) + &
572 f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
573 f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
574 END DO
575 END IF
576
577 ! *** Shift of angular momentum component x from a to b ***
578
579 IF (ax == 0) THEN
580 DO by = 0, lb - 1
581 bz = lb - 1 - by
582 v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
583 rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
584 rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1)
585 END DO
586 DO bx = 2, lb
587 f3 = f2*real(bx - 1, dp)
588 DO by = 0, lb - bx
589 bz = lb - bx - by
590 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
591 rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
592 rpw(1)*v(coset(ax, ay, az), &
593 coset(bx - 1, by, bz), 1, n + 1) + &
594 f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
595 f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
596 END DO
597 END DO
598 ELSE
599 DO by = 0, lb - 1
600 bz = lb - 1 - by
601 v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
602 rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
603 rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1) + &
604 fx*(v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n) + &
605 f4*v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n + 1))
606 END DO
607 DO bx = 2, lb
608 f3 = f2*real(bx - 1, dp)
609 DO by = 0, lb - bx
610 bz = lb - bx - by
611 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
612 rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
613 rpw(1)*v(coset(ax, ay, az), &
614 coset(bx - 1, by, bz), 1, n + 1) + &
615 fx*(v(coset(ax - 1, ay, az), &
616 coset(bx - 1, by, bz), 1, n) + &
617 f4*v(coset(ax - 1, ay, az), &
618 coset(bx - 1, by, bz), 1, n + 1)) + &
619 f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
620 f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
621 END DO
622 END DO
623 END IF
624
625 END DO
626 END DO
627 END DO
628
629 END DO
630
631 END IF
632
633 ELSE
634
635 IF (lb_max > 0) THEN
636
637 ! *** Vertical recurrence steps: [ss||s] -> [sb||s] ***
638
639 rbp(:) = rap(:) - rab(:)
640
641 ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
642 ! *** (Wi - Pi)*[ss||s]{n+1} ***
643
644 DO n = 1, nmax - 1
645 v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
646 v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
647 v(1, 4, 1, n) = rbp(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
648 END DO
649
650 ! *** [sb||s]{n} = (Pi - Bi)*[s(b-1i)||s]{n} + ***
651 ! *** (Wi - Pi)*[s(b-1i)||s]{n+1} + ***
652 ! *** f2*Ni(b-1i)*( [s(b-2i)||s]{n} + ***
653 ! *** f4*[s(b-2i)||s]{n+1}) ***
654
655 DO lb = 2, lb_max
656
657 DO n = 1, nmax - lb
658
659 ! *** Increase the angular momentum component z of b ***
660
661 v(1, coset(0, 0, lb), 1, n) = &
662 rbp(3)*v(1, coset(0, 0, lb - 1), 1, n) + &
663 rpw(3)*v(1, coset(0, 0, lb - 1), 1, n + 1) + &
664 f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), 1, n) + &
665 f4*v(1, coset(0, 0, lb - 2), 1, n + 1))
666
667 ! *** Increase the angular momentum component y of b ***
668
669 bz = lb - 1
670 v(1, coset(0, 1, bz), 1, n) = &
671 rbp(2)*v(1, coset(0, 0, bz), 1, n) + &
672 rpw(2)*v(1, coset(0, 0, bz), 1, n + 1)
673
674 DO by = 2, lb
675 bz = lb - by
676 v(1, coset(0, by, bz), 1, n) = &
677 rbp(2)*v(1, coset(0, by - 1, bz), 1, n) + &
678 rpw(2)*v(1, coset(0, by - 1, bz), 1, n + 1) + &
679 f2*real(by - 1, dp)*(v(1, coset(0, by - 2, bz), 1, n) + &
680 f4*v(1, coset(0, by - 2, bz), 1, n + 1))
681 END DO
682
683 ! *** Increase the angular momentum component x of b ***
684
685 DO by = 0, lb - 1
686 bz = lb - 1 - by
687 v(1, coset(1, by, bz), 1, n) = &
688 rbp(1)*v(1, coset(0, by, bz), 1, n) + &
689 rpw(1)*v(1, coset(0, by, bz), 1, n + 1)
690 END DO
691
692 DO bx = 2, lb
693 f3 = f2*real(bx - 1, dp)
694 DO by = 0, lb - bx
695 bz = lb - bx - by
696 v(1, coset(bx, by, bz), 1, n) = &
697 rbp(1)*v(1, coset(bx - 1, by, bz), 1, n) + &
698 rpw(1)*v(1, coset(bx - 1, by, bz), 1, n + 1) + &
699 f3*(v(1, coset(bx - 2, by, bz), 1, n) + &
700 f4*v(1, coset(bx - 2, by, bz), 1, n + 1))
701 END DO
702 END DO
703
704 END DO
705
706 END DO
707
708 END IF
709
710 END IF
711
712 ! *** Recurrence steps: [ab||s] -> [ab||c] ***
713
714 IF (lc_max > 0) THEN
715
716 ! *** Vertical recurrence steps: [ss||s] -> [ss||c] ***
717
718 f5 = -zetw/zetp
719 f6 = 0.5_dp*zetw
720 f7 = 0.5_dp*zetq
721
722 rcw(:) = rcp(:) + rpw(:)
723
724 ! *** [ss||p]{n} = (Wi - Ci)*[ss||s]{n+1} (i = x,y,z) ***
725
726 DO n = 1, nmax - 1
727 v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
728 v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
729 v(1, 1, 4, n) = rcw(3)*v(1, 1, 1, n + 1)
730 END DO
731
732 ! *** [ss||c]{n} = (Wi - Ci)*[ss||c-1i]{n+1} + ***
733 ! *** f7*Ni(c-1i)*[ss||c-2i]{n} + ***
734 ! *** f5*[ss||c-2i]{n+1} ***
735
736 DO lc = 2, lc_max
737
738 DO n = 1, nmax - lc
739
740 ! *** Increase the angular momentum component z of c ***
741
742 v(1, 1, coset(0, 0, lc), n) = &
743 rcw(3)*v(1, 1, coset(0, 0, lc - 1), n + 1) + &
744 f7*real(lc - 1, dp)*(v(1, 1, coset(0, 0, lc - 2), n) + &
745 f5*v(1, 1, coset(0, 0, lc - 2), n + 1))
746
747 ! *** Increase the angular momentum component y of c ***
748
749 cz = lc - 1
750 v(1, 1, coset(0, 1, cz), n) = rcw(2)*v(1, 1, coset(0, 0, cz), n + 1)
751
752 DO cy = 2, lc
753 cz = lc - cy
754 v(1, 1, coset(0, cy, cz), n) = &
755 rcw(2)*v(1, 1, coset(0, cy - 1, cz), n + 1) + &
756 f7*real(cy - 1, dp)*(v(1, 1, coset(0, cy - 2, cz), n) + &
757 f5*v(1, 1, coset(0, cy - 2, cz), n + 1))
758 END DO
759
760 ! *** Increase the angular momentum component x of c ***
761
762 DO cy = 0, lc - 1
763 cz = lc - 1 - cy
764 v(1, 1, coset(1, cy, cz), n) = rcw(1)*v(1, 1, coset(0, cy, cz), n + 1)
765 END DO
766
767 DO cx = 2, lc
768 DO cy = 0, lc - cx
769 cz = lc - cx - cy
770 v(1, 1, coset(cx, cy, cz), n) = &
771 rcw(1)*v(1, 1, coset(cx - 1, cy, cz), n + 1) + &
772 f7*real(cx - 1, dp)*(v(1, 1, coset(cx - 2, cy, cz), n) + &
773 f5*v(1, 1, coset(cx - 2, cy, cz), n + 1))
774 END DO
775 END DO
776
777 END DO
778
779 END DO
780
781 ! *** Recurrence steps: [ss||c] -> [ab||c] ***
782
783 DO lc = 1, lc_max
784
785 DO cx = 0, lc
786 DO cy = 0, lc - cx
787 cz = lc - cx - cy
788
789 coc = coset(cx, cy, cz)
790 cocx = coset(max(0, cx - 1), cy, cz)
791 cocy = coset(cx, max(0, cy - 1), cz)
792 cocz = coset(cx, cy, max(0, cz - 1))
793
794 fcx = f6*real(cx, dp)
795 fcy = f6*real(cy, dp)
796 fcz = f6*real(cz, dp)
797
798 ! *** Recurrence steps: [ss||c] -> [as||c] ***
799
800 IF (la_max > 0) THEN
801
802 ! *** Vertical recurrence steps: [ss||c] -> [as||c] ***
803
804 ! *** [ps||c]{n} = (Pi - Ai)*[ss||c]{n} + ***
805 ! *** (Wi - Pi)*[ss||c]{n+1} + ***
806 ! *** f6*Ni(c)*[ss||c-1i]{n+1} (i = x,y,z) ***
807
808 DO n = 1, nmax - 1 - lc
809 v(2, 1, coc, n) = rap(1)*v(1, 1, coc, n) + &
810 rpw(1)*v(1, 1, coc, n + 1) + &
811 fcx*v(1, 1, cocx, n + 1)
812 v(3, 1, coc, n) = rap(2)*v(1, 1, coc, n) + &
813 rpw(2)*v(1, 1, coc, n + 1) + &
814 fcy*v(1, 1, cocy, n + 1)
815 v(4, 1, coc, n) = rap(3)*v(1, 1, coc, n) + &
816 rpw(3)*v(1, 1, coc, n + 1) + &
817 fcz*v(1, 1, cocz, n + 1)
818 END DO
819
820 ! *** [as||c]{n} = (Pi - Ai)*[(a-1i)s||c]{n} + ***
821 ! *** (Wi - Pi)*[(a-1i)s||c]{n+1} + ***
822 ! *** f2*Ni(a-1i)*( [(a-2i)s||c]{n} + ***
823 ! *** f4*[(a-2i)s||c]{n+1}) + ***
824 ! *** f6*Ni(c)*[(a-1i)s||c-1i]{n+1} ***
825
826 DO la = 2, la_max
827
828 DO n = 1, nmax - la - lc
829
830 ! *** Increase the angular momentum component z of a ***
831
832 v(coset(0, 0, la), 1, coc, n) = &
833 rap(3)*v(coset(0, 0, la - 1), 1, coc, n) + &
834 rpw(3)*v(coset(0, 0, la - 1), 1, coc, n + 1) + &
835 f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, coc, n) + &
836 f4*v(coset(0, 0, la - 2), 1, coc, n + 1)) + &
837 fcz*v(coset(0, 0, la - 1), 1, cocz, n + 1)
838
839 ! *** Increase the angular momentum component y of a ***
840
841 az = la - 1
842 v(coset(0, 1, az), 1, coc, n) = &
843 rap(2)*v(coset(0, 0, az), 1, coc, n) + &
844 rpw(2)*v(coset(0, 0, az), 1, coc, n + 1) + &
845 fcy*v(coset(0, 0, az), 1, cocy, n + 1)
846
847 DO ay = 2, la
848 f3 = f2*real(ay - 1, dp)
849 az = la - ay
850 v(coset(0, ay, az), 1, coc, n) = &
851 rap(2)*v(coset(0, ay - 1, az), 1, coc, n) + &
852 rpw(2)*v(coset(0, ay - 1, az), 1, coc, n + 1) + &
853 f3*(v(coset(0, ay - 2, az), 1, coc, n) + &
854 f4*v(coset(0, ay - 2, az), 1, coc, n + 1)) + &
855 fcy*v(coset(0, ay - 1, az), 1, cocy, n + 1)
856 END DO
857
858 ! *** Increase the angular momentum component x of a ***
859
860 DO ay = 0, la - 1
861 az = la - 1 - ay
862 v(coset(1, ay, az), 1, coc, n) = &
863 rap(1)*v(coset(0, ay, az), 1, coc, n) + &
864 rpw(1)*v(coset(0, ay, az), 1, coc, n + 1) + &
865 fcx*v(coset(0, ay, az), 1, cocx, n + 1)
866 END DO
867
868 DO ax = 2, la
869 f3 = f2*real(ax - 1, dp)
870 DO ay = 0, la - ax
871 az = la - ax - ay
872 v(coset(ax, ay, az), 1, coc, n) = &
873 rap(1)*v(coset(ax - 1, ay, az), 1, coc, n) + &
874 rpw(1)*v(coset(ax - 1, ay, az), 1, coc, n + 1) + &
875 f3*(v(coset(ax - 2, ay, az), 1, coc, n) + &
876 f4*v(coset(ax - 2, ay, az), 1, coc, n + 1)) + &
877 fcx*v(coset(ax - 1, ay, az), 1, cocx, n + 1)
878 END DO
879 END DO
880
881 END DO
882
883 END DO
884
885 ! *** Recurrence steps: [as||c] -> [ab||c] ***
886
887 IF (lb_max > 0) THEN
888
889 ! *** Horizontal recurrence steps ***
890
891 ! *** [ap||c]{n} = [(a+1i)s||c]{n} - (Bi - Ai)*[as||c]{n} ***
892
893 la_start = max(0, la_min - 1)
894
895 DO la = la_start, la_max - 1
896 DO n = 1, nmax - la - 1 - lc
897 DO ax = 0, la
898 DO ay = 0, la - ax
899 az = la - ax - ay
900 v(coset(ax, ay, az), 2, coc, n) = &
901 v(coset(ax + 1, ay, az), 1, coc, n) - &
902 rab(1)*v(coset(ax, ay, az), 1, coc, n)
903 v(coset(ax, ay, az), 3, coc, n) = &
904 v(coset(ax, ay + 1, az), 1, coc, n) - &
905 rab(2)*v(coset(ax, ay, az), 1, coc, n)
906 v(coset(ax, ay, az), 4, coc, n) = &
907 v(coset(ax, ay, az + 1), 1, coc, n) - &
908 rab(3)*v(coset(ax, ay, az), 1, coc, n)
909 END DO
910 END DO
911 END DO
912 END DO
913
914 ! *** Vertical recurrence step ***
915
916 ! *** [ap||c]{n} = (Pi - Bi)*[as||c]{n} + ***
917 ! *** (Wi - Pi)*[as||c]{n+1} + ***
918 ! *** f2*Ni(a)*( [(a-1i)s||c]{n} + ***
919 ! *** f4*[(a-1i)s||c]{n+1}) + ***
920 ! *** f6*Ni(c)*[(as||c-1i]{n+1}) ***
921
922 DO n = 1, nmax - la_max - 1 - lc
923 DO ax = 0, la_max
924 fx = f2*real(ax, dp)
925 DO ay = 0, la_max - ax
926 fy = f2*real(ay, dp)
927 az = la_max - ax - ay
928 fz = f2*real(az, dp)
929
930 IF (ax == 0) THEN
931 v(coset(ax, ay, az), 2, coc, n) = &
932 rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
933 rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
934 fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
935 ELSE
936 v(coset(ax, ay, az), 2, coc, n) = &
937 rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
938 rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
939 fx*(v(coset(ax - 1, ay, az), 1, coc, n) + &
940 f4*v(coset(ax - 1, ay, az), 1, coc, n + 1)) + &
941 fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
942 END IF
943
944 IF (ay == 0) THEN
945 v(coset(ax, ay, az), 3, coc, n) = &
946 rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
947 rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
948 fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
949 ELSE
950 v(coset(ax, ay, az), 3, coc, n) = &
951 rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
952 rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
953 fy*(v(coset(ax, ay - 1, az), 1, coc, n) + &
954 f4*v(coset(ax, ay - 1, az), 1, coc, n + 1)) + &
955 fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
956 END IF
957
958 IF (az == 0) THEN
959 v(coset(ax, ay, az), 4, coc, n) = &
960 rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
961 rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
962 fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
963 ELSE
964 v(coset(ax, ay, az), 4, coc, n) = &
965 rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
966 rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
967 fz*(v(coset(ax, ay, az - 1), 1, coc, n) + &
968 f4*v(coset(ax, ay, az - 1), 1, coc, n + 1)) + &
969 fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
970 END IF
971
972 END DO
973 END DO
974 END DO
975
976 ! *** Recurrence steps: [ap||c] -> [ab||c] ***
977
978 DO lb = 2, lb_max
979
980 ! *** Horizontal recurrence steps ***
981
982 ! *** [ab||c]{n} = [(a+1i)(b-1i)||c]{n} - ***
983 ! *** (Bi - Ai)*[a(b-1i)||c]{n} ***
984
985 la_start = max(0, la_min - 1)
986
987 DO la = la_start, la_max - 1
988 DO n = 1, nmax - la - lb - lc
989 DO ax = 0, la
990 DO ay = 0, la - ax
991 az = la - ax - ay
992
993 ! *** Shift of angular momentum component z ***
994
995 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
996 v(coset(ax, ay, az + 1), &
997 coset(0, 0, lb - 1), coc, n) - &
998 rab(3)*v(coset(ax, ay, az), &
999 coset(0, 0, lb - 1), coc, n)
1000
1001 ! *** Shift of angular momentum component y ***
1002
1003 DO by = 1, lb
1004 bz = lb - by
1005 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1006 v(coset(ax, ay + 1, az), &
1007 coset(0, by - 1, bz), coc, n) - &
1008 rab(2)*v(coset(ax, ay, az), &
1009 coset(0, by - 1, bz), coc, n)
1010 END DO
1011
1012 ! *** Shift of angular momentum component x ***
1013
1014 DO bx = 1, lb
1015 DO by = 0, lb - bx
1016 bz = lb - bx - by
1017 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1018 v(coset(ax + 1, ay, az), &
1019 coset(bx - 1, by, bz), coc, n) - &
1020 rab(1)*v(coset(ax, ay, az), &
1021 coset(bx - 1, by, bz), coc, n)
1022 END DO
1023 END DO
1024
1025 END DO
1026 END DO
1027 END DO
1028 END DO
1029
1030 ! *** Vertical recurrence step ***
1031
1032 ! *** [ab||c]{n} = (Pi - Bi)*[a(b-1i)||c]{n} + ***
1033 ! *** (Wi - Pi)*[a(b-1i)||c]{n+1} + ***
1034 ! *** f2*Ni(a)*( [(a-1i)(b-1i)||c]{n} + ***
1035 ! *** f4*[(a-1i)(b-1i)||c]{n+1}) ***
1036 ! *** f2*Ni(b-1i)*( [a(b-2i)||c]{n} + ***
1037 ! *** f4*[a(b-2i)||c]{n+1}) + ***
1038 ! *** f6*Ni(c)*[a(b-1i)||c-1i]{n+1}) ***
1039
1040 DO n = 1, nmax - la_max - lb - lc
1041 DO ax = 0, la_max
1042 fx = f2*real(ax, dp)
1043 DO ay = 0, la_max - ax
1044 fy = f2*real(ay, dp)
1045 az = la_max - ax - ay
1046 fz = f2*real(az, dp)
1047
1048 ! *** Shift of angular momentum component z from a to b ***
1049
1050 f3 = f2*real(lb - 1, dp)
1051
1052 IF (az == 0) THEN
1053 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1054 rbp(3)*v(coset(ax, ay, az), &
1055 coset(0, 0, lb - 1), coc, n) + &
1056 rpw(3)*v(coset(ax, ay, az), &
1057 coset(0, 0, lb - 1), coc, n + 1) + &
1058 f3*(v(coset(ax, ay, az), &
1059 coset(0, 0, lb - 2), coc, n) + &
1060 f4*v(coset(ax, ay, az), &
1061 coset(0, 0, lb - 2), coc, n + 1)) + &
1062 fcz*v(coset(ax, ay, az), &
1063 coset(0, 0, lb - 1), cocz, n + 1)
1064 ELSE
1065 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1066 rbp(3)*v(coset(ax, ay, az), &
1067 coset(0, 0, lb - 1), coc, n) + &
1068 rpw(3)*v(coset(ax, ay, az), &
1069 coset(0, 0, lb - 1), coc, n + 1) + &
1070 fz*(v(coset(ax, ay, az - 1), &
1071 coset(0, 0, lb - 1), coc, n) + &
1072 f4*v(coset(ax, ay, az - 1), &
1073 coset(0, 0, lb - 1), coc, n + 1)) + &
1074 f3*(v(coset(ax, ay, az), &
1075 coset(0, 0, lb - 2), coc, n) + &
1076 f4*v(coset(ax, ay, az), &
1077 coset(0, 0, lb - 2), coc, n + 1)) + &
1078 fcz*v(coset(ax, ay, az), &
1079 coset(0, 0, lb - 1), cocz, n + 1)
1080 END IF
1081
1082 ! *** Shift of angular momentum component y from a to b ***
1083
1084 IF (ay == 0) THEN
1085 bz = lb - 1
1086 v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1087 rbp(2)*v(coset(ax, ay, az), &
1088 coset(0, 0, bz), coc, n) + &
1089 rpw(2)*v(coset(ax, ay, az), &
1090 coset(0, 0, bz), coc, n + 1) + &
1091 fcy*v(coset(ax, ay, az), &
1092 coset(0, 0, bz), cocy, n + 1)
1093 DO by = 2, lb
1094 bz = lb - by
1095 f3 = f2*real(by - 1, dp)
1096 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1097 rbp(2)*v(coset(ax, ay, az), &
1098 coset(0, by - 1, bz), coc, n) + &
1099 rpw(2)*v(coset(ax, ay, az), &
1100 coset(0, by - 1, bz), coc, n + 1) + &
1101 f3*(v(coset(ax, ay, az), &
1102 coset(0, by - 2, bz), coc, n) + &
1103 f4*v(coset(ax, ay, az), &
1104 coset(0, by - 2, bz), coc, n + 1)) + &
1105 fcy*v(coset(ax, ay, az), &
1106 coset(0, by - 1, bz), cocy, n + 1)
1107 END DO
1108 ELSE
1109 bz = lb - 1
1110 v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1111 rbp(2)*v(coset(ax, ay, az), &
1112 coset(0, 0, bz), coc, n) + &
1113 rpw(2)*v(coset(ax, ay, az), &
1114 coset(0, 0, bz), coc, n + 1) + &
1115 fy*(v(coset(ax, ay - 1, az), &
1116 coset(0, 0, bz), coc, n) + &
1117 f4*v(coset(ax, ay - 1, az), &
1118 coset(0, 0, bz), coc, n + 1)) + &
1119 fcy*v(coset(ax, ay, az), &
1120 coset(0, 0, bz), cocy, n + 1)
1121 DO by = 2, lb
1122 bz = lb - by
1123 f3 = f2*real(by - 1, dp)
1124 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1125 rbp(2)*v(coset(ax, ay, az), &
1126 coset(0, by - 1, bz), coc, n) + &
1127 rpw(2)*v(coset(ax, ay, az), &
1128 coset(0, by - 1, bz), coc, n + 1) + &
1129 fy*(v(coset(ax, ay - 1, az), &
1130 coset(0, by - 1, bz), coc, n) + &
1131 f4*v(coset(ax, ay - 1, az), &
1132 coset(0, by - 1, bz), coc, n + 1)) + &
1133 f3*(v(coset(ax, ay, az), &
1134 coset(0, by - 2, bz), coc, n) + &
1135 f4*v(coset(ax, ay, az), &
1136 coset(0, by - 2, bz), coc, n + 1)) + &
1137 fcy*v(coset(ax, ay, az), &
1138 coset(0, by - 1, bz), cocy, n + 1)
1139 END DO
1140 END IF
1141
1142 ! *** Shift of angular momentum component x from a to b ***
1143
1144 IF (ax == 0) THEN
1145 DO by = 0, lb - 1
1146 bz = lb - 1 - by
1147 v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1148 rbp(1)*v(coset(ax, ay, az), &
1149 coset(0, by, bz), coc, n) + &
1150 rpw(1)*v(coset(ax, ay, az), &
1151 coset(0, by, bz), coc, n + 1) + &
1152 fcx*v(coset(ax, ay, az), &
1153 coset(0, by, bz), cocx, n + 1)
1154 END DO
1155 DO bx = 2, lb
1156 f3 = f2*real(bx - 1, dp)
1157 DO by = 0, lb - bx
1158 bz = lb - bx - by
1159 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1160 rbp(1)*v(coset(ax, ay, az), &
1161 coset(bx - 1, by, bz), coc, n) + &
1162 rpw(1)*v(coset(ax, ay, az), &
1163 coset(bx - 1, by, bz), coc, n + 1) + &
1164 f3*(v(coset(ax, ay, az), &
1165 coset(bx - 2, by, bz), coc, n) + &
1166 f4*v(coset(ax, ay, az), &
1167 coset(bx - 2, by, bz), coc, n + 1)) + &
1168 fcx*v(coset(ax, ay, az), &
1169 coset(bx - 1, by, bz), cocx, n + 1)
1170 END DO
1171 END DO
1172 ELSE
1173 DO by = 0, lb - 1
1174 bz = lb - 1 - by
1175 v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1176 rbp(1)*v(coset(ax, ay, az), &
1177 coset(0, by, bz), coc, n) + &
1178 rpw(1)*v(coset(ax, ay, az), &
1179 coset(0, by, bz), coc, n + 1) + &
1180 fx*(v(coset(ax - 1, ay, az), &
1181 coset(0, by, bz), coc, n) + &
1182 f4*v(coset(ax - 1, ay, az), &
1183 coset(0, by, bz), coc, n + 1)) + &
1184 fcx*v(coset(ax, ay, az), &
1185 coset(0, by, bz), cocx, n + 1)
1186 END DO
1187 DO bx = 2, lb
1188 f3 = f2*real(bx - 1, dp)
1189 DO by = 0, lb - bx
1190 bz = lb - bx - by
1191 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1192 rbp(1)*v(coset(ax, ay, az), &
1193 coset(bx - 1, by, bz), coc, n) + &
1194 rpw(1)*v(coset(ax, ay, az), &
1195 coset(bx - 1, by, bz), coc, n + 1) + &
1196 fx*(v(coset(ax - 1, ay, az), &
1197 coset(bx - 1, by, bz), coc, n) + &
1198 f4*v(coset(ax - 1, ay, az), &
1199 coset(bx - 1, by, bz), coc, n + 1)) + &
1200 f3*(v(coset(ax, ay, az), &
1201 coset(bx - 2, by, bz), coc, n) + &
1202 f4*v(coset(ax, ay, az), &
1203 coset(bx - 2, by, bz), coc, n + 1)) + &
1204 fcx*v(coset(ax, ay, az), &
1205 coset(bx - 1, by, bz), cocx, n + 1)
1206 END DO
1207 END DO
1208 END IF
1209
1210 END DO
1211 END DO
1212 END DO
1213
1214 END DO
1215 END IF
1216
1217 ELSE
1218
1219 IF (lb_max > 0) THEN
1220
1221 ! *** Vertical recurrence steps: [ss||c] -> [sb||c] ***
1222
1223 ! *** [sp||c]{n} = (Pi - Bi)*[ss||c]{n} + ***
1224 ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1225 ! *** f6*Ni(c)**[ss||c-1i]{n+1} ***
1226
1227 DO n = 1, nmax - 1 - lc
1228 v(1, 2, coc, n) = rbp(1)*v(1, 1, coc, n) + &
1229 rpw(1)*v(1, 1, coc, n + 1) + &
1230 fcx*v(1, 1, cocx, n + 1)
1231 v(1, 3, coc, n) = rbp(2)*v(1, 1, coc, n) + &
1232 rpw(2)*v(1, 1, coc, n + 1) + &
1233 fcy*v(1, 1, cocy, n + 1)
1234 v(1, 4, coc, n) = rbp(3)*v(1, 1, coc, n) + &
1235 rpw(3)*v(1, 1, coc, n + 1) + &
1236 fcz*v(1, 1, cocz, n + 1)
1237 END DO
1238
1239 ! *** [sb||c]{n} = (Pi - Bi)*[s(b-1i)||c]{n} + ***
1240 ! *** (Wi - Pi)*[s(b-1i)||c]{n+1} + ***
1241 ! *** f2*Ni(b-1i)*( [s(b-2i)||c]{n} + ***
1242 ! *** f4*[s(b-2i)||c]{n+1}) + ***
1243 ! *** f6*Ni(c)**[s(b-1i)||c-1i]{n+1} ***
1244
1245 DO lb = 2, lb_max
1246
1247 DO n = 1, nmax - lb - lc
1248
1249 ! *** Increase the angular momentum component z of b ***
1250
1251 v(1, coset(0, 0, lb), coc, n) = &
1252 rbp(3)*v(1, coset(0, 0, lb - 1), coc, n) + &
1253 rpw(3)*v(1, coset(0, 0, lb - 1), coc, n + 1) + &
1254 f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), coc, n) + &
1255 f4*v(1, coset(0, 0, lb - 2), coc, n + 1)) + &
1256 fcz*v(1, coset(0, 0, lb - 1), cocz, n + 1)
1257
1258 ! *** Increase the angular momentum component y of b ***
1259
1260 bz = lb - 1
1261 v(1, coset(0, 1, bz), coc, n) = &
1262 rbp(2)*v(1, coset(0, 0, bz), coc, n) + &
1263 rpw(2)*v(1, coset(0, 0, bz), coc, n + 1) + &
1264 fcy*v(1, coset(0, 0, bz), cocy, n + 1)
1265
1266 DO by = 2, lb
1267 f3 = f2*real(by - 1, dp)
1268 bz = lb - by
1269 v(1, coset(0, by, bz), coc, n) = &
1270 rbp(2)*v(1, coset(0, by - 1, bz), coc, n) + &
1271 rpw(2)*v(1, coset(0, by - 1, bz), coc, n + 1) + &
1272 f3*(v(1, coset(0, by - 2, bz), coc, n) + &
1273 f4*v(1, coset(0, by - 2, bz), coc, n + 1)) + &
1274 fcy*v(1, coset(0, by - 1, bz), cocy, n + 1)
1275 END DO
1276
1277 ! *** Increase the angular momentum component x of b ***
1278
1279 DO by = 0, lb - 1
1280 bz = lb - 1 - by
1281 v(1, coset(1, by, bz), coc, n) = &
1282 rbp(1)*v(1, coset(0, by, bz), coc, n) + &
1283 rpw(1)*v(1, coset(0, by, bz), coc, n + 1) + &
1284 fcx*v(1, coset(0, by, bz), cocx, n + 1)
1285 END DO
1286
1287 DO bx = 2, lb
1288 f3 = f2*real(bx - 1, dp)
1289 DO by = 0, lb - bx
1290 bz = lb - bx - by
1291 v(1, coset(bx, by, bz), coc, n) = &
1292 rbp(1)*v(1, coset(bx - 1, by, bz), coc, n) + &
1293 rpw(1)*v(1, coset(bx - 1, by, bz), coc, n + 1) + &
1294 f3*(v(1, coset(bx - 2, by, bz), coc, n) + &
1295 f4*v(1, coset(bx - 2, by, bz), coc, n + 1)) + &
1296 fcx*v(1, coset(bx - 1, by, bz), cocx, n + 1)
1297 END DO
1298 END DO
1299
1300 END DO
1301
1302 END DO
1303
1304 END IF
1305
1306 END IF
1307
1308 END DO
1309 END DO
1310
1311 END DO
1312
1313 END IF
1314
1315 ! *** Add the contribution of the current pair ***
1316 ! *** of primitive Gaussian-type functions ***
1317
1318 DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1319 kk = k - ncoset(lc_min - 1)
1320 DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1321 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
1322 vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1323 int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1324 END DO
1325 END DO
1326 END DO
1327
1328 IF (PRESENT(maxder)) THEN
1329 DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1330 kk = k - ncoset(lc_min - 1)
1331 DO j = 1, ncoset(lb_max)
1332 DO i = 1, ncoset(la_max)
1333 vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1334 END DO
1335 END DO
1336 END DO
1337 END IF
1338
1339 nb = nb + ncoset(lb_max)
1340
1341 END DO
1342
1343 na = na + ncoset(la_max - maxder_local)
1344 nap = nap + ncoset(la_max)
1345
1346 END DO
1347
1348 END SUBROUTINE coulomb3
1349
1350END MODULE ai_coulomb
ai_coulomb
Calculation of Coulomb integrals over Cartesian Gaussian-type functions (electron repulsion integrals...
Definition ai_coulomb.F:41
ai_coulomb::coulomb2
subroutine, public coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, rac, rac2, vac, v, f, screening)
Calculation of the primitive two-center Coulomb integrals over Cartesian Gaussian-type functions.
Definition ai_coulomb.F:86
ai_coulomb::coulomb3
subroutine, public coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, v, f, maxder, vabc_plus)
Calculation of the primitive three-center Coulomb integrals over Cartesian Gaussian-type functions (e...
Definition ai_coulomb.F:195
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::operator2_recurrence
subroutine, public operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta_a, zeta_c, zetp, zetq, zetw, rho, rac, vac, v, na, nc, nap, ncp, maxder, vac_plus)
Apply the common two-center OS recurrence to one primitive pair. The caller initializes v(1,...
Definition ai_operators_r12.F:198
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
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