Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of 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 : ! **************************************************************************************************
41 : MODULE 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 :
58 : CONTAINS
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 300 : SUBROUTINE coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, &
85 450 : 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 5148738 : v = 0.0_dp
107 :
108 150 : maxder_local = 0
109 150 : IF (PRESENT(maxder)) THEN
110 0 : maxder_local = maxder
111 0 : vac_plus = 0.0_dp
112 : END IF
113 :
114 150 : nmax = la_max + lc_max + 1
115 :
116 : ! *** Calculate the distance of the centers a and c ***
117 :
118 150 : dac = SQRT(rac2)
119 :
120 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
121 :
122 150 : na = 0
123 150 : nap = 0
124 :
125 750 : DO ipgf = 1, npgfa
126 :
127 600 : nc = 0
128 :
129 3000 : DO jpgf = 1, npgfc
130 :
131 : ! *** Screening ***
132 :
133 2400 : IF (rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
134 0 : DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
135 0 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max - maxder_local)
136 0 : vac(i, j) = 0.0_dp
137 : END DO
138 : END DO
139 0 : nc = nc + ncoset(lc_max)
140 0 : CYCLE
141 : END IF
142 :
143 : ! *** Calculate some prefactors ***
144 :
145 2400 : zetp = 1.0_dp/zeta(ipgf)
146 2400 : zetq = 1.0_dp/zetc(jpgf)
147 2400 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
148 :
149 2400 : rho = zeta(ipgf)*zetc(jpgf)*zetw
150 :
151 2400 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
152 :
153 : ! *** Calculate the incomplete Gamma function ***
154 :
155 2400 : t = rho*rac2
156 :
157 2400 : CALL fgamma(nmax - 1, t, f)
158 :
159 : ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
160 :
161 8096 : DO n = 1, nmax
162 8096 : v(1, 1, n) = f0*f(n - 1)
163 : END DO
164 :
165 : ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
166 :
167 2400 : IF (lc_max > 0) THEN
168 :
169 800 : f1 = 0.5_dp*zetq
170 800 : f2 = -rho*zetq
171 :
172 3200 : rcw(:) = -zeta(ipgf)*zetw*rac(:)
173 :
174 : ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
175 :
176 4096 : DO n = 1, nmax - 1
177 3296 : v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
178 3296 : v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
179 4096 : 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 1648 : DO lc = 2, lc_max
187 :
188 6592 : 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 4944 : f2*v(1, coset(0, 0, lc - 2), n + 1))
196 :
197 : ! *** Increase the angular momentum component y of c ***
198 :
199 4944 : cz = lc - 1
200 4944 : v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
201 :
202 15136 : DO cy = 2, lc
203 10192 : 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 15136 : 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 20080 : DO cy = 0, lc - 1
213 15136 : cz = lc - 1 - cy
214 20080 : v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
215 : END DO
216 :
217 15984 : DO cx = 2, lc
218 10192 : f6 = f1*REAL(cx - 1, dp)
219 34096 : DO cy = 0, lc - cx
220 18960 : 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 29152 : 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 2400 : IF (la_max > 0) THEN
237 :
238 800 : f3 = 0.5_dp*zetp
239 800 : f4 = -rho*zetp
240 800 : f5 = 0.5_dp*zetw
241 :
242 3200 : raw(:) = zetc(jpgf)*zetw*rac(:)
243 :
244 : ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
245 :
246 4096 : DO n = 1, nmax - 1
247 3296 : v(2, 1, n) = raw(1)*v(1, 1, n + 1)
248 3296 : v(3, 1, n) = raw(2)*v(1, 1, n + 1)
249 4096 : 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 1648 : DO la = 2, la_max
257 :
258 6592 : 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 4944 : f4*v(coset(0, 0, la - 2), 1, n + 1))
266 :
267 : ! *** Increase the angular momentum component y of a ***
268 :
269 4944 : az = la - 1
270 4944 : v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
271 :
272 15136 : DO ay = 2, la
273 10192 : 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 15136 : 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 20080 : DO ay = 0, la - 1
283 15136 : az = la - 1 - ay
284 20080 : v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
285 : END DO
286 :
287 15984 : DO ax = 2, la
288 10192 : f6 = f3*REAL(ax - 1, dp)
289 34096 : DO ay = 0, la - ax
290 18960 : 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 29152 : 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 2448 : DO lc = 1, lc_max
303 :
304 7392 : DO cx = 0, lc
305 17792 : DO cy = 0, lc - cx
306 11200 : cz = lc - cx - cy
307 :
308 11200 : coc = coset(cx, cy, cz)
309 11200 : cocx = coset(MAX(0, cx - 1), cy, cz)
310 11200 : cocy = coset(cx, MAX(0, cy - 1), cz)
311 11200 : cocz = coset(cx, cy, MAX(0, cz - 1))
312 :
313 11200 : fcx = f5*REAL(cx, dp)
314 11200 : fcy = f5*REAL(cy, dp)
315 11200 : 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 64064 : DO n = 1, nmax - 1 - lc
321 52864 : v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
322 52864 : v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
323 64064 : 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 48224 : DO la = 2, la_max
332 :
333 157136 : 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 113856 : fcz*v(coset(0, 0, la - 1), cocz, n + 1)
342 :
343 : ! *** Increase the angular momentum component y of a ***
344 :
345 113856 : az = la - 1
346 : v(coset(0, 1, az), coc, n) = &
347 : raw(2)*v(coset(0, 0, az), coc, n + 1) + &
348 113856 : fcy*v(coset(0, 0, az), cocy, n + 1)
349 :
350 366944 : DO ay = 2, la
351 253088 : 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 366944 : 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 480800 : DO ay = 0, la - 1
362 366944 : az = la - 1 - ay
363 : v(coset(1, ay, az), coc, n) = &
364 : raw(1)*v(coset(0, ay, az), coc, n + 1) + &
365 480800 : fcx*v(coset(0, ay, az), cocx, n + 1)
366 : END DO
367 :
368 399024 : DO ax = 2, la
369 253088 : f6 = f3*REAL(ax - 1, dp)
370 858784 : DO ay = 0, la - ax
371 491840 : 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 744928 : 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 9744 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
392 84480 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
393 82080 : vac(na + i, nc + j) = v(i, j, 1)
394 : END DO
395 : END DO
396 :
397 2400 : IF (PRESENT(maxder)) THEN
398 0 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
399 0 : DO i = 1, ncoset(la_max)
400 0 : vac_plus(nap + i, nc + j) = v(i, j, 1)
401 : END DO
402 : END DO
403 : END IF
404 :
405 3000 : nc = nc + ncoset(lc_max)
406 :
407 : END DO
408 :
409 600 : na = na + ncoset(la_max - maxder_local)
410 750 : nap = nap + ncoset(la_max)
411 :
412 : END DO
413 :
414 150 : 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 62448 : SUBROUTINE coulomb2_new(la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
440 124896 : 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 9567511 : v = 0.0_dp
462 :
463 62448 : maxder_local = 0
464 62448 : IF (PRESENT(maxder)) THEN
465 0 : maxder_local = maxder
466 0 : vac_plus = 0.0_dp
467 : END IF
468 :
469 62448 : nmax = la_max + lc_max + 1
470 :
471 : ! *** Calculate the distance of the centers a and c ***
472 :
473 62448 : dac = SQRT(rac2)
474 :
475 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
476 :
477 62448 : na = 0
478 62448 : nap = 0
479 :
480 124896 : DO ipgf = 1, npgfa
481 :
482 62448 : nc = 0
483 62448 : ncp = 0
484 :
485 124896 : DO jpgf = 1, npgfc
486 :
487 : ! *** Calculate some prefactors ***
488 :
489 62448 : zetp = 1.0_dp/zeta(ipgf)
490 62448 : zetq = 1.0_dp/zetc(jpgf)
491 62448 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
492 :
493 62448 : rho = zeta(ipgf)*zetc(jpgf)*zetw
494 :
495 62448 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
496 :
497 : ! *** Calculate the incomplete Gamma function ***
498 :
499 62448 : t = rho*rac2
500 :
501 62448 : CALL fgamma(nmax - 1, t, f)
502 :
503 : ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
504 :
505 242916 : DO n = 1, nmax
506 242916 : v(1, 1, n) = f0*f(n - 1)
507 : END DO
508 :
509 : ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
510 :
511 62448 : IF (lc_max > 0) THEN
512 :
513 37117 : f1 = 0.5_dp*zetq
514 37117 : f2 = -rho*zetq
515 :
516 148468 : rcw(:) = -zeta(ipgf)*zetw*rac(:)
517 :
518 : ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
519 :
520 131558 : DO n = 1, nmax - 1
521 94441 : v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
522 94441 : v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
523 131558 : 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 59800 : DO lc = 2, lc_max
531 :
532 108529 : 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 48729 : f2*v(1, coset(0, 0, lc - 2), n + 1))
540 :
541 : ! *** Increase the angular momentum component y of c ***
542 :
543 48729 : cz = lc - 1
544 48729 : v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
545 :
546 106800 : DO cy = 2, lc
547 58071 : 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 106800 : 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 155529 : DO cy = 0, lc - 1
557 106800 : cz = lc - 1 - cy
558 155529 : v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
559 : END DO
560 :
561 129483 : DO cx = 2, lc
562 58071 : f6 = f1*REAL(cx - 1, dp)
563 174213 : DO cy = 0, lc - cx
564 67413 : 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 125484 : 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 62448 : IF (la_max > 0) THEN
581 :
582 36634 : f3 = 0.5_dp*zetp
583 36634 : f4 = -rho*zetp
584 36634 : f5 = 0.5_dp*zetw
585 :
586 146536 : raw(:) = zetc(jpgf)*zetw*rac(:)
587 :
588 : ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
589 :
590 129974 : DO n = 1, nmax - 1
591 93340 : v(2, 1, n) = raw(1)*v(1, 1, n + 1)
592 93340 : v(3, 1, n) = raw(2)*v(1, 1, n + 1)
593 129974 : 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 58220 : DO la = 2, la_max
601 :
602 104887 : 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 46667 : f4*v(coset(0, 0, la - 2), 1, n + 1))
610 :
611 : ! *** Increase the angular momentum component y of a ***
612 :
613 46667 : az = la - 1
614 46667 : v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
615 :
616 101838 : DO ay = 2, la
617 55171 : 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 101838 : 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 148505 : DO ay = 0, la - 1
627 101838 : az = la - 1 - ay
628 148505 : v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
629 : END DO
630 :
631 123424 : DO ax = 2, la
632 55171 : f6 = f3*REAL(ax - 1, dp)
633 165513 : DO ay = 0, la - ax
634 63675 : 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 118846 : 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 71754 : DO lc = 1, lc_max
647 :
648 158165 : DO cx = 0, lc
649 278241 : DO cy = 0, lc - cx
650 156710 : cz = lc - cx - cy
651 :
652 156710 : coc = coset(cx, cy, cz)
653 156710 : cocx = coset(MAX(0, cx - 1), cy, cz)
654 156710 : cocy = coset(cx, MAX(0, cy - 1), cz)
655 156710 : cocz = coset(cx, cy, MAX(0, cz - 1))
656 :
657 156710 : fcx = f5*REAL(cx, dp)
658 156710 : fcy = f5*REAL(cy, dp)
659 156710 : 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 463161 : DO n = 1, nmax - 1 - lc
665 306451 : v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
666 306451 : v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
667 463161 : 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 335838 : DO la = 2, la_max
676 :
677 394670 : 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 145243 : fcz*v(coset(0, 0, la - 1), cocz, n + 1)
686 :
687 : ! *** Increase the angular momentum component y of a ***
688 :
689 145243 : az = la - 1
690 : v(coset(0, 1, az), coc, n) = &
691 : raw(2)*v(coset(0, 0, az), coc, n + 1) + &
692 145243 : fcy*v(coset(0, 0, az), cocy, n + 1)
693 :
694 316018 : DO ay = 2, la
695 170775 : 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 316018 : 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 461261 : DO ay = 0, la - 1
706 316018 : az = la - 1 - ay
707 : v(coset(1, ay, az), coc, n) = &
708 : raw(1)*v(coset(0, ay, az), coc, n + 1) + &
709 461261 : fcx*v(coset(0, ay, az), cocx, n + 1)
710 : END DO
711 :
712 408735 : DO ax = 2, la
713 170775 : f6 = f3*REAL(ax - 1, dp)
714 512325 : DO ay = 0, la - ax
715 196307 : 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 367082 : 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 271993 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
736 960289 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
737 897841 : vac(na + i, nc + j) = v(i, j, 1)
738 : END DO
739 : END DO
740 :
741 62448 : IF (PRESENT(maxder)) THEN
742 0 : DO j = 1, ncoset(lc_max)
743 0 : DO i = 1, ncoset(la_max)
744 0 : vac_plus(nap + i, ncp + j) = v(i, j, 1)
745 : END DO
746 : END DO
747 : END IF
748 :
749 62448 : nc = nc + ncoset(lc_max - maxder_local)
750 124896 : ncp = ncp + ncoset(lc_max)
751 :
752 : END DO
753 :
754 62448 : na = na + ncoset(la_max - maxder_local)
755 124896 : nap = nap + ncoset(la_max)
756 :
757 : END DO
758 :
759 62448 : 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 23706 : SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
796 47412 : lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
797 71118 : 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 20791632 : v = 0.0_dp
827 :
828 23706 : maxder_local = 0
829 23706 : IF (PRESENT(maxder)) THEN
830 0 : maxder_local = maxder
831 : END IF
832 :
833 23706 : nmax = la_max + lb_max + lc_max + 1
834 :
835 : ! *** Calculate the distances of the centers a, b and c ***
836 :
837 23706 : dab = SQRT(rab2)
838 23706 : dac = SQRT(rac2)
839 23706 : dbc = SQRT(rbc2)
840 :
841 : ! *** Initialize integrals array
842 2513361 : int_abc = 0.0_dp
843 :
844 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
845 :
846 23706 : na = 0
847 23706 : nap = 0
848 :
849 61956 : DO ipgf = 1, npgfa
850 :
851 : ! *** Screening ***
852 38250 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
853 0 : na = na + ncoset(la_max - maxder_local)
854 0 : nap = nap + ncoset(la_max)
855 0 : CYCLE
856 : END IF
857 :
858 38250 : nb = 0
859 :
860 100260 : DO jpgf = 1, npgfb
861 :
862 : ! *** Screening ***
863 : IF ( &
864 62010 : (rpgfb(jpgf) + rpgfc < dbc) .OR. &
865 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
866 0 : nb = nb + ncoset(lb_max)
867 0 : CYCLE
868 : END IF
869 :
870 : ! *** Calculate some prefactors ***
871 :
872 62010 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
873 62010 : zetq = 1.0_dp/zetc
874 62010 : zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
875 :
876 62010 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
877 62010 : f1 = zetb(jpgf)*zetp
878 62010 : f2 = 0.5_dp*zetp
879 62010 : f4 = -zetc*zetw
880 :
881 62010 : f0 = f0*EXP(-zeta(ipgf)*f1*rab2)
882 :
883 248040 : rap(:) = f1*rab(:)
884 248040 : rcp(:) = rap(:) - rac(:)
885 248040 : rpw(:) = f4*rcp(:)
886 :
887 : ! *** Calculate the incomplete Gamma function ***
888 :
889 62010 : t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
890 :
891 62010 : CALL fgamma(nmax - 1, t, f)
892 :
893 : ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
894 :
895 255660 : DO n = 1, nmax
896 255660 : v(1, 1, 1, n) = f0*f(n - 1)
897 : END DO
898 :
899 : ! *** Recurrence steps: [ss||s] -> [as||s] ***
900 :
901 62010 : 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 112306 : DO n = 1, nmax - 1
909 82966 : v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
910 82966 : v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
911 112306 : 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 34140 : DO la = 2, la_max
920 :
921 47188 : 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 13048 : f4*v(coset(0, 0, la - 2), 1, 1, n + 1))
930 :
931 : ! *** Increase the angular momentum component y of a ***
932 :
933 13048 : az = la - 1
934 : v(coset(0, 1, az), 1, 1, n) = &
935 : rap(2)*v(coset(0, 0, az), 1, 1, n) + &
936 13048 : rpw(2)*v(coset(0, 0, az), 1, 1, n + 1)
937 :
938 26096 : DO ay = 2, la
939 13048 : 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 26096 : 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 39144 : DO ay = 0, la - 1
950 26096 : 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 39144 : rpw(1)*v(coset(0, ay, az), 1, 1, n + 1)
954 : END DO
955 :
956 30896 : DO ax = 2, la
957 13048 : f3 = f2*REAL(ax - 1, dp)
958 39144 : DO ay = 0, la - ax
959 13048 : 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 26096 : 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 29340 : IF (lb_max > 0) THEN
975 :
976 : ! *** Horizontal recurrence steps ***
977 :
978 65952 : rbp(:) = rap(:) - rab(:)
979 :
980 : ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
981 :
982 16488 : la_start = MAX(0, la_min - 1)
983 :
984 32976 : DO la = la_start, la_max - 1
985 85430 : DO n = 1, nmax - la - 1
986 130604 : DO ax = 0, la
987 184986 : DO ay = 0, la - ax
988 70870 : az = la - ax - ay
989 : v(coset(ax, ay, az), 2, 1, n) = &
990 : v(coset(ax + 1, ay, az), 1, 1, n) - &
991 70870 : 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 70870 : 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 132532 : 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 52454 : DO n = 1, nmax - la_max - 1
1011 130714 : DO ax = 0, la_max
1012 78260 : fx = f2*REAL(ax, dp)
1013 241108 : DO ay = 0, la_max - ax
1014 126882 : fy = f2*REAL(ay, dp)
1015 126882 : az = la_max - ax - ay
1016 126882 : fz = f2*REAL(az, dp)
1017 :
1018 126882 : 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 78260 : 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 48622 : f4*v(coset(ax - 1, ay, az), 1, 1, n + 1))
1028 : END IF
1029 :
1030 126882 : 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 78260 : 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 48622 : f4*v(coset(ax, ay - 1, az), 1, 1, n + 1))
1040 : END IF
1041 :
1042 205142 : 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 78260 : 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 48622 : 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 19416 : 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 5856 : la_start = MAX(0, la_min - 1)
1068 :
1069 5856 : DO la = la_start, la_max - 1
1070 14652 : DO n = 1, nmax - la - lb
1071 22196 : DO ax = 0, la
1072 31416 : DO ay = 0, la - ax
1073 12148 : 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 12148 : 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 36444 : DO by = 1, lb
1084 24296 : 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 36444 : 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 46916 : DO bx = 1, lb
1093 72888 : DO by = 0, lb - bx
1094 36444 : 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 60740 : 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 25284 : DO n = 1, nmax - la_max - lb
1116 21650 : DO ax = 0, la_max
1117 12854 : fx = f2*REAL(ax, dp)
1118 39680 : DO ay = 0, la_max - ax
1119 20958 : fy = f2*REAL(ay, dp)
1120 20958 : az = la_max - ax - ay
1121 20958 : fz = f2*REAL(az, dp)
1122 :
1123 : ! *** Shift of angular momentum component z from a to b ***
1124 :
1125 20958 : f3 = f2*REAL(lb - 1, dp)
1126 :
1127 20958 : 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 12854 : 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 8104 : 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 20958 : IF (ay == 0) THEN
1146 12854 : 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 12854 : rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1)
1150 25708 : DO by = 2, lb
1151 12854 : bz = lb - by
1152 12854 : 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 25708 : f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
1158 : END DO
1159 : ELSE
1160 8104 : 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 8104 : f4*v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n + 1))
1166 16208 : DO by = 2, lb
1167 8104 : bz = lb - by
1168 8104 : 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 16208 : 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 33812 : IF (ax == 0) THEN
1183 38562 : DO by = 0, lb - 1
1184 25708 : 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 38562 : rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1)
1188 : END DO
1189 25708 : DO bx = 2, lb
1190 12854 : f3 = f2*REAL(bx - 1, dp)
1191 38562 : DO by = 0, lb - bx
1192 12854 : 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 25708 : f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
1199 : END DO
1200 : END DO
1201 : ELSE
1202 24312 : DO by = 0, lb - 1
1203 16208 : 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 24312 : f4*v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n + 1))
1209 : END DO
1210 16208 : DO bx = 2, lb
1211 8104 : f3 = f2*REAL(bx - 1, dp)
1212 24312 : DO by = 0, lb - bx
1213 8104 : 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 16208 : 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 32670 : IF (lb_max > 0) THEN
1239 :
1240 : ! *** Vertical recurrence steps: [ss||s] -> [sb||s] ***
1241 :
1242 55248 : rbp(:) = rap(:) - rab(:)
1243 :
1244 : ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
1245 : ! *** (Wi - Pi)*[ss||s]{n+1} ***
1246 :
1247 43596 : DO n = 1, nmax - 1
1248 29784 : v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
1249 29784 : v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
1250 43596 : 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 15924 : DO lb = 2, lb_max
1259 :
1260 20156 : 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 4232 : f4*v(1, coset(0, 0, lb - 2), 1, n + 1))
1269 :
1270 : ! *** Increase the angular momentum component y of b ***
1271 :
1272 4232 : bz = lb - 1
1273 : v(1, coset(0, 1, bz), 1, n) = &
1274 : rbp(2)*v(1, coset(0, 0, bz), 1, n) + &
1275 4232 : rpw(2)*v(1, coset(0, 0, bz), 1, n + 1)
1276 :
1277 8464 : DO by = 2, lb
1278 4232 : 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 8464 : 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 12696 : DO by = 0, lb - 1
1289 8464 : 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 12696 : rpw(1)*v(1, coset(0, by, bz), 1, n + 1)
1293 : END DO
1294 :
1295 10576 : DO bx = 2, lb
1296 4232 : f3 = f2*REAL(bx - 1, dp)
1297 12696 : DO by = 0, lb - bx
1298 4232 : 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 8464 : 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 62010 : IF (lc_max > 0) THEN
1318 :
1319 : ! *** Vertical recurrence steps: [ss||s] -> [ss||c] ***
1320 :
1321 37296 : f5 = -zetw/zetp
1322 37296 : f6 = 0.5_dp*zetw
1323 37296 : f7 = 0.5_dp*zetq
1324 :
1325 149184 : rcw(:) = rcp(:) + rpw(:)
1326 :
1327 : ! *** [ss||p]{n} = (Wi - Ci)*[ss||s]{n+1} (i = x,y,z) ***
1328 :
1329 141270 : DO n = 1, nmax - 1
1330 103974 : v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
1331 103974 : v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
1332 141270 : 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 62160 : DO lc = 2, lc_max
1340 :
1341 121041 : 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 58881 : f5*v(1, 1, coset(0, 0, lc - 2), n + 1))
1349 :
1350 : ! *** Increase the angular momentum component y of c ***
1351 :
1352 58881 : cz = lc - 1
1353 58881 : v(1, 1, coset(0, 1, cz), n) = rcw(2)*v(1, 1, coset(0, 0, cz), n + 1)
1354 :
1355 130767 : DO cy = 2, lc
1356 71886 : 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 130767 : 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 189648 : DO cy = 0, lc - 1
1366 130767 : cz = lc - 1 - cy
1367 189648 : v(1, 1, coset(1, cy, cz), n) = rcw(1)*v(1, 1, coset(0, cy, cz), n + 1)
1368 : END DO
1369 :
1370 155631 : DO cx = 2, lc
1371 215658 : DO cy = 0, lc - cx
1372 84891 : 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 156777 : 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 99456 : DO lc = 1, lc_max
1387 :
1388 254781 : DO cx = 0, lc
1389 503121 : DO cy = 0, lc - cx
1390 285636 : cz = lc - cx - cy
1391 :
1392 285636 : coc = coset(cx, cy, cz)
1393 285636 : cocx = coset(MAX(0, cx - 1), cy, cz)
1394 285636 : cocy = coset(cx, MAX(0, cy - 1), cz)
1395 285636 : cocz = coset(cx, cy, MAX(0, cz - 1))
1396 :
1397 285636 : fcx = f6*REAL(cx, dp)
1398 285636 : fcy = f6*REAL(cy, dp)
1399 285636 : fcz = f6*REAL(cz, dp)
1400 :
1401 : ! *** Recurrence steps: [ss||c] -> [as||c] ***
1402 :
1403 440961 : 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 434632 : 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 299486 : 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 299486 : 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 434632 : 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 157252 : DO la = 2, la_max
1430 :
1431 203832 : 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 46580 : fcz*v(coset(0, 0, la - 1), 1, cocz, n + 1)
1441 :
1442 : ! *** Increase the angular momentum component y of a ***
1443 :
1444 46580 : 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 46580 : fcy*v(coset(0, 0, az), 1, cocy, n + 1)
1449 :
1450 93160 : DO ay = 2, la
1451 46580 : f3 = f2*REAL(ay - 1, dp)
1452 46580 : 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 93160 : 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 139740 : DO ay = 0, la - 1
1464 93160 : 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 139740 : fcx*v(coset(0, ay, az), 1, cocx, n + 1)
1469 : END DO
1470 :
1471 115266 : DO ax = 2, la
1472 46580 : f3 = f2*REAL(ax - 1, dp)
1473 139740 : DO ay = 0, la - ax
1474 46580 : 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 93160 : 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 135146 : 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 76006 : la_start = MAX(0, la_min - 1)
1497 :
1498 152012 : DO la = la_start, la_max - 1
1499 347128 : DO n = 1, nmax - la - 1 - lc
1500 500530 : DO ax = 0, la
1501 688224 : DO ay = 0, la - ax
1502 263700 : az = la - ax - ay
1503 : v(coset(ax, ay, az), 2, coc, n) = &
1504 : v(coset(ax + 1, ay, az), 1, coc, n) - &
1505 263700 : 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 263700 : 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 493108 : 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 195116 : DO n = 1, nmax - la_max - 1 - lc
1526 454354 : DO ax = 0, la_max
1527 259238 : fx = f2*REAL(ax, dp)
1528 798732 : DO ay = 0, la_max - ax
1529 420384 : fy = f2*REAL(ay, dp)
1530 420384 : az = la_max - ax - ay
1531 420384 : fz = f2*REAL(az, dp)
1532 :
1533 420384 : 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 259238 : 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 161146 : fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1545 : END IF
1546 :
1547 420384 : 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 259238 : 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 161146 : fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1559 : END IF
1560 :
1561 679622 : 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 259238 : 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 161146 : 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 89506 : 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 27000 : la_start = MAX(0, la_min - 1)
1589 :
1590 27000 : DO la = la_start, la_max - 1
1591 59256 : DO n = 1, nmax - la - lb - lc
1592 84158 : DO ax = 0, la
1593 115206 : DO ay = 0, la - ax
1594 44548 : 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 44548 : coset(0, 0, lb - 1), coc, n)
1603 :
1604 : ! *** Shift of angular momentum component y ***
1605 :
1606 133644 : DO by = 1, lb
1607 89096 : 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 133644 : coset(0, by - 1, bz), coc, n)
1613 : END DO
1614 :
1615 : ! *** Shift of angular momentum component x ***
1616 :
1617 172046 : DO bx = 1, lb
1618 267288 : DO by = 0, lb - bx
1619 133644 : 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 222740 : 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 108262 : DO n = 1, nmax - la_max - lb - lc
1644 73342 : DO ax = 0, la_max
1645 41086 : fx = f2*REAL(ax, dp)
1646 126832 : DO ay = 0, la_max - ax
1647 66990 : fy = f2*REAL(ay, dp)
1648 66990 : az = la_max - ax - ay
1649 66990 : fz = f2*REAL(az, dp)
1650 :
1651 : ! *** Shift of angular momentum component z from a to b ***
1652 :
1653 66990 : f3 = f2*REAL(lb - 1, dp)
1654 :
1655 66990 : 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 41086 : 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 25904 : 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 66990 : IF (ay == 0) THEN
1688 41086 : 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 41086 : coset(0, 0, bz), cocy, n + 1)
1696 82172 : DO by = 2, lb
1697 41086 : bz = lb - by
1698 41086 : 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 82172 : coset(0, by - 1, bz), cocy, n + 1)
1710 : END DO
1711 : ELSE
1712 25904 : 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 25904 : coset(0, 0, bz), cocy, n + 1)
1724 51808 : DO by = 2, lb
1725 25904 : bz = lb - by
1726 25904 : 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 51808 : 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 108076 : IF (ax == 0) THEN
1748 123258 : DO by = 0, lb - 1
1749 82172 : 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 123258 : coset(0, by, bz), cocx, n + 1)
1757 : END DO
1758 82172 : DO bx = 2, lb
1759 41086 : f3 = f2*REAL(bx - 1, dp)
1760 123258 : DO by = 0, lb - bx
1761 41086 : 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 82172 : coset(bx - 1, by, bz), cocx, n + 1)
1773 : END DO
1774 : END DO
1775 : ELSE
1776 77712 : DO by = 0, lb - 1
1777 51808 : 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 77712 : coset(0, by, bz), cocx, n + 1)
1789 : END DO
1790 51808 : DO bx = 2, lb
1791 25904 : f3 = f2*REAL(bx - 1, dp)
1792 77712 : DO by = 0, lb - bx
1793 25904 : 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 51808 : 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 150490 : 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 161860 : 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 98200 : 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 98200 : 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 161860 : 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 73396 : DO lb = 2, lb_max
1849 :
1850 86924 : 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 13528 : fcz*v(1, coset(0, 0, lb - 1), cocz, n + 1)
1860 :
1861 : ! *** Increase the angular momentum component y of b ***
1862 :
1863 13528 : 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 13528 : fcy*v(1, coset(0, 0, bz), cocy, n + 1)
1868 :
1869 27056 : DO by = 2, lb
1870 13528 : f3 = f2*REAL(by - 1, dp)
1871 13528 : 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 27056 : 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 40584 : DO by = 0, lb - 1
1883 27056 : 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 40584 : fcx*v(1, coset(0, by, bz), cocx, n + 1)
1888 : END DO
1889 :
1890 36792 : DO bx = 2, lb
1891 13528 : f3 = f2*REAL(bx - 1, dp)
1892 40584 : DO by = 0, lb - bx
1893 13528 : 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 27056 : 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 279345 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1922 217335 : kk = k - ncoset(lc_min - 1)
1923 820800 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1924 2181894 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
1925 1423104 : vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1926 1964559 : int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1927 : END DO
1928 : END DO
1929 : END DO
1930 :
1931 62010 : IF (PRESENT(maxder)) THEN
1932 0 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1933 0 : kk = k - ncoset(lc_min - 1)
1934 0 : DO j = 1, ncoset(lb_max)
1935 0 : DO i = 1, ncoset(la_max)
1936 0 : 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 100260 : nb = nb + ncoset(lb_max)
1943 :
1944 : END DO
1945 :
1946 38250 : na = na + ncoset(la_max - maxder_local)
1947 61956 : nap = nap + ncoset(la_max)
1948 :
1949 : END DO
1950 :
1951 23706 : END SUBROUTINE coulomb3
1952 :
1953 : END MODULE ai_coulomb
|
