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