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 Build up the nuclear potential part of the core Hamiltonian matrix
10 : !> in the case of an allelectron calculation by computing the nuclear
11 : !> attraction integral matrix <a|-Z/r|b> and the integral matrix
12 : !> <a|Z*erf(r)/r|b>. The integrals <a|Z*erf(r)/r|b> can be rewritten
13 : !> as the three-center Coulomb integrals <ab||c> with a primitive
14 : !> s-type Gaussian function c multiplied by a factor N.
15 : !>
16 : !> erfc(r)
17 : !> <a|V|b> = -Z*<a|---------|b>
18 : !> r
19 : !>
20 : !> 1 - erf(r)
21 : !> = -Z*<a|------------|b>
22 : !> r
23 : !>
24 : !> 1 erf(r)
25 : !> = -Z*(<a|---|b> - <a|--------|b>)
26 : !> r r
27 : !>
28 : !> 1
29 : !> = -Z*(<a|---|b> - N*<ab||c>)
30 : !> r
31 : !>
32 : !> 1
33 : !> = -Z*<a|---|b> + Z*N*<ab||c>
34 : !> r
35 : !> \_______/ \_____/
36 : !> | |
37 : !> nuclear coulomb
38 : !> \par Literature
39 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
40 : !> \par Parameters
41 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
42 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
43 : !> - coset : Cartesian orbital set pointer.
44 : !> - dab : Distance between the atomic centers a and b.
45 : !> - dac : Distance between the atomic centers a and c.
46 : !> - dbc : Distance between the atomic centers b and c.
47 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
48 : !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
49 : !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
50 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
51 : !> - npgf{a,b} : Degree of contraction of shell a or b.
52 : !> - rab : Distance vector between the atomic centers a and b.
53 : !> - rab2 : Square of the distance between the atomic centers a and b.
54 : !> - rac : Distance vector between the atomic centers a and c.
55 : !> - rac2 : Square of the distance between the atomic centers a and c.
56 : !> - rbc2 : Square of the distance between the atomic centers b and c.
57 : !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
58 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
59 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
60 : !> - zetw : Reciprocal of the sum of the exponents of orbital a, b and c.
61 : !> Indices for pVp matrices
62 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
63 : !> - n{x1,y1,z1}{x2,y2,z2} : indices for the pVp matrix elements
64 : !> <p_{x1,y1,z1}a|V|p_{x2,y2,z2}>
65 : !> - co{a,b}{m,p}{x,y,z} : coefficients of the pVp matrix elements: co == coset(lx,ly,lz)
66 : !> m == li - 1; p == li +1; li= x,y,z (which l quantum number is modified); l = a,b
67 : !> coset ( ... -1 ...) = 1 , should be zero, lmi = 1.0 for li >=1 and lmi = 0.0 for li = 0
68 : !> guarantees this as a prefactor
69 : !> - f{a,b}{x,y,z},ft{a,b},z : f(t)li: prefactors as a result of derivations of l with respect to i
70 : !> - pVpout : pVp matrix elements
71 : !> - pVp : local pVp matrix elements
72 : !> - lamin,lbmin : minimal angular quantum number used in pVp matrices
73 : !> _ vnucp : potential of the nuclei
74 : !> - vnuc : potential of the electrons
75 : !> _ pVpa, pVpb : pVpl=coset(lx,ly,lz)
76 : !> - na_pgf,nb_pgf : indices for primitive gaussian functions
77 : !> \par History
78 : !> 10.2008 added pVp matrix elements (jens)
79 : !> \author Matthias Krack (04.10.2000)
80 : ! **************************************************************************************************
81 :
82 : MODULE ai_verfc
83 :
84 : USE gamma, ONLY: fgamma => fgamma_0
85 : USE kinds, ONLY: dp
86 : USE mathconstants, ONLY: pi
87 : USE orbital_pointers, ONLY: coset,&
88 : ncoset
89 : #include "../base/base_uses.f90"
90 :
91 : IMPLICIT NONE
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_verfc'
93 : PRIVATE
94 :
95 : ! *** Public subroutines ***
96 :
97 : PUBLIC :: verfc
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief Calculation of the primitive three-center nuclear potential
103 : !> integrals <a|Z*erfc(r)/r|b> over Cartesian Gaussian-type
104 : !> functions.
105 : !> \param la_max1 ...
106 : !> \param npgfa ...
107 : !> \param zeta ...
108 : !> \param rpgfa ...
109 : !> \param la_min1 ...
110 : !> \param lb_max1 ...
111 : !> \param npgfb ...
112 : !> \param zetb ...
113 : !> \param rpgfb ...
114 : !> \param lb_min1 ...
115 : !> \param zetc ...
116 : !> \param rpgfc ...
117 : !> \param zc ...
118 : !> \param cerf ...
119 : !> \param rab ...
120 : !> \param rab2 ...
121 : !> \param rac ...
122 : !> \param rac2 ...
123 : !> \param rbc2 ...
124 : !> \param vabc ...
125 : !> \param verf ...
126 : !> \param vnuc ...
127 : !> \param f ...
128 : !> \param maxder ...
129 : !> \param vabc_plus ...
130 : !> \param vnabc ...
131 : !> \param pVp_sum ...
132 : !> \param pVp ...
133 : !> \param dkh_erfc ...
134 : !> \date 04.10.2000
135 : !> \author Matthias Krack
136 : !> \version 1.0
137 : ! **************************************************************************************************
138 1603879 : SUBROUTINE verfc(la_max1, npgfa, zeta, rpgfa, la_min1, lb_max1, npgfb, zetb, rpgfb, lb_min1, &
139 4811637 : zetc, rpgfc, zc, cerf, rab, rab2, rac, rac2, rbc2, vabc, verf, vnuc, f, &
140 4811637 : maxder, vabc_plus, vnabc, pVp_sum, pVp, dkh_erfc) !JT
141 : INTEGER, INTENT(IN) :: la_max1, npgfa
142 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
143 : INTEGER, INTENT(IN) :: la_min1, lb_max1, npgfb
144 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
145 : INTEGER, INTENT(IN) :: lb_min1
146 : REAL(KIND=dp), INTENT(IN) :: zetc, rpgfc, zc, cerf
147 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
148 : REAL(KIND=dp), INTENT(IN) :: rab2
149 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
150 : REAL(KIND=dp), INTENT(IN) :: rac2, rbc2
151 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
152 : REAL(KIND=dp), DIMENSION(:, :, :) :: verf, vnuc
153 : REAL(KIND=dp), DIMENSION(0:) :: f
154 : INTEGER, INTENT(IN), OPTIONAL :: maxder
155 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus, vnabc, pVp_sum, pVp
156 : LOGICAL, OPTIONAL :: dkh_erfc
157 :
158 : INTEGER :: ax, ay, az, bx, by, bz, coamx, coamy, coamz, coapx, coapy, coapz, cobmx, cobmy, &
159 : cobmz, cobpx, cobpy, cobpz, i, ipgf, j, jpgf, la, la_max, la_min, la_start, lb, lb_max, &
160 : lb_min, maxder_local, n, na, nap, nb, nmax, nxx, nxy, nxz, nyx, nyy, nyz, nzx, nzy, nzz, &
161 : pVpa, pVpb
162 : LOGICAL :: do_dkh
163 : REAL(KIND=dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, fax, &
164 : fay, faz, fbx, fby, fbz, ferf, fnuc, &
165 : ftaz, ftbz, fx, fy, fz, rcp2, t, zetp, &
166 : zetq, zetw
167 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp, rcp, rpw
168 :
169 : ! ---------------------------------------------------------------------------
170 :
171 1603879 : IF (PRESENT(maxder)) THEN
172 20310969 : verf(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
173 20310969 : vnuc(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
174 : END IF
175 :
176 1603879 : do_dkh = .FALSE.
177 :
178 1603879 : IF (PRESENT(pVp)) THEN
179 0 : do_dkh = .TRUE.
180 0 : pVp = 0.0_dp
181 : ! pVp matrices requires angular momentum quantum numbers l+1 and l-1
182 :
183 0 : la_max = la_max1 + 1
184 0 : IF (PRESENT(maxder)) THEN
185 0 : IF (maxder > 0) THEN
186 0 : la_max = la_max1
187 : END IF
188 : END IF
189 0 : lb_max = lb_max1 + 1
190 :
191 0 : la_min = MAX(0, la_min1 - 1)
192 0 : lb_min = MAX(0, lb_min1 - 1)
193 : ELSE
194 1603879 : do_dkh = .FALSE.
195 1603879 : la_max = la_max1
196 1603879 : lb_max = lb_max1
197 1603879 : la_min = la_min1
198 1603879 : lb_min = lb_min1
199 : END IF
200 :
201 1603879 : nmax = la_max + lb_max + 1
202 :
203 1603879 : maxder_local = 0
204 1603879 : IF (PRESENT(maxder)) maxder_local = maxder
205 :
206 : ! *** Calculate the distances of the centers a, b and c ***
207 :
208 1603879 : dab = SQRT(rab2)
209 1603879 : dac = SQRT(rac2)
210 1603879 : dbc = SQRT(rbc2)
211 :
212 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
213 :
214 1603879 : na = 0
215 1603879 : nap = 0
216 :
217 4027201 : DO ipgf = 1, npgfa
218 :
219 : ! *** Screening ***
220 :
221 2423322 : IF (do_dkh) THEN
222 0 : IF (dkh_erfc) THEN
223 0 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
224 0 : na = na + ncoset(la_max1 - maxder_local) !JT
225 0 : nap = nap + ncoset(la_max1) !JT
226 0 : CYCLE
227 : END IF
228 : END IF
229 : ELSE
230 2423322 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
231 522866 : na = na + ncoset(la_max1 - maxder_local) !JT
232 522866 : nap = nap + ncoset(la_max1) !JT
233 522866 : CYCLE
234 : END IF
235 : END IF
236 :
237 1900456 : nb = 0
238 :
239 4960039 : DO jpgf = 1, npgfb
240 :
241 : ! *** Screening ***
242 3059583 : IF (do_dkh) THEN
243 0 : IF (dkh_erfc) THEN
244 0 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
245 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
246 0 : nb = nb + ncoset(lb_max1) !JT
247 0 : CYCLE
248 : END IF
249 : END IF
250 : ELSE
251 3059583 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
252 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
253 781877 : nb = nb + ncoset(lb_max1) !JT
254 781877 : CYCLE
255 : END IF
256 : END IF
257 : ! *** Calculate some prefactors ***
258 :
259 2277706 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
260 2277706 : zetq = 1.0_dp/zetc
261 2277706 : zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
262 :
263 2277706 : f1 = zetb(jpgf)*zetp
264 2277706 : f2 = 0.5_dp*zetp
265 2277706 : f4 = -zetc*zetw
266 :
267 2277706 : f0 = EXP(-zeta(ipgf)*f1*rab2)
268 :
269 2277706 : fnuc = 2.0_dp*pi*zetp*f0
270 2277706 : ferf = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq*f0
271 :
272 9110824 : rap(:) = f1*rab(:)
273 9110824 : rcp(:) = rap(:) - rac(:)
274 9110824 : rpw(:) = f4*rcp(:)
275 :
276 : ! *** Calculate the incomplete Gamma function values ***
277 :
278 2277706 : rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3)
279 :
280 2277706 : t = rcp2/zetp
281 :
282 2277706 : CALL fgamma(nmax - 1, t, f)
283 :
284 : ! *** Calculate the basic nuclear attraction integrals [s|A(0)|s]{n} ***
285 :
286 8307599 : DO n = 1, nmax
287 8307599 : vnuc(1, 1, n) = fnuc*f(n - 1)
288 : END DO
289 :
290 : ! *** Calculate the incomplete Gamma function values ***
291 :
292 2277706 : t = -f4*rcp2/zetp
293 :
294 2277706 : CALL fgamma(nmax - 1, t, f)
295 :
296 : ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
297 :
298 8307599 : DO n = 1, nmax
299 8307599 : verf(1, 1, n) = ferf*f(n - 1)
300 : END DO
301 :
302 : ! *** Recurrence steps: [s|A(0)|s] -> [a|A(0)|b] ***
303 : ! *** [ss||s] -> [ab||s] ***
304 :
305 2277706 : IF (la_max > 0) THEN
306 :
307 : ! *** Vertical recurrence steps: [s|A(0)|s] -> [a|A(0)|s] ***
308 : ! *** [ss||s] -> [as||s] ***
309 :
310 : ! *** [p|A(0)|s]{n} = (Pi - Ai)*[s|A(0)|s]{n} - ***
311 : ! *** (Pi - Ci)*[s|A(0)|s]{n+1} ***
312 : ! *** [ps||s]{n} = (Pi - Ai)*[ss||s]{n} + ***
313 : ! *** (Wi - Pi)*[ss||s]{n+1} (i = x,y,z) ***
314 :
315 4988572 : DO n = 1, nmax - 1
316 3388518 : vnuc(2, 1, n) = rap(1)*vnuc(1, 1, n) - rcp(1)*vnuc(1, 1, n + 1)
317 3388518 : verf(2, 1, n) = rap(1)*verf(1, 1, n) + rpw(1)*verf(1, 1, n + 1)
318 3388518 : vnuc(3, 1, n) = rap(2)*vnuc(1, 1, n) - rcp(2)*vnuc(1, 1, n + 1)
319 3388518 : verf(3, 1, n) = rap(2)*verf(1, 1, n) + rpw(2)*verf(1, 1, n + 1)
320 3388518 : vnuc(4, 1, n) = rap(3)*vnuc(1, 1, n) - rcp(3)*vnuc(1, 1, n + 1)
321 4988572 : verf(4, 1, n) = rap(3)*verf(1, 1, n) + rpw(3)*verf(1, 1, n + 1)
322 : END DO
323 :
324 : ! *** [a|A(0)|s]{n} = (Pi - Ai)*[a-1i|A(0)|s]{n} - ***
325 : ! *** (Pi - Ci)*[a-1i|A(0)|s]{n+1} + ***
326 : ! *** f2*Ni(a-1i)*([a-2i|A(0)|s]{n} - ***
327 : ! *** [a-2i|A(0)|s]{n+1} ***
328 : ! *** [as||s]{n} = (Pi - Ai)*[(a-1i)s||s]{n} + ***
329 : ! *** (Wi - Pi)*[(a-1i)s||s]{n+1} + ***
330 : ! *** f2*Ni(a-1i)*( [(a-2i)s||s]{n} + ***
331 : ! *** f4*[(a-2i)s||s]{n+1}) ***
332 :
333 2266312 : DO la = 2, la_max
334 :
335 3486749 : DO n = 1, nmax - la
336 :
337 : ! *** Increase the angular momentum component z of function a ***
338 :
339 : vnuc(coset(0, 0, la), 1, n) = &
340 : rap(3)*vnuc(coset(0, 0, la - 1), 1, n) - &
341 : rcp(3)*vnuc(coset(0, 0, la - 1), 1, n + 1) + &
342 : f2*REAL(la - 1, dp)*(vnuc(coset(0, 0, la - 2), 1, n) - &
343 1220437 : vnuc(coset(0, 0, la - 2), 1, n + 1))
344 : verf(coset(0, 0, la), 1, n) = &
345 : rap(3)*verf(coset(0, 0, la - 1), 1, n) + &
346 : rpw(3)*verf(coset(0, 0, la - 1), 1, n + 1) + &
347 : f2*REAL(la - 1, dp)*(verf(coset(0, 0, la - 2), 1, n) + &
348 1220437 : f4*verf(coset(0, 0, la - 2), 1, n + 1))
349 :
350 : ! *** Increase the angular momentum component y of function a ***
351 :
352 1220437 : az = la - 1
353 : vnuc(coset(0, 1, az), 1, n) = &
354 : rap(2)*vnuc(coset(0, 0, az), 1, n) - &
355 1220437 : rcp(2)*vnuc(coset(0, 0, az), 1, n + 1)
356 : verf(coset(0, 1, az), 1, n) = &
357 : rap(2)*verf(coset(0, 0, az), 1, n) + &
358 1220437 : rpw(2)*verf(coset(0, 0, az), 1, n + 1)
359 :
360 2558249 : DO ay = 2, la
361 1337812 : az = la - ay
362 : vnuc(coset(0, ay, az), 1, n) = &
363 : rap(2)*vnuc(coset(0, ay - 1, az), 1, n) - &
364 : rcp(2)*vnuc(coset(0, ay - 1, az), 1, n + 1) + &
365 : f2*REAL(ay - 1, dp)*(vnuc(coset(0, ay - 2, az), 1, n) - &
366 1337812 : vnuc(coset(0, ay - 2, az), 1, n + 1))
367 : verf(coset(0, ay, az), 1, n) = &
368 : rap(2)*verf(coset(0, ay - 1, az), 1, n) + &
369 : rpw(2)*verf(coset(0, ay - 1, az), 1, n + 1) + &
370 : f2*REAL(ay - 1, dp)*(verf(coset(0, ay - 2, az), 1, n) + &
371 2558249 : f4*verf(coset(0, ay - 2, az), 1, n + 1))
372 : END DO
373 :
374 : ! *** Increase the angular momentum component x of function a ***
375 :
376 3778686 : DO ay = 0, la - 1
377 2558249 : az = la - 1 - ay
378 : vnuc(coset(1, ay, az), 1, n) = &
379 : rap(1)*vnuc(coset(0, ay, az), 1, n) - &
380 2558249 : rcp(1)*vnuc(coset(0, ay, az), 1, n + 1)
381 : verf(coset(1, ay, az), 1, n) = &
382 : rap(1)*verf(coset(0, ay, az), 1, n) + &
383 3778686 : rpw(1)*verf(coset(0, ay, az), 1, n + 1)
384 : END DO
385 :
386 3224507 : DO ax = 2, la
387 1337812 : f3 = f2*REAL(ax - 1, dp)
388 4014807 : DO ay = 0, la - ax
389 1456558 : az = la - ax - ay
390 : vnuc(coset(ax, ay, az), 1, n) = &
391 : rap(1)*vnuc(coset(ax - 1, ay, az), 1, n) - &
392 : rcp(1)*vnuc(coset(ax - 1, ay, az), 1, n + 1) + &
393 : f3*(vnuc(coset(ax - 2, ay, az), 1, n) - &
394 1456558 : vnuc(coset(ax - 2, ay, az), 1, n + 1))
395 : verf(coset(ax, ay, az), 1, n) = &
396 : rap(1)*verf(coset(ax - 1, ay, az), 1, n) + &
397 : rpw(1)*verf(coset(ax - 1, ay, az), 1, n + 1) + &
398 : f3*(verf(coset(ax - 2, ay, az), 1, n) + &
399 2794370 : f4*verf(coset(ax - 2, ay, az), 1, n + 1))
400 : END DO
401 : END DO
402 :
403 : END DO
404 :
405 : END DO
406 :
407 : ! *** Recurrence steps: [a|A(0)|s] -> [a|A(0)|b] ***
408 : ! *** [as||s] -> [ab||s] ***
409 :
410 1600054 : IF (lb_max > 0) THEN
411 :
412 : ! *** Horizontal recurrence steps ***
413 :
414 3928988 : rbp(:) = rap(:) - rab(:)
415 :
416 : ! *** [a||A(0)|p]{n} = [a+1i|A(0)|s]{n} - (Bi - Ai)*[a|A(0)|s]{n} ***
417 : ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
418 :
419 982247 : la_start = MAX(0, la_min - 1)
420 :
421 2282825 : DO la = la_start, la_max - 1
422 5393929 : DO n = 1, nmax - la - 1
423 8557105 : DO ax = 0, la
424 12528548 : DO ay = 0, la - ax
425 5272021 : az = la - ax - ay
426 : vnuc(coset(ax, ay, az), 2, n) = &
427 : vnuc(coset(ax + 1, ay, az), 1, n) - &
428 5272021 : rab(1)*vnuc(coset(ax, ay, az), 1, n)
429 : verf(coset(ax, ay, az), 2, n) = &
430 : verf(coset(ax + 1, ay, az), 1, n) - &
431 5272021 : rab(1)*verf(coset(ax, ay, az), 1, n)
432 : vnuc(coset(ax, ay, az), 3, n) = &
433 : vnuc(coset(ax, ay + 1, az), 1, n) - &
434 5272021 : rab(2)*vnuc(coset(ax, ay, az), 1, n)
435 : verf(coset(ax, ay, az), 3, n) = &
436 : verf(coset(ax, ay + 1, az), 1, n) - &
437 5272021 : rab(2)*verf(coset(ax, ay, az), 1, n)
438 : vnuc(coset(ax, ay, az), 4, n) = &
439 : vnuc(coset(ax, ay, az + 1), 1, n) - &
440 5272021 : rab(3)*vnuc(coset(ax, ay, az), 1, n)
441 : verf(coset(ax, ay, az), 4, n) = &
442 : verf(coset(ax, ay, az + 1), 1, n) - &
443 9417444 : rab(3)*verf(coset(ax, ay, az), 1, n)
444 : END DO
445 : END DO
446 : END DO
447 : END DO
448 :
449 : ! *** Vertical recurrence step ***
450 :
451 : ! *** [a|A(0)|p]{n} = (Pi - Bi)*[a|A(0)|s]{n} - ***
452 : ! *** (Pi - Ci)*[a|A(0)|s]{n+1} + ***
453 : ! *** f2*Ni(a)*([a-1i|A(0)|s]{n} - ***
454 : ! *** [a-1i|A(0)|s]{n+1}) ***
455 : ! *** [ap||s]{n} = (Pi - Bi)*[as||s]{n} + ***
456 : ! *** (Wi - Pi)*[as||s]{n+1} + ***
457 : ! *** f2*Ni(a)*( [(a-1i)s||s]{n} + ***
458 : ! *** f4*[(a-1i)s||s]{n+1}) ***
459 :
460 2104453 : DO n = 1, nmax - la_max - 1
461 4836249 : DO ax = 0, la_max
462 2731796 : fx = f2*REAL(ax, dp)
463 8732781 : DO ay = 0, la_max - ax
464 4878779 : fy = f2*REAL(ay, dp)
465 4878779 : az = la_max - ax - ay
466 4878779 : fz = f2*REAL(az, dp)
467 :
468 4878779 : IF (ax == 0) THEN
469 : vnuc(coset(ax, ay, az), 2, n) = &
470 : rbp(1)*vnuc(coset(ax, ay, az), 1, n) - &
471 2731796 : rcp(1)*vnuc(coset(ax, ay, az), 1, n + 1)
472 : verf(coset(ax, ay, az), 2, n) = &
473 : rbp(1)*verf(coset(ax, ay, az), 1, n) + &
474 2731796 : rpw(1)*verf(coset(ax, ay, az), 1, n + 1)
475 : ELSE
476 : vnuc(coset(ax, ay, az), 2, n) = &
477 : rbp(1)*vnuc(coset(ax, ay, az), 1, n) - &
478 : rcp(1)*vnuc(coset(ax, ay, az), 1, n + 1) + &
479 : fx*(vnuc(coset(ax - 1, ay, az), 1, n) - &
480 2146983 : vnuc(coset(ax - 1, ay, az), 1, n + 1))
481 : verf(coset(ax, ay, az), 2, n) = &
482 : rbp(1)*verf(coset(ax, ay, az), 1, n) + &
483 : rpw(1)*verf(coset(ax, ay, az), 1, n + 1) + &
484 : fx*(verf(coset(ax - 1, ay, az), 1, n) + &
485 2146983 : f4*verf(coset(ax - 1, ay, az), 1, n + 1))
486 : END IF
487 :
488 4878779 : IF (ay == 0) THEN
489 : vnuc(coset(ax, ay, az), 3, n) = &
490 : rbp(2)*vnuc(coset(ax, ay, az), 1, n) - &
491 2731796 : rcp(2)*vnuc(coset(ax, ay, az), 1, n + 1)
492 : verf(coset(ax, ay, az), 3, n) = &
493 : rbp(2)*verf(coset(ax, ay, az), 1, n) + &
494 2731796 : rpw(2)*verf(coset(ax, ay, az), 1, n + 1)
495 : ELSE
496 : vnuc(coset(ax, ay, az), 3, n) = &
497 : rbp(2)*vnuc(coset(ax, ay, az), 1, n) - &
498 : rcp(2)*vnuc(coset(ax, ay, az), 1, n + 1) + &
499 : fy*(vnuc(coset(ax, ay - 1, az), 1, n) - &
500 2146983 : vnuc(coset(ax, ay - 1, az), 1, n + 1))
501 : verf(coset(ax, ay, az), 3, n) = &
502 : rbp(2)*verf(coset(ax, ay, az), 1, n) + &
503 : rpw(2)*verf(coset(ax, ay, az), 1, n + 1) + &
504 : fy*(verf(coset(ax, ay - 1, az), 1, n) + &
505 2146983 : f4*verf(coset(ax, ay - 1, az), 1, n + 1))
506 : END IF
507 :
508 7610575 : IF (az == 0) THEN
509 : vnuc(coset(ax, ay, az), 4, n) = &
510 : rbp(3)*vnuc(coset(ax, ay, az), 1, n) - &
511 2731796 : rcp(3)*vnuc(coset(ax, ay, az), 1, n + 1)
512 : verf(coset(ax, ay, az), 4, n) = &
513 : rbp(3)*verf(coset(ax, ay, az), 1, n) + &
514 2731796 : rpw(3)*verf(coset(ax, ay, az), 1, n + 1)
515 : ELSE
516 : vnuc(coset(ax, ay, az), 4, n) = &
517 : rbp(3)*vnuc(coset(ax, ay, az), 1, n) - &
518 : rcp(3)*vnuc(coset(ax, ay, az), 1, n + 1) + &
519 : fz*(vnuc(coset(ax, ay, az - 1), 1, n) - &
520 2146983 : vnuc(coset(ax, ay, az - 1), 1, n + 1))
521 : verf(coset(ax, ay, az), 4, n) = &
522 : rbp(3)*verf(coset(ax, ay, az), 1, n) + &
523 : rpw(3)*verf(coset(ax, ay, az), 1, n + 1) + &
524 : fz*(verf(coset(ax, ay, az - 1), 1, n) + &
525 2146983 : f4*verf(coset(ax, ay, az - 1), 1, n + 1))
526 : END IF
527 :
528 : END DO
529 : END DO
530 : END DO
531 :
532 : ! *** Recurrence steps: [a|A(0)|p] -> [a|A(0)|b] ***
533 : ! *** [ap||s] -> [ab||s] ***
534 :
535 1122206 : DO lb = 2, lb_max
536 :
537 : ! *** Horizontal recurrence steps ***
538 :
539 : ! *** [a||A(0)|b]{n} = [a+1i|A(0)|b-1i]{n} - ***
540 : ! *** (Bi - Ai)*[a|A(0)|b-1i]{n} ***
541 : ! *** [ab||s]{n} = [(a+1i)(b-1i)||s]{n} - ***
542 : ! *** (Bi - Ai)*[a(b-1i)||s]{n} ***
543 :
544 330526 : la_start = MAX(0, la_min - 1)
545 :
546 330526 : DO la = la_start, la_max - 1
547 765085 : DO n = 1, nmax - la - lb
548 1214414 : DO ax = 0, la
549 1784317 : DO ay = 0, la - ax
550 760470 : az = la - ax - ay
551 :
552 : ! *** Shift of angular momentum component z from a to b ***
553 :
554 : vnuc(coset(ax, ay, az), coset(0, 0, lb), n) = &
555 : vnuc(coset(ax, ay, az + 1), coset(0, 0, lb - 1), n) - &
556 760470 : rab(3)*vnuc(coset(ax, ay, az), coset(0, 0, lb - 1), n)
557 : verf(coset(ax, ay, az), coset(0, 0, lb), n) = &
558 : verf(coset(ax, ay, az + 1), coset(0, 0, lb - 1), n) - &
559 760470 : rab(3)*verf(coset(ax, ay, az), coset(0, 0, lb - 1), n)
560 :
561 : ! *** Shift of angular momentum component y from a to b ***
562 :
563 2295605 : DO by = 1, lb
564 1535135 : bz = lb - by
565 : vnuc(coset(ax, ay, az), coset(0, by, bz), n) = &
566 : vnuc(coset(ax, ay + 1, az), coset(0, by - 1, bz), n) - &
567 1535135 : rab(2)*vnuc(coset(ax, ay, az), coset(0, by - 1, bz), n)
568 : verf(coset(ax, ay, az), coset(0, by, bz), n) = &
569 : verf(coset(ax, ay + 1, az), coset(0, by - 1, bz), n) - &
570 2295605 : rab(2)*verf(coset(ax, ay, az), coset(0, by - 1, bz), n)
571 : END DO
572 :
573 : ! *** Shift of angular momentum component x from a to b ***
574 :
575 2884893 : DO bx = 1, lb
576 4622654 : DO by = 0, lb - bx
577 2327049 : bz = lb - bx - by
578 : vnuc(coset(ax, ay, az), coset(bx, by, bz), n) = &
579 : vnuc(coset(ax + 1, ay, az), coset(bx - 1, by, bz), n) - &
580 2327049 : rab(1)*vnuc(coset(ax, ay, az), coset(bx - 1, by, bz), n)
581 : verf(coset(ax, ay, az), coset(bx, by, bz), n) = &
582 : verf(coset(ax + 1, ay, az), coset(bx - 1, by, bz), n) - &
583 3862184 : rab(1)*verf(coset(ax, ay, az), coset(bx - 1, by, bz), n)
584 : END DO
585 : END DO
586 :
587 : END DO
588 : END DO
589 : END DO
590 : END DO
591 :
592 : ! *** Vertical recurrence step ***
593 :
594 : ! *** [a|A(0)|b]{n} = (Pi - Bi)*[a|A(0)|b-1i]{n} - ***
595 : ! *** (Pi - Ci)*[a|A(0)|b-1i]{n+1} + ***
596 : ! *** f2*Ni(a)*([a-1i|A(0)|b-1i]{n} - ***
597 : ! *** [a-1i|A(0)|b-1i]{n+1}) + ***
598 : ! *** f2*Ni(b-1i)*([a|A(0)|b-2i]{n} - ***
599 : ! *** [a|A(0)|b-2i]{n+1}) ***
600 : ! *** [ab||s]{n} = (Pi - Bi)*[a(b-1i)||s]{n} + ***
601 : ! *** (Wi - Pi)*[a(b-1i)||s]{n+1} + ***
602 : ! *** f2*Ni(a)*( [(a-1i)(b-1i)||s]{n} + ***
603 : ! *** f4*[(a-1i)(b-1i)||s]{n+1}) ***
604 : ! *** f2*Ni(b-1i)*( [a(b-2i)||s]{n} + ***
605 : ! *** f4*[a(b-2i)||s]{n+1}) ***
606 :
607 1264100 : DO n = 1, nmax - la_max - lb
608 634872 : DO ax = 0, la_max
609 353019 : fx = f2*REAL(ax, dp)
610 1137002 : DO ay = 0, la_max - ax
611 642089 : fy = f2*REAL(ay, dp)
612 642089 : az = la_max - ax - ay
613 642089 : fz = f2*REAL(az, dp)
614 :
615 : ! *** Shift of angular momentum component z from a to b ***
616 :
617 642089 : f3 = f2*REAL(lb - 1, dp)
618 :
619 642089 : IF (az == 0) THEN
620 : vnuc(coset(ax, ay, az), coset(0, 0, lb), n) = &
621 : rbp(3)*vnuc(coset(ax, ay, az), coset(0, 0, lb - 1), n) - &
622 : rcp(3)*vnuc(coset(ax, ay, az), coset(0, 0, lb - 1), n + 1) + &
623 : f3*(vnuc(coset(ax, ay, az), coset(0, 0, lb - 2), n) - &
624 353019 : vnuc(coset(ax, ay, az), coset(0, 0, lb - 2), n + 1))
625 : verf(coset(ax, ay, az), coset(0, 0, lb), n) = &
626 : rbp(3)*verf(coset(ax, ay, az), coset(0, 0, lb - 1), n) + &
627 : rpw(3)*verf(coset(ax, ay, az), coset(0, 0, lb - 1), n + 1) + &
628 : f3*(verf(coset(ax, ay, az), coset(0, 0, lb - 2), n) + &
629 353019 : f4*verf(coset(ax, ay, az), coset(0, 0, lb - 2), n + 1))
630 : ELSE
631 : vnuc(coset(ax, ay, az), coset(0, 0, lb), n) = &
632 : rbp(3)*vnuc(coset(ax, ay, az), coset(0, 0, lb - 1), n) - &
633 : rcp(3)*vnuc(coset(ax, ay, az), coset(0, 0, lb - 1), n + 1) + &
634 : fz*(vnuc(coset(ax, ay, az - 1), coset(0, 0, lb - 1), n) - &
635 : vnuc(coset(ax, ay, az - 1), coset(0, 0, lb - 1), n + 1)) + &
636 : f3*(vnuc(coset(ax, ay, az), coset(0, 0, lb - 2), n) - &
637 289070 : vnuc(coset(ax, ay, az), coset(0, 0, lb - 2), n + 1))
638 : verf(coset(ax, ay, az), coset(0, 0, lb), n) = &
639 : rbp(3)*verf(coset(ax, ay, az), coset(0, 0, lb - 1), n) + &
640 : rpw(3)*verf(coset(ax, ay, az), coset(0, 0, lb - 1), n + 1) + &
641 : fz*(verf(coset(ax, ay, az - 1), coset(0, 0, lb - 1), n) + &
642 : f4*verf(coset(ax, ay, az - 1), coset(0, 0, lb - 1), n + 1)) + &
643 : f3*(verf(coset(ax, ay, az), coset(0, 0, lb - 2), n) + &
644 289070 : f4*verf(coset(ax, ay, az), coset(0, 0, lb - 2), n + 1))
645 : END IF
646 :
647 : ! *** Shift of angular momentum component y from a to b ***
648 :
649 642089 : IF (ay == 0) THEN
650 353019 : bz = lb - 1
651 : vnuc(coset(ax, ay, az), coset(0, 1, bz), n) = &
652 : rbp(2)*vnuc(coset(ax, ay, az), coset(0, 0, bz), n) - &
653 353019 : rcp(2)*vnuc(coset(ax, ay, az), coset(0, 0, bz), n + 1)
654 : verf(coset(ax, ay, az), coset(0, 1, bz), n) = &
655 : rbp(2)*verf(coset(ax, ay, az), coset(0, 0, bz), n) + &
656 353019 : rpw(2)*verf(coset(ax, ay, az), coset(0, 0, bz), n + 1)
657 712705 : DO by = 2, lb
658 359686 : bz = lb - by
659 359686 : f3 = f2*REAL(by - 1, dp)
660 : vnuc(coset(ax, ay, az), coset(0, by, bz), n) = &
661 : rbp(2)*vnuc(coset(ax, ay, az), coset(0, by - 1, bz), n) - &
662 : rcp(2)*vnuc(coset(ax, ay, az), coset(0, by - 1, bz), n + 1) + &
663 : f3*(vnuc(coset(ax, ay, az), coset(0, by - 2, bz), n) - &
664 359686 : vnuc(coset(ax, ay, az), coset(0, by - 2, bz), n + 1))
665 : verf(coset(ax, ay, az), coset(0, by, bz), n) = &
666 : rbp(2)*verf(coset(ax, ay, az), coset(0, by - 1, bz), n) + &
667 : rpw(2)*verf(coset(ax, ay, az), coset(0, by - 1, bz), n + 1) + &
668 : f3*(verf(coset(ax, ay, az), coset(0, by - 2, bz), n) + &
669 712705 : f4*verf(coset(ax, ay, az), coset(0, by - 2, bz), n + 1))
670 : END DO
671 : ELSE
672 289070 : bz = lb - 1
673 : vnuc(coset(ax, ay, az), coset(0, 1, bz), n) = &
674 : rbp(2)*vnuc(coset(ax, ay, az), coset(0, 0, bz), n) - &
675 : rcp(2)*vnuc(coset(ax, ay, az), coset(0, 0, bz), n + 1) + &
676 : fy*(vnuc(coset(ax, ay - 1, az), coset(0, 0, bz), n) - &
677 289070 : vnuc(coset(ax, ay - 1, az), coset(0, 0, bz), n + 1))
678 : verf(coset(ax, ay, az), coset(0, 1, bz), n) = &
679 : rbp(2)*verf(coset(ax, ay, az), coset(0, 0, bz), n) + &
680 : rpw(2)*verf(coset(ax, ay, az), coset(0, 0, bz), n + 1) + &
681 : fy*(verf(coset(ax, ay - 1, az), coset(0, 0, bz), n) + &
682 289070 : f4*verf(coset(ax, ay - 1, az), coset(0, 0, bz), n + 1))
683 585947 : DO by = 2, lb
684 296877 : bz = lb - by
685 296877 : f3 = f2*REAL(by - 1, dp)
686 : vnuc(coset(ax, ay, az), coset(0, by, bz), n) = &
687 : rbp(2)*vnuc(coset(ax, ay, az), coset(0, by - 1, bz), n) - &
688 : rcp(2)*vnuc(coset(ax, ay, az), coset(0, by - 1, bz), n + 1) + &
689 : fy*(vnuc(coset(ax, ay - 1, az), coset(0, by - 1, bz), n) - &
690 : vnuc(coset(ax, ay - 1, az), coset(0, by - 1, bz), n + 1)) + &
691 : f3*(vnuc(coset(ax, ay, az), coset(0, by - 2, bz), n) - &
692 296877 : vnuc(coset(ax, ay, az), coset(0, by - 2, bz), n + 1))
693 : verf(coset(ax, ay, az), coset(0, by, bz), n) = &
694 : rbp(2)*verf(coset(ax, ay, az), coset(0, by - 1, bz), n) + &
695 : rpw(2)*verf(coset(ax, ay, az), coset(0, by - 1, bz), n + 1) + &
696 : fy*(verf(coset(ax, ay - 1, az), coset(0, by - 1, bz), n) + &
697 : f4*verf(coset(ax, ay - 1, az), &
698 : coset(0, by - 1, bz), n + 1)) + &
699 : f3*(verf(coset(ax, ay, az), coset(0, by - 2, bz), n) + &
700 585947 : f4*verf(coset(ax, ay, az), coset(0, by - 2, bz), n + 1))
701 : END DO
702 : END IF
703 :
704 : ! *** Shift of angular momentum component x from a to b ***
705 :
706 995108 : IF (ax == 0) THEN
707 1065724 : DO by = 0, lb - 1
708 712705 : bz = lb - 1 - by
709 : vnuc(coset(ax, ay, az), coset(1, by, bz), n) = &
710 : rbp(1)*vnuc(coset(ax, ay, az), coset(0, by, bz), n) - &
711 712705 : rcp(1)*vnuc(coset(ax, ay, az), coset(0, by, bz), n + 1)
712 : verf(coset(ax, ay, az), coset(1, by, bz), n) = &
713 : rbp(1)*verf(coset(ax, ay, az), coset(0, by, bz), n) + &
714 1065724 : rpw(1)*verf(coset(ax, ay, az), coset(0, by, bz), n + 1)
715 : END DO
716 712705 : DO bx = 2, lb
717 359686 : f3 = f2*REAL(bx - 1, dp)
718 1080268 : DO by = 0, lb - bx
719 367563 : bz = lb - bx - by
720 : vnuc(coset(ax, ay, az), coset(bx, by, bz), n) = &
721 : rbp(1)*vnuc(coset(ax, ay, az), coset(bx - 1, by, bz), n) - &
722 : rcp(1)*vnuc(coset(ax, ay, az), &
723 : coset(bx - 1, by, bz), n + 1) + &
724 : f3*(vnuc(coset(ax, ay, az), coset(bx - 2, by, bz), n) - &
725 367563 : vnuc(coset(ax, ay, az), coset(bx - 2, by, bz), n + 1))
726 : verf(coset(ax, ay, az), coset(bx, by, bz), n) = &
727 : rbp(1)*verf(coset(ax, ay, az), coset(bx - 1, by, bz), n) + &
728 : rpw(1)*verf(coset(ax, ay, az), &
729 : coset(bx - 1, by, bz), n + 1) + &
730 : f3*(verf(coset(ax, ay, az), coset(bx - 2, by, bz), n) + &
731 727249 : f4*verf(coset(ax, ay, az), coset(bx - 2, by, bz), n + 1))
732 : END DO
733 : END DO
734 : ELSE
735 875017 : DO by = 0, lb - 1
736 585947 : bz = lb - 1 - by
737 : vnuc(coset(ax, ay, az), coset(1, by, bz), n) = &
738 : rbp(1)*vnuc(coset(ax, ay, az), coset(0, by, bz), n) - &
739 : rcp(1)*vnuc(coset(ax, ay, az), coset(0, by, bz), n + 1) + &
740 : fx*(vnuc(coset(ax - 1, ay, az), coset(0, by, bz), n) - &
741 585947 : vnuc(coset(ax - 1, ay, az), coset(0, by, bz), n + 1))
742 : verf(coset(ax, ay, az), coset(1, by, bz), n) = &
743 : rbp(1)*verf(coset(ax, ay, az), coset(0, by, bz), n) + &
744 : rpw(1)*verf(coset(ax, ay, az), coset(0, by, bz), n + 1) + &
745 : fx*(verf(coset(ax - 1, ay, az), coset(0, by, bz), n) + &
746 875017 : f4*verf(coset(ax - 1, ay, az), coset(0, by, bz), n + 1))
747 : END DO
748 585947 : DO bx = 2, lb
749 296877 : f3 = f2*REAL(bx - 1, dp)
750 892266 : DO by = 0, lb - bx
751 306319 : bz = lb - bx - by
752 : vnuc(coset(ax, ay, az), coset(bx, by, bz), n) = &
753 : rbp(1)*vnuc(coset(ax, ay, az), coset(bx - 1, by, bz), n) - &
754 : rcp(1)*vnuc(coset(ax, ay, az), &
755 : coset(bx - 1, by, bz), n + 1) + &
756 : fx*(vnuc(coset(ax - 1, ay, az), coset(bx - 1, by, bz), n) - &
757 : vnuc(coset(ax - 1, ay, az), &
758 : coset(bx - 1, by, bz), n + 1)) + &
759 : f3*(vnuc(coset(ax, ay, az), coset(bx - 2, by, bz), n) - &
760 306319 : vnuc(coset(ax, ay, az), coset(bx - 2, by, bz), n + 1))
761 : verf(coset(ax, ay, az), coset(bx, by, bz), n) = &
762 : rbp(1)*verf(coset(ax, ay, az), coset(bx - 1, by, bz), n) + &
763 : rpw(1)*verf(coset(ax, ay, az), &
764 : coset(bx - 1, by, bz), n + 1) + &
765 : fx*(verf(coset(ax - 1, ay, az), &
766 : coset(bx - 1, by, bz), n) + &
767 : f4*verf(coset(ax - 1, ay, az), &
768 : coset(bx - 1, by, bz), n + 1)) + &
769 : f3*(verf(coset(ax, ay, az), coset(bx - 2, by, bz), n) + &
770 603196 : f4*verf(coset(ax, ay, az), coset(bx - 2, by, bz), n + 1))
771 : END DO
772 : END DO
773 : END IF
774 :
775 : END DO
776 : END DO
777 : END DO
778 :
779 : END DO
780 :
781 : END IF
782 :
783 : ELSE
784 :
785 677652 : IF (lb_max > 0) THEN
786 :
787 : ! *** Vertical recurrence steps: [s|A(0)|s] -> [s|A(0)|b] ***
788 : ! *** [ss||s] -> [sb||s] ***
789 :
790 1274008 : rbp(:) = rap(:) - rab(:)
791 :
792 : ! *** [s|A(0)|p]{n} = (Pi - Bi)*[s|A(0)|s]{n} - ***
793 : ! *** (Pi - Ci)*[s|A(0)|s]{n+1} ***
794 : ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
795 : ! *** (Wi - Pi)*[ss||s]{n+1} ***
796 :
797 682171 : DO n = 1, nmax - 1
798 363669 : vnuc(1, 2, n) = rbp(1)*vnuc(1, 1, n) - rcp(1)*vnuc(1, 1, n + 1)
799 363669 : verf(1, 2, n) = rbp(1)*verf(1, 1, n) + rpw(1)*verf(1, 1, n + 1)
800 363669 : vnuc(1, 3, n) = rbp(2)*vnuc(1, 1, n) - rcp(2)*vnuc(1, 1, n + 1)
801 363669 : verf(1, 3, n) = rbp(2)*verf(1, 1, n) + rpw(2)*verf(1, 1, n + 1)
802 363669 : vnuc(1, 4, n) = rbp(3)*vnuc(1, 1, n) - rcp(3)*vnuc(1, 1, n + 1)
803 682171 : verf(1, 4, n) = rbp(3)*verf(1, 1, n) + rpw(3)*verf(1, 1, n + 1)
804 : END DO
805 :
806 : ! *** [s|A(0)|b]{n} = (Pi - Bi)*[s|A(0)|b-1i]{n} - ***
807 : ! *** (Pi - Ci)*[s|A(0)|b-1i]{n+1} + ***
808 : ! *** f2*Ni(b-1i)*([s|A(0)|b-2i]{n} - ***
809 : ! *** [s|A(0)|b-2i]{n+1} ***
810 : ! *** [sb||s]{n} = (Pi - Bi)*[s(b-1i)||s]{n} + ***
811 : ! *** (Wi - Pi)*[s(b-1i)||s]{n+1} + ***
812 : ! *** f2*Ni(b-1i)*( [s(b-2i)||s]{n} + ***
813 : ! *** f4*[s(b-2i)||s]{n+1}) ***
814 :
815 363669 : DO lb = 2, lb_max
816 :
817 410229 : DO n = 1, nmax - lb
818 :
819 : ! *** Increase the angular momentum component z of function b ***
820 :
821 : vnuc(1, coset(0, 0, lb), n) = &
822 : rbp(3)*vnuc(1, coset(0, 0, lb - 1), n) - &
823 : rcp(3)*vnuc(1, coset(0, 0, lb - 1), n + 1) + &
824 : f2*REAL(lb - 1, dp)*(vnuc(1, coset(0, 0, lb - 2), n) - &
825 46560 : vnuc(1, coset(0, 0, lb - 2), n + 1))
826 : verf(1, coset(0, 0, lb), n) = &
827 : rbp(3)*verf(1, coset(0, 0, lb - 1), n) + &
828 : rpw(3)*verf(1, coset(0, 0, lb - 1), n + 1) + &
829 : f2*REAL(lb - 1, dp)*(verf(1, coset(0, 0, lb - 2), n) + &
830 46560 : f4*verf(1, coset(0, 0, lb - 2), n + 1))
831 :
832 : ! *** Increase the angular momentum component y of function b ***
833 :
834 46560 : bz = lb - 1
835 : vnuc(1, coset(0, 1, bz), n) = &
836 : rbp(2)*vnuc(1, coset(0, 0, bz), n) - &
837 46560 : rcp(2)*vnuc(1, coset(0, 0, bz), n + 1)
838 : verf(1, coset(0, 1, bz), n) = &
839 : rbp(2)*verf(1, coset(0, 0, bz), n) + &
840 46560 : rpw(2)*verf(1, coset(0, 0, bz), n + 1)
841 :
842 94750 : DO by = 2, lb
843 48190 : bz = lb - by
844 : vnuc(1, coset(0, by, bz), n) = &
845 : rbp(2)*vnuc(1, coset(0, by - 1, bz), n) - &
846 : rcp(2)*vnuc(1, coset(0, by - 1, bz), n + 1) + &
847 : f2*REAL(by - 1, dp)*(vnuc(1, coset(0, by - 2, bz), n) - &
848 48190 : vnuc(1, coset(0, by - 2, bz), n + 1))
849 : verf(1, coset(0, by, bz), n) = &
850 : rbp(2)*verf(1, coset(0, by - 1, bz), n) + &
851 : rpw(2)*verf(1, coset(0, by - 1, bz), n + 1) + &
852 : f2*REAL(by - 1, dp)*(verf(1, coset(0, by - 2, bz), n) + &
853 94750 : f4*verf(1, coset(0, by - 2, bz), n + 1))
854 : END DO
855 :
856 : ! *** Increase the angular momentum component x of function b ***
857 :
858 141310 : DO by = 0, lb - 1
859 94750 : bz = lb - 1 - by
860 : vnuc(1, coset(1, by, bz), n) = &
861 : rbp(1)*vnuc(1, coset(0, by, bz), n) - &
862 94750 : rcp(1)*vnuc(1, coset(0, by, bz), n + 1)
863 : verf(1, coset(1, by, bz), n) = &
864 : rbp(1)*verf(1, coset(0, by, bz), n) + &
865 141310 : rpw(1)*verf(1, coset(0, by, bz), n + 1)
866 : END DO
867 :
868 139917 : DO bx = 2, lb
869 48190 : f3 = f2*REAL(bx - 1, dp)
870 144837 : DO by = 0, lb - bx
871 50087 : bz = lb - bx - by
872 : vnuc(1, coset(bx, by, bz), n) = &
873 : rbp(1)*vnuc(1, coset(bx - 1, by, bz), n) - &
874 : rcp(1)*vnuc(1, coset(bx - 1, by, bz), n + 1) + &
875 : f3*(vnuc(1, coset(bx - 2, by, bz), n) - &
876 50087 : vnuc(1, coset(bx - 2, by, bz), n + 1))
877 : verf(1, coset(bx, by, bz), n) = &
878 : rbp(1)*verf(1, coset(bx - 1, by, bz), n) + &
879 : rpw(1)*verf(1, coset(bx - 1, by, bz), n + 1) + &
880 : f3*(verf(1, coset(bx - 2, by, bz), n) + &
881 98277 : f4*verf(1, coset(bx - 2, by, bz), n + 1))
882 : END DO
883 : END DO
884 :
885 : END DO
886 :
887 : END DO
888 :
889 : END IF
890 :
891 : END IF
892 :
893 : ! *** Add the contribution of the current pair ***
894 : ! *** of primitive Gaussian-type functions ***
895 :
896 : !JT
897 2277706 : IF (do_dkh) THEN
898 0 : DO j = 1, ncoset(lb_max1)
899 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
900 0 : vnabc(na + i, nb + j) = vnabc(na + i, nb + j) + cerf*verf(i, j, 1)
901 : END DO
902 : END DO
903 0 : DO j = 1, ncoset(lb_max1)
904 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
905 0 : vabc(na + i, nb + j) = vabc(na + i, nb + j) - zc*vnuc(i, j, 1)
906 : END DO
907 : END DO
908 : ELSE
909 9582298 : DO j = 1, ncoset(lb_max1) !JT
910 29023037 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local) !JT
911 : vabc(na + i, nb + j) = vabc(na + i, nb + j) - &
912 26745331 : zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
913 : END DO
914 : END DO
915 : END IF
916 : !JT
917 :
918 2277706 : IF (PRESENT(maxder)) THEN
919 3409539 : DO j = 1, ncoset(lb_max1) !JT
920 25909253 : DO i = 1, ncoset(la_max1) !JT
921 : vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) - &
922 25121480 : zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
923 : END DO
924 : END DO
925 : END IF
926 : !JTs
927 2277706 : IF (do_dkh) THEN
928 : ! *********** Calculate pVp matrix **************
929 : ! [pa|V|pb]=4*zeta*zetb*[a+i|V|b+i]-2*zetb*Ni(a)*[a-i|V|b+i]-2*zeta*Ni(b)*[a+i|V|b-i]+Ni(a)*Ni(b)*[a-i|V|b-i]
930 : ! Every integral has been calculated before (except [-1|V|-1],[a|V|-1] and [-1|V|b],
931 : ! which do not contribute, because their prefactor Ni(a) or Ni(b) is zero)
932 : ! and is given by verf/vnuc(coset(ax,ay,az),coset(bx,by,bz),1)
933 : ! ***********************************************
934 0 : ftaz = 2.0_dp*zeta(ipgf)
935 0 : ftbz = 2.0_dp*zetb(jpgf)
936 0 : nxx = 1
937 0 : nyy = 2
938 0 : nzz = 3
939 0 : nxy = 4
940 0 : nxz = 5
941 0 : nyx = 6
942 0 : nyz = 7
943 0 : nzx = 8
944 0 : nzy = 9
945 0 : DO la = 0, la_max1
946 0 : DO ax = 0, la
947 0 : fax = REAL(ax, dp)
948 0 : DO ay = 0, la - ax
949 0 : fay = REAL(ay, dp)
950 0 : az = la - ax - ay
951 0 : faz = REAL(az, dp)
952 0 : pVpa = coset(ax, ay, az)
953 0 : coamx = coset(ax - 1, ay, az)
954 0 : coamy = coset(ax, ay - 1, az)
955 0 : coamz = coset(ax, ay, az - 1)
956 0 : coapx = coset(ax + 1, ay, az)
957 0 : coapy = coset(ax, ay + 1, az)
958 0 : coapz = coset(ax, ay, az + 1)
959 0 : DO lb = 0, lb_max1
960 0 : DO bx = 0, lb
961 0 : fbx = REAL(bx, dp)
962 0 : DO by = 0, lb - bx
963 0 : fby = REAL(by, dp)
964 0 : bz = lb - bx - by
965 0 : fbz = REAL(bz, dp)
966 0 : pVpb = coset(bx, by, bz)
967 0 : cobmx = coset(bx - 1, by, bz)
968 0 : cobmy = coset(bx, by - 1, bz)
969 0 : cobmz = coset(bx, by, bz - 1)
970 0 : cobpx = coset(bx + 1, by, bz)
971 0 : cobpy = coset(bx, by + 1, bz)
972 0 : cobpz = coset(bx, by, bz + 1)
973 0 : IF (dkh_erfc) THEN
974 : pVp(pVpa, pVpb) = ftaz*ftbz*(-zc*vnuc(coapx, cobpx, 1) + cerf*verf(coapx, cobpx, 1)) - &
975 : ftaz*fbx*(-zc*vnuc(coapx, cobmx, 1) + cerf*verf(coapx, cobmx, 1)) - &
976 : ftbz*fax*(-zc*vnuc(coamx, cobpx, 1) + cerf*verf(coamx, cobpx, 1)) + &
977 : fax*fbx*(-zc*vnuc(coamx, cobmx, 1) + cerf*verf(coamx, cobmx, 1)) + &
978 : ftaz*ftbz*(-zc*vnuc(coapy, cobpy, 1) + cerf*verf(coapy, cobpy, 1)) - &
979 : ftaz*fby*(-zc*vnuc(coapy, cobmy, 1) + cerf*verf(coapy, cobmy, 1)) - &
980 : ftbz*fay*(-zc*vnuc(coamy, cobpy, 1) + cerf*verf(coamy, cobpy, 1)) + &
981 : fay*fby*(-zc*vnuc(coamy, cobmy, 1) + cerf*verf(coamy, cobmy, 1)) + &
982 : ftaz*ftbz*(-zc*vnuc(coapz, cobpz, 1) + cerf*verf(coapz, cobpz, 1)) - &
983 : ftaz*fbz*(-zc*vnuc(coapz, cobmz, 1) + cerf*verf(coapz, cobmz, 1)) - &
984 : ftbz*faz*(-zc*vnuc(coamz, cobpz, 1) + cerf*verf(coamz, cobpz, 1)) + &
985 0 : faz*fbz*(-zc*vnuc(coamz, cobmz, 1) + cerf*verf(coamz, cobmz, 1))
986 : ELSE
987 : pVp(pVpa, pVpb) = ftaz*ftbz*(-zc*vnuc(coapx, cobpx, 1)) - &
988 : ftaz*fbx*(-zc*vnuc(coapx, cobmx, 1)) - &
989 : ftbz*fax*(-zc*vnuc(coamx, cobpx, 1)) + &
990 : fax*fbx*(-zc*vnuc(coamx, cobmx, 1)) + &
991 : ftaz*ftbz*(-zc*vnuc(coapy, cobpy, 1)) - &
992 : ftaz*fby*(-zc*vnuc(coapy, cobmy, 1)) - &
993 : ftbz*fay*(-zc*vnuc(coamy, cobpy, 1)) + &
994 : fay*fby*(-zc*vnuc(coamy, cobmy, 1)) + &
995 : ftaz*ftbz*(-zc*vnuc(coapz, cobpz, 1)) - &
996 : ftaz*fbz*(-zc*vnuc(coapz, cobmz, 1)) - &
997 : ftbz*faz*(-zc*vnuc(coamz, cobpz, 1)) + &
998 0 : faz*fbz*(-zc*vnuc(coamz, cobmz, 1))
999 : END IF
1000 : END DO
1001 : END DO
1002 : END DO
1003 : END DO
1004 : END DO
1005 : END DO
1006 :
1007 0 : DO j = 1, ncoset(lb_max1)
1008 0 : DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
1009 0 : pVp_sum(na + i, nb + j) = pVp_sum(na + i, nb + j) + pVp(i, j)
1010 : END DO
1011 : END DO
1012 : END IF
1013 : !JTe
1014 4178162 : nb = nb + ncoset(lb_max1) !JT
1015 :
1016 : END DO
1017 :
1018 1900456 : na = na + ncoset(la_max1 - maxder_local) !JT
1019 3504335 : nap = nap + ncoset(la_max1) !JT
1020 :
1021 : END DO
1022 :
1023 1603879 : END SUBROUTINE verfc
1024 :
1025 : END MODULE ai_verfc
|
