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