Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : !!****** cp2k/ai_overlap3 [1.0] *
8 : !!
9 : !! NAME
10 : !! ai_overlap3
11 : !!
12 : !! FUNCTION
13 : !! Calculation of three-center overlap integrals over Cartesian
14 : !! Gaussian-type functions.
15 : !!
16 : !! AUTHOR
17 : !! Matthias Krack (26.06.2001)
18 : !!
19 : !! LITERATURE
20 : !! S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
21 : !!
22 : !******************************************************************************
23 :
24 : MODULE ai_overlap3
25 :
26 : ! **************************************************************************************************
27 :
28 : ! ax,ay,az : Angular momentum index numbers of orbital a.
29 : ! bx,by,bz : Angular momentum index numbers of orbital b.
30 : ! coset : Cartesian orbital set pointer.
31 : ! dab : Distance between the atomic centers a and b.
32 : ! dac : Distance between the atomic centers a and c.
33 : ! dbc : Distance between the atomic centers b and c.
34 : ! l{a,b,c} : Angular momentum quantum number of shell a, b or c.
35 : ! l{a,b}_max : Maximum angular momentum quantum number of shell a, b or c.
36 : ! ncoset : Number of Cartesian orbitals up to l.
37 : ! rab : Distance vector between the atomic centers a and b.
38 : ! rac : Distance vector between the atomic centers a and c.
39 : ! rbc : Distance vector between the atomic centers b and c.
40 : ! rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b.
41 : ! zet{a,b,c} : Exponents of the Gaussian-type functions a or b.
42 : ! zetg : Reciprocal of the sum of the exponents of orbital a, b and c.
43 : ! zetp : Reciprocal of the sum of the exponents of orbital a and b.
44 :
45 : ! **************************************************************************************************
46 :
47 : USE kinds, ONLY: dp
48 : USE mathconstants, ONLY: pi
49 : USE orbital_pointers, ONLY: coset,&
50 : ncoset
51 : #include "../base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3'
58 :
59 : ! *** Public subroutines ***
60 :
61 : PUBLIC :: overlap3
62 :
63 : !!***
64 : ! **************************************************************************************************
65 :
66 : CONTAINS
67 :
68 : ! ***************************************************************************************************
69 : !> \brief Calculation of three-center overlap integrals [a|b|c] over primitive
70 : !> Cartesian Gaussian functions
71 : !> \param la_max_set ...
72 : !> \param npgfa ...
73 : !> \param zeta ...
74 : !> \param rpgfa ...
75 : !> \param la_min_set ...
76 : !> \param lb_max_set ...
77 : !> \param npgfb ...
78 : !> \param zetb ...
79 : !> \param rpgfb ...
80 : !> \param lb_min_set ...
81 : !> \param lc_max_set ...
82 : !> \param npgfc ...
83 : !> \param zetc ...
84 : !> \param rpgfc ...
85 : !> \param lc_min_set ...
86 : !> \param rab ...
87 : !> \param dab ...
88 : !> \param rac ...
89 : !> \param dac ...
90 : !> \param rbc ...
91 : !> \param dbc ...
92 : !> \param sabc integrals [a|b|c]
93 : !> \param sdabc derivative [da/dAi|b|c]
94 : !> \param sabdc derivative [a|b|dc/dCi]
95 : !> \param int_abc_ext the extremal value of sabc, i.e., MAXVAL(ABS(sabc))
96 : !> \par History
97 : !> 05.2014 created (Dorothea Golze)
98 : !> \author Dorothea Golze
99 : !> \note overlap3 essentially uses the setup of overlap3_old
100 : ! **************************************************************************************************
101 :
102 24 : SUBROUTINE overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
103 48 : lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
104 48 : lc_max_set, npgfc, zetc, rpgfc, lc_min_set, &
105 24 : rab, dab, rac, dac, rbc, dbc, sabc, &
106 24 : sdabc, sabdc, int_abc_ext)
107 :
108 : INTEGER, INTENT(IN) :: la_max_set, npgfa
109 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
110 : INTEGER, INTENT(IN) :: la_min_set, lb_max_set, npgfb
111 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
112 : INTEGER, INTENT(IN) :: lb_min_set, lc_max_set, npgfc
113 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
114 : INTEGER, INTENT(IN) :: lc_min_set
115 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
116 : REAL(KIND=dp), INTENT(IN) :: dab
117 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
118 : REAL(KIND=dp), INTENT(IN) :: dac
119 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
120 : REAL(KIND=dp), INTENT(IN) :: dbc
121 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabc
122 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
123 : INTENT(INOUT), OPTIONAL :: sdabc, sabdc
124 : REAL(dp), INTENT(OUT), OPTIONAL :: int_abc_ext
125 :
126 : CHARACTER(len=*), PARAMETER :: routineN = 'overlap3'
127 :
128 : INTEGER :: ax, ay, az, bx, by, bz, coa, coax, coay, coaz, coc, cocx, cocy, cocz, cx, cy, cz, &
129 : handle, i, ipgf, j, jpgf, k, kpgf, l, la, la_max, la_min, la_start, lai, lb, lb_max, &
130 : lb_min, lc, lc_max, lc_min, lci, na, nb, nc, nda, ndc
131 : REAL(KIND=dp) :: f0, f1, f2, f3, fcx, fcy, fcz, fx, fy, &
132 : fz, rcp2, zetg, zetp
133 : REAL(KIND=dp), DIMENSION(3) :: rag, rbg, rcg, rcp
134 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: s
135 24 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc
136 :
137 : ! ---------------------------------------------------------------------------
138 :
139 24 : CALL timeset(routineN, handle)
140 :
141 24 : NULLIFY (s, sda, sdc)
142 :
143 24 : lai = 0
144 24 : lci = 0
145 :
146 24 : IF (PRESENT(sdabc)) lai = 1
147 24 : IF (PRESENT(sabdc)) lci = 1
148 :
149 24 : la_max = la_max_set + lai
150 24 : la_min = MAX(0, la_min_set - lai)
151 24 : lb_max = lb_max_set
152 24 : lb_min = lb_min_set
153 24 : lc_max = lc_max_set + lci
154 24 : lc_min = MAX(0, lc_min_set - lci)
155 :
156 120 : ALLOCATE (s(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
157 7890 : s = 0._dp
158 24 : IF (PRESENT(sdabc)) THEN
159 72 : ALLOCATE (sda(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
160 19488 : sda = 0._dp
161 : END IF
162 24 : IF (PRESENT(sabdc)) THEN
163 72 : ALLOCATE (sdc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
164 19488 : sdc = 0._dp
165 : END IF
166 24 : IF (PRESENT(int_abc_ext)) THEN
167 0 : int_abc_ext = 0.0_dp
168 : END IF
169 :
170 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
171 :
172 24 : na = 0
173 24 : nda = 0
174 144 : DO ipgf = 1, npgfa
175 :
176 120 : nb = 0
177 720 : DO jpgf = 1, npgfb
178 :
179 : ! *** Screening ***
180 600 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
181 : sabc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
182 0 : 1:npgfc*ncoset(lc_max_set)) = 0.0_dp
183 0 : IF (PRESENT(sdabc)) THEN
184 : sdabc(nda + 1:nda + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
185 0 : 1:npgfc*ncoset(lc_max_set), 1:3) = 0.0_dp
186 : END IF
187 0 : IF (PRESENT(sabdc)) THEN
188 : sabdc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
189 0 : 1:npgfc*ncoset(lc_max_set), 1:3) = 0.0_dp
190 : END IF
191 0 : nb = nb + ncoset(lb_max_set)
192 0 : CYCLE
193 : END IF
194 :
195 600 : nc = 0
196 600 : ndc = 0
197 2700 : DO kpgf = 1, npgfc
198 :
199 : ! *** Screening ***
200 2100 : IF ((rpgfb(jpgf) + rpgfc(kpgf) < dbc) .OR. &
201 : (rpgfa(ipgf) + rpgfc(kpgf) < dac)) THEN
202 : sabc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
203 0 : nc + 1:nc + ncoset(lc_max_set)) = 0.0_dp
204 0 : IF (PRESENT(sdabc)) THEN
205 : sdabc(nda + 1:nda + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
206 0 : nc + 1:nc + ncoset(lc_max_set), 1:3) = 0.0_dp
207 : END IF
208 0 : IF (PRESENT(sabdc)) THEN
209 : sabdc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
210 0 : ndc + 1:ndc + ncoset(lc_max_set), 1:3) = 0.0_dp
211 : END IF
212 0 : nc = nc + ncoset(lc_max_set)
213 0 : ndc = ndc + ncoset(lc_max_set)
214 0 : CYCLE
215 : END IF
216 :
217 : ! *** Calculate some prefactors ***
218 2100 : zetg = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc(kpgf))
219 2100 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
220 2100 : f0 = (pi*zetg)**1.5_dp
221 2100 : f1 = zetb(jpgf)*zetp
222 2100 : f2 = 0.5_dp*zetg
223 8400 : rcp(:) = f1*rab(:) - rac(:)
224 2100 : rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3)
225 :
226 : ! *** Calculate the basic three-center overlap integral [s|s|s] ***
227 2100 : s(1, 1, 1) = f0*EXP(-(zeta(ipgf)*f1*dab*dab + zetc(kpgf)*zetg*rcp2/zetp))
228 :
229 : ! *** Recurrence steps: [s|s|s] -> [a|s|s] ***
230 :
231 2100 : IF (la_max > 0) THEN
232 :
233 : ! *** Vertical recurrence steps: [s|s|s] -> [a|s|s] ***
234 :
235 8400 : rag(:) = zetg*(zetb(jpgf)*rab(:) + zetc(kpgf)*rac(:))
236 :
237 : ! *** [p|s|s] = (Gi - Ai)*[s|s|s] (i = x,y,z) ***
238 :
239 2100 : s(2, 1, 1) = rag(1)*s(1, 1, 1)
240 2100 : s(3, 1, 1) = rag(2)*s(1, 1, 1)
241 2100 : s(4, 1, 1) = rag(3)*s(1, 1, 1)
242 :
243 : ! *** [a|s|s] = (Gi - Ai)*[a-1i|s|s] + f2*Ni(a-1i)*[a-2i|s|s] ***
244 :
245 3150 : DO la = 2, la_max
246 :
247 : ! *** Increase the angular momentum component z of function a ***
248 :
249 : s(coset(0, 0, la), 1, 1) = rag(3)*s(coset(0, 0, la - 1), 1, 1) + &
250 1050 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, 1)
251 :
252 : ! *** Increase the angular momentum component y of function a ***
253 :
254 1050 : az = la - 1
255 1050 : s(coset(0, 1, az), 1, 1) = rag(2)*s(coset(0, 0, az), 1, 1)
256 :
257 2100 : DO ay = 2, la
258 1050 : az = la - ay
259 : s(coset(0, ay, az), 1, 1) = rag(2)*s(coset(0, ay - 1, az), 1, 1) + &
260 2100 : f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, 1)
261 : END DO
262 :
263 : ! *** Increase the angular momentum component x of function a ***
264 :
265 3150 : DO ay = 0, la - 1
266 2100 : az = la - 1 - ay
267 3150 : s(coset(1, ay, az), 1, 1) = rag(1)*s(coset(0, ay, az), 1, 1)
268 : END DO
269 :
270 4200 : DO ax = 2, la
271 1050 : f3 = f2*REAL(ax - 1, dp)
272 3150 : DO ay = 0, la - ax
273 1050 : az = la - ax - ay
274 : s(coset(ax, ay, az), 1, 1) = rag(1)*s(coset(ax - 1, ay, az), 1, 1) + &
275 2100 : f3*s(coset(ax - 2, ay, az), 1, 1)
276 : END DO
277 : END DO
278 :
279 : END DO
280 :
281 : ! *** Recurrence steps: [a|s|s] -> [a|s|b] ***
282 :
283 2100 : IF (lb_max > 0) THEN
284 :
285 : ! *** Horizontal recurrence steps ***
286 :
287 8400 : rbg(:) = rag(:) - rab(:)
288 :
289 : ! *** [a|s|p] = [a+1i|s|s] - (Bi - Ai)*[a|s|s] ***
290 :
291 2100 : IF (lb_max == 1) THEN
292 : la_start = la_min
293 : ELSE
294 0 : la_start = MAX(0, la_min - 1)
295 : END IF
296 :
297 5250 : DO la = la_start, la_max - 1
298 9450 : DO ax = 0, la
299 12600 : DO ay = 0, la - ax
300 5250 : az = la - ax - ay
301 5250 : coa = coset(ax, ay, az)
302 5250 : coax = coset(ax + 1, ay, az)
303 5250 : coay = coset(ax, ay + 1, az)
304 5250 : coaz = coset(ax, ay, az + 1)
305 5250 : s(coset(ax, ay, az), 2, 1) = s(coax, 1, 1) - rab(1)*s(coa, 1, 1)
306 5250 : s(coset(ax, ay, az), 3, 1) = s(coay, 1, 1) - rab(2)*s(coa, 1, 1)
307 9450 : s(coset(ax, ay, az), 4, 1) = s(coaz, 1, 1) - rab(3)*s(coa, 1, 1)
308 : END DO
309 : END DO
310 : END DO
311 :
312 : ! *** Vertical recurrence step ***
313 :
314 : ! *** [a|s|p] = (Gi - Bi)*[a|s|s] + f2*Ni(a)*[a-1i|s|s] ***
315 :
316 7350 : DO ax = 0, la_max
317 5250 : fx = f2*REAL(ax, dp)
318 16800 : DO ay = 0, la_max - ax
319 9450 : fy = f2*REAL(ay, dp)
320 9450 : az = la_max - ax - ay
321 9450 : fz = f2*REAL(az, dp)
322 9450 : coa = coset(ax, ay, az)
323 9450 : IF (ax == 0) THEN
324 5250 : s(coa, 2, 1) = rbg(1)*s(coa, 1, 1)
325 : ELSE
326 4200 : s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) + fx*s(coset(ax - 1, ay, az), 1, 1)
327 : END IF
328 9450 : IF (ay == 0) THEN
329 5250 : s(coa, 3, 1) = rbg(2)*s(coa, 1, 1)
330 : ELSE
331 4200 : s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) + fy*s(coset(ax, ay - 1, az), 1, 1)
332 : END IF
333 14700 : IF (az == 0) THEN
334 5250 : s(coa, 4, 1) = rbg(3)*s(coa, 1, 1)
335 : ELSE
336 4200 : s(coa, 4, 1) = rbg(3)*s(coa, 1, 1) + fz*s(coset(ax, ay, az - 1), 1, 1)
337 : END IF
338 : END DO
339 : END DO
340 :
341 : ! *** Recurrence steps: [a|s|p] -> [a|s|b] ***
342 :
343 2100 : DO lb = 2, lb_max
344 :
345 : ! *** Horizontal recurrence steps ***
346 :
347 : ! *** [a|s|b] = [a+1i|s|b-1i] - (Bi - Ai)*[a|s|b-1i] ***
348 :
349 0 : IF (lb == lb_max) THEN
350 : la_start = la_min
351 : ELSE
352 0 : la_start = MAX(0, la_min - 1)
353 : END IF
354 :
355 0 : DO la = la_start, la_max - 1
356 0 : DO ax = 0, la
357 0 : DO ay = 0, la - ax
358 0 : az = la - ax - ay
359 :
360 0 : coa = coset(ax, ay, az)
361 0 : coax = coset(ax + 1, ay, az)
362 0 : coay = coset(ax, ay + 1, az)
363 0 : coaz = coset(ax, ay, az + 1)
364 :
365 : ! *** Shift of angular momentum component z from a to b ***
366 :
367 : s(coa, coset(0, 0, lb), 1) = &
368 : s(coaz, coset(0, 0, lb - 1), 1) - &
369 0 : rab(3)*s(coa, coset(0, 0, lb - 1), 1)
370 :
371 : ! *** Shift of angular momentum component y from a to b ***
372 :
373 0 : DO by = 1, lb
374 0 : bz = lb - by
375 : s(coa, coset(0, by, bz), 1) = &
376 : s(coay, coset(0, by - 1, bz), 1) - &
377 0 : rab(2)*s(coa, coset(0, by - 1, bz), 1)
378 : END DO
379 :
380 : ! *** Shift of angular momentum component x from a to b ***
381 :
382 0 : DO bx = 1, lb
383 0 : DO by = 0, lb - bx
384 0 : bz = lb - bx - by
385 : s(coa, coset(bx, by, bz), 1) = &
386 : s(coax, coset(bx - 1, by, bz), 1) - &
387 0 : rab(1)*s(coa, coset(bx - 1, by, bz), 1)
388 : END DO
389 : END DO
390 :
391 : END DO
392 : END DO
393 : END DO
394 :
395 : ! *** Vertical recurrence step ***
396 :
397 : ! *** [a|s|b] = (Gi - Bi)*[a|s|b-1i] + ***
398 : ! *** f2*Ni(a)*[a-1i|s|b-1i] + ***
399 : ! *** f2*Ni(b-1i)*[a|s|b-2i] ***
400 :
401 2100 : DO ax = 0, la_max
402 0 : fx = f2*REAL(ax, dp)
403 0 : DO ay = 0, la_max - ax
404 0 : fy = f2*REAL(ay, dp)
405 0 : az = la_max - ax - ay
406 0 : fz = f2*REAL(az, dp)
407 :
408 0 : coa = coset(ax, ay, az)
409 :
410 0 : f3 = f2*REAL(lb - 1, dp)
411 :
412 : ! *** Shift of angular momentum component z from a to b ***
413 :
414 0 : IF (az == 0) THEN
415 : s(coa, coset(0, 0, lb), 1) = &
416 : rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
417 0 : f3*s(coa, coset(0, 0, lb - 2), 1)
418 : ELSE
419 0 : coaz = coset(ax, ay, az - 1)
420 : s(coa, coset(0, 0, lb), 1) = &
421 : rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
422 : fz*s(coaz, coset(0, 0, lb - 1), 1) + &
423 0 : f3*s(coa, coset(0, 0, lb - 2), 1)
424 : END IF
425 :
426 : ! *** Shift of angular momentum component y from a to b ***
427 :
428 0 : IF (ay == 0) THEN
429 0 : bz = lb - 1
430 : s(coa, coset(0, 1, bz), 1) = &
431 0 : rbg(2)*s(coa, coset(0, 0, bz), 1)
432 0 : DO by = 2, lb
433 0 : bz = lb - by
434 0 : f3 = f2*REAL(by - 1, dp)
435 : s(coa, coset(0, by, bz), 1) = &
436 : rbg(2)*s(coa, coset(0, by - 1, bz), 1) + &
437 0 : f3*s(coa, coset(0, by - 2, bz), 1)
438 : END DO
439 : ELSE
440 0 : coay = coset(ax, ay - 1, az)
441 0 : bz = lb - 1
442 : s(coa, coset(0, 1, bz), 1) = &
443 : rbg(2)*s(coa, coset(0, 0, bz), 1) + &
444 0 : fy*s(coay, coset(0, 0, bz), 1)
445 0 : DO by = 2, lb
446 0 : bz = lb - by
447 0 : f3 = f2*REAL(by - 1, dp)
448 : s(coa, coset(0, by, bz), 1) = &
449 : rbg(2)*s(coa, coset(0, by - 1, bz), 1) + &
450 : fy*s(coay, coset(0, by - 1, bz), 1) + &
451 0 : f3*s(coa, coset(0, by - 2, bz), 1)
452 : END DO
453 : END IF
454 :
455 : ! *** Shift of angular momentum component x from a to b ***
456 :
457 0 : IF (ax == 0) THEN
458 0 : DO by = 0, lb - 1
459 0 : bz = lb - 1 - by
460 : s(coa, coset(1, by, bz), 1) = &
461 0 : rbg(1)*s(coa, coset(0, by, bz), 1)
462 : END DO
463 0 : DO bx = 2, lb
464 0 : f3 = f2*REAL(bx - 1, dp)
465 0 : DO by = 0, lb - bx
466 0 : bz = lb - bx - by
467 : s(coa, coset(bx, by, bz), 1) = &
468 : rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
469 0 : f3*s(coa, coset(bx - 2, by, bz), 1)
470 : END DO
471 : END DO
472 : ELSE
473 0 : coax = coset(ax - 1, ay, az)
474 0 : DO by = 0, lb - 1
475 0 : bz = lb - 1 - by
476 : s(coa, coset(1, by, bz), 1) = &
477 : rbg(1)*s(coa, coset(0, by, bz), 1) + &
478 0 : fx*s(coax, coset(0, by, bz), 1)
479 : END DO
480 0 : DO bx = 2, lb
481 0 : f3 = f2*REAL(bx - 1, dp)
482 0 : DO by = 0, lb - bx
483 0 : bz = lb - bx - by
484 : s(coa, coset(bx, by, bz), 1) = &
485 : rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
486 : fx*s(coax, coset(bx - 1, by, bz), 1) + &
487 0 : f3*s(coa, coset(bx - 2, by, bz), 1)
488 : END DO
489 : END DO
490 : END IF
491 :
492 : END DO
493 : END DO
494 :
495 : END DO
496 :
497 : END IF
498 :
499 : ELSE
500 :
501 0 : IF (lb_max > 0) THEN
502 :
503 : ! *** Vertical recurrence steps: [s|s|s] -> [s|s|b] ***
504 :
505 0 : rbg(:) = -zetg*(zeta(ipgf)*rab(:) - zetc(kpgf)*rbc(:))
506 :
507 : ! *** [s|s|p] = (Gi - Bi)*[s|s|s] ***
508 :
509 0 : s(1, 2, 1) = rbg(1)*s(1, 1, 1)
510 0 : s(1, 3, 1) = rbg(2)*s(1, 1, 1)
511 0 : s(1, 4, 1) = rbg(3)*s(1, 1, 1)
512 :
513 : ! *** [s|s|b] = (Gi - Bi)*[s|s|b-1i] + f2*Ni(b-1i)*[s|s|b-2i] ***
514 :
515 0 : DO lb = 2, lb_max
516 :
517 : ! *** Increase the angular momentum component z of function b ***
518 :
519 : s(1, coset(0, 0, lb), 1) = rbg(3)*s(1, coset(0, 0, lb - 1), 1) + &
520 0 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), 1)
521 :
522 : ! *** Increase the angular momentum component y of function b ***
523 :
524 0 : bz = lb - 1
525 0 : s(1, coset(0, 1, bz), 1) = rbg(2)*s(1, coset(0, 0, bz), 1)
526 :
527 0 : DO by = 2, lb
528 0 : bz = lb - by
529 : s(1, coset(0, by, bz), 1) = &
530 : rbg(2)*s(1, coset(0, by - 1, bz), 1) + &
531 0 : f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), 1)
532 : END DO
533 :
534 : ! *** Increase the angular momentum component x of function b ***
535 :
536 0 : DO by = 0, lb - 1
537 0 : bz = lb - 1 - by
538 0 : s(1, coset(1, by, bz), 1) = rbg(1)*s(1, coset(0, by, bz), 1)
539 : END DO
540 :
541 0 : DO bx = 2, lb
542 0 : f3 = f2*REAL(bx - 1, dp)
543 0 : DO by = 0, lb - bx
544 0 : bz = lb - bx - by
545 : s(1, coset(bx, by, bz), 1) = rbg(1)*s(1, coset(bx - 1, by, bz), 1) + &
546 0 : f3*s(1, coset(bx - 2, by, bz), 1)
547 : END DO
548 : END DO
549 :
550 : END DO
551 :
552 : END IF
553 :
554 : END IF
555 :
556 : ! *** Recurrence steps: [a|s|b] -> [a|c|b] ***
557 :
558 2100 : IF (lc_max > 0) THEN
559 :
560 : ! *** Vertical recurrence steps: [s|s|s] -> [s|c|s] ***
561 :
562 4800 : rcg(:) = -zetg*(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))
563 :
564 : ! *** [s|p|s] = (Gi - Ci)*[s|s|s] (i = x,y,z) ***
565 :
566 1200 : s(1, 1, 2) = rcg(1)*s(1, 1, 1)
567 1200 : s(1, 1, 3) = rcg(2)*s(1, 1, 1)
568 1200 : s(1, 1, 4) = rcg(3)*s(1, 1, 1)
569 :
570 : ! *** [s|c|s] = (Gi - Ci)*[s|c-1i|s] + f2*Ni(c-1i)*[s|c-2i|s] ***
571 :
572 1650 : DO lc = 2, lc_max
573 :
574 : ! *** Increase the angular momentum component z of function c ***
575 :
576 : s(1, 1, coset(0, 0, lc)) = rcg(3)*s(1, 1, coset(0, 0, lc - 1)) + &
577 450 : f2*REAL(lc - 1, dp)*s(1, 1, coset(0, 0, lc - 2))
578 :
579 : ! *** Increase the angular momentum component y of function c ***
580 :
581 450 : cz = lc - 1
582 450 : s(1, 1, coset(0, 1, cz)) = rcg(2)*s(1, 1, coset(0, 0, cz))
583 :
584 1050 : DO cy = 2, lc
585 600 : cz = lc - cy
586 : s(1, 1, coset(0, cy, cz)) = rcg(2)*s(1, 1, coset(0, cy - 1, cz)) + &
587 1050 : f2*REAL(cy - 1, dp)*s(1, 1, coset(0, cy - 2, cz))
588 : END DO
589 :
590 : ! *** Increase the angular momentum component x of function c ***
591 :
592 1500 : DO cy = 0, lc - 1
593 1050 : cz = lc - 1 - cy
594 1500 : s(1, 1, coset(1, cy, cz)) = rcg(1)*s(1, 1, coset(0, cy, cz))
595 : END DO
596 :
597 2250 : DO cx = 2, lc
598 600 : f3 = f2*REAL(cx - 1, dp)
599 1800 : DO cy = 0, lc - cx
600 750 : cz = lc - cx - cy
601 : s(1, 1, coset(cx, cy, cz)) = rcg(1)*s(1, 1, coset(cx - 1, cy, cz)) + &
602 1350 : f3*s(1, 1, coset(cx - 2, cy, cz))
603 : END DO
604 : END DO
605 :
606 : END DO
607 :
608 : ! *** Recurrence steps: [s|c|s] -> [a|c|b] ***
609 :
610 2850 : DO lc = 1, lc_max
611 :
612 6750 : DO cx = 0, lc
613 12450 : DO cy = 0, lc - cx
614 6900 : cz = lc - cx - cy
615 :
616 6900 : coc = coset(cx, cy, cz)
617 6900 : cocx = coset(MAX(0, cx - 1), cy, cz)
618 6900 : cocy = coset(cx, MAX(0, cy - 1), cz)
619 6900 : cocz = coset(cx, cy, MAX(0, cz - 1))
620 :
621 6900 : fcx = f2*REAL(cx, dp)
622 6900 : fcy = f2*REAL(cy, dp)
623 6900 : fcz = f2*REAL(cz, dp)
624 :
625 : ! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
626 :
627 10800 : IF (la_max > 0) THEN
628 :
629 : ! *** Vertical recurrence steps: [s|c|s] -> [a|c|s] ***
630 :
631 27600 : rag(:) = rcg(:) + rac(:)
632 :
633 : ! *** [p|c|s] = (Gi - Ai)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
634 :
635 6900 : s(2, 1, coc) = rag(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
636 6900 : s(3, 1, coc) = rag(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
637 6900 : s(4, 1, coc) = rag(3)*s(1, 1, coc) + fcz*s(1, 1, cocz)
638 :
639 : ! *** [a|c|s] = (Gi - Ai)*[a-1i|c|s] + ***
640 : ! *** f2*Ni(a-1i)*[a-2i|c|s] + ***
641 : ! *** f2*Ni(c)*[a-1i|c-1i|s] ***
642 :
643 12450 : DO la = 2, la_max
644 :
645 : ! *** Increase the angular momentum component z of a ***
646 :
647 : s(coset(0, 0, la), 1, coc) = &
648 : rag(3)*s(coset(0, 0, la - 1), 1, coc) + &
649 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, coc) + &
650 5550 : fcz*s(coset(0, 0, la - 1), 1, cocz)
651 :
652 : ! *** Increase the angular momentum component y of a ***
653 :
654 5550 : az = la - 1
655 : s(coset(0, 1, az), 1, coc) = &
656 : rag(2)*s(coset(0, 0, az), 1, coc) + &
657 5550 : fcy*s(coset(0, 0, az), 1, cocy)
658 :
659 11100 : DO ay = 2, la
660 5550 : az = la - ay
661 : s(coset(0, ay, az), 1, coc) = &
662 : rag(2)*s(coset(0, ay - 1, az), 1, coc) + &
663 : f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, coc) + &
664 11100 : fcy*s(coset(0, ay - 1, az), 1, cocy)
665 : END DO
666 :
667 : ! *** Increase the angular momentum component x of a ***
668 :
669 16650 : DO ay = 0, la - 1
670 11100 : az = la - 1 - ay
671 : s(coset(1, ay, az), 1, coc) = &
672 : rag(1)*s(coset(0, ay, az), 1, coc) + &
673 16650 : fcx*s(coset(0, ay, az), 1, cocx)
674 : END DO
675 :
676 18000 : DO ax = 2, la
677 5550 : f3 = f2*REAL(ax - 1, dp)
678 16650 : DO ay = 0, la - ax
679 5550 : az = la - ax - ay
680 : s(coset(ax, ay, az), 1, coc) = &
681 : rag(1)*s(coset(ax - 1, ay, az), 1, coc) + &
682 : f3*s(coset(ax - 2, ay, az), 1, coc) + &
683 11100 : fcx*s(coset(ax - 1, ay, az), 1, cocx)
684 : END DO
685 : END DO
686 :
687 : END DO
688 :
689 : ! *** Recurrence steps: [a|c|s] -> [a|c|b] ***
690 :
691 6900 : IF (lb_max > 0) THEN
692 :
693 : ! *** Horizontal recurrence steps ***
694 :
695 27600 : rbg(:) = rag(:) - rab(:)
696 :
697 : ! *** [a|c|p] = [a+1i|c|s] - (Bi - Ai)*[a|c|s] ***
698 :
699 6900 : IF (lb_max == 1) THEN
700 : la_start = la_min
701 : ELSE
702 0 : la_start = MAX(0, la_min - 1)
703 : END IF
704 :
705 19350 : DO la = la_start, la_max - 1
706 37350 : DO ax = 0, la
707 54000 : DO ay = 0, la - ax
708 23550 : az = la - ax - ay
709 23550 : coa = coset(ax, ay, az)
710 23550 : coax = coset(ax + 1, ay, az)
711 23550 : coay = coset(ax, ay + 1, az)
712 23550 : coaz = coset(ax, ay, az + 1)
713 23550 : s(coa, 2, coc) = s(coax, 1, coc) - rab(1)*s(coa, 1, coc)
714 23550 : s(coa, 3, coc) = s(coay, 1, coc) - rab(2)*s(coa, 1, coc)
715 41550 : s(coa, 4, coc) = s(coaz, 1, coc) - rab(3)*s(coa, 1, coc)
716 : END DO
717 : END DO
718 : END DO
719 :
720 : ! *** Vertical recurrence step ***
721 :
722 : ! *** [a|c|p] = (Gi - Bi)*[a|c|s] + ***
723 : ! f2*Ni(a)*[a-1i|c|s] + ***
724 : ! f2*Ni(c)*[a|c-1i|s] ***
725 :
726 26250 : DO ax = 0, la_max
727 19350 : fx = f2*REAL(ax, dp)
728 63600 : DO ay = 0, la_max - ax
729 37350 : fy = f2*REAL(ay, dp)
730 37350 : az = la_max - ax - ay
731 37350 : fz = f2*REAL(az, dp)
732 37350 : coa = coset(ax, ay, az)
733 37350 : IF (ax == 0) THEN
734 : s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
735 19350 : fcx*s(coa, 1, cocx)
736 : ELSE
737 : s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
738 : fx*s(coset(ax - 1, ay, az), 1, coc) + &
739 18000 : fcx*s(coa, 1, cocx)
740 : END IF
741 37350 : IF (ay == 0) THEN
742 : s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
743 19350 : fcy*s(coa, 1, cocy)
744 : ELSE
745 : s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
746 : fy*s(coset(ax, ay - 1, az), 1, coc) + &
747 18000 : fcy*s(coa, 1, cocy)
748 : END IF
749 56700 : IF (az == 0) THEN
750 : s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
751 19350 : fcz*s(coa, 1, cocz)
752 : ELSE
753 : s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
754 : fz*s(coset(ax, ay, az - 1), 1, coc) + &
755 18000 : fcz*s(coa, 1, cocz)
756 : END IF
757 : END DO
758 : END DO
759 :
760 : ! *** Recurrence steps: [a|c|p] -> [a|c|b] ***
761 :
762 6900 : DO lb = 2, lb_max
763 :
764 : ! *** Horizontal recurrence steps ***
765 :
766 : ! *** [a|c|b] = [a+1i|c|b-1i] - (Bi - Ai)*[a|c|b-1i] ***
767 :
768 0 : IF (lb == lb_max) THEN
769 : la_start = la_min
770 : ELSE
771 0 : la_start = MAX(0, la_min - 1)
772 : END IF
773 :
774 0 : DO la = la_start, la_max - 1
775 0 : DO ax = 0, la
776 0 : DO ay = 0, la - ax
777 0 : az = la - ax - ay
778 :
779 0 : coa = coset(ax, ay, az)
780 0 : coax = coset(ax + 1, ay, az)
781 0 : coay = coset(ax, ay + 1, az)
782 0 : coaz = coset(ax, ay, az + 1)
783 :
784 : ! *** Shift of angular momentum ***
785 : ! *** component z from a to b ***
786 :
787 : s(coa, coset(0, 0, lb), coc) = &
788 : s(coaz, coset(0, 0, lb - 1), coc) - &
789 0 : rab(3)*s(coa, coset(0, 0, lb - 1), coc)
790 :
791 : ! *** Shift of angular momentum ***
792 : ! *** component y from a to b ***
793 :
794 0 : DO by = 1, lb
795 0 : bz = lb - by
796 : s(coa, coset(0, by, bz), coc) = &
797 : s(coay, coset(0, by - 1, bz), coc) - &
798 0 : rab(2)*s(coa, coset(0, by - 1, bz), coc)
799 : END DO
800 :
801 : ! *** Shift of angular momentum ***
802 : ! *** component x from a to b ***
803 :
804 0 : DO bx = 1, lb
805 0 : DO by = 0, lb - bx
806 0 : bz = lb - bx - by
807 : s(coa, coset(bx, by, bz), coc) = &
808 : s(coax, coset(bx - 1, by, bz), coc) - &
809 0 : rab(1)*s(coa, coset(bx - 1, by, bz), coc)
810 : END DO
811 : END DO
812 :
813 : END DO
814 : END DO
815 : END DO
816 :
817 : ! *** Vertical recurrence step ***
818 :
819 : ! *** [a|c|b] = (Gi - Bi)*[a|c|b-1i] + ***
820 : ! *** f2*Ni(a)*[a-1i|c|b-1i] + ***
821 : ! *** f2*Ni(b-1i)*[a|c|b-2i] + ***
822 : ! *** f2*Ni(c)*[a|c-1i|b-1i] ***
823 :
824 6900 : DO ax = 0, la_max
825 0 : fx = f2*REAL(ax, dp)
826 0 : DO ay = 0, la_max - ax
827 0 : fy = f2*REAL(ay, dp)
828 0 : az = la_max - ax - ay
829 0 : fz = f2*REAL(az, dp)
830 :
831 0 : coa = coset(ax, ay, az)
832 0 : coax = coset(MAX(0, ax - 1), ay, az)
833 0 : coay = coset(ax, MAX(0, ay - 1), az)
834 0 : coaz = coset(ax, ay, MAX(0, az - 1))
835 :
836 0 : f3 = f2*REAL(lb - 1, dp)
837 :
838 : ! *** Shift of angular momentum ***
839 : ! *** component z from a to b ***
840 :
841 0 : IF (az == 0) THEN
842 : s(coa, coset(0, 0, lb), coc) = &
843 : rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + &
844 : f3*s(coa, coset(0, 0, lb - 2), coc) + &
845 0 : fcz*s(coa, coset(0, 0, lb - 1), cocz)
846 : ELSE
847 : s(coa, coset(0, 0, lb), coc) = &
848 : rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + &
849 : fz*s(coaz, coset(0, 0, lb - 1), coc) + &
850 : f3*s(coa, coset(0, 0, lb - 2), coc) + &
851 0 : fcz*s(coa, coset(0, 0, lb - 1), cocz)
852 : END IF
853 :
854 : ! *** Shift of angular momentum ***
855 : ! *** component y from a to b ***
856 :
857 0 : IF (ay == 0) THEN
858 0 : bz = lb - 1
859 : s(coa, coset(0, 1, bz), coc) = &
860 : rbg(2)*s(coa, coset(0, 0, bz), coc) + &
861 0 : fcy*s(coa, coset(0, 0, bz), cocy)
862 0 : DO by = 2, lb
863 0 : bz = lb - by
864 0 : f3 = f2*REAL(by - 1, dp)
865 : s(coa, coset(0, by, bz), coc) = &
866 : rbg(2)*s(coa, coset(0, by - 1, bz), coc) + &
867 : f3*s(coa, coset(0, by - 2, bz), coc) + &
868 0 : fcy*s(coa, coset(0, by - 1, bz), cocy)
869 : END DO
870 : ELSE
871 0 : bz = lb - 1
872 : s(coa, coset(0, 1, bz), coc) = &
873 : rbg(2)*s(coa, coset(0, 0, bz), coc) + &
874 : fy*s(coay, coset(0, 0, bz), coc) + &
875 0 : fcy*s(coa, coset(0, 0, bz), cocy)
876 0 : DO by = 2, lb
877 0 : bz = lb - by
878 0 : f3 = f2*REAL(by - 1, dp)
879 : s(coa, coset(0, by, bz), coc) = &
880 : rbg(2)*s(coa, coset(0, by - 1, bz), coc) + &
881 : fy*s(coay, coset(0, by - 1, bz), coc) + &
882 : f3*s(coa, coset(0, by - 2, bz), coc) + &
883 0 : fcy*s(coa, coset(0, by - 1, bz), cocy)
884 : END DO
885 : END IF
886 :
887 : ! *** Shift of angular momentum ***
888 : ! *** component x from a to b ***
889 :
890 0 : IF (ax == 0) THEN
891 0 : DO by = 0, lb - 1
892 0 : bz = lb - 1 - by
893 : s(coa, coset(1, by, bz), coc) = &
894 : rbg(1)*s(coa, coset(0, by, bz), coc) + &
895 0 : fcx*s(coa, coset(0, by, bz), cocx)
896 : END DO
897 0 : DO bx = 2, lb
898 0 : f3 = f2*REAL(bx - 1, dp)
899 0 : DO by = 0, lb - bx
900 0 : bz = lb - bx - by
901 : s(coa, coset(bx, by, bz), coc) = &
902 : rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + &
903 : f3*s(coa, coset(bx - 2, by, bz), coc) + &
904 0 : fcx*s(coa, coset(bx - 1, by, bz), cocx)
905 : END DO
906 : END DO
907 : ELSE
908 0 : DO by = 0, lb - 1
909 0 : bz = lb - 1 - by
910 : s(coa, coset(1, by, bz), coc) = &
911 : rbg(1)*s(coa, coset(0, by, bz), coc) + &
912 : fx*s(coax, coset(0, by, bz), coc) + &
913 0 : fcx*s(coa, coset(0, by, bz), cocx)
914 : END DO
915 0 : DO bx = 2, lb
916 0 : f3 = f2*REAL(bx - 1, dp)
917 0 : DO by = 0, lb - bx
918 0 : bz = lb - bx - by
919 : s(coa, coset(bx, by, bz), coc) = &
920 : rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + &
921 : fx*s(coax, coset(bx - 1, by, bz), coc) + &
922 : f3*s(coa, coset(bx - 2, by, bz), coc) + &
923 0 : fcx*s(coa, coset(bx - 1, by, bz), cocx)
924 : END DO
925 : END DO
926 : END IF
927 :
928 : END DO
929 : END DO
930 :
931 : END DO
932 :
933 : END IF
934 :
935 : ELSE
936 :
937 0 : IF (lb_max > 0) THEN
938 :
939 : ! *** Vertical recurrence steps: [s|c|s] -> [s|c|b] ***
940 :
941 0 : rbg(:) = rcg(:) + rbc(:)
942 :
943 : ! *** [s|c|p] = (Gi - Bi)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
944 :
945 0 : s(1, 2, coc) = rbg(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
946 0 : s(1, 3, coc) = rbg(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
947 0 : s(1, 4, coc) = rbg(3)*s(1, 1, coc) + fcz*s(1, 1, cocz)
948 :
949 : ! *** [s|c|b] = (Gi - Bi)*[s|c|b-1i] + ***
950 : ! *** f2*Ni(b-1i)*[s|c|b-2i] ***
951 : ! *** f2*Ni(c)*[s|c-1i|b-1i] ***
952 :
953 0 : DO lb = 2, lb_max
954 :
955 : ! *** Increase the angular momentum component z of b ***
956 :
957 : s(1, coset(0, 0, lb), coc) = &
958 : rbg(3)*s(1, coset(0, 0, lb - 1), coc) + &
959 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), coc) + &
960 0 : fcz*s(1, coset(0, 0, lb - 1), cocz)
961 :
962 : ! *** Increase the angular momentum component y of b ***
963 :
964 0 : bz = lb - 1
965 : s(1, coset(0, 1, bz), coc) = &
966 : rbg(2)*s(1, coset(0, 0, bz), coc) + &
967 0 : fcy*s(1, coset(0, 0, bz), cocy)
968 :
969 0 : DO by = 2, lb
970 0 : bz = lb - by
971 : s(1, coset(0, by, bz), coc) = &
972 : rbg(2)*s(1, coset(0, by - 1, bz), coc) + &
973 : f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), coc) + &
974 0 : fcy*s(1, coset(0, by - 1, bz), cocy)
975 : END DO
976 :
977 : ! *** Increase the angular momentum component x of b ***
978 :
979 0 : DO by = 0, lb - 1
980 0 : bz = lb - 1 - by
981 : s(1, coset(1, by, bz), coc) = &
982 : rbg(1)*s(1, coset(0, by, bz), coc) + &
983 0 : fcx*s(1, coset(0, by, bz), cocx)
984 : END DO
985 :
986 0 : DO bx = 2, lb
987 0 : f3 = f2*REAL(bx - 1, dp)
988 0 : DO by = 0, lb - bx
989 0 : bz = lb - bx - by
990 : s(1, coset(bx, by, bz), coc) = &
991 : rbg(1)*s(1, coset(bx - 1, by, bz), coc) + &
992 : f3*s(1, coset(bx - 2, by, bz), coc) + &
993 0 : fcx*s(1, coset(bx - 1, by, bz), cocx)
994 : END DO
995 : END DO
996 :
997 : END DO
998 :
999 : END IF
1000 :
1001 : END IF
1002 :
1003 : END DO
1004 : END DO
1005 :
1006 : END DO
1007 :
1008 : END IF
1009 :
1010 : ! *** Store integrals
1011 :
1012 2100 : IF (PRESENT(int_abc_ext)) THEN
1013 0 : DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1014 0 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1015 0 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1016 0 : sabc(na + i, nb + j, nc + k) = s(i, j, k)
1017 0 : int_abc_ext = MAX(int_abc_ext, ABS(s(i, j, k)))
1018 : END DO
1019 : END DO
1020 : END DO
1021 : ELSE
1022 6600 : DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1023 24600 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1024 94500 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1025 90000 : sabc(na + i, nb + j, nc + k) = s(i, j, k)
1026 : END DO
1027 : END DO
1028 : END DO
1029 : END IF
1030 :
1031 : ! *** Calculate the requested derivatives with respect to ***
1032 : ! *** the nuclear coordinates of the atomic center a and c ***
1033 :
1034 2100 : IF (PRESENT(sdabc) .OR. PRESENT(sabdc)) THEN
1035 : CALL derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, &
1036 : lc_max_set, lc_min_set, zeta(ipgf), zetc(kpgf), &
1037 1050 : s, sda, sdc)
1038 : END IF
1039 :
1040 : ! *** Store the first derivatives of the primitive overlap integrals ***
1041 :
1042 1050 : IF (PRESENT(sdabc)) THEN
1043 4200 : DO k = 1, 3
1044 11400 : DO l = 1, ncoset(lc_max_set)
1045 39150 : DO j = 1, ncoset(lb_max_set)
1046 151200 : DO i = 1, ncoset(la_max_set)
1047 144000 : sdabc(nda + i, nb + j, nc + l, k) = sda(i, j, l, k)
1048 : END DO
1049 : END DO
1050 : END DO
1051 : END DO
1052 : END IF
1053 :
1054 2100 : IF (PRESENT(sabdc)) THEN
1055 4200 : DO k = 1, 3
1056 11400 : DO l = 1, ncoset(lc_max_set)
1057 39150 : DO j = 1, ncoset(lb_max_set)
1058 151200 : DO i = 1, ncoset(la_max_set)
1059 144000 : sabdc(na + i, nb + j, ndc + l, k) = sdc(i, j, l, k)
1060 : END DO
1061 : END DO
1062 : END DO
1063 : END DO
1064 : END IF
1065 :
1066 2100 : nc = nc + ncoset(lc_max_set)
1067 2700 : ndc = ndc + ncoset(lc_max_set)
1068 : END DO
1069 :
1070 720 : nb = nb + ncoset(lb_max)
1071 : END DO
1072 :
1073 120 : na = na + ncoset(la_max_set)
1074 144 : nda = nda + ncoset(la_max_set)
1075 : END DO
1076 :
1077 24 : DEALLOCATE (s)
1078 24 : IF (PRESENT(sdabc)) THEN
1079 12 : DEALLOCATE (sda)
1080 : END IF
1081 24 : IF (PRESENT(sabdc)) THEN
1082 12 : DEALLOCATE (sdc)
1083 : END IF
1084 :
1085 24 : CALL timestop(handle)
1086 :
1087 24 : END SUBROUTINE overlap3
1088 :
1089 : ! **************************************************************************************************
1090 : !> \brief Calculates the derivatives of the three-center overlap integral [a|b|c]
1091 : !> with respect to the nuclear coordinates of the atomic center a and c
1092 : !> \param la_max_set ...
1093 : !> \param la_min_set ...
1094 : !> \param lb_max_set ...
1095 : !> \param lb_min_set ...
1096 : !> \param lc_max_set ...
1097 : !> \param lc_min_set ...
1098 : !> \param zeta ...
1099 : !> \param zetc ...
1100 : !> \param s integrals [a|b|c]
1101 : !> \param sda derivative [da/dAi|b|c]
1102 : !> \param sdc derivative [a|b|dc/dCi]
1103 : ! **************************************************************************************************
1104 1050 : SUBROUTINE derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, &
1105 : lc_max_set, lc_min_set, zeta, zetc, s, sda, sdc)
1106 :
1107 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, lb_max_set, &
1108 : lb_min_set, lc_max_set, lc_min_set
1109 : REAL(KIND=dp), INTENT(IN) :: zeta, zetc
1110 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: s
1111 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc
1112 :
1113 : CHARACTER(len=*), PARAMETER :: routineN = 'derivatives_overlap3'
1114 :
1115 : INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, coamy, coamz, coapx, coapy, coapz, cob, coc, &
1116 : cocmx, cocmy, cocmz, cocpx, cocpy, cocpz, cx, cy, cz, devx, devy, devz, handle, la, lb, lc
1117 : REAL(KIND=dp) :: fax, fay, faz, fcx, fcy, fcz, fexpa, &
1118 : fexpc
1119 :
1120 1050 : CALL timeset(routineN, handle)
1121 :
1122 1050 : fexpa = 2.0_dp*zeta
1123 1050 : fexpc = 2.0_dp*zetc
1124 :
1125 : ! derivative with respec to x,y,z
1126 :
1127 1050 : devx = 1
1128 1050 : devy = 2
1129 1050 : devz = 3
1130 :
1131 : ! *** [da/dAi|b|c] = 2*zeta*[a+1i|b|c] - Ni(a)[a-1i|b|c] ***
1132 : ! *** [a|b|dc/dCi] = 2*zetc*[a|b|c+1i] - Ni(c)[a|b|c-1i] ***
1133 :
1134 3150 : DO la = la_min_set, la_max_set
1135 6300 : DO ax = 0, la
1136 3150 : fax = REAL(ax, dp)
1137 9450 : DO ay = 0, la - ax
1138 4200 : fay = REAL(ay, dp)
1139 4200 : az = la - ax - ay
1140 4200 : faz = REAL(az, dp)
1141 4200 : coa = coset(ax, ay, az)
1142 4200 : coamx = coset(ax - 1, ay, az)
1143 4200 : coamy = coset(ax, ay - 1, az)
1144 4200 : coamz = coset(ax, ay, az - 1)
1145 4200 : coapx = coset(ax + 1, ay, az)
1146 4200 : coapy = coset(ax, ay + 1, az)
1147 4200 : coapz = coset(ax, ay, az + 1)
1148 15750 : DO lb = lb_min_set, lb_max_set
1149 25200 : DO bx = 0, lb
1150 37800 : DO by = 0, lb - bx
1151 16800 : bz = lb - bx - by
1152 16800 : cob = coset(bx, by, bz)
1153 48600 : DO lc = lc_min_set, lc_max_set
1154 62400 : DO cx = 0, lc
1155 26400 : fcx = REAL(cx, dp)
1156 81600 : DO cy = 0, lc - cx
1157 36000 : fcy = REAL(cy, dp)
1158 36000 : cz = lc - cx - cy
1159 36000 : fcz = REAL(cz, dp)
1160 36000 : coc = coset(cx, cy, cz)
1161 36000 : cocmx = coset(cx - 1, cy, cz)
1162 36000 : cocmy = coset(cx, cy - 1, cz)
1163 36000 : cocmz = coset(cx, cy, cz - 1)
1164 36000 : cocpx = coset(cx + 1, cy, cz)
1165 36000 : cocpy = coset(cx, cy + 1, cz)
1166 36000 : cocpz = coset(cx, cy, cz + 1)
1167 36000 : IF (ASSOCIATED(sda)) THEN
1168 : sda(coa, cob, coc, devx) = fexpa*s(coapx, cob, coc) - &
1169 36000 : fax*s(coamx, cob, coc)
1170 : sda(coa, cob, coc, devy) = fexpa*s(coapy, cob, coc) - &
1171 36000 : fay*s(coamy, cob, coc)
1172 : sda(coa, cob, coc, devz) = fexpa*s(coapz, cob, coc) - &
1173 36000 : faz*s(coamz, cob, coc)
1174 : END IF
1175 62400 : IF (ASSOCIATED(sdc)) THEN
1176 : sdc(coa, cob, coc, devx) = fexpc*s(coa, cob, cocpx) - &
1177 36000 : fcx*s(coa, cob, cocmx)
1178 : sdc(coa, cob, coc, devy) = fexpc*s(coa, cob, cocpy) - &
1179 36000 : fcy*s(coa, cob, cocmy)
1180 : sdc(coa, cob, coc, devz) = fexpc*s(coa, cob, cocpz) - &
1181 36000 : fcz*s(coa, cob, cocmz)
1182 : END IF
1183 : END DO
1184 : END DO
1185 : END DO
1186 : END DO
1187 : END DO
1188 : END DO
1189 : END DO
1190 : END DO
1191 : END DO
1192 :
1193 1050 : CALL timestop(handle)
1194 :
1195 1050 : END SUBROUTINE derivatives_overlap3
1196 :
1197 : END MODULE ai_overlap3
|
