(git:298f7c4)
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-2025 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 gamma, ONLY: fgamma => fgamma_0
44 USE kinds, ONLY: dp
45 USE mathconstants, ONLY: pi
46 USE orbital_pointers, ONLY: coset,&
47 ncoset
48#include "../base/base_uses.f90"
49
50 IMPLICIT NONE
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_coulomb'
52 PRIVATE
53
54 ! *** Public subroutines ***
55
56 PUBLIC :: coulomb2, coulomb2_new, coulomb3
57
58CONTAINS
59
60! **************************************************************************************************
61!> \brief Calculation of the primitive two-center Coulomb integrals over
62!> Cartesian Gaussian-type functions.
63!> \param la_max ...
64!> \param npgfa ...
65!> \param zeta ...
66!> \param rpgfa ...
67!> \param la_min ...
68!> \param lc_max ...
69!> \param npgfc ...
70!> \param zetc ...
71!> \param rpgfc ...
72!> \param lc_min ...
73!> \param rac ...
74!> \param rac2 ...
75!> \param vac ...
76!> \param v ...
77!> \param f ...
78!> \param maxder ...
79!> \param vac_plus ...
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, maxder, vac_plus)
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 INTEGER, INTENT(IN), OPTIONAL :: maxder
97 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
98
99 INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
100 cy, cz, i, ipgf, j, jpgf, la, lc, &
101 maxder_local, n, na, nap, nc, nmax
102 REAL(kind=dp) :: dac, f0, f1, f2, f3, f4, f5, f6, fcx, &
103 fcy, fcz, rho, t, zetp, zetq, zetw
104 REAL(kind=dp), DIMENSION(3) :: raw, rcw
105
106 v = 0.0_dp
107
108 maxder_local = 0
109 IF (PRESENT(maxder)) THEN
110 maxder_local = maxder
111 vac_plus = 0.0_dp
112 END IF
113
114 nmax = la_max + lc_max + 1
115
116 ! *** Calculate the distance of the centers a and c ***
117
118 dac = sqrt(rac2)
119
120 ! *** Loop over all pairs of primitive Gaussian-type functions ***
121
122 na = 0
123 nap = 0
124
125 DO ipgf = 1, npgfa
126
127 nc = 0
128
129 DO jpgf = 1, npgfc
130
131 ! *** Screening ***
132
133 IF (rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
134 DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
135 DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max - maxder_local)
136 vac(i, j) = 0.0_dp
137 END DO
138 END DO
139 nc = nc + ncoset(lc_max)
140 cycle
141 END IF
142
143 ! *** Calculate some prefactors ***
144
145 zetp = 1.0_dp/zeta(ipgf)
146 zetq = 1.0_dp/zetc(jpgf)
147 zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
148
149 rho = zeta(ipgf)*zetc(jpgf)*zetw
150
151 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
152
153 ! *** Calculate the incomplete Gamma function ***
154
155 t = rho*rac2
156
157 CALL fgamma(nmax - 1, t, f)
158
159 ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
160
161 DO n = 1, nmax
162 v(1, 1, n) = f0*f(n - 1)
163 END DO
164
165 ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
166
167 IF (lc_max > 0) THEN
168
169 f1 = 0.5_dp*zetq
170 f2 = -rho*zetq
171
172 rcw(:) = -zeta(ipgf)*zetw*rac(:)
173
174 ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
175
176 DO n = 1, nmax - 1
177 v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
178 v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
179 v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
180 END DO
181
182 ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
183 ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
184 ! ** f2*[s||c-2i]{n+1} ***
185
186 DO lc = 2, lc_max
187
188 DO n = 1, nmax - lc
189
190 ! **** Increase the angular momentum component z of c ***
191
192 v(1, coset(0, 0, lc), n) = &
193 rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
194 f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
195 f2*v(1, coset(0, 0, lc - 2), n + 1))
196
197 ! *** Increase the angular momentum component y of c ***
198
199 cz = lc - 1
200 v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
201
202 DO cy = 2, lc
203 cz = lc - cy
204 v(1, coset(0, cy, cz), n) = &
205 rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
206 f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
207 f2*v(1, coset(0, cy - 2, cz), n + 1))
208 END DO
209
210 ! *** Increase the angular momentum component x of c ***
211
212 DO cy = 0, lc - 1
213 cz = lc - 1 - cy
214 v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
215 END DO
216
217 DO cx = 2, lc
218 f6 = f1*real(cx - 1, dp)
219 DO cy = 0, lc - cx
220 cz = lc - cx - cy
221 v(1, coset(cx, cy, cz), n) = &
222 rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
223 f6*(v(1, coset(cx - 2, cy, cz), n) + &
224 f2*v(1, coset(cx - 2, cy, cz), n + 1))
225 END DO
226 END DO
227
228 END DO
229
230 END DO
231
232 END IF
233
234 ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
235
236 IF (la_max > 0) THEN
237
238 f3 = 0.5_dp*zetp
239 f4 = -rho*zetp
240 f5 = 0.5_dp*zetw
241
242 raw(:) = zetc(jpgf)*zetw*rac(:)
243
244 ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
245
246 DO n = 1, nmax - 1
247 v(2, 1, n) = raw(1)*v(1, 1, n + 1)
248 v(3, 1, n) = raw(2)*v(1, 1, n + 1)
249 v(4, 1, n) = raw(3)*v(1, 1, n + 1)
250 END DO
251
252 ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
253 ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
254 ! *** f4*[a-2i||s]{n+1}) ***
255
256 DO la = 2, la_max
257
258 DO n = 1, nmax - la
259
260 ! *** Increase the angular momentum component z of a ***
261
262 v(coset(0, 0, la), 1, n) = &
263 raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
264 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
265 f4*v(coset(0, 0, la - 2), 1, n + 1))
266
267 ! *** Increase the angular momentum component y of a ***
268
269 az = la - 1
270 v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
271
272 DO ay = 2, la
273 az = la - ay
274 v(coset(0, ay, az), 1, n) = &
275 raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
276 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
277 f4*v(coset(0, ay - 2, az), 1, n + 1))
278 END DO
279
280 ! *** Increase the angular momentum component x of a ***
281
282 DO ay = 0, la - 1
283 az = la - 1 - ay
284 v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
285 END DO
286
287 DO ax = 2, la
288 f6 = f3*real(ax - 1, dp)
289 DO ay = 0, la - ax
290 az = la - ax - ay
291 v(coset(ax, ay, az), 1, n) = &
292 raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
293 f6*(v(coset(ax - 2, ay, az), 1, n) + &
294 f4*v(coset(ax - 2, ay, az), 1, n + 1))
295 END DO
296 END DO
297
298 END DO
299
300 END DO
301
302 DO lc = 1, lc_max
303
304 DO cx = 0, lc
305 DO cy = 0, lc - cx
306 cz = lc - cx - cy
307
308 coc = coset(cx, cy, cz)
309 cocx = coset(max(0, cx - 1), cy, cz)
310 cocy = coset(cx, max(0, cy - 1), cz)
311 cocz = coset(cx, cy, max(0, cz - 1))
312
313 fcx = f5*real(cx, dp)
314 fcy = f5*real(cy, dp)
315 fcz = f5*real(cz, dp)
316
317 ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
318 ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
319
320 DO n = 1, nmax - 1 - lc
321 v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
322 v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
323 v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
324 END DO
325
326 ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
327 ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
328 ! *** f4*[a-2i||c]{n+1}) + ***
329 ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
330
331 DO la = 2, la_max
332
333 DO n = 1, nmax - la - lc
334
335 ! *** Increase the angular momentum component z of a ***
336
337 v(coset(0, 0, la), coc, n) = &
338 raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
339 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
340 f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
341 fcz*v(coset(0, 0, la - 1), cocz, n + 1)
342
343 ! *** Increase the angular momentum component y of a ***
344
345 az = la - 1
346 v(coset(0, 1, az), coc, n) = &
347 raw(2)*v(coset(0, 0, az), coc, n + 1) + &
348 fcy*v(coset(0, 0, az), cocy, n + 1)
349
350 DO ay = 2, la
351 az = la - ay
352 v(coset(0, ay, az), coc, n) = &
353 raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
354 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
355 f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
356 fcy*v(coset(0, ay - 1, az), cocy, n + 1)
357 END DO
358
359 ! *** Increase the angular momentum component x of a ***
360
361 DO ay = 0, la - 1
362 az = la - 1 - ay
363 v(coset(1, ay, az), coc, n) = &
364 raw(1)*v(coset(0, ay, az), coc, n + 1) + &
365 fcx*v(coset(0, ay, az), cocx, n + 1)
366 END DO
367
368 DO ax = 2, la
369 f6 = f3*real(ax - 1, dp)
370 DO ay = 0, la - ax
371 az = la - ax - ay
372 v(coset(ax, ay, az), coc, n) = &
373 raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
374 f6*(v(coset(ax - 2, ay, az), coc, n) + &
375 f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
376 fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
377 END DO
378 END DO
379
380 END DO
381
382 END DO
383
384 END DO
385 END DO
386
387 END DO
388
389 END IF
390
391 DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
392 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
393 vac(na + i, nc + j) = v(i, j, 1)
394 END DO
395 END DO
396
397 IF (PRESENT(maxder)) THEN
398 DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
399 DO i = 1, ncoset(la_max)
400 vac_plus(nap + i, nc + j) = v(i, j, 1)
401 END DO
402 END DO
403 END IF
404
405 nc = nc + ncoset(lc_max)
406
407 END DO
408
409 na = na + ncoset(la_max - maxder_local)
410 nap = nap + ncoset(la_max)
411
412 END DO
413
414 END SUBROUTINE coulomb2
415! **************************************************************************************************
416!> \brief Calculation of the primitive two-center Coulomb integrals over
417!> Cartesian Gaussian-type functions. Same as coulomb2 treats derivative
418!> different (now vac_plus is symmetric)
419!> \param la_max ...
420!> \param npgfa ...
421!> \param zeta ...
422!> \param la_min ...
423!> \param lc_max ...
424!> \param npgfc ...
425!> \param zetc ...
426!> \param lc_min ...
427!> \param rac ...
428!> \param rac2 ...
429!> \param vac ...
430!> \param v ...
431!> \param f ...
432!> \param maxder ...
433!> \param vac_plus ...
434!> \date 6.02.2008
435!> \author CJM
436!> \version 1.0
437! **************************************************************************************************
438
439 SUBROUTINE coulomb2_new(la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
440 rac, rac2, vac, v, f, maxder, vac_plus)
441 INTEGER, INTENT(IN) :: la_max, npgfa
442 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
443 INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
444 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc
445 INTEGER, INTENT(IN) :: lc_min
446 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
447 REAL(kind=dp), INTENT(IN) :: rac2
448 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
449 REAL(kind=dp), DIMENSION(:, :, :) :: v
450 REAL(kind=dp), DIMENSION(0:) :: f
451 INTEGER, INTENT(IN), OPTIONAL :: maxder
452 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
453
454 INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
455 cy, cz, i, ipgf, j, jpgf, la, lc, &
456 maxder_local, n, na, nap, nc, ncp, nmax
457 REAL(kind=dp) :: dac, f0, f1, f2, f3, f4, f5, f6, fcx, &
458 fcy, fcz, rho, t, zetp, zetq, zetw
459 REAL(kind=dp), DIMENSION(3) :: raw, rcw
460
461 v = 0.0_dp
462
463 maxder_local = 0
464 IF (PRESENT(maxder)) THEN
465 maxder_local = maxder
466 vac_plus = 0.0_dp
467 END IF
468
469 nmax = la_max + lc_max + 1
470
471 ! *** Calculate the distance of the centers a and c ***
472
473 dac = sqrt(rac2)
474
475 ! *** Loop over all pairs of primitive Gaussian-type functions ***
476
477 na = 0
478 nap = 0
479
480 DO ipgf = 1, npgfa
481
482 nc = 0
483 ncp = 0
484
485 DO jpgf = 1, npgfc
486
487 ! *** Calculate some prefactors ***
488
489 zetp = 1.0_dp/zeta(ipgf)
490 zetq = 1.0_dp/zetc(jpgf)
491 zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
492
493 rho = zeta(ipgf)*zetc(jpgf)*zetw
494
495 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
496
497 ! *** Calculate the incomplete Gamma function ***
498
499 t = rho*rac2
500
501 CALL fgamma(nmax - 1, t, f)
502
503 ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
504
505 DO n = 1, nmax
506 v(1, 1, n) = f0*f(n - 1)
507 END DO
508
509 ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
510
511 IF (lc_max > 0) THEN
512
513 f1 = 0.5_dp*zetq
514 f2 = -rho*zetq
515
516 rcw(:) = -zeta(ipgf)*zetw*rac(:)
517
518 ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
519
520 DO n = 1, nmax - 1
521 v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
522 v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
523 v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
524 END DO
525
526 ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
527 ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
528 ! ** f2*[s||c-2i]{n+1} ***
529
530 DO lc = 2, lc_max
531
532 DO n = 1, nmax - lc
533
534 ! **** Increase the angular momentum component z of c ***
535
536 v(1, coset(0, 0, lc), n) = &
537 rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
538 f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
539 f2*v(1, coset(0, 0, lc - 2), n + 1))
540
541 ! *** Increase the angular momentum component y of c ***
542
543 cz = lc - 1
544 v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
545
546 DO cy = 2, lc
547 cz = lc - cy
548 v(1, coset(0, cy, cz), n) = &
549 rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
550 f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
551 f2*v(1, coset(0, cy - 2, cz), n + 1))
552 END DO
553
554 ! *** Increase the angular momentum component x of c ***
555
556 DO cy = 0, lc - 1
557 cz = lc - 1 - cy
558 v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
559 END DO
560
561 DO cx = 2, lc
562 f6 = f1*real(cx - 1, dp)
563 DO cy = 0, lc - cx
564 cz = lc - cx - cy
565 v(1, coset(cx, cy, cz), n) = &
566 rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
567 f6*(v(1, coset(cx - 2, cy, cz), n) + &
568 f2*v(1, coset(cx - 2, cy, cz), n + 1))
569 END DO
570 END DO
571
572 END DO
573
574 END DO
575
576 END IF
577
578 ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
579
580 IF (la_max > 0) THEN
581
582 f3 = 0.5_dp*zetp
583 f4 = -rho*zetp
584 f5 = 0.5_dp*zetw
585
586 raw(:) = zetc(jpgf)*zetw*rac(:)
587
588 ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
589
590 DO n = 1, nmax - 1
591 v(2, 1, n) = raw(1)*v(1, 1, n + 1)
592 v(3, 1, n) = raw(2)*v(1, 1, n + 1)
593 v(4, 1, n) = raw(3)*v(1, 1, n + 1)
594 END DO
595
596 ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
597 ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
598 ! *** f4*[a-2i||s]{n+1}) ***
599
600 DO la = 2, la_max
601
602 DO n = 1, nmax - la
603
604 ! *** Increase the angular momentum component z of a ***
605
606 v(coset(0, 0, la), 1, n) = &
607 raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
608 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
609 f4*v(coset(0, 0, la - 2), 1, n + 1))
610
611 ! *** Increase the angular momentum component y of a ***
612
613 az = la - 1
614 v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
615
616 DO ay = 2, la
617 az = la - ay
618 v(coset(0, ay, az), 1, n) = &
619 raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
620 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
621 f4*v(coset(0, ay - 2, az), 1, n + 1))
622 END DO
623
624 ! *** Increase the angular momentum component x of a ***
625
626 DO ay = 0, la - 1
627 az = la - 1 - ay
628 v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
629 END DO
630
631 DO ax = 2, la
632 f6 = f3*real(ax - 1, dp)
633 DO ay = 0, la - ax
634 az = la - ax - ay
635 v(coset(ax, ay, az), 1, n) = &
636 raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
637 f6*(v(coset(ax - 2, ay, az), 1, n) + &
638 f4*v(coset(ax - 2, ay, az), 1, n + 1))
639 END DO
640 END DO
641
642 END DO
643
644 END DO
645
646 DO lc = 1, lc_max
647
648 DO cx = 0, lc
649 DO cy = 0, lc - cx
650 cz = lc - cx - cy
651
652 coc = coset(cx, cy, cz)
653 cocx = coset(max(0, cx - 1), cy, cz)
654 cocy = coset(cx, max(0, cy - 1), cz)
655 cocz = coset(cx, cy, max(0, cz - 1))
656
657 fcx = f5*real(cx, dp)
658 fcy = f5*real(cy, dp)
659 fcz = f5*real(cz, dp)
660
661 ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
662 ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
663
664 DO n = 1, nmax - 1 - lc
665 v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
666 v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
667 v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
668 END DO
669
670 ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
671 ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
672 ! *** f4*[a-2i||c]{n+1}) + ***
673 ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
674
675 DO la = 2, la_max
676
677 DO n = 1, nmax - la - lc
678
679 ! *** Increase the angular momentum component z of a ***
680
681 v(coset(0, 0, la), coc, n) = &
682 raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
683 f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
684 f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
685 fcz*v(coset(0, 0, la - 1), cocz, n + 1)
686
687 ! *** Increase the angular momentum component y of a ***
688
689 az = la - 1
690 v(coset(0, 1, az), coc, n) = &
691 raw(2)*v(coset(0, 0, az), coc, n + 1) + &
692 fcy*v(coset(0, 0, az), cocy, n + 1)
693
694 DO ay = 2, la
695 az = la - ay
696 v(coset(0, ay, az), coc, n) = &
697 raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
698 f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
699 f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
700 fcy*v(coset(0, ay - 1, az), cocy, n + 1)
701 END DO
702
703 ! *** Increase the angular momentum component x of a ***
704
705 DO ay = 0, la - 1
706 az = la - 1 - ay
707 v(coset(1, ay, az), coc, n) = &
708 raw(1)*v(coset(0, ay, az), coc, n + 1) + &
709 fcx*v(coset(0, ay, az), cocx, n + 1)
710 END DO
711
712 DO ax = 2, la
713 f6 = f3*real(ax - 1, dp)
714 DO ay = 0, la - ax
715 az = la - ax - ay
716 v(coset(ax, ay, az), coc, n) = &
717 raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
718 f6*(v(coset(ax - 2, ay, az), coc, n) + &
719 f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
720 fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
721 END DO
722 END DO
723
724 END DO
725
726 END DO
727
728 END DO
729 END DO
730
731 END DO
732
733 END IF
734
735 DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
736 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
737 vac(na + i, nc + j) = v(i, j, 1)
738 END DO
739 END DO
740
741 IF (PRESENT(maxder)) THEN
742 DO j = 1, ncoset(lc_max)
743 DO i = 1, ncoset(la_max)
744 vac_plus(nap + i, ncp + j) = v(i, j, 1)
745 END DO
746 END DO
747 END IF
748
749 nc = nc + ncoset(lc_max - maxder_local)
750 ncp = ncp + ncoset(lc_max)
751
752 END DO
753
754 na = na + ncoset(la_max - maxder_local)
755 nap = nap + ncoset(la_max)
756
757 END DO
758
759 END SUBROUTINE coulomb2_new
760
761! **************************************************************************************************
762!> \brief Calculation of the primitive three-center Coulomb integrals over
763!> Cartesian Gaussian-type functions (electron repulsion integrals,
764!> ERIs).
765!> \param la_max ...
766!> \param npgfa ...
767!> \param zeta ...
768!> \param rpgfa ...
769!> \param la_min ...
770!> \param lb_max ...
771!> \param npgfb ...
772!> \param zetb ...
773!> \param rpgfb ...
774!> \param lb_min ...
775!> \param lc_max ...
776!> \param zetc ...
777!> \param rpgfc ...
778!> \param lc_min ...
779!> \param gccc ...
780!> \param rab ...
781!> \param rab2 ...
782!> \param rac ...
783!> \param rac2 ...
784!> \param rbc2 ...
785!> \param vabc ...
786!> \param int_abc ...
787!> \param v ...
788!> \param f ...
789!> \param maxder ...
790!> \param vabc_plus ...
791!> \date 06.11.2000
792!> \author Matthias Krack
793!> \version 1.0
794! **************************************************************************************************
795 SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
796 lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
797 v, f, maxder, vabc_plus)
798 INTEGER, INTENT(IN) :: la_max, npgfa
799 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
800 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
801 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
802 INTEGER, INTENT(IN) :: lb_min, lc_max
803 REAL(kind=dp), INTENT(IN) :: zetc, rpgfc
804 INTEGER, INTENT(IN) :: lc_min
805 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: gccc
806 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
807 REAL(kind=dp), INTENT(IN) :: rab2
808 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
809 REAL(kind=dp), INTENT(IN) :: rac2, rbc2
810 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
811 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: int_abc
812 REAL(kind=dp), DIMENSION(:, :, :, :) :: v
813 REAL(kind=dp), DIMENSION(0:) :: f
814 INTEGER, INTENT(IN), OPTIONAL :: maxder
815 REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus
816
817 INTEGER :: ax, ay, az, bx, by, bz, coc, cocx, cocy, &
818 cocz, cx, cy, cz, i, ipgf, j, jpgf, k, &
819 kk, la, la_start, lb, lc, &
820 maxder_local, n, na, nap, nb, nmax
821 REAL(kind=dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, f5, &
822 f6, f7, fcx, fcy, fcz, fx, fy, fz, t, &
823 zetp, zetq, zetw
824 REAL(kind=dp), DIMENSION(3) :: rap, rbp, rcp, rcw, rpw
825
826 v = 0.0_dp
827
828 maxder_local = 0
829 IF (PRESENT(maxder)) THEN
830 maxder_local = maxder
831 END IF
832
833 nmax = la_max + lb_max + lc_max + 1
834
835 ! *** Calculate the distances of the centers a, b and c ***
836
837 dab = sqrt(rab2)
838 dac = sqrt(rac2)
839 dbc = sqrt(rbc2)
840
841 ! *** Initialize integrals array
842 int_abc = 0.0_dp
843
844 ! *** Loop over all pairs of primitive Gaussian-type functions ***
845
846 na = 0
847 nap = 0
848
849 DO ipgf = 1, npgfa
850
851 ! *** Screening ***
852 IF (rpgfa(ipgf) + rpgfc < dac) THEN
853 na = na + ncoset(la_max - maxder_local)
854 nap = nap + ncoset(la_max)
855 cycle
856 END IF
857
858 nb = 0
859
860 DO jpgf = 1, npgfb
861
862 ! *** Screening ***
863 IF ( &
864 (rpgfb(jpgf) + rpgfc < dbc) .OR. &
865 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
866 nb = nb + ncoset(lb_max)
867 cycle
868 END IF
869
870 ! *** Calculate some prefactors ***
871
872 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
873 zetq = 1.0_dp/zetc
874 zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
875
876 f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
877 f1 = zetb(jpgf)*zetp
878 f2 = 0.5_dp*zetp
879 f4 = -zetc*zetw
880
881 f0 = f0*exp(-zeta(ipgf)*f1*rab2)
882
883 rap(:) = f1*rab(:)
884 rcp(:) = rap(:) - rac(:)
885 rpw(:) = f4*rcp(:)
886
887 ! *** Calculate the incomplete Gamma function ***
888
889 t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
890
891 CALL fgamma(nmax - 1, t, f)
892
893 ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
894
895 DO n = 1, nmax
896 v(1, 1, 1, n) = f0*f(n - 1)
897 END DO
898
899 ! *** Recurrence steps: [ss||s] -> [as||s] ***
900
901 IF (la_max > 0) THEN
902
903 ! *** Vertical recurrence steps: [ss||s] -> [as||s] ***
904
905 ! *** [ps||s]{n} = (Pi - Ai)*[ss||s]{n} + ***
906 ! *** (Wi - Pi)*[ss||s]{n+1} (i = x,y,z) ***
907
908 DO n = 1, nmax - 1
909 v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
910 v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
911 v(4, 1, 1, n) = rap(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
912 END DO
913
914 ! *** [as||s]{n} = (Pi - Ai)*[(a-1i)s||s]{n} + ***
915 ! *** (Wi - Pi)*[(a-1i)s||s]{n+1} + ***
916 ! *** f2*Ni(a-1i)*( [(a-2i)s||s]{n} + ***
917 ! *** f4*[(a-2i)s||s]{n+1}) ***
918
919 DO la = 2, la_max
920
921 DO n = 1, nmax - la
922
923 ! *** Increase the angular momentum component z of a ***
924
925 v(coset(0, 0, la), 1, 1, n) = &
926 rap(3)*v(coset(0, 0, la - 1), 1, 1, n) + &
927 rpw(3)*v(coset(0, 0, la - 1), 1, 1, n + 1) + &
928 f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, 1, n) + &
929 f4*v(coset(0, 0, la - 2), 1, 1, n + 1))
930
931 ! *** Increase the angular momentum component y of a ***
932
933 az = la - 1
934 v(coset(0, 1, az), 1, 1, n) = &
935 rap(2)*v(coset(0, 0, az), 1, 1, n) + &
936 rpw(2)*v(coset(0, 0, az), 1, 1, n + 1)
937
938 DO ay = 2, la
939 az = la - ay
940 v(coset(0, ay, az), 1, 1, n) = &
941 rap(2)*v(coset(0, ay - 1, az), 1, 1, n) + &
942 rpw(2)*v(coset(0, ay - 1, az), 1, 1, n + 1) + &
943 f2*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, 1, n) + &
944 f4*v(coset(0, ay - 2, az), 1, 1, n + 1))
945 END DO
946
947 ! *** Increase the angular momentum component x of a ***
948
949 DO ay = 0, la - 1
950 az = la - 1 - ay
951 v(coset(1, ay, az), 1, 1, n) = &
952 rap(1)*v(coset(0, ay, az), 1, 1, n) + &
953 rpw(1)*v(coset(0, ay, az), 1, 1, n + 1)
954 END DO
955
956 DO ax = 2, la
957 f3 = f2*real(ax - 1, dp)
958 DO ay = 0, la - ax
959 az = la - ax - ay
960 v(coset(ax, ay, az), 1, 1, n) = &
961 rap(1)*v(coset(ax - 1, ay, az), 1, 1, n) + &
962 rpw(1)*v(coset(ax - 1, ay, az), 1, 1, n + 1) + &
963 f3*(v(coset(ax - 2, ay, az), 1, 1, n) + &
964 f4*v(coset(ax - 2, ay, az), 1, 1, n + 1))
965 END DO
966 END DO
967
968 END DO
969
970 END DO
971
972 ! *** Recurrence steps: [as||s] -> [ab||s] ***
973
974 IF (lb_max > 0) THEN
975
976 ! *** Horizontal recurrence steps ***
977
978 rbp(:) = rap(:) - rab(:)
979
980 ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
981
982 la_start = max(0, la_min - 1)
983
984 DO la = la_start, la_max - 1
985 DO n = 1, nmax - la - 1
986 DO ax = 0, la
987 DO ay = 0, la - ax
988 az = la - ax - ay
989 v(coset(ax, ay, az), 2, 1, n) = &
990 v(coset(ax + 1, ay, az), 1, 1, n) - &
991 rab(1)*v(coset(ax, ay, az), 1, 1, n)
992 v(coset(ax, ay, az), 3, 1, n) = &
993 v(coset(ax, ay + 1, az), 1, 1, n) - &
994 rab(2)*v(coset(ax, ay, az), 1, 1, n)
995 v(coset(ax, ay, az), 4, 1, n) = &
996 v(coset(ax, ay, az + 1), 1, 1, n) - &
997 rab(3)*v(coset(ax, ay, az), 1, 1, n)
998 END DO
999 END DO
1000 END DO
1001 END DO
1002
1003 ! *** Vertical recurrence step ***
1004
1005 ! *** [ap||s]{n} = (Pi - Bi)*[as||s]{n} + ***
1006 ! *** (Wi - Pi)*[as||s]{n+1} + ***
1007 ! *** f2*Ni(a)*( [(a-1i)s||s]{n} + ***
1008 ! *** f4*[(a-1i)s||s]{n+1}) ***
1009
1010 DO n = 1, nmax - la_max - 1
1011 DO ax = 0, la_max
1012 fx = f2*real(ax, dp)
1013 DO ay = 0, la_max - ax
1014 fy = f2*real(ay, dp)
1015 az = la_max - ax - ay
1016 fz = f2*real(az, dp)
1017
1018 IF (ax == 0) THEN
1019 v(coset(ax, ay, az), 2, 1, n) = &
1020 rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
1021 rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1)
1022 ELSE
1023 v(coset(ax, ay, az), 2, 1, n) = &
1024 rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
1025 rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1026 fx*(v(coset(ax - 1, ay, az), 1, 1, n) + &
1027 f4*v(coset(ax - 1, ay, az), 1, 1, n + 1))
1028 END IF
1029
1030 IF (ay == 0) THEN
1031 v(coset(ax, ay, az), 3, 1, n) = &
1032 rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
1033 rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1)
1034 ELSE
1035 v(coset(ax, ay, az), 3, 1, n) = &
1036 rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
1037 rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1038 fy*(v(coset(ax, ay - 1, az), 1, 1, n) + &
1039 f4*v(coset(ax, ay - 1, az), 1, 1, n + 1))
1040 END IF
1041
1042 IF (az == 0) THEN
1043 v(coset(ax, ay, az), 4, 1, n) = &
1044 rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
1045 rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1)
1046 ELSE
1047 v(coset(ax, ay, az), 4, 1, n) = &
1048 rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
1049 rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1050 fz*(v(coset(ax, ay, az - 1), 1, 1, n) + &
1051 f4*v(coset(ax, ay, az - 1), 1, 1, n + 1))
1052 END IF
1053
1054 END DO
1055 END DO
1056 END DO
1057
1058 ! *** Recurrence steps: [ap||s] -> [ab||s] ***
1059
1060 DO lb = 2, lb_max
1061
1062 ! *** Horizontal recurrence steps ***
1063
1064 ! *** [ab||s]{n} = [(a+1i)(b-1i)||s]{n} - ***
1065 ! *** (Bi - Ai)*[a(b-1i)||s]{n} ***
1066
1067 la_start = max(0, la_min - 1)
1068
1069 DO la = la_start, la_max - 1
1070 DO n = 1, nmax - la - lb
1071 DO ax = 0, la
1072 DO ay = 0, la - ax
1073 az = la - ax - ay
1074
1075 ! *** Shift of angular momentum component z from a to b ***
1076
1077 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1078 v(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1, n) - &
1079 rab(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n)
1080
1081 ! *** Shift of angular momentum component y from a to b ***
1082
1083 DO by = 1, lb
1084 bz = lb - by
1085 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1086 v(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1, n) - &
1087 rab(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n)
1088 END DO
1089
1090 ! *** Shift of angular momentum component x from a to b ***
1091
1092 DO bx = 1, lb
1093 DO by = 0, lb - bx
1094 bz = lb - bx - by
1095 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1096 v(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1, n) - &
1097 rab(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n)
1098 END DO
1099 END DO
1100
1101 END DO
1102 END DO
1103 END DO
1104 END DO
1105
1106 ! *** Vertical recurrence step ***
1107
1108 ! *** [ab||s]{n} = (Pi - Bi)*[a(b-1i)||s]{n} + ***
1109 ! *** (Wi - Pi)*[a(b-1i)||s]{n+1} + ***
1110 ! *** f2*Ni(a)*( [(a-1i)(b-1i)||s]{n} + ***
1111 ! *** f4*[(a-1i)(b-1i)||s]{n+1}) ***
1112 ! *** f2*Ni(b-1i)*( [a(b-2i)||s]{n} + ***
1113 ! *** f4*[a(b-2i)||s]{n+1}) ***
1114
1115 DO n = 1, nmax - la_max - lb
1116 DO ax = 0, la_max
1117 fx = f2*real(ax, dp)
1118 DO ay = 0, la_max - ax
1119 fy = f2*real(ay, dp)
1120 az = la_max - ax - ay
1121 fz = f2*real(az, dp)
1122
1123 ! *** Shift of angular momentum component z from a to b ***
1124
1125 f3 = f2*real(lb - 1, dp)
1126
1127 IF (az == 0) THEN
1128 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1129 rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
1130 rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
1131 f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
1132 f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
1133 ELSE
1134 v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1135 rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
1136 rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
1137 fz*(v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n) + &
1138 f4*v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n + 1)) + &
1139 f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
1140 f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
1141 END IF
1142
1143 ! *** Shift of angular momentum component y from a to b ***
1144
1145 IF (ay == 0) THEN
1146 bz = lb - 1
1147 v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
1148 rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
1149 rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1)
1150 DO by = 2, lb
1151 bz = lb - by
1152 f3 = f2*real(by - 1, dp)
1153 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1154 rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
1155 rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
1156 f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
1157 f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
1158 END DO
1159 ELSE
1160 bz = lb - 1
1161 v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
1162 rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
1163 rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1) + &
1164 fy*(v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n) + &
1165 f4*v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n + 1))
1166 DO by = 2, lb
1167 bz = lb - by
1168 f3 = f2*real(by - 1, dp)
1169 v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1170 rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
1171 rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
1172 fy*(v(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1, n) + &
1173 f4*v(coset(ax, ay - 1, az), &
1174 coset(0, by - 1, bz), 1, n + 1)) + &
1175 f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
1176 f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
1177 END DO
1178 END IF
1179
1180 ! *** Shift of angular momentum component x from a to b ***
1181
1182 IF (ax == 0) THEN
1183 DO by = 0, lb - 1
1184 bz = lb - 1 - by
1185 v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
1186 rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
1187 rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1)
1188 END DO
1189 DO bx = 2, lb
1190 f3 = f2*real(bx - 1, dp)
1191 DO by = 0, lb - bx
1192 bz = lb - bx - by
1193 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1194 rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
1195 rpw(1)*v(coset(ax, ay, az), &
1196 coset(bx - 1, by, bz), 1, n + 1) + &
1197 f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
1198 f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
1199 END DO
1200 END DO
1201 ELSE
1202 DO by = 0, lb - 1
1203 bz = lb - 1 - by
1204 v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
1205 rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
1206 rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1) + &
1207 fx*(v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n) + &
1208 f4*v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n + 1))
1209 END DO
1210 DO bx = 2, lb
1211 f3 = f2*real(bx - 1, dp)
1212 DO by = 0, lb - bx
1213 bz = lb - bx - by
1214 v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1215 rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
1216 rpw(1)*v(coset(ax, ay, az), &
1217 coset(bx - 1, by, bz), 1, n + 1) + &
1218 fx*(v(coset(ax - 1, ay, az), &
1219 coset(bx - 1, by, bz), 1, n) + &
1220 f4*v(coset(ax - 1, ay, az), &
1221 coset(bx - 1, by, bz), 1, n + 1)) + &
1222 f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
1223 f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
1224 END DO
1225 END DO
1226 END IF
1227
1228 END DO
1229 END DO
1230 END DO
1231
1232 END DO
1233
1234 END IF
1235
1236 ELSE
1237
1238 IF (lb_max > 0) THEN
1239
1240 ! *** Vertical recurrence steps: [ss||s] -> [sb||s] ***
1241
1242 rbp(:) = rap(:) - rab(:)
1243
1244 ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
1245 ! *** (Wi - Pi)*[ss||s]{n+1} ***
1246
1247 DO n = 1, nmax - 1
1248 v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
1249 v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
1250 v(1, 4, 1, n) = rbp(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
1251 END DO
1252
1253 ! *** [sb||s]{n} = (Pi - Bi)*[s(b-1i)||s]{n} + ***
1254 ! *** (Wi - Pi)*[s(b-1i)||s]{n+1} + ***
1255 ! *** f2*Ni(b-1i)*( [s(b-2i)||s]{n} + ***
1256 ! *** f4*[s(b-2i)||s]{n+1}) ***
1257
1258 DO lb = 2, lb_max
1259
1260 DO n = 1, nmax - lb
1261
1262 ! *** Increase the angular momentum component z of b ***
1263
1264 v(1, coset(0, 0, lb), 1, n) = &
1265 rbp(3)*v(1, coset(0, 0, lb - 1), 1, n) + &
1266 rpw(3)*v(1, coset(0, 0, lb - 1), 1, n + 1) + &
1267 f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), 1, n) + &
1268 f4*v(1, coset(0, 0, lb - 2), 1, n + 1))
1269
1270 ! *** Increase the angular momentum component y of b ***
1271
1272 bz = lb - 1
1273 v(1, coset(0, 1, bz), 1, n) = &
1274 rbp(2)*v(1, coset(0, 0, bz), 1, n) + &
1275 rpw(2)*v(1, coset(0, 0, bz), 1, n + 1)
1276
1277 DO by = 2, lb
1278 bz = lb - by
1279 v(1, coset(0, by, bz), 1, n) = &
1280 rbp(2)*v(1, coset(0, by - 1, bz), 1, n) + &
1281 rpw(2)*v(1, coset(0, by - 1, bz), 1, n + 1) + &
1282 f2*real(by - 1, dp)*(v(1, coset(0, by - 2, bz), 1, n) + &
1283 f4*v(1, coset(0, by - 2, bz), 1, n + 1))
1284 END DO
1285
1286 ! *** Increase the angular momentum component x of b ***
1287
1288 DO by = 0, lb - 1
1289 bz = lb - 1 - by
1290 v(1, coset(1, by, bz), 1, n) = &
1291 rbp(1)*v(1, coset(0, by, bz), 1, n) + &
1292 rpw(1)*v(1, coset(0, by, bz), 1, n + 1)
1293 END DO
1294
1295 DO bx = 2, lb
1296 f3 = f2*real(bx - 1, dp)
1297 DO by = 0, lb - bx
1298 bz = lb - bx - by
1299 v(1, coset(bx, by, bz), 1, n) = &
1300 rbp(1)*v(1, coset(bx - 1, by, bz), 1, n) + &
1301 rpw(1)*v(1, coset(bx - 1, by, bz), 1, n + 1) + &
1302 f3*(v(1, coset(bx - 2, by, bz), 1, n) + &
1303 f4*v(1, coset(bx - 2, by, bz), 1, n + 1))
1304 END DO
1305 END DO
1306
1307 END DO
1308
1309 END DO
1310
1311 END IF
1312
1313 END IF
1314
1315 ! *** Recurrence steps: [ab||s] -> [ab||c] ***
1316
1317 IF (lc_max > 0) THEN
1318
1319 ! *** Vertical recurrence steps: [ss||s] -> [ss||c] ***
1320
1321 f5 = -zetw/zetp
1322 f6 = 0.5_dp*zetw
1323 f7 = 0.5_dp*zetq
1324
1325 rcw(:) = rcp(:) + rpw(:)
1326
1327 ! *** [ss||p]{n} = (Wi - Ci)*[ss||s]{n+1} (i = x,y,z) ***
1328
1329 DO n = 1, nmax - 1
1330 v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
1331 v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
1332 v(1, 1, 4, n) = rcw(3)*v(1, 1, 1, n + 1)
1333 END DO
1334
1335 ! *** [ss||c]{n} = (Wi - Ci)*[ss||c-1i]{n+1} + ***
1336 ! *** f7*Ni(c-1i)*[ss||c-2i]{n} + ***
1337 ! *** f5*[ss||c-2i]{n+1} ***
1338
1339 DO lc = 2, lc_max
1340
1341 DO n = 1, nmax - lc
1342
1343 ! *** Increase the angular momentum component z of c ***
1344
1345 v(1, 1, coset(0, 0, lc), n) = &
1346 rcw(3)*v(1, 1, coset(0, 0, lc - 1), n + 1) + &
1347 f7*real(lc - 1, dp)*(v(1, 1, coset(0, 0, lc - 2), n) + &
1348 f5*v(1, 1, coset(0, 0, lc - 2), n + 1))
1349
1350 ! *** Increase the angular momentum component y of c ***
1351
1352 cz = lc - 1
1353 v(1, 1, coset(0, 1, cz), n) = rcw(2)*v(1, 1, coset(0, 0, cz), n + 1)
1354
1355 DO cy = 2, lc
1356 cz = lc - cy
1357 v(1, 1, coset(0, cy, cz), n) = &
1358 rcw(2)*v(1, 1, coset(0, cy - 1, cz), n + 1) + &
1359 f7*real(cy - 1, dp)*(v(1, 1, coset(0, cy - 2, cz), n) + &
1360 f5*v(1, 1, coset(0, cy - 2, cz), n + 1))
1361 END DO
1362
1363 ! *** Increase the angular momentum component x of c ***
1364
1365 DO cy = 0, lc - 1
1366 cz = lc - 1 - cy
1367 v(1, 1, coset(1, cy, cz), n) = rcw(1)*v(1, 1, coset(0, cy, cz), n + 1)
1368 END DO
1369
1370 DO cx = 2, lc
1371 DO cy = 0, lc - cx
1372 cz = lc - cx - cy
1373 v(1, 1, coset(cx, cy, cz), n) = &
1374 rcw(1)*v(1, 1, coset(cx - 1, cy, cz), n + 1) + &
1375 f7*real(cx - 1, dp)*(v(1, 1, coset(cx - 2, cy, cz), n) + &
1376 f5*v(1, 1, coset(cx - 2, cy, cz), n + 1))
1377 END DO
1378 END DO
1379
1380 END DO
1381
1382 END DO
1383
1384 ! *** Recurrence steps: [ss||c] -> [ab||c] ***
1385
1386 DO lc = 1, lc_max
1387
1388 DO cx = 0, lc
1389 DO cy = 0, lc - cx
1390 cz = lc - cx - cy
1391
1392 coc = coset(cx, cy, cz)
1393 cocx = coset(max(0, cx - 1), cy, cz)
1394 cocy = coset(cx, max(0, cy - 1), cz)
1395 cocz = coset(cx, cy, max(0, cz - 1))
1396
1397 fcx = f6*real(cx, dp)
1398 fcy = f6*real(cy, dp)
1399 fcz = f6*real(cz, dp)
1400
1401 ! *** Recurrence steps: [ss||c] -> [as||c] ***
1402
1403 IF (la_max > 0) THEN
1404
1405 ! *** Vertical recurrence steps: [ss||c] -> [as||c] ***
1406
1407 ! *** [ps||c]{n} = (Pi - Ai)*[ss||c]{n} + ***
1408 ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1409 ! *** f6*Ni(c)*[ss||c-1i]{n+1} (i = x,y,z) ***
1410
1411 DO n = 1, nmax - 1 - lc
1412 v(2, 1, coc, n) = rap(1)*v(1, 1, coc, n) + &
1413 rpw(1)*v(1, 1, coc, n + 1) + &
1414 fcx*v(1, 1, cocx, n + 1)
1415 v(3, 1, coc, n) = rap(2)*v(1, 1, coc, n) + &
1416 rpw(2)*v(1, 1, coc, n + 1) + &
1417 fcy*v(1, 1, cocy, n + 1)
1418 v(4, 1, coc, n) = rap(3)*v(1, 1, coc, n) + &
1419 rpw(3)*v(1, 1, coc, n + 1) + &
1420 fcz*v(1, 1, cocz, n + 1)
1421 END DO
1422
1423 ! *** [as||c]{n} = (Pi - Ai)*[(a-1i)s||c]{n} + ***
1424 ! *** (Wi - Pi)*[(a-1i)s||c]{n+1} + ***
1425 ! *** f2*Ni(a-1i)*( [(a-2i)s||c]{n} + ***
1426 ! *** f4*[(a-2i)s||c]{n+1}) + ***
1427 ! *** f6*Ni(c)*[(a-1i)s||c-1i]{n+1} ***
1428
1429 DO la = 2, la_max
1430
1431 DO n = 1, nmax - la - lc
1432
1433 ! *** Increase the angular momentum component z of a ***
1434
1435 v(coset(0, 0, la), 1, coc, n) = &
1436 rap(3)*v(coset(0, 0, la - 1), 1, coc, n) + &
1437 rpw(3)*v(coset(0, 0, la - 1), 1, coc, n + 1) + &
1438 f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, coc, n) + &
1439 f4*v(coset(0, 0, la - 2), 1, coc, n + 1)) + &
1440 fcz*v(coset(0, 0, la - 1), 1, cocz, n + 1)
1441
1442 ! *** Increase the angular momentum component y of a ***
1443
1444 az = la - 1
1445 v(coset(0, 1, az), 1, coc, n) = &
1446 rap(2)*v(coset(0, 0, az), 1, coc, n) + &
1447 rpw(2)*v(coset(0, 0, az), 1, coc, n + 1) + &
1448 fcy*v(coset(0, 0, az), 1, cocy, n + 1)
1449
1450 DO ay = 2, la
1451 f3 = f2*real(ay - 1, dp)
1452 az = la - ay
1453 v(coset(0, ay, az), 1, coc, n) = &
1454 rap(2)*v(coset(0, ay - 1, az), 1, coc, n) + &
1455 rpw(2)*v(coset(0, ay - 1, az), 1, coc, n + 1) + &
1456 f3*(v(coset(0, ay - 2, az), 1, coc, n) + &
1457 f4*v(coset(0, ay - 2, az), 1, coc, n + 1)) + &
1458 fcy*v(coset(0, ay - 1, az), 1, cocy, n + 1)
1459 END DO
1460
1461 ! *** Increase the angular momentum component x of a ***
1462
1463 DO ay = 0, la - 1
1464 az = la - 1 - ay
1465 v(coset(1, ay, az), 1, coc, n) = &
1466 rap(1)*v(coset(0, ay, az), 1, coc, n) + &
1467 rpw(1)*v(coset(0, ay, az), 1, coc, n + 1) + &
1468 fcx*v(coset(0, ay, az), 1, cocx, n + 1)
1469 END DO
1470
1471 DO ax = 2, la
1472 f3 = f2*real(ax - 1, dp)
1473 DO ay = 0, la - ax
1474 az = la - ax - ay
1475 v(coset(ax, ay, az), 1, coc, n) = &
1476 rap(1)*v(coset(ax - 1, ay, az), 1, coc, n) + &
1477 rpw(1)*v(coset(ax - 1, ay, az), 1, coc, n + 1) + &
1478 f3*(v(coset(ax - 2, ay, az), 1, coc, n) + &
1479 f4*v(coset(ax - 2, ay, az), 1, coc, n + 1)) + &
1480 fcx*v(coset(ax - 1, ay, az), 1, cocx, n + 1)
1481 END DO
1482 END DO
1483
1484 END DO
1485
1486 END DO
1487
1488 ! *** Recurrence steps: [as||c] -> [ab||c] ***
1489
1490 IF (lb_max > 0) THEN
1491
1492 ! *** Horizontal recurrence steps ***
1493
1494 ! *** [ap||c]{n} = [(a+1i)s||c]{n} - (Bi - Ai)*[as||c]{n} ***
1495
1496 la_start = max(0, la_min - 1)
1497
1498 DO la = la_start, la_max - 1
1499 DO n = 1, nmax - la - 1 - lc
1500 DO ax = 0, la
1501 DO ay = 0, la - ax
1502 az = la - ax - ay
1503 v(coset(ax, ay, az), 2, coc, n) = &
1504 v(coset(ax + 1, ay, az), 1, coc, n) - &
1505 rab(1)*v(coset(ax, ay, az), 1, coc, n)
1506 v(coset(ax, ay, az), 3, coc, n) = &
1507 v(coset(ax, ay + 1, az), 1, coc, n) - &
1508 rab(2)*v(coset(ax, ay, az), 1, coc, n)
1509 v(coset(ax, ay, az), 4, coc, n) = &
1510 v(coset(ax, ay, az + 1), 1, coc, n) - &
1511 rab(3)*v(coset(ax, ay, az), 1, coc, n)
1512 END DO
1513 END DO
1514 END DO
1515 END DO
1516
1517 ! *** Vertical recurrence step ***
1518
1519 ! *** [ap||c]{n} = (Pi - Bi)*[as||c]{n} + ***
1520 ! *** (Wi - Pi)*[as||c]{n+1} + ***
1521 ! *** f2*Ni(a)*( [(a-1i)s||c]{n} + ***
1522 ! *** f4*[(a-1i)s||c]{n+1}) + ***
1523 ! *** f6*Ni(c)*[(as||c-1i]{n+1}) ***
1524
1525 DO n = 1, nmax - la_max - 1 - lc
1526 DO ax = 0, la_max
1527 fx = f2*real(ax, dp)
1528 DO ay = 0, la_max - ax
1529 fy = f2*real(ay, dp)
1530 az = la_max - ax - ay
1531 fz = f2*real(az, dp)
1532
1533 IF (ax == 0) THEN
1534 v(coset(ax, ay, az), 2, coc, n) = &
1535 rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1536 rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1537 fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1538 ELSE
1539 v(coset(ax, ay, az), 2, coc, n) = &
1540 rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1541 rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1542 fx*(v(coset(ax - 1, ay, az), 1, coc, n) + &
1543 f4*v(coset(ax - 1, ay, az), 1, coc, n + 1)) + &
1544 fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1545 END IF
1546
1547 IF (ay == 0) THEN
1548 v(coset(ax, ay, az), 3, coc, n) = &
1549 rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1550 rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1551 fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1552 ELSE
1553 v(coset(ax, ay, az), 3, coc, n) = &
1554 rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1555 rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1556 fy*(v(coset(ax, ay - 1, az), 1, coc, n) + &
1557 f4*v(coset(ax, ay - 1, az), 1, coc, n + 1)) + &
1558 fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1559 END IF
1560
1561 IF (az == 0) THEN
1562 v(coset(ax, ay, az), 4, coc, n) = &
1563 rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1564 rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1565 fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1566 ELSE
1567 v(coset(ax, ay, az), 4, coc, n) = &
1568 rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1569 rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1570 fz*(v(coset(ax, ay, az - 1), 1, coc, n) + &
1571 f4*v(coset(ax, ay, az - 1), 1, coc, n + 1)) + &
1572 fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1573 END IF
1574
1575 END DO
1576 END DO
1577 END DO
1578
1579 ! *** Recurrence steps: [ap||c] -> [ab||c] ***
1580
1581 DO lb = 2, lb_max
1582
1583 ! *** Horizontal recurrence steps ***
1584
1585 ! *** [ab||c]{n} = [(a+1i)(b-1i)||c]{n} - ***
1586 ! *** (Bi - Ai)*[a(b-1i)||c]{n} ***
1587
1588 la_start = max(0, la_min - 1)
1589
1590 DO la = la_start, la_max - 1
1591 DO n = 1, nmax - la - lb - lc
1592 DO ax = 0, la
1593 DO ay = 0, la - ax
1594 az = la - ax - ay
1595
1596 ! *** Shift of angular momentum component z ***
1597
1598 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1599 v(coset(ax, ay, az + 1), &
1600 coset(0, 0, lb - 1), coc, n) - &
1601 rab(3)*v(coset(ax, ay, az), &
1602 coset(0, 0, lb - 1), coc, n)
1603
1604 ! *** Shift of angular momentum component y ***
1605
1606 DO by = 1, lb
1607 bz = lb - by
1608 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1609 v(coset(ax, ay + 1, az), &
1610 coset(0, by - 1, bz), coc, n) - &
1611 rab(2)*v(coset(ax, ay, az), &
1612 coset(0, by - 1, bz), coc, n)
1613 END DO
1614
1615 ! *** Shift of angular momentum component x ***
1616
1617 DO bx = 1, lb
1618 DO by = 0, lb - bx
1619 bz = lb - bx - by
1620 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1621 v(coset(ax + 1, ay, az), &
1622 coset(bx - 1, by, bz), coc, n) - &
1623 rab(1)*v(coset(ax, ay, az), &
1624 coset(bx - 1, by, bz), coc, n)
1625 END DO
1626 END DO
1627
1628 END DO
1629 END DO
1630 END DO
1631 END DO
1632
1633 ! *** Vertical recurrence step ***
1634
1635 ! *** [ab||c]{n} = (Pi - Bi)*[a(b-1i)||c]{n} + ***
1636 ! *** (Wi - Pi)*[a(b-1i)||c]{n+1} + ***
1637 ! *** f2*Ni(a)*( [(a-1i)(b-1i)||c]{n} + ***
1638 ! *** f4*[(a-1i)(b-1i)||c]{n+1}) ***
1639 ! *** f2*Ni(b-1i)*( [a(b-2i)||c]{n} + ***
1640 ! *** f4*[a(b-2i)||c]{n+1}) + ***
1641 ! *** f6*Ni(c)*[a(b-1i)||c-1i]{n+1}) ***
1642
1643 DO n = 1, nmax - la_max - lb - lc
1644 DO ax = 0, la_max
1645 fx = f2*real(ax, dp)
1646 DO ay = 0, la_max - ax
1647 fy = f2*real(ay, dp)
1648 az = la_max - ax - ay
1649 fz = f2*real(az, dp)
1650
1651 ! *** Shift of angular momentum component z from a to b ***
1652
1653 f3 = f2*real(lb - 1, dp)
1654
1655 IF (az == 0) THEN
1656 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1657 rbp(3)*v(coset(ax, ay, az), &
1658 coset(0, 0, lb - 1), coc, n) + &
1659 rpw(3)*v(coset(ax, ay, az), &
1660 coset(0, 0, lb - 1), coc, n + 1) + &
1661 f3*(v(coset(ax, ay, az), &
1662 coset(0, 0, lb - 2), coc, n) + &
1663 f4*v(coset(ax, ay, az), &
1664 coset(0, 0, lb - 2), coc, n + 1)) + &
1665 fcz*v(coset(ax, ay, az), &
1666 coset(0, 0, lb - 1), cocz, n + 1)
1667 ELSE
1668 v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1669 rbp(3)*v(coset(ax, ay, az), &
1670 coset(0, 0, lb - 1), coc, n) + &
1671 rpw(3)*v(coset(ax, ay, az), &
1672 coset(0, 0, lb - 1), coc, n + 1) + &
1673 fz*(v(coset(ax, ay, az - 1), &
1674 coset(0, 0, lb - 1), coc, n) + &
1675 f4*v(coset(ax, ay, az - 1), &
1676 coset(0, 0, lb - 1), coc, n + 1)) + &
1677 f3*(v(coset(ax, ay, az), &
1678 coset(0, 0, lb - 2), coc, n) + &
1679 f4*v(coset(ax, ay, az), &
1680 coset(0, 0, lb - 2), coc, n + 1)) + &
1681 fcz*v(coset(ax, ay, az), &
1682 coset(0, 0, lb - 1), cocz, n + 1)
1683 END IF
1684
1685 ! *** Shift of angular momentum component y from a to b ***
1686
1687 IF (ay == 0) THEN
1688 bz = lb - 1
1689 v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1690 rbp(2)*v(coset(ax, ay, az), &
1691 coset(0, 0, bz), coc, n) + &
1692 rpw(2)*v(coset(ax, ay, az), &
1693 coset(0, 0, bz), coc, n + 1) + &
1694 fcy*v(coset(ax, ay, az), &
1695 coset(0, 0, bz), cocy, n + 1)
1696 DO by = 2, lb
1697 bz = lb - by
1698 f3 = f2*real(by - 1, dp)
1699 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1700 rbp(2)*v(coset(ax, ay, az), &
1701 coset(0, by - 1, bz), coc, n) + &
1702 rpw(2)*v(coset(ax, ay, az), &
1703 coset(0, by - 1, bz), coc, n + 1) + &
1704 f3*(v(coset(ax, ay, az), &
1705 coset(0, by - 2, bz), coc, n) + &
1706 f4*v(coset(ax, ay, az), &
1707 coset(0, by - 2, bz), coc, n + 1)) + &
1708 fcy*v(coset(ax, ay, az), &
1709 coset(0, by - 1, bz), cocy, n + 1)
1710 END DO
1711 ELSE
1712 bz = lb - 1
1713 v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1714 rbp(2)*v(coset(ax, ay, az), &
1715 coset(0, 0, bz), coc, n) + &
1716 rpw(2)*v(coset(ax, ay, az), &
1717 coset(0, 0, bz), coc, n + 1) + &
1718 fy*(v(coset(ax, ay - 1, az), &
1719 coset(0, 0, bz), coc, n) + &
1720 f4*v(coset(ax, ay - 1, az), &
1721 coset(0, 0, bz), coc, n + 1)) + &
1722 fcy*v(coset(ax, ay, az), &
1723 coset(0, 0, bz), cocy, n + 1)
1724 DO by = 2, lb
1725 bz = lb - by
1726 f3 = f2*real(by - 1, dp)
1727 v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1728 rbp(2)*v(coset(ax, ay, az), &
1729 coset(0, by - 1, bz), coc, n) + &
1730 rpw(2)*v(coset(ax, ay, az), &
1731 coset(0, by - 1, bz), coc, n + 1) + &
1732 fy*(v(coset(ax, ay - 1, az), &
1733 coset(0, by - 1, bz), coc, n) + &
1734 f4*v(coset(ax, ay - 1, az), &
1735 coset(0, by - 1, bz), coc, n + 1)) + &
1736 f3*(v(coset(ax, ay, az), &
1737 coset(0, by - 2, bz), coc, n) + &
1738 f4*v(coset(ax, ay, az), &
1739 coset(0, by - 2, bz), coc, n + 1)) + &
1740 fcy*v(coset(ax, ay, az), &
1741 coset(0, by - 1, bz), cocy, n + 1)
1742 END DO
1743 END IF
1744
1745 ! *** Shift of angular momentum component x from a to b ***
1746
1747 IF (ax == 0) THEN
1748 DO by = 0, lb - 1
1749 bz = lb - 1 - by
1750 v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1751 rbp(1)*v(coset(ax, ay, az), &
1752 coset(0, by, bz), coc, n) + &
1753 rpw(1)*v(coset(ax, ay, az), &
1754 coset(0, by, bz), coc, n + 1) + &
1755 fcx*v(coset(ax, ay, az), &
1756 coset(0, by, bz), cocx, n + 1)
1757 END DO
1758 DO bx = 2, lb
1759 f3 = f2*real(bx - 1, dp)
1760 DO by = 0, lb - bx
1761 bz = lb - bx - by
1762 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1763 rbp(1)*v(coset(ax, ay, az), &
1764 coset(bx - 1, by, bz), coc, n) + &
1765 rpw(1)*v(coset(ax, ay, az), &
1766 coset(bx - 1, by, bz), coc, n + 1) + &
1767 f3*(v(coset(ax, ay, az), &
1768 coset(bx - 2, by, bz), coc, n) + &
1769 f4*v(coset(ax, ay, az), &
1770 coset(bx - 2, by, bz), coc, n + 1)) + &
1771 fcx*v(coset(ax, ay, az), &
1772 coset(bx - 1, by, bz), cocx, n + 1)
1773 END DO
1774 END DO
1775 ELSE
1776 DO by = 0, lb - 1
1777 bz = lb - 1 - by
1778 v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1779 rbp(1)*v(coset(ax, ay, az), &
1780 coset(0, by, bz), coc, n) + &
1781 rpw(1)*v(coset(ax, ay, az), &
1782 coset(0, by, bz), coc, n + 1) + &
1783 fx*(v(coset(ax - 1, ay, az), &
1784 coset(0, by, bz), coc, n) + &
1785 f4*v(coset(ax - 1, ay, az), &
1786 coset(0, by, bz), coc, n + 1)) + &
1787 fcx*v(coset(ax, ay, az), &
1788 coset(0, by, bz), cocx, n + 1)
1789 END DO
1790 DO bx = 2, lb
1791 f3 = f2*real(bx - 1, dp)
1792 DO by = 0, lb - bx
1793 bz = lb - bx - by
1794 v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1795 rbp(1)*v(coset(ax, ay, az), &
1796 coset(bx - 1, by, bz), coc, n) + &
1797 rpw(1)*v(coset(ax, ay, az), &
1798 coset(bx - 1, by, bz), coc, n + 1) + &
1799 fx*(v(coset(ax - 1, ay, az), &
1800 coset(bx - 1, by, bz), coc, n) + &
1801 f4*v(coset(ax - 1, ay, az), &
1802 coset(bx - 1, by, bz), coc, n + 1)) + &
1803 f3*(v(coset(ax, ay, az), &
1804 coset(bx - 2, by, bz), coc, n) + &
1805 f4*v(coset(ax, ay, az), &
1806 coset(bx - 2, by, bz), coc, n + 1)) + &
1807 fcx*v(coset(ax, ay, az), &
1808 coset(bx - 1, by, bz), cocx, n + 1)
1809 END DO
1810 END DO
1811 END IF
1812
1813 END DO
1814 END DO
1815 END DO
1816
1817 END DO
1818 END IF
1819
1820 ELSE
1821
1822 IF (lb_max > 0) THEN
1823
1824 ! *** Vertical recurrence steps: [ss||c] -> [sb||c] ***
1825
1826 ! *** [sp||c]{n} = (Pi - Bi)*[ss||c]{n} + ***
1827 ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1828 ! *** f6*Ni(c)**[ss||c-1i]{n+1} ***
1829
1830 DO n = 1, nmax - 1 - lc
1831 v(1, 2, coc, n) = rbp(1)*v(1, 1, coc, n) + &
1832 rpw(1)*v(1, 1, coc, n + 1) + &
1833 fcx*v(1, 1, cocx, n + 1)
1834 v(1, 3, coc, n) = rbp(2)*v(1, 1, coc, n) + &
1835 rpw(2)*v(1, 1, coc, n + 1) + &
1836 fcy*v(1, 1, cocy, n + 1)
1837 v(1, 4, coc, n) = rbp(3)*v(1, 1, coc, n) + &
1838 rpw(3)*v(1, 1, coc, n + 1) + &
1839 fcz*v(1, 1, cocz, n + 1)
1840 END DO
1841
1842 ! *** [sb||c]{n} = (Pi - Bi)*[s(b-1i)||c]{n} + ***
1843 ! *** (Wi - Pi)*[s(b-1i)||c]{n+1} + ***
1844 ! *** f2*Ni(b-1i)*( [s(b-2i)||c]{n} + ***
1845 ! *** f4*[s(b-2i)||c]{n+1}) + ***
1846 ! *** f6*Ni(c)**[s(b-1i)||c-1i]{n+1} ***
1847
1848 DO lb = 2, lb_max
1849
1850 DO n = 1, nmax - lb - lc
1851
1852 ! *** Increase the angular momentum component z of b ***
1853
1854 v(1, coset(0, 0, lb), coc, n) = &
1855 rbp(3)*v(1, coset(0, 0, lb - 1), coc, n) + &
1856 rpw(3)*v(1, coset(0, 0, lb - 1), coc, n + 1) + &
1857 f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), coc, n) + &
1858 f4*v(1, coset(0, 0, lb - 2), coc, n + 1)) + &
1859 fcz*v(1, coset(0, 0, lb - 1), cocz, n + 1)
1860
1861 ! *** Increase the angular momentum component y of b ***
1862
1863 bz = lb - 1
1864 v(1, coset(0, 1, bz), coc, n) = &
1865 rbp(2)*v(1, coset(0, 0, bz), coc, n) + &
1866 rpw(2)*v(1, coset(0, 0, bz), coc, n + 1) + &
1867 fcy*v(1, coset(0, 0, bz), cocy, n + 1)
1868
1869 DO by = 2, lb
1870 f3 = f2*real(by - 1, dp)
1871 bz = lb - by
1872 v(1, coset(0, by, bz), coc, n) = &
1873 rbp(2)*v(1, coset(0, by - 1, bz), coc, n) + &
1874 rpw(2)*v(1, coset(0, by - 1, bz), coc, n + 1) + &
1875 f3*(v(1, coset(0, by - 2, bz), coc, n) + &
1876 f4*v(1, coset(0, by - 2, bz), coc, n + 1)) + &
1877 fcy*v(1, coset(0, by - 1, bz), cocy, n + 1)
1878 END DO
1879
1880 ! *** Increase the angular momentum component x of b ***
1881
1882 DO by = 0, lb - 1
1883 bz = lb - 1 - by
1884 v(1, coset(1, by, bz), coc, n) = &
1885 rbp(1)*v(1, coset(0, by, bz), coc, n) + &
1886 rpw(1)*v(1, coset(0, by, bz), coc, n + 1) + &
1887 fcx*v(1, coset(0, by, bz), cocx, n + 1)
1888 END DO
1889
1890 DO bx = 2, lb
1891 f3 = f2*real(bx - 1, dp)
1892 DO by = 0, lb - bx
1893 bz = lb - bx - by
1894 v(1, coset(bx, by, bz), coc, n) = &
1895 rbp(1)*v(1, coset(bx - 1, by, bz), coc, n) + &
1896 rpw(1)*v(1, coset(bx - 1, by, bz), coc, n + 1) + &
1897 f3*(v(1, coset(bx - 2, by, bz), coc, n) + &
1898 f4*v(1, coset(bx - 2, by, bz), coc, n + 1)) + &
1899 fcx*v(1, coset(bx - 1, by, bz), cocx, n + 1)
1900 END DO
1901 END DO
1902
1903 END DO
1904
1905 END DO
1906
1907 END IF
1908
1909 END IF
1910
1911 END DO
1912 END DO
1913
1914 END DO
1915
1916 END IF
1917
1918 ! *** Add the contribution of the current pair ***
1919 ! *** of primitive Gaussian-type functions ***
1920
1921 DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1922 kk = k - ncoset(lc_min - 1)
1923 DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1924 DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
1925 vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1926 int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1927 END DO
1928 END DO
1929 END DO
1930
1931 IF (PRESENT(maxder)) THEN
1932 DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1933 kk = k - ncoset(lc_min - 1)
1934 DO j = 1, ncoset(lb_max)
1935 DO i = 1, ncoset(la_max)
1936 vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1937 END DO
1938 END DO
1939 END DO
1940 END IF
1941
1942 nb = nb + ncoset(lb_max)
1943
1944 END DO
1945
1946 na = na + ncoset(la_max - maxder_local)
1947 nap = nap + ncoset(la_max)
1948
1949 END DO
1950
1951 END SUBROUTINE coulomb3
1952
1953END 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, maxder, vac_plus)
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:798
ai_coulomb::coulomb2_new
subroutine, public coulomb2_new(la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, rac, rac2, vac, v, f, maxder, vac_plus)
Calculation of the primitive two-center Coulomb integrals over Cartesian Gaussian-type functions....
Definition ai_coulomb.F:441
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