(git:d2a9ebd)
Loading...
Searching...
No Matches
ai_overlap3.F
Go to the documentation of this file.
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
24MODULE 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
66CONTAINS
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, 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 lci = 0
145
146 IF (PRESENT(sdabc)) lai = 1
147 IF (PRESENT(sabdc)) lci = 1
148
149 la_max = la_max_set + lai
150 la_min = max(0, la_min_set - lai)
151 lb_max = lb_max_set
152 lb_min = lb_min_set
153 lc_max = lc_max_set + lci
154 lc_min = max(0, lc_min_set - lci)
155
156 ALLOCATE (s(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
157 s = 0._dp
158 IF (PRESENT(sdabc)) THEN
159 ALLOCATE (sda(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
160 sda = 0._dp
161 END IF
162 IF (PRESENT(sabdc)) THEN
163 ALLOCATE (sdc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3))
164 sdc = 0._dp
165 END IF
166 IF (PRESENT(int_abc_ext)) THEN
167 int_abc_ext = 0.0_dp
168 END IF
169
170! *** Loop over all pairs of primitive Gaussian-type functions ***
171
172 na = 0
173 nda = 0
174 DO ipgf = 1, npgfa
175
176 nb = 0
177 DO jpgf = 1, npgfb
178
179 ! *** Screening ***
180 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
181 sabc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
182 1:npgfc*ncoset(lc_max_set)) = 0.0_dp
183 IF (PRESENT(sdabc)) THEN
184 sdabc(nda + 1:nda + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
185 1:npgfc*ncoset(lc_max_set), 1:3) = 0.0_dp
186 END IF
187 IF (PRESENT(sabdc)) THEN
188 sabdc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
189 1:npgfc*ncoset(lc_max_set), 1:3) = 0.0_dp
190 END IF
191 nb = nb + ncoset(lb_max_set)
192 cycle
193 END IF
194
195 nc = 0
196 ndc = 0
197 DO kpgf = 1, npgfc
198
199 ! *** Screening ***
200 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 nc + 1:nc + ncoset(lc_max_set)) = 0.0_dp
204 IF (PRESENT(sdabc)) THEN
205 sdabc(nda + 1:nda + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
206 nc + 1:nc + ncoset(lc_max_set), 1:3) = 0.0_dp
207 END IF
208 IF (PRESENT(sabdc)) THEN
209 sabdc(na + 1:na + ncoset(la_max_set), nb + 1:nb + ncoset(lb_max_set), &
210 ndc + 1:ndc + ncoset(lc_max_set), 1:3) = 0.0_dp
211 END IF
212 nc = nc + ncoset(lc_max_set)
213 ndc = ndc + ncoset(lc_max_set)
214 cycle
215 END IF
216
217 ! *** Calculate some prefactors ***
218 zetg = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc(kpgf))
219 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
220 f0 = (pi*zetg)**1.5_dp
221 f1 = zetb(jpgf)*zetp
222 f2 = 0.5_dp*zetg
223 rcp(:) = f1*rab(:) - rac(:)
224 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 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 IF (la_max > 0) THEN
232
233! *** Vertical recurrence steps: [s|s|s] -> [a|s|s] ***
234
235 rag(:) = zetg*(zetb(jpgf)*rab(:) + zetc(kpgf)*rac(:))
236
237! *** [p|s|s] = (Gi - Ai)*[s|s|s] (i = x,y,z) ***
238
239 s(2, 1, 1) = rag(1)*s(1, 1, 1)
240 s(3, 1, 1) = rag(2)*s(1, 1, 1)
241 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 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 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 az = la - 1
255 s(coset(0, 1, az), 1, 1) = rag(2)*s(coset(0, 0, az), 1, 1)
256
257 DO ay = 2, la
258 az = la - ay
259 s(coset(0, ay, az), 1, 1) = rag(2)*s(coset(0, ay - 1, az), 1, 1) + &
260 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 DO ay = 0, la - 1
266 az = la - 1 - ay
267 s(coset(1, ay, az), 1, 1) = rag(1)*s(coset(0, ay, az), 1, 1)
268 END DO
269
270 DO ax = 2, la
271 f3 = f2*real(ax - 1, dp)
272 DO ay = 0, la - ax
273 az = la - ax - ay
274 s(coset(ax, ay, az), 1, 1) = rag(1)*s(coset(ax - 1, ay, az), 1, 1) + &
275 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 IF (lb_max > 0) THEN
284
285! *** Horizontal recurrence steps ***
286
287 rbg(:) = rag(:) - rab(:)
288
289! *** [a|s|p] = [a+1i|s|s] - (Bi - Ai)*[a|s|s] ***
290
291 IF (lb_max == 1) THEN
292 la_start = la_min
293 ELSE
294 la_start = max(0, la_min - 1)
295 END IF
296
297 DO la = la_start, la_max - 1
298 DO ax = 0, la
299 DO ay = 0, la - ax
300 az = la - ax - ay
301 coa = coset(ax, ay, az)
302 coax = coset(ax + 1, ay, az)
303 coay = coset(ax, ay + 1, az)
304 coaz = coset(ax, ay, az + 1)
305 s(coset(ax, ay, az), 2, 1) = s(coax, 1, 1) - rab(1)*s(coa, 1, 1)
306 s(coset(ax, ay, az), 3, 1) = s(coay, 1, 1) - rab(2)*s(coa, 1, 1)
307 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 DO ax = 0, la_max
317 fx = f2*real(ax, dp)
318 DO ay = 0, la_max - ax
319 fy = f2*real(ay, dp)
320 az = la_max - ax - ay
321 fz = f2*real(az, dp)
322 coa = coset(ax, ay, az)
323 IF (ax == 0) THEN
324 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1)
325 ELSE
326 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) + fx*s(coset(ax - 1, ay, az), 1, 1)
327 END IF
328 IF (ay == 0) THEN
329 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1)
330 ELSE
331 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) + fy*s(coset(ax, ay - 1, az), 1, 1)
332 END IF
333 IF (az == 0) THEN
334 s(coa, 4, 1) = rbg(3)*s(coa, 1, 1)
335 ELSE
336 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 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 IF (lb == lb_max) THEN
350 la_start = la_min
351 ELSE
352 la_start = max(0, la_min - 1)
353 END IF
354
355 DO la = la_start, la_max - 1
356 DO ax = 0, la
357 DO ay = 0, la - ax
358 az = la - ax - ay
359
360 coa = coset(ax, ay, az)
361 coax = coset(ax + 1, ay, az)
362 coay = coset(ax, ay + 1, az)
363 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 rab(3)*s(coa, coset(0, 0, lb - 1), 1)
370
371! *** Shift of angular momentum component y from a to b ***
372
373 DO by = 1, lb
374 bz = lb - by
375 s(coa, coset(0, by, bz), 1) = &
376 s(coay, coset(0, by - 1, bz), 1) - &
377 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 DO bx = 1, lb
383 DO by = 0, lb - bx
384 bz = lb - bx - by
385 s(coa, coset(bx, by, bz), 1) = &
386 s(coax, coset(bx - 1, by, bz), 1) - &
387 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 DO ax = 0, la_max
402 fx = f2*real(ax, dp)
403 DO ay = 0, la_max - ax
404 fy = f2*real(ay, dp)
405 az = la_max - ax - ay
406 fz = f2*real(az, dp)
407
408 coa = coset(ax, ay, az)
409
410 f3 = f2*real(lb - 1, dp)
411
412! *** Shift of angular momentum component z from a to b ***
413
414 IF (az == 0) THEN
415 s(coa, coset(0, 0, lb), 1) = &
416 rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + &
417 f3*s(coa, coset(0, 0, lb - 2), 1)
418 ELSE
419 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 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 IF (ay == 0) THEN
429 bz = lb - 1
430 s(coa, coset(0, 1, bz), 1) = &
431 rbg(2)*s(coa, coset(0, 0, bz), 1)
432 DO by = 2, lb
433 bz = lb - by
434 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 f3*s(coa, coset(0, by - 2, bz), 1)
438 END DO
439 ELSE
440 coay = coset(ax, ay - 1, az)
441 bz = lb - 1
442 s(coa, coset(0, 1, bz), 1) = &
443 rbg(2)*s(coa, coset(0, 0, bz), 1) + &
444 fy*s(coay, coset(0, 0, bz), 1)
445 DO by = 2, lb
446 bz = lb - by
447 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 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 IF (ax == 0) THEN
458 DO by = 0, lb - 1
459 bz = lb - 1 - by
460 s(coa, coset(1, by, bz), 1) = &
461 rbg(1)*s(coa, coset(0, by, bz), 1)
462 END DO
463 DO bx = 2, lb
464 f3 = f2*real(bx - 1, dp)
465 DO by = 0, lb - bx
466 bz = lb - bx - by
467 s(coa, coset(bx, by, bz), 1) = &
468 rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + &
469 f3*s(coa, coset(bx - 2, by, bz), 1)
470 END DO
471 END DO
472 ELSE
473 coax = coset(ax - 1, ay, az)
474 DO by = 0, lb - 1
475 bz = lb - 1 - by
476 s(coa, coset(1, by, bz), 1) = &
477 rbg(1)*s(coa, coset(0, by, bz), 1) + &
478 fx*s(coax, coset(0, by, bz), 1)
479 END DO
480 DO bx = 2, lb
481 f3 = f2*real(bx - 1, dp)
482 DO by = 0, lb - bx
483 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 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 IF (lb_max > 0) THEN
502
503! *** Vertical recurrence steps: [s|s|s] -> [s|s|b] ***
504
505 rbg(:) = -zetg*(zeta(ipgf)*rab(:) - zetc(kpgf)*rbc(:))
506
507! *** [s|s|p] = (Gi - Bi)*[s|s|s] ***
508
509 s(1, 2, 1) = rbg(1)*s(1, 1, 1)
510 s(1, 3, 1) = rbg(2)*s(1, 1, 1)
511 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 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 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 bz = lb - 1
525 s(1, coset(0, 1, bz), 1) = rbg(2)*s(1, coset(0, 0, bz), 1)
526
527 DO by = 2, lb
528 bz = lb - by
529 s(1, coset(0, by, bz), 1) = &
530 rbg(2)*s(1, coset(0, by - 1, bz), 1) + &
531 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 DO by = 0, lb - 1
537 bz = lb - 1 - by
538 s(1, coset(1, by, bz), 1) = rbg(1)*s(1, coset(0, by, bz), 1)
539 END DO
540
541 DO bx = 2, lb
542 f3 = f2*real(bx - 1, dp)
543 DO by = 0, lb - bx
544 bz = lb - bx - by
545 s(1, coset(bx, by, bz), 1) = rbg(1)*s(1, coset(bx - 1, by, bz), 1) + &
546 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 IF (lc_max > 0) THEN
559
560! *** Vertical recurrence steps: [s|s|s] -> [s|c|s] ***
561
562 rcg(:) = -zetg*(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))
563
564! *** [s|p|s] = (Gi - Ci)*[s|s|s] (i = x,y,z) ***
565
566 s(1, 1, 2) = rcg(1)*s(1, 1, 1)
567 s(1, 1, 3) = rcg(2)*s(1, 1, 1)
568 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 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 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 cz = lc - 1
582 s(1, 1, coset(0, 1, cz)) = rcg(2)*s(1, 1, coset(0, 0, cz))
583
584 DO cy = 2, lc
585 cz = lc - cy
586 s(1, 1, coset(0, cy, cz)) = rcg(2)*s(1, 1, coset(0, cy - 1, cz)) + &
587 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 DO cy = 0, lc - 1
593 cz = lc - 1 - cy
594 s(1, 1, coset(1, cy, cz)) = rcg(1)*s(1, 1, coset(0, cy, cz))
595 END DO
596
597 DO cx = 2, lc
598 f3 = f2*real(cx - 1, dp)
599 DO cy = 0, lc - cx
600 cz = lc - cx - cy
601 s(1, 1, coset(cx, cy, cz)) = rcg(1)*s(1, 1, coset(cx - 1, cy, cz)) + &
602 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 DO lc = 1, lc_max
611
612 DO cx = 0, lc
613 DO cy = 0, lc - cx
614 cz = lc - cx - cy
615
616 coc = coset(cx, cy, cz)
617 cocx = coset(max(0, cx - 1), cy, cz)
618 cocy = coset(cx, max(0, cy - 1), cz)
619 cocz = coset(cx, cy, max(0, cz - 1))
620
621 fcx = f2*real(cx, dp)
622 fcy = f2*real(cy, dp)
623 fcz = f2*real(cz, dp)
624
625! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
626
627 IF (la_max > 0) THEN
628
629! *** Vertical recurrence steps: [s|c|s] -> [a|c|s] ***
630
631 rag(:) = rcg(:) + rac(:)
632
633! *** [p|c|s] = (Gi - Ai)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
634
635 s(2, 1, coc) = rag(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
636 s(3, 1, coc) = rag(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
637 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 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 fcz*s(coset(0, 0, la - 1), 1, cocz)
651
652! *** Increase the angular momentum component y of a ***
653
654 az = la - 1
655 s(coset(0, 1, az), 1, coc) = &
656 rag(2)*s(coset(0, 0, az), 1, coc) + &
657 fcy*s(coset(0, 0, az), 1, cocy)
658
659 DO ay = 2, la
660 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 fcy*s(coset(0, ay - 1, az), 1, cocy)
665 END DO
666
667! *** Increase the angular momentum component x of a ***
668
669 DO ay = 0, la - 1
670 az = la - 1 - ay
671 s(coset(1, ay, az), 1, coc) = &
672 rag(1)*s(coset(0, ay, az), 1, coc) + &
673 fcx*s(coset(0, ay, az), 1, cocx)
674 END DO
675
676 DO ax = 2, la
677 f3 = f2*real(ax - 1, dp)
678 DO ay = 0, la - ax
679 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 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 IF (lb_max > 0) THEN
692
693! *** Horizontal recurrence steps ***
694
695 rbg(:) = rag(:) - rab(:)
696
697! *** [a|c|p] = [a+1i|c|s] - (Bi - Ai)*[a|c|s] ***
698
699 IF (lb_max == 1) THEN
700 la_start = la_min
701 ELSE
702 la_start = max(0, la_min - 1)
703 END IF
704
705 DO la = la_start, la_max - 1
706 DO ax = 0, la
707 DO ay = 0, la - ax
708 az = la - ax - ay
709 coa = coset(ax, ay, az)
710 coax = coset(ax + 1, ay, az)
711 coay = coset(ax, ay + 1, az)
712 coaz = coset(ax, ay, az + 1)
713 s(coa, 2, coc) = s(coax, 1, coc) - rab(1)*s(coa, 1, coc)
714 s(coa, 3, coc) = s(coay, 1, coc) - rab(2)*s(coa, 1, coc)
715 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 DO ax = 0, la_max
727 fx = f2*real(ax, dp)
728 DO ay = 0, la_max - ax
729 fy = f2*real(ay, dp)
730 az = la_max - ax - ay
731 fz = f2*real(az, dp)
732 coa = coset(ax, ay, az)
733 IF (ax == 0) THEN
734 s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + &
735 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 fcx*s(coa, 1, cocx)
740 END IF
741 IF (ay == 0) THEN
742 s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + &
743 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 fcy*s(coa, 1, cocy)
748 END IF
749 IF (az == 0) THEN
750 s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + &
751 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 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 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 IF (lb == lb_max) THEN
769 la_start = la_min
770 ELSE
771 la_start = max(0, la_min - 1)
772 END IF
773
774 DO la = la_start, la_max - 1
775 DO ax = 0, la
776 DO ay = 0, la - ax
777 az = la - ax - ay
778
779 coa = coset(ax, ay, az)
780 coax = coset(ax + 1, ay, az)
781 coay = coset(ax, ay + 1, az)
782 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 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 DO by = 1, lb
795 bz = lb - by
796 s(coa, coset(0, by, bz), coc) = &
797 s(coay, coset(0, by - 1, bz), coc) - &
798 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 DO bx = 1, lb
805 DO by = 0, lb - bx
806 bz = lb - bx - by
807 s(coa, coset(bx, by, bz), coc) = &
808 s(coax, coset(bx - 1, by, bz), coc) - &
809 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 DO ax = 0, la_max
825 fx = f2*real(ax, dp)
826 DO ay = 0, la_max - ax
827 fy = f2*real(ay, dp)
828 az = la_max - ax - ay
829 fz = f2*real(az, dp)
830
831 coa = coset(ax, ay, az)
832 coax = coset(max(0, ax - 1), ay, az)
833 coay = coset(ax, max(0, ay - 1), az)
834 coaz = coset(ax, ay, max(0, az - 1))
835
836 f3 = f2*real(lb - 1, dp)
837
838! *** Shift of angular momentum ***
839! *** component z from a to b ***
840
841 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 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 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 IF (ay == 0) THEN
858 bz = lb - 1
859 s(coa, coset(0, 1, bz), coc) = &
860 rbg(2)*s(coa, coset(0, 0, bz), coc) + &
861 fcy*s(coa, coset(0, 0, bz), cocy)
862 DO by = 2, lb
863 bz = lb - by
864 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 fcy*s(coa, coset(0, by - 1, bz), cocy)
869 END DO
870 ELSE
871 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 fcy*s(coa, coset(0, 0, bz), cocy)
876 DO by = 2, lb
877 bz = lb - by
878 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 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 IF (ax == 0) THEN
891 DO by = 0, lb - 1
892 bz = lb - 1 - by
893 s(coa, coset(1, by, bz), coc) = &
894 rbg(1)*s(coa, coset(0, by, bz), coc) + &
895 fcx*s(coa, coset(0, by, bz), cocx)
896 END DO
897 DO bx = 2, lb
898 f3 = f2*real(bx - 1, dp)
899 DO by = 0, lb - bx
900 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 fcx*s(coa, coset(bx - 1, by, bz), cocx)
905 END DO
906 END DO
907 ELSE
908 DO by = 0, lb - 1
909 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 fcx*s(coa, coset(0, by, bz), cocx)
914 END DO
915 DO bx = 2, lb
916 f3 = f2*real(bx - 1, dp)
917 DO by = 0, lb - bx
918 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 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 IF (lb_max > 0) THEN
938
939! *** Vertical recurrence steps: [s|c|s] -> [s|c|b] ***
940
941 rbg(:) = rcg(:) + rbc(:)
942
943! *** [s|c|p] = (Gi - Bi)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] ***
944
945 s(1, 2, coc) = rbg(1)*s(1, 1, coc) + fcx*s(1, 1, cocx)
946 s(1, 3, coc) = rbg(2)*s(1, 1, coc) + fcy*s(1, 1, cocy)
947 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 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 fcz*s(1, coset(0, 0, lb - 1), cocz)
961
962! *** Increase the angular momentum component y of b ***
963
964 bz = lb - 1
965 s(1, coset(0, 1, bz), coc) = &
966 rbg(2)*s(1, coset(0, 0, bz), coc) + &
967 fcy*s(1, coset(0, 0, bz), cocy)
968
969 DO by = 2, lb
970 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 fcy*s(1, coset(0, by - 1, bz), cocy)
975 END DO
976
977! *** Increase the angular momentum component x of b ***
978
979 DO by = 0, lb - 1
980 bz = lb - 1 - by
981 s(1, coset(1, by, bz), coc) = &
982 rbg(1)*s(1, coset(0, by, bz), coc) + &
983 fcx*s(1, coset(0, by, bz), cocx)
984 END DO
985
986 DO bx = 2, lb
987 f3 = f2*real(bx - 1, dp)
988 DO by = 0, lb - bx
989 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 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 IF (PRESENT(int_abc_ext)) THEN
1013 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1014 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1015 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1016 sabc(na + i, nb + j, nc + k) = s(i, j, k)
1017 int_abc_ext = max(int_abc_ext, abs(s(i, j, k)))
1018 END DO
1019 END DO
1020 END DO
1021 ELSE
1022 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set)
1023 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
1024 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
1025 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 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 s, sda, sdc)
1038 END IF
1039
1040! *** Store the first derivatives of the primitive overlap integrals ***
1041
1042 IF (PRESENT(sdabc)) THEN
1043 DO k = 1, 3
1044 DO l = 1, ncoset(lc_max_set)
1045 DO j = 1, ncoset(lb_max_set)
1046 DO i = 1, ncoset(la_max_set)
1047 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 IF (PRESENT(sabdc)) THEN
1055 DO k = 1, 3
1056 DO l = 1, ncoset(lc_max_set)
1057 DO j = 1, ncoset(lb_max_set)
1058 DO i = 1, ncoset(la_max_set)
1059 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 nc = nc + ncoset(lc_max_set)
1067 ndc = ndc + ncoset(lc_max_set)
1068 END DO
1069
1070 nb = nb + ncoset(lb_max)
1071 END DO
1072
1073 na = na + ncoset(la_max_set)
1074 nda = nda + ncoset(la_max_set)
1075 END DO
1076
1077 DEALLOCATE (s)
1078 IF (PRESENT(sdabc)) THEN
1079 DEALLOCATE (sda)
1080 END IF
1081 IF (PRESENT(sabdc)) THEN
1082 DEALLOCATE (sdc)
1083 END IF
1084
1085 CALL timestop(handle)
1086
1087 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 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 CALL timeset(routinen, handle)
1121
1122 fexpa = 2.0_dp*zeta
1123 fexpc = 2.0_dp*zetc
1124
1125! derivative with respec to x,y,z
1126
1127 devx = 1
1128 devy = 2
1129 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 DO la = la_min_set, la_max_set
1135 DO ax = 0, la
1136 fax = real(ax, dp)
1137 DO ay = 0, la - ax
1138 fay = real(ay, dp)
1139 az = la - ax - ay
1140 faz = real(az, dp)
1141 coa = coset(ax, ay, az)
1142 coamx = coset(ax - 1, ay, az)
1143 coamy = coset(ax, ay - 1, az)
1144 coamz = coset(ax, ay, az - 1)
1145 coapx = coset(ax + 1, ay, az)
1146 coapy = coset(ax, ay + 1, az)
1147 coapz = coset(ax, ay, az + 1)
1148 DO lb = lb_min_set, lb_max_set
1149 DO bx = 0, lb
1150 DO by = 0, lb - bx
1151 bz = lb - bx - by
1152 cob = coset(bx, by, bz)
1153 DO lc = lc_min_set, lc_max_set
1154 DO cx = 0, lc
1155 fcx = real(cx, dp)
1156 DO cy = 0, lc - cx
1157 fcy = real(cy, dp)
1158 cz = lc - cx - cy
1159 fcz = real(cz, dp)
1160 coc = coset(cx, cy, cz)
1161 cocmx = coset(cx - 1, cy, cz)
1162 cocmy = coset(cx, cy - 1, cz)
1163 cocmz = coset(cx, cy, cz - 1)
1164 cocpx = coset(cx + 1, cy, cz)
1165 cocpy = coset(cx, cy + 1, cz)
1166 cocpz = coset(cx, cy, cz + 1)
1167 IF (ASSOCIATED(sda)) THEN
1168 sda(coa, cob, coc, devx) = fexpa*s(coapx, cob, coc) - &
1169 fax*s(coamx, cob, coc)
1170 sda(coa, cob, coc, devy) = fexpa*s(coapy, cob, coc) - &
1171 fay*s(coamy, cob, coc)
1172 sda(coa, cob, coc, devz) = fexpa*s(coapz, cob, coc) - &
1173 faz*s(coamz, cob, coc)
1174 END IF
1175 IF (ASSOCIATED(sdc)) THEN
1176 sdc(coa, cob, coc, devx) = fexpc*s(coa, cob, cocpx) - &
1177 fcx*s(coa, cob, cocmx)
1178 sdc(coa, cob, coc, devy) = fexpc*s(coa, cob, cocpy) - &
1179 fcy*s(coa, cob, cocmy)
1180 sdc(coa, cob, coc, devz) = fexpc*s(coa, cob, cocpz) - &
1181 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 CALL timestop(handle)
1194
1195 END SUBROUTINE derivatives_overlap3
1196
1197END 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