(git:b279b6b)
ai_verfc.F
Go to the documentation of this file.
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  SUBROUTINE verfc(la_max1, npgfa, zeta, rpgfa, la_min1, lb_max1, npgfb, zetb, rpgfb, lb_min1, &
139  zetc, rpgfc, zc, cerf, rab, rab2, rac, rac2, rbc2, vabc, verf, vnuc, f, &
140  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  IF (PRESENT(maxder)) THEN
172  verf(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
173  vnuc(1:ncoset(la_max1), 1:ncoset(lb_max1), 1) = 0._dp
174  END IF
175 
176  do_dkh = .false.
177 
178  IF (PRESENT(pvp)) THEN
179  do_dkh = .true.
180  pvp = 0.0_dp
181  ! pVp matrices requires angular momentum quantum numbers l+1 and l-1
182 
183  la_max = la_max1 + 1
184  IF (PRESENT(maxder)) THEN
185  IF (maxder > 0) THEN
186  la_max = la_max1
187  END IF
188  END IF
189  lb_max = lb_max1 + 1
190 
191  la_min = max(0, la_min1 - 1)
192  lb_min = max(0, lb_min1 - 1)
193  ELSE
194  do_dkh = .false.
195  la_max = la_max1
196  lb_max = lb_max1
197  la_min = la_min1
198  lb_min = lb_min1
199  END IF
200 
201  nmax = la_max + lb_max + 1
202 
203  maxder_local = 0
204  IF (PRESENT(maxder)) maxder_local = maxder
205 
206 ! *** Calculate the distances of the centers a, b and c ***
207 
208  dab = sqrt(rab2)
209  dac = sqrt(rac2)
210  dbc = sqrt(rbc2)
211 
212 ! *** Loop over all pairs of primitive Gaussian-type functions ***
213 
214  na = 0
215  nap = 0
216 
217  DO ipgf = 1, npgfa
218 
219 ! *** Screening ***
220 
221  IF (do_dkh) THEN
222  IF (dkh_erfc) THEN
223  IF (rpgfa(ipgf) + rpgfc < dac) THEN
224  na = na + ncoset(la_max1 - maxder_local) !JT
225  nap = nap + ncoset(la_max1) !JT
226  cycle
227  END IF
228  END IF
229  ELSE
230  IF (rpgfa(ipgf) + rpgfc < dac) THEN
231  na = na + ncoset(la_max1 - maxder_local) !JT
232  nap = nap + ncoset(la_max1) !JT
233  cycle
234  END IF
235  END IF
236 
237  nb = 0
238 
239  DO jpgf = 1, npgfb
240 
241 ! *** Screening ***
242  IF (do_dkh) THEN
243  IF (dkh_erfc) THEN
244  IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
245  (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
246  nb = nb + ncoset(lb_max1) !JT
247  cycle
248  END IF
249  END IF
250  ELSE
251  IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
252  (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
253  nb = nb + ncoset(lb_max1) !JT
254  cycle
255  END IF
256  END IF
257 ! *** Calculate some prefactors ***
258 
259  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
260  zetq = 1.0_dp/zetc
261  zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
262 
263  f1 = zetb(jpgf)*zetp
264  f2 = 0.5_dp*zetp
265  f4 = -zetc*zetw
266 
267  f0 = exp(-zeta(ipgf)*f1*rab2)
268 
269  fnuc = 2.0_dp*pi*zetp*f0
270  ferf = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq*f0
271 
272  rap(:) = f1*rab(:)
273  rcp(:) = rap(:) - rac(:)
274  rpw(:) = f4*rcp(:)
275 
276 ! *** Calculate the incomplete Gamma function values ***
277 
278  rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3)
279 
280  t = rcp2/zetp
281 
282  CALL fgamma(nmax - 1, t, f)
283 
284 ! *** Calculate the basic nuclear attraction integrals [s|A(0)|s]{n} ***
285 
286  DO n = 1, nmax
287  vnuc(1, 1, n) = fnuc*f(n - 1)
288  END DO
289 
290 ! *** Calculate the incomplete Gamma function values ***
291 
292  t = -f4*rcp2/zetp
293 
294  CALL fgamma(nmax - 1, t, f)
295 
296 ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
297 
298  DO n = 1, nmax
299  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  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  DO n = 1, nmax - 1
316  vnuc(2, 1, n) = rap(1)*vnuc(1, 1, n) - rcp(1)*vnuc(1, 1, n + 1)
317  verf(2, 1, n) = rap(1)*verf(1, 1, n) + rpw(1)*verf(1, 1, n + 1)
318  vnuc(3, 1, n) = rap(2)*vnuc(1, 1, n) - rcp(2)*vnuc(1, 1, n + 1)
319  verf(3, 1, n) = rap(2)*verf(1, 1, n) + rpw(2)*verf(1, 1, n + 1)
320  vnuc(4, 1, n) = rap(3)*vnuc(1, 1, n) - rcp(3)*vnuc(1, 1, n + 1)
321  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  DO la = 2, la_max
334 
335  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  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  f4*verf(coset(0, 0, la - 2), 1, n + 1))
349 
350 ! *** Increase the angular momentum component y of function a ***
351 
352  az = la - 1
353  vnuc(coset(0, 1, az), 1, n) = &
354  rap(2)*vnuc(coset(0, 0, az), 1, n) - &
355  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  rpw(2)*verf(coset(0, 0, az), 1, n + 1)
359 
360  DO ay = 2, la
361  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  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  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  DO ay = 0, la - 1
377  az = la - 1 - ay
378  vnuc(coset(1, ay, az), 1, n) = &
379  rap(1)*vnuc(coset(0, ay, az), 1, n) - &
380  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  rpw(1)*verf(coset(0, ay, az), 1, n + 1)
384  END DO
385 
386  DO ax = 2, la
387  f3 = f2*real(ax - 1, dp)
388  DO ay = 0, la - ax
389  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  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  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  IF (lb_max > 0) THEN
411 
412 ! *** Horizontal recurrence steps ***
413 
414  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  la_start = max(0, la_min - 1)
420 
421  DO la = la_start, la_max - 1
422  DO n = 1, nmax - la - 1
423  DO ax = 0, la
424  DO ay = 0, la - ax
425  az = la - ax - ay
426  vnuc(coset(ax, ay, az), 2, n) = &
427  vnuc(coset(ax + 1, ay, az), 1, n) - &
428  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  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  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  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  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  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  DO n = 1, nmax - la_max - 1
461  DO ax = 0, la_max
462  fx = f2*real(ax, dp)
463  DO ay = 0, la_max - ax
464  fy = f2*real(ay, dp)
465  az = la_max - ax - ay
466  fz = f2*real(az, dp)
467 
468  IF (ax == 0) THEN
469  vnuc(coset(ax, ay, az), 2, n) = &
470  rbp(1)*vnuc(coset(ax, ay, az), 1, n) - &
471  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  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  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  f4*verf(coset(ax - 1, ay, az), 1, n + 1))
486  END IF
487 
488  IF (ay == 0) THEN
489  vnuc(coset(ax, ay, az), 3, n) = &
490  rbp(2)*vnuc(coset(ax, ay, az), 1, n) - &
491  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  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  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  f4*verf(coset(ax, ay - 1, az), 1, n + 1))
506  END IF
507 
508  IF (az == 0) THEN
509  vnuc(coset(ax, ay, az), 4, n) = &
510  rbp(3)*vnuc(coset(ax, ay, az), 1, n) - &
511  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  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  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  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  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  la_start = max(0, la_min - 1)
545 
546  DO la = la_start, la_max - 1
547  DO n = 1, nmax - la - lb
548  DO ax = 0, la
549  DO ay = 0, la - ax
550  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  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  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  DO by = 1, lb
564  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  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  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  DO bx = 1, lb
576  DO by = 0, lb - bx
577  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  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  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  DO n = 1, nmax - la_max - lb
608  DO ax = 0, la_max
609  fx = f2*real(ax, dp)
610  DO ay = 0, la_max - ax
611  fy = f2*real(ay, dp)
612  az = la_max - ax - ay
613  fz = f2*real(az, dp)
614 
615 ! *** Shift of angular momentum component z from a to b ***
616 
617  f3 = f2*real(lb - 1, dp)
618 
619  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  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  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  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  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  IF (ay == 0) THEN
650  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  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  rpw(2)*verf(coset(ax, ay, az), coset(0, 0, bz), n + 1)
657  DO by = 2, lb
658  bz = lb - by
659  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  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  f4*verf(coset(ax, ay, az), coset(0, by - 2, bz), n + 1))
670  END DO
671  ELSE
672  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  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  f4*verf(coset(ax, ay - 1, az), coset(0, 0, bz), n + 1))
683  DO by = 2, lb
684  bz = lb - by
685  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  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  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  IF (ax == 0) THEN
707  DO by = 0, lb - 1
708  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  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  rpw(1)*verf(coset(ax, ay, az), coset(0, by, bz), n + 1)
715  END DO
716  DO bx = 2, lb
717  f3 = f2*real(bx - 1, dp)
718  DO by = 0, lb - bx
719  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  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  f4*verf(coset(ax, ay, az), coset(bx - 2, by, bz), n + 1))
732  END DO
733  END DO
734  ELSE
735  DO by = 0, lb - 1
736  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  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  f4*verf(coset(ax - 1, ay, az), coset(0, by, bz), n + 1))
747  END DO
748  DO bx = 2, lb
749  f3 = f2*real(bx - 1, dp)
750  DO by = 0, lb - bx
751  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  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  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  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  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  DO n = 1, nmax - 1
798  vnuc(1, 2, n) = rbp(1)*vnuc(1, 1, n) - rcp(1)*vnuc(1, 1, n + 1)
799  verf(1, 2, n) = rbp(1)*verf(1, 1, n) + rpw(1)*verf(1, 1, n + 1)
800  vnuc(1, 3, n) = rbp(2)*vnuc(1, 1, n) - rcp(2)*vnuc(1, 1, n + 1)
801  verf(1, 3, n) = rbp(2)*verf(1, 1, n) + rpw(2)*verf(1, 1, n + 1)
802  vnuc(1, 4, n) = rbp(3)*vnuc(1, 1, n) - rcp(3)*vnuc(1, 1, n + 1)
803  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  DO lb = 2, lb_max
816 
817  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  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  f4*verf(1, coset(0, 0, lb - 2), n + 1))
831 
832 ! *** Increase the angular momentum component y of function b ***
833 
834  bz = lb - 1
835  vnuc(1, coset(0, 1, bz), n) = &
836  rbp(2)*vnuc(1, coset(0, 0, bz), n) - &
837  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  rpw(2)*verf(1, coset(0, 0, bz), n + 1)
841 
842  DO by = 2, lb
843  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  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  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  DO by = 0, lb - 1
859  bz = lb - 1 - by
860  vnuc(1, coset(1, by, bz), n) = &
861  rbp(1)*vnuc(1, coset(0, by, bz), n) - &
862  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  rpw(1)*verf(1, coset(0, by, bz), n + 1)
866  END DO
867 
868  DO bx = 2, lb
869  f3 = f2*real(bx - 1, dp)
870  DO by = 0, lb - bx
871  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  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  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  IF (do_dkh) THEN
898  DO j = 1, ncoset(lb_max1)
899  DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
900  vnabc(na + i, nb + j) = vnabc(na + i, nb + j) + cerf*verf(i, j, 1)
901  END DO
902  END DO
903  DO j = 1, ncoset(lb_max1)
904  DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
905  vabc(na + i, nb + j) = vabc(na + i, nb + j) - zc*vnuc(i, j, 1)
906  END DO
907  END DO
908  ELSE
909  DO j = 1, ncoset(lb_max1) !JT
910  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  zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
913  END DO
914  END DO
915  END IF
916 !JT
917 
918  IF (PRESENT(maxder)) THEN
919  DO j = 1, ncoset(lb_max1) !JT
920  DO i = 1, ncoset(la_max1) !JT
921  vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) - &
922  zc*vnuc(i, j, 1) + cerf*verf(i, j, 1)
923  END DO
924  END DO
925  END IF
926 !JTs
927  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  ftaz = 2.0_dp*zeta(ipgf)
935  ftbz = 2.0_dp*zetb(jpgf)
936  nxx = 1
937  nyy = 2
938  nzz = 3
939  nxy = 4
940  nxz = 5
941  nyx = 6
942  nyz = 7
943  nzx = 8
944  nzy = 9
945  DO la = 0, la_max1
946  DO ax = 0, la
947  fax = real(ax, dp)
948  DO ay = 0, la - ax
949  fay = real(ay, dp)
950  az = la - ax - ay
951  faz = real(az, dp)
952  pvpa = coset(ax, ay, az)
953  coamx = coset(ax - 1, ay, az)
954  coamy = coset(ax, ay - 1, az)
955  coamz = coset(ax, ay, az - 1)
956  coapx = coset(ax + 1, ay, az)
957  coapy = coset(ax, ay + 1, az)
958  coapz = coset(ax, ay, az + 1)
959  DO lb = 0, lb_max1
960  DO bx = 0, lb
961  fbx = real(bx, dp)
962  DO by = 0, lb - bx
963  fby = real(by, dp)
964  bz = lb - bx - by
965  fbz = real(bz, dp)
966  pvpb = coset(bx, by, bz)
967  cobmx = coset(bx - 1, by, bz)
968  cobmy = coset(bx, by - 1, bz)
969  cobmz = coset(bx, by, bz - 1)
970  cobpx = coset(bx + 1, by, bz)
971  cobpy = coset(bx, by + 1, bz)
972  cobpz = coset(bx, by, bz + 1)
973  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  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  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  DO j = 1, ncoset(lb_max1)
1008  DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1 - maxder_local)
1009  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  nb = nb + ncoset(lb_max1) !JT
1015 
1016  END DO
1017 
1018  na = na + ncoset(la_max1 - maxder_local) !JT
1019  nap = nap + ncoset(la_max1) !JT
1020 
1021  END DO
1022 
1023  END SUBROUTINE verfc
1024 
1025 END MODULE ai_verfc
ai_verfc
Build up the nuclear potential part of the core Hamiltonian matrix in the case of an allelectron calc...
Definition: ai_verfc.F:82
ai_verfc::verfc
subroutine, public verfc(la_max1, npgfa, zeta, rpgfa, la_min1, lb_max1, npgfb, zetb, rpgfb, lb_min1, zetc, rpgfc, zc, cerf, rab, rab2, rac, rac2, rbc2, vabc, verf, vnuc, f, maxder, vabc_plus, vnabc, pVp_sum, pVp, dkh_erfc)
Calculation of the primitive three-center nuclear potential integrals <a|Z*erfc(r)/r|b> over Cartesia...
Definition: ai_verfc.F:141
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