(git:e7e05ae)
ai_overlap3.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 !!****** 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  SUBROUTINE overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
103  lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
104  lc_max_set, npgfc, zetc, rpgfc, lc_min_set, &
105  rab, dab, rac, dac, rbc, dbc, sabc, &
106  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  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc
136 
137 ! ---------------------------------------------------------------------------
138 
139  CALL timeset(routinen, handle)
140 
141  NULLIFY (s, sda, sdc)
142 
143  lai = 0
144  lbi = 0
145  lci = 0
146 
147  IF (PRESENT(sdabc)) lai = 1
148  IF (PRESENT(sabdc)) lci = 1
149 
150  la_max = la_max_set + lai
151  la_min = max(0, la_min_set - lai)
152  lb_max = lb_max_set
153  lb_min = lb_min_set
154  lc_max = lc_max_set + lci
155  lc_min = max(0, lc_min_set - lci)
156 
157  ALLOCATE (s(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
158  s = 0._dp
159  IF (PRESENT(sdabc)) THEN
160  ALLOCATE (sda(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
161  sda = 0._dp
162  END IF
163  IF (PRESENT(sabdc)) THEN
164  ALLOCATE (sdc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
165  sdc = 0._dp
166  END IF
167  IF (PRESENT(int_abc_ext)) THEN
168  int_abc_ext = 0.0_dp
169  END IF
170 
171 ! *** Loop over all pairs of primitive Gaussian-type functions ***
172 
173  na = 0
174  nda = 0
175  DO ipgf = 1, npgfa
176 
177  nb = 0
178  DO jpgf = 1, npgfb
179 
180  ! *** Screening ***
181  IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
182  nb = nb + ncoset(lb_max_set)
183  cycle
184  END IF
185 
186  nc = 0
187  ndc = 0
188  DO kpgf = 1, npgfc
189 
190  ! *** Screening ***
191  IF ((rpgfb(jpgf) + rpgfc(kpgf) < dbc) .OR. &
192  (rpgfa(ipgf) + rpgfc(kpgf) < dac)) THEN
193  nc = nc + ncoset(lc_max_set)
194  ndc = ndc + ncoset(lc_max_set)
195  cycle
196  END IF
197 
198  ! *** Calculate some prefactors ***
199  zetg = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc(kpgf))
200  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
201  f0 = (pi*zetg)**1.5_dp
202  f1 = zetb(jpgf)*zetp
203  f2 = 0.5_dp*zetg
204  rcp(:) = f1*rab(:) - rac(:)
205  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  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  IF (la_max > 0) THEN
213 
214 ! *** Vertical recurrence steps: [s|s|s] -> [a|s|s] ***
215 
216  rag(:) = zetg*(zetb(jpgf)*rab(:) + zetc(kpgf)*rac(:))
217 
218 ! *** [p|s|s] = (Gi - Ai)*[s|s|s] (i = x,y,z) ***
219 
220  s(2, 1, 1) = rag(1)*s(1, 1, 1)
221  s(3, 1, 1) = rag(2)*s(1, 1, 1)
222  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  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  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  az = la - 1
236  s(coset(0, 1, az), 1, 1) = rag(2)*s(coset(0, 0, az), 1, 1)
237 
238  DO ay = 2, la
239  az = la - ay
240  s(coset(0, ay, az), 1, 1) = rag(2)*s(coset(0, ay - 1, az), 1, 1) + &
241  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  DO ay = 0, la - 1
247  az = la - 1 - ay
248  s(coset(1, ay, az), 1, 1) = rag(1)*s(coset(0, ay, az), 1, 1)
249  END DO
250 
251  DO ax = 2, la
252  f3 = f2*real(ax - 1, dp)
253  DO ay = 0, la - ax
254  az = la - ax - ay
255  s(coset(ax, ay, az), 1, 1) = rag(1)*s(coset(ax - 1, ay, az), 1, 1) + &
256  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  IF (lb_max > 0) THEN
265 
266 ! *** Horizontal recurrence steps ***
267 
268  rbg(:) = rag(:) - rab(:)
269 
270 ! *** [a|s|p] = [a+1i|s|s] - (Bi - Ai)*[a|s|s] ***
271 
272  IF (lb_max == 1) THEN
273  la_start = la_min
274  ELSE
275  la_start = max(0, la_min - 1)
276  END IF
277 
278  DO la = la_start, la_max - 1
279  DO ax = 0, la
280  DO ay = 0, la - ax
281  az = la - ax - ay
282  coa = coset(ax, ay, az)
283  coax = coset(ax + 1, ay, az)
284  coay = coset(ax, ay + 1, az)
285  coaz = coset(ax, ay, az + 1)
286  s(coset(ax, ay, az), 2, 1) = s(coax, 1, 1) - rab(1)*s(coa, 1, 1)
287  s(coset(ax, ay, az), 3, 1) = s(coay, 1, 1) - rab(2)*s(coa, 1, 1)
288  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  DO ax = 0, la_max
298  fx = f2*real(ax, dp)
299  DO ay = 0, la_max - ax
300  fy = f2*real(ay, dp)
301  az = la_max - ax - ay
302  fz = f2*real(az, dp)
303  coa = coset(ax, ay, az)
304  IF (ax == 0) THEN
305  s(coa, 2, 1) = rbg(1)*s(coa, 1, 1)
306  ELSE
307  s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) + fx*s(coset(ax - 1, ay, az), 1, 1)
308  END IF
309  IF (ay == 0) THEN
310  s(coa, 3, 1) = rbg(2)*s(coa, 1, 1)
311  ELSE
312  s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) + fy*s(coset(ax, ay - 1, az), 1, 1)
313  END IF
314  IF (az == 0) THEN
315  s(coa, 4, 1) = rbg(3)*s(coa, 1, 1)
316  ELSE
317  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  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  IF (lb == lb_max) THEN
331  la_start = la_min
332  ELSE
333  la_start = max(0, la_min - 1)
334  END IF
335 
336  DO la = la_start, la_max - 1
337  DO ax = 0, la
338  DO ay = 0, la - ax
339  az = la - ax - ay
340 
341  coa = coset(ax, ay, az)
342  coax = coset(ax + 1, ay, az)
343  coay = coset(ax, ay + 1, az)
344  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  rab(3)*s(coa, coset(0, 0, lb - 1), 1)
351 
352 ! *** Shift of angular momentum component y from a to b ***
353 
354  DO by = 1, lb
355  bz = lb - by
356  s(coa, coset(0, by, bz), 1) = &
357  s(coay, coset(0, by - 1, bz), 1) - &
358  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  DO bx = 1, lb
364  DO by = 0, lb - bx
365  bz = lb - bx - by
366  s(coa, coset(bx, by, bz), 1) = &
367  s(coax, coset(bx - 1, by, bz), 1) - &
368  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  DO ax = 0, la_max
383  fx = f2*real(ax, dp)
384  DO ay = 0, la_max - ax
385  fy = f2*real(ay, dp)
386  az = la_max - ax - ay
387  fz = f2*real(az, dp)
388 
389  coa = coset(ax, ay, az)
390 
391  f3 = f2*real(lb - 1, dp)
392 
393 ! *** Shift of angular momentum component z from a to b ***
394 
395  IF (az == 0) THEN
396  s(coa, coset(0, 0, lb), 1) = &
397  rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
398  f3*s(coa, coset(0, 0, lb - 2), 1)
399  ELSE
400  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  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  IF (ay == 0) THEN
410  bz = lb - 1
411  s(coa, coset(0, 1, bz), 1) = &
412  rbg(2)*s(coa, coset(0, 0, bz), 1)
413  DO by = 2, lb
414  bz = lb - by
415  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  f3*s(coa, coset(0, by - 2, bz), 1)
419  END DO
420  ELSE
421  coay = coset(ax, ay - 1, az)
422  bz = lb - 1
423  s(coa, coset(0, 1, bz), 1) = &
424  rbg(2)*s(coa, coset(0, 0, bz), 1) + &
425  fy*s(coay, coset(0, 0, bz), 1)
426  DO by = 2, lb
427  bz = lb - by
428  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  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  IF (ax == 0) THEN
439  DO by = 0, lb - 1
440  bz = lb - 1 - by
441  s(coa, coset(1, by, bz), 1) = &
442  rbg(1)*s(coa, coset(0, by, bz), 1)
443  END DO
444  DO bx = 2, lb
445  f3 = f2*real(bx - 1, dp)
446  DO by = 0, lb - bx
447  bz = lb - bx - by
448  s(coa, coset(bx, by, bz), 1) = &
449  rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
450  f3*s(coa, coset(bx - 2, by, bz), 1)
451  END DO
452  END DO
453  ELSE
454  coax = coset(ax - 1, ay, az)
455  DO by = 0, lb - 1
456  bz = lb - 1 - by
457  s(coa, coset(1, by, bz), 1) = &
458  rbg(1)*s(coa, coset(0, by, bz), 1) + &
459  fx*s(coax, coset(0, by, bz), 1)
460  END DO
461  DO bx = 2, lb
462  f3 = f2*real(bx - 1, dp)
463  DO by = 0, lb - bx
464  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  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  IF (lb_max > 0) THEN
483 
484 ! *** Vertical recurrence steps: [s|s|s] -> [s|s|b] ***
485 
486  rbg(:) = -zetg*(zeta(ipgf)*rab(:) - zetc(kpgf)*rbc(:))
487 
488 ! *** [s|s|p] = (Gi - Bi)*[s|s|s] ***
489 
490  s(1, 2, 1) = rbg(1)*s(1, 1, 1)
491  s(1, 3, 1) = rbg(2)*s(1, 1, 1)
492  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  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  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  bz = lb - 1
506  s(1, coset(0, 1, bz), 1) = rbg(2)*s(1, coset(0, 0, bz), 1)
507 
508  DO by = 2, lb
509  bz = lb - by
510  s(1, coset(0, by, bz), 1) = &
511  rbg(2)*s(1, coset(0, by - 1, bz), 1) + &
512  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  DO by = 0, lb - 1
518  bz = lb - 1 - by
519  s(1, coset(1, by, bz), 1) = rbg(1)*s(1, coset(0, by, bz), 1)
520  END DO
521 
522  DO bx = 2, lb
523  f3 = f2*real(bx - 1, dp)
524  DO by = 0, lb - bx
525  bz = lb - bx - by
526  s(1, coset(bx, by, bz), 1) = rbg(1)*s(1, coset(bx - 1, by, bz), 1) + &
527  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  IF (lc_max > 0) THEN
540 
541 ! *** Vertical recurrence steps: [s|s|s] -> [s|c|s] ***
542 
543  rcg(:) = -zetg*(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))
544 
545 ! *** [s|p|s] = (Gi - Ci)*[s|s|s] (i = x,y,z) ***
546 
547  s(1, 1, 2) = rcg(1)*s(1, 1, 1)
548  s(1, 1, 3) = rcg(2)*s(1, 1, 1)
549  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  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  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  cz = lc - 1
563  s(1, 1, coset(0, 1, cz)) = rcg(2)*s(1, 1, coset(0, 0, cz))
564 
565  DO cy = 2, lc
566  cz = lc - cy
567  s(1, 1, coset(0, cy, cz)) = rcg(2)*s(1, 1, coset(0, cy - 1, cz)) + &
568  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  DO cy = 0, lc - 1
574  cz = lc - 1 - cy
575  s(1, 1, coset(1, cy, cz)) = rcg(1)*s(1, 1, coset(0, cy, cz))
576  END DO
577 
578  DO cx = 2, lc
579  f3 = f2*real(cx - 1, dp)
580  DO cy = 0, lc - cx
581  cz = lc - cx - cy
582  s(1, 1, coset(cx, cy, cz)) = rcg(1)*s(1, 1, coset(cx - 1, cy, cz)) + &
583  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  DO lc = 1, lc_max
592 
593  DO cx = 0, lc
594  DO cy = 0, lc - cx
595  cz = lc - cx - cy
596 
597  coc = coset(cx, cy, cz)
598  cocx = coset(max(0, cx - 1), cy, cz)
599  cocy = coset(cx, max(0, cy - 1), cz)
600  cocz = coset(cx, cy, max(0, cz - 1))
601 
602  fcx = f2*real(cx, dp)
603  fcy = f2*real(cy, dp)
604  fcz = f2*real(cz, dp)
605 
606 ! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
607 
608  IF (la_max > 0) THEN
609 
610 ! *** Vertical recurrence steps: [s|c|s] -> [a|c|s] ***
611 
612  rag(:) = rcg(:) + rac(:)
613 
614 ! *** [p|c|s] = (Gi - Ai)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
615 
616  s(2, 1, coc) = rag(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
617  s(3, 1, coc) = rag(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
618  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  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  fcz*s(coset(0, 0, la - 1), 1, cocz)
632 
633 ! *** Increase the angular momentum component y of a ***
634 
635  az = la - 1
636  s(coset(0, 1, az), 1, coc) = &
637  rag(2)*s(coset(0, 0, az), 1, coc) + &
638  fcy*s(coset(0, 0, az), 1, cocy)
639 
640  DO ay = 2, la
641  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  fcy*s(coset(0, ay - 1, az), 1, cocy)
646  END DO
647 
648 ! *** Increase the angular momentum component x of a ***
649 
650  DO ay = 0, la - 1
651  az = la - 1 - ay
652  s(coset(1, ay, az), 1, coc) = &
653  rag(1)*s(coset(0, ay, az), 1, coc) + &
654  fcx*s(coset(0, ay, az), 1, cocx)
655  END DO
656 
657  DO ax = 2, la
658  f3 = f2*real(ax - 1, dp)
659  DO ay = 0, la - ax
660  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  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  IF (lb_max > 0) THEN
673 
674 ! *** Horizontal recurrence steps ***
675 
676  rbg(:) = rag(:) - rab(:)
677 
678 ! *** [a|c|p] = [a+1i|c|s] - (Bi - Ai)*[a|c|s] ***
679 
680  IF (lb_max == 1) THEN
681  la_start = la_min
682  ELSE
683  la_start = max(0, la_min - 1)
684  END IF
685 
686  DO la = la_start, la_max - 1
687  DO ax = 0, la
688  DO ay = 0, la - ax
689  az = la - ax - ay
690  coa = coset(ax, ay, az)
691  coax = coset(ax + 1, ay, az)
692  coay = coset(ax, ay + 1, az)
693  coaz = coset(ax, ay, az + 1)
694  s(coa, 2, coc) = s(coax, 1, coc) - rab(1)*s(coa, 1, coc)
695  s(coa, 3, coc) = s(coay, 1, coc) - rab(2)*s(coa, 1, coc)
696  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  DO ax = 0, la_max
708  fx = f2*real(ax, dp)
709  DO ay = 0, la_max - ax
710  fy = f2*real(ay, dp)
711  az = la_max - ax - ay
712  fz = f2*real(az, dp)
713  coa = coset(ax, ay, az)
714  IF (ax == 0) THEN
715  s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
716  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  fcx*s(coa, 1, cocx)
721  END IF
722  IF (ay == 0) THEN
723  s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
724  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  fcy*s(coa, 1, cocy)
729  END IF
730  IF (az == 0) THEN
731  s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
732  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  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  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  IF (lb == lb_max) THEN
750  la_start = la_min
751  ELSE
752  la_start = max(0, la_min - 1)
753  END IF
754 
755  DO la = la_start, la_max - 1
756  DO ax = 0, la
757  DO ay = 0, la - ax
758  az = la - ax - ay
759 
760  coa = coset(ax, ay, az)
761  coax = coset(ax + 1, ay, az)
762  coay = coset(ax, ay + 1, az)
763  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  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  DO by = 1, lb
776  bz = lb - by
777  s(coa, coset(0, by, bz), coc) = &
778  s(coay, coset(0, by - 1, bz), coc) - &
779  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  DO bx = 1, lb
786  DO by = 0, lb - bx
787  bz = lb - bx - by
788  s(coa, coset(bx, by, bz), coc) = &
789  s(coax, coset(bx - 1, by, bz), coc) - &
790  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  DO ax = 0, la_max
806  fx = f2*real(ax, dp)
807  DO ay = 0, la_max - ax
808  fy = f2*real(ay, dp)
809  az = la_max - ax - ay
810  fz = f2*real(az, dp)
811 
812  coa = coset(ax, ay, az)
813  coax = coset(max(0, ax - 1), ay, az)
814  coay = coset(ax, max(0, ay - 1), az)
815  coaz = coset(ax, ay, max(0, az - 1))
816 
817  f3 = f2*real(lb - 1, dp)
818 
819 ! *** Shift of angular momentum ***
820 ! *** component z from a to b ***
821 
822  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  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  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  IF (ay == 0) THEN
839  bz = lb - 1
840  s(coa, coset(0, 1, bz), coc) = &
841  rbg(2)*s(coa, coset(0, 0, bz), coc) + &
842  fcy*s(coa, coset(0, 0, bz), cocy)
843  DO by = 2, lb
844  bz = lb - by
845  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  fcy*s(coa, coset(0, by - 1, bz), cocy)
850  END DO
851  ELSE
852  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  fcy*s(coa, coset(0, 0, bz), cocy)
857  DO by = 2, lb
858  bz = lb - by
859  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  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  IF (ax == 0) THEN
872  DO by = 0, lb - 1
873  bz = lb - 1 - by
874  s(coa, coset(1, by, bz), coc) = &
875  rbg(1)*s(coa, coset(0, by, bz), coc) + &
876  fcx*s(coa, coset(0, by, bz), cocx)
877  END DO
878  DO bx = 2, lb
879  f3 = f2*real(bx - 1, dp)
880  DO by = 0, lb - bx
881  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  fcx*s(coa, coset(bx - 1, by, bz), cocx)
886  END DO
887  END DO
888  ELSE
889  DO by = 0, lb - 1
890  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  fcx*s(coa, coset(0, by, bz), cocx)
895  END DO
896  DO bx = 2, lb
897  f3 = f2*real(bx - 1, dp)
898  DO by = 0, lb - bx
899  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  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  IF (lb_max > 0) THEN
919 
920 ! *** Vertical recurrence steps: [s|c|s] -> [s|c|b] ***
921 
922  rbg(:) = rcg(:) + rbc(:)
923 
924 ! *** [s|c|p] = (Gi - Bi)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
925 
926  s(1, 2, coc) = rbg(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
927  s(1, 3, coc) = rbg(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
928  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  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  fcz*s(1, coset(0, 0, lb - 1), cocz)
942 
943 ! *** Increase the angular momentum component y of b ***
944 
945  bz = lb - 1
946  s(1, coset(0, 1, bz), coc) = &
947  rbg(2)*s(1, coset(0, 0, bz), coc) + &
948  fcy*s(1, coset(0, 0, bz), cocy)
949 
950  DO by = 2, lb
951  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  fcy*s(1, coset(0, by - 1, bz), cocy)
956  END DO
957 
958 ! *** Increase the angular momentum component x of b ***
959 
960  DO by = 0, lb - 1
961  bz = lb - 1 - by
962  s(1, coset(1, by, bz), coc) = &
963  rbg(1)*s(1, coset(0, by, bz), coc) + &
964  fcx*s(1, coset(0, by, bz), cocx)
965  END DO
966 
967  DO bx = 2, lb
968  f3 = f2*real(bx - 1, dp)
969  DO by = 0, lb - bx
970  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  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  IF (PRESENT(int_abc_ext)) THEN
994  DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
995  DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
996  DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
997  sabc(na + i, nb + j, nc + k) = s(i, j, k)
998  int_abc_ext = max(int_abc_ext, abs(s(i, j, k)))
999  END DO
1000  END DO
1001  END DO
1002  ELSE
1003  DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1004  DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1005  DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1006  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  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  s, sda, sdc)
1019  END IF
1020 
1021 ! *** Store the first derivatives of the primitive overlap integrals ***
1022 
1023  IF (PRESENT(sdabc)) THEN
1024  DO k = 1, 3
1025  DO l = 1, ncoset(lc_max_set)
1026  DO j = 1, ncoset(lb_max_set)
1027  DO i = 1, ncoset(la_max_set)
1028  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  IF (PRESENT(sabdc)) THEN
1036  DO k = 1, 3
1037  DO l = 1, ncoset(lc_max_set)
1038  DO j = 1, ncoset(lb_max_set)
1039  DO i = 1, ncoset(la_max_set)
1040  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  nc = nc + ncoset(lc_max_set)
1048  ndc = ndc + ncoset(lc_max_set)
1049  END DO
1050 
1051  nb = nb + ncoset(lb_max)
1052  END DO
1053 
1054  na = na + ncoset(la_max_set)
1055  nda = nda + ncoset(la_max_set)
1056  END DO
1057 
1058  DEALLOCATE (s)
1059  IF (PRESENT(sdabc)) THEN
1060  DEALLOCATE (sda)
1061  END IF
1062  IF (PRESENT(sabdc)) THEN
1063  DEALLOCATE (sdc)
1064  END IF
1065 
1066  CALL timestop(handle)
1067 
1068  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  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  CALL timeset(routinen, handle)
1102 
1103  fexpa = 2.0_dp*zeta
1104  fexpc = 2.0_dp*zetc
1105 
1106 ! derivative with respec to x,y,z
1107 
1108  devx = 1
1109  devy = 2
1110  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  DO la = la_min_set, la_max_set
1116  DO ax = 0, la
1117  fax = real(ax, dp)
1118  DO ay = 0, la - ax
1119  fay = real(ay, dp)
1120  az = la - ax - ay
1121  faz = real(az, dp)
1122  coa = coset(ax, ay, az)
1123  coamx = coset(ax - 1, ay, az)
1124  coamy = coset(ax, ay - 1, az)
1125  coamz = coset(ax, ay, az - 1)
1126  coapx = coset(ax + 1, ay, az)
1127  coapy = coset(ax, ay + 1, az)
1128  coapz = coset(ax, ay, az + 1)
1129  DO lb = lb_min_set, lb_max_set
1130  DO bx = 0, lb
1131  DO by = 0, lb - bx
1132  bz = lb - bx - by
1133  cob = coset(bx, by, bz)
1134  DO lc = lc_min_set, lc_max_set
1135  DO cx = 0, lc
1136  fcx = real(cx, dp)
1137  DO cy = 0, lc - cx
1138  fcy = real(cy, dp)
1139  cz = lc - cx - cy
1140  fcz = real(cz, dp)
1141  coc = coset(cx, cy, cz)
1142  cocmx = coset(cx - 1, cy, cz)
1143  cocmy = coset(cx, cy - 1, cz)
1144  cocmz = coset(cx, cy, cz - 1)
1145  cocpx = coset(cx + 1, cy, cz)
1146  cocpy = coset(cx, cy + 1, cz)
1147  cocpz = coset(cx, cy, cz + 1)
1148  IF (ASSOCIATED(sda)) THEN
1149  sda(coa, cob, coc, devx) = fexpa*s(coapx, cob, coc) - &
1150  fax*s(coamx, cob, coc)
1151  sda(coa, cob, coc, devy) = fexpa*s(coapy, cob, coc) - &
1152  fay*s(coamy, cob, coc)
1153  sda(coa, cob, coc, devz) = fexpa*s(coapz, cob, coc) - &
1154  faz*s(coamz, cob, coc)
1155  END IF
1156  IF (ASSOCIATED(sdc)) THEN
1157  sdc(coa, cob, coc, devx) = fexpc*s(coa, cob, cocpx) - &
1158  fcx*s(coa, cob, cocmx)
1159  sdc(coa, cob, coc, devy) = fexpc*s(coa, cob, cocpy) - &
1160  fcy*s(coa, cob, cocmy)
1161  sdc(coa, cob, coc, devz) = fexpc*s(coa, cob, cocpz) - &
1162  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  CALL timestop(handle)
1175 
1176  END SUBROUTINE derivatives_overlap3
1177 
1178 END MODULE ai_overlap3
ai_overlap3
Definition: ai_overlap3.F:24
ai_overlap3::overlap3
subroutine, public overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, lc_max_set, npgfc, zetc, rpgfc, lc_min_set, rab, dab, rac, dac, rbc, dbc, sabc, sdabc, sabdc, int_abc_ext)
Calculation of three-center overlap integrals [a|b|c] over primitive Cartesian Gaussian functions.
Definition: ai_overlap3.F:107
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