(git:374b731)
Loading...
Searching...
No Matches
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
82MODULE 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
99CONTAINS
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
1025END 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