(git:34ef472)
ai_overlap.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculation of the overlap integrals over Cartesian Gaussian-type
10 !> functions.
11 !> \par Literature
12 !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 !> \par History
14 !> - Derivatives added (02.05.2002,MK)
15 !> - New OS routine with simpler logic (11.07.2014, JGH)
16 !> \author Matthias Krack (08.10.1999)
17 ! **************************************************************************************************
18 MODULE ai_overlap
19  USE ai_os_rr, ONLY: os_rr_ovlp
20  USE kinds, ONLY: dp
21  USE mathconstants, ONLY: pi,&
22  twopi
23  USE orbital_pointers, ONLY: coset,&
24  nco,&
25  ncoset,&
26  nso
27  USE orbital_transformation_matrices, ONLY: orbtramat
28 #include "../base/base_uses.f90"
29 
30  IMPLICIT NONE
31 
32  PRIVATE
33 
34  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap'
35 
36 ! *** Public subroutines ***
37  PUBLIC :: overlap, overlap_ab, overlap_aab, overlap_ab_s, overlap_ab_sp, &
38  overlap_abb
39 
40 CONTAINS
41 
42 ! **************************************************************************************************
43 !> \brief Purpose: Calculation of the two-center overlap integrals [a|b] over
44 !> Cartesian Gaussian-type functions.
45 !> \param la_max_set Max L on center A
46 !> \param la_min_set Min L on center A
47 !> \param npgfa Number of primitives on center A
48 !> \param rpgfa Range of functions on A, used for screening
49 !> \param zeta Exponents on center A
50 !> \param lb_max_set Max L on center B
51 !> \param lb_min_set Min L on center B
52 !> \param npgfb Number of primitives on center B
53 !> \param rpgfb Range of functions on B, used for screening
54 !> \param zetb Exponents on center B
55 !> \param rab Distance vector A-B
56 !> \param dab Distance A-B
57 !> \param sab Final Integrals, basic and derivatives
58 !> \param da_max_set Some additional derivative information
59 !> \param return_derivatives Return integral derivatives
60 !> \param s Work space
61 !> \param lds Leading dimension of s
62 !> \param sdab Return additional derivative integrals
63 !> \param pab Density matrix block, used to calculate forces
64 !> \param force_a Force vector [da/dR|b]
65 !> \date 19.09.2000
66 !> \author MK
67 !> \version 1.0
68 ! **************************************************************************************************
69  SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
70  lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
71  rab, dab, sab, da_max_set, return_derivatives, s, lds, &
72  sdab, pab, force_a)
73  INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
74  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
75  INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
76  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
77  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
78  REAL(kind=dp), INTENT(IN) :: dab
79  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
80  INTEGER, INTENT(IN) :: da_max_set
81  LOGICAL, INTENT(IN) :: return_derivatives
82  INTEGER, INTENT(IN) :: lds
83  REAL(kind=dp), DIMENSION(lds, lds, *), &
84  INTENT(INOUT) :: s
85  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
86  OPTIONAL :: sdab
87  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
88  OPTIONAL :: pab
89  REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a
90 
91  INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
92  coapy, coapz, cob, cobm2x, cobm2y, cobm2z, cobmx, cobmy, cobmz, da, da_max, dax, day, &
93  daz, i, ipgf, j, jk, jpgf, jstart, k, la, la_max, la_min, la_start, lb, lb_max, lb_min, &
94  lb_start, na, nb, nda
95  LOGICAL :: calculate_force_a
96  REAL(kind=dp) :: f0, f1, f2, f3, f4, fax, fay, faz, ftz, &
97  zetp
98  REAL(kind=dp), DIMENSION(3) :: rap, rbp
99 
100  IF (PRESENT(pab) .AND. PRESENT(force_a)) THEN
101  calculate_force_a = .true.
102  force_a(:) = 0.0_dp
103  ELSE
104  calculate_force_a = .false.
105  END IF
106 
107  IF (PRESENT(sdab) .OR. calculate_force_a) THEN
108  IF (da_max_set == 0) THEN
109  da_max = 1
110  la_max = la_max_set + 1
111  la_min = max(0, la_min_set - 1)
112  ELSE
113  da_max = da_max_set
114  la_max = la_max_set + da_max_set + 1
115  la_min = max(0, la_min_set - da_max_set - 1)
116  END IF
117  ELSE
118  da_max = da_max_set
119  la_max = la_max_set + da_max_set
120  la_min = max(0, la_min_set - da_max_set)
121  END IF
122 
123  lb_max = lb_max_set
124  lb_min = lb_min_set
125 
126 ! *** Loop over all pairs of primitive Gaussian-type functions ***
127 
128  na = 0
129  nda = 0
130 
131  DO ipgf = 1, npgfa
132 
133  nb = 0
134 
135  DO jpgf = 1, npgfb
136 
137 ! *** Screening ***
138 
139  IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
140  DO j = nb + 1, nb + ncoset(lb_max_set)
141  DO i = na + 1, na + ncoset(la_max_set)
142  sab(i, j) = 0.0_dp
143  END DO
144  END DO
145  IF (return_derivatives) THEN
146  DO k = 2, ncoset(da_max_set)
147  jstart = (k - 1)*SIZE(sab, 1)
148  DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
149  DO i = na + 1, na + ncoset(la_max_set)
150  sab(i, j) = 0.0_dp
151  END DO
152  END DO
153  END DO
154  END IF
155  nb = nb + ncoset(lb_max_set)
156  cycle
157  END IF
158 
159 ! *** Calculate some prefactors ***
160 
161  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
162 
163  f0 = sqrt((pi*zetp)**3)
164  f1 = zetb(jpgf)*zetp
165  f2 = 0.5_dp*zetp
166 
167 ! *** Calculate the basic two-center overlap integral [s|s] ***
168 
169  s(1, 1, 1) = f0*exp(-zeta(ipgf)*f1*dab*dab)
170 
171 ! *** Recurrence steps: [s|s] -> [a|b] ***
172 
173  IF (la_max > 0) THEN
174 
175 ! *** Vertical recurrence steps: [s|s] -> [a|s] ***
176 
177  rap(:) = f1*rab(:)
178 
179 ! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) ***
180 
181  s(2, 1, 1) = rap(1)*s(1, 1, 1) ! [px|s]
182  s(3, 1, 1) = rap(2)*s(1, 1, 1) ! [py|s]
183  s(4, 1, 1) = rap(3)*s(1, 1, 1) ! [pz|s]
184 
185  IF (la_max > 1) THEN
186 
187 ! *** [d|s] ***
188 
189  f3 = f2*s(1, 1, 1)
190 
191  s(5, 1, 1) = rap(1)*s(2, 1, 1) + f3 ! [dx2|s]
192  s(6, 1, 1) = rap(1)*s(3, 1, 1) ! [dxy|s]
193  s(7, 1, 1) = rap(1)*s(4, 1, 1) ! [dxz|s]
194  s(8, 1, 1) = rap(2)*s(3, 1, 1) + f3 ! [dy2|s]
195  s(9, 1, 1) = rap(2)*s(4, 1, 1) ! [dyz|s]
196  s(10, 1, 1) = rap(3)*s(4, 1, 1) + f3 ! [dz2|s]
197 
198  IF (la_max > 2) THEN
199 
200 ! *** [f|s] ***
201 
202  f3 = 2.0_dp*f2
203 
204  s(11, 1, 1) = rap(1)*s(5, 1, 1) + f3*s(2, 1, 1) ! [fx3 |s]
205  s(12, 1, 1) = rap(1)*s(6, 1, 1) + f2*s(3, 1, 1) ! [fx2y|s]
206  s(13, 1, 1) = rap(1)*s(7, 1, 1) + f2*s(4, 1, 1) ! [fx2z|s]
207  s(14, 1, 1) = rap(1)*s(8, 1, 1) ! [fxy2|s]
208  s(15, 1, 1) = rap(1)*s(9, 1, 1) ! [fxyz|s]
209  s(16, 1, 1) = rap(1)*s(10, 1, 1) ! [fxz2|s]
210  s(17, 1, 1) = rap(2)*s(8, 1, 1) + f3*s(3, 1, 1) ! [fy3 |s]
211  s(18, 1, 1) = rap(2)*s(9, 1, 1) + f2*s(4, 1, 1) ! [fy2z|s]
212  s(19, 1, 1) = rap(2)*s(10, 1, 1) ! [fyz2|s]
213  s(20, 1, 1) = rap(3)*s(10, 1, 1) + f3*s(4, 1, 1) ! [fz3 |s]
214 
215  IF (la_max > 3) THEN
216 
217 ! *** [g|s] ***
218 
219  f4 = 3.0_dp*f2
220 
221  s(21, 1, 1) = rap(1)*s(11, 1, 1) + f4*s(5, 1, 1) ! [gx4 |s]
222  s(22, 1, 1) = rap(1)*s(12, 1, 1) + f3*s(6, 1, 1) ! [gx3y |s]
223  s(23, 1, 1) = rap(1)*s(13, 1, 1) + f3*s(7, 1, 1) ! [gx3z |s]
224  s(24, 1, 1) = rap(1)*s(14, 1, 1) + f2*s(8, 1, 1) ! [gx2y2|s]
225  s(25, 1, 1) = rap(1)*s(15, 1, 1) + f2*s(9, 1, 1) ! [gx2yz|s]
226  s(26, 1, 1) = rap(1)*s(16, 1, 1) + f2*s(10, 1, 1) ! [gx2z2|s]
227  s(27, 1, 1) = rap(1)*s(17, 1, 1) ! [gxy3 |s]
228  s(28, 1, 1) = rap(1)*s(18, 1, 1) ! [gxy2z|s]
229  s(29, 1, 1) = rap(1)*s(19, 1, 1) ! [gxyz2|s]
230  s(30, 1, 1) = rap(1)*s(20, 1, 1) ! [gxz3 |s]
231  s(31, 1, 1) = rap(2)*s(17, 1, 1) + f4*s(8, 1, 1) ! [gy4 |s]
232  s(32, 1, 1) = rap(2)*s(18, 1, 1) + f3*s(9, 1, 1) ! [gy3z |s]
233  s(33, 1, 1) = rap(2)*s(19, 1, 1) + f2*s(10, 1, 1) ! [gy2z2|s]
234  s(34, 1, 1) = rap(2)*s(20, 1, 1) ! [gyz3 |s]
235  s(35, 1, 1) = rap(3)*s(20, 1, 1) + f4*s(10, 1, 1) ! [gz4 |s]
236 
237 ! *** [a|s] = (Pi - Ai)*[a-1i|s] + f2*Ni(a-1i)*[a-2i|s] ***
238 
239  DO la = 5, la_max
240 
241 ! *** Increase the angular momentum component z of a ***
242 
243  s(coset(0, 0, la), 1, 1) = &
244  rap(3)*s(coset(0, 0, la - 1), 1, 1) + &
245  f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1, 1)
246 
247 ! *** Increase the angular momentum component y of a ***
248 
249  az = la - 1
250  s(coset(0, 1, az), 1, 1) = rap(2)*s(coset(0, 0, az), 1, 1)
251  DO ay = 2, la
252  az = la - ay
253  s(coset(0, ay, az), 1, 1) = &
254  rap(2)*s(coset(0, ay - 1, az), 1, 1) + &
255  f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1, 1)
256  END DO
257 
258 ! *** Increase the angular momentum component x of a ***
259 
260  DO ay = 0, la - 1
261  az = la - 1 - ay
262  s(coset(1, ay, az), 1, 1) = rap(1)*s(coset(0, ay, az), 1, 1)
263  END DO
264  DO ax = 2, la
265  f3 = f2*real(ax - 1, dp)
266  DO ay = 0, la - ax
267  az = la - ax - ay
268  s(coset(ax, ay, az), 1, 1) = &
269  rap(1)*s(coset(ax - 1, ay, az), 1, 1) + &
270  f3*s(coset(ax - 2, ay, az), 1, 1)
271  END DO
272  END DO
273 
274  END DO
275 
276  END IF
277 
278  END IF
279 
280  END IF
281 
282 ! *** Recurrence steps: [a|s] -> [a|b] ***
283 
284  IF (lb_max > 0) THEN
285 
286  DO j = 2, ncoset(lb_max)
287  DO i = 1, ncoset(la_min)
288  s(i, j, 1) = 0.0_dp
289  END DO
290  END DO
291 
292 ! *** Horizontal recurrence steps ***
293 
294  rbp(:) = rap(:) - rab(:)
295 
296 ! *** [a|p] = [a+1i|s] - (Bi - Ai)*[a|s] ***
297 
298  IF (lb_max == 1) THEN
299  la_start = la_min
300  ELSE
301  la_start = max(0, la_min - 1)
302  END IF
303 
304  DO la = la_start, la_max - 1
305  DO ax = 0, la
306  DO ay = 0, la - ax
307  az = la - ax - ay
308  coa = coset(ax, ay, az)
309  coapx = coset(ax + 1, ay, az)
310  coapy = coset(ax, ay + 1, az)
311  coapz = coset(ax, ay, az + 1)
312  s(coa, 2, 1) = s(coapx, 1, 1) - rab(1)*s(coa, 1, 1)
313  s(coa, 3, 1) = s(coapy, 1, 1) - rab(2)*s(coa, 1, 1)
314  s(coa, 4, 1) = s(coapz, 1, 1) - rab(3)*s(coa, 1, 1)
315  END DO
316  END DO
317  END DO
318 
319 ! *** Vertical recurrence step ***
320 
321 ! *** [a|p] = (Pi - Bi)*[a|s] + f2*Ni(a)*[a-1i|s] ***
322 
323  DO ax = 0, la_max
324  fax = f2*real(ax, dp)
325  DO ay = 0, la_max - ax
326  fay = f2*real(ay, dp)
327  az = la_max - ax - ay
328  faz = f2*real(az, dp)
329  coa = coset(ax, ay, az)
330  coamx = coset(ax - 1, ay, az)
331  coamy = coset(ax, ay - 1, az)
332  coamz = coset(ax, ay, az - 1)
333  s(coa, 2, 1) = rbp(1)*s(coa, 1, 1) + fax*s(coamx, 1, 1)
334  s(coa, 3, 1) = rbp(2)*s(coa, 1, 1) + fay*s(coamy, 1, 1)
335  s(coa, 4, 1) = rbp(3)*s(coa, 1, 1) + faz*s(coamz, 1, 1)
336  END DO
337  END DO
338 
339 ! *** Recurrence steps: [a|p] -> [a|b] ***
340 
341  DO lb = 2, lb_max
342 
343 ! *** Horizontal recurrence steps ***
344 
345 ! *** [a|b] = [a+1i|b-1i] - (Bi - Ai)*[a|b-1i] ***
346 
347  IF (lb == lb_max) THEN
348  la_start = la_min
349  ELSE
350  la_start = max(0, la_min - 1)
351  END IF
352 
353  DO la = la_start, la_max - 1
354  DO ax = 0, la
355  DO ay = 0, la - ax
356  az = la - ax - ay
357  coa = coset(ax, ay, az)
358  coapx = coset(ax + 1, ay, az)
359  coapy = coset(ax, ay + 1, az)
360  coapz = coset(ax, ay, az + 1)
361 
362 ! *** Shift of angular momentum component z from a to b ***
363 
364  cob = coset(0, 0, lb)
365  cobmz = coset(0, 0, lb - 1)
366  s(coa, cob, 1) = s(coapz, cobmz, 1) - rab(3)*s(coa, cobmz, 1)
367 
368 ! *** Shift of angular momentum component y from a to b ***
369 
370  DO by = 1, lb
371  bz = lb - by
372  cob = coset(0, by, bz)
373  cobmy = coset(0, by - 1, bz)
374  s(coa, cob, 1) = s(coapy, cobmy, 1) - rab(2)*s(coa, cobmy, 1)
375  END DO
376 
377 ! *** Shift of angular momentum component x from a to b ***
378 
379  DO bx = 1, lb
380  DO by = 0, lb - bx
381  bz = lb - bx - by
382  cob = coset(bx, by, bz)
383  cobmx = coset(bx - 1, by, bz)
384  s(coa, cob, 1) = s(coapx, cobmx, 1) - rab(1)*s(coa, cobmx, 1)
385  END DO
386  END DO
387 
388  END DO
389  END DO
390  END DO
391 
392 ! *** Vertical recurrence step ***
393 
394 ! *** [a|b] = (Pi - Bi)*[a|b-1i] + f2*Ni(a)*[a-1i|b-1i] + ***
395 ! *** f2*Ni(b-1i)*[a|b-2i] ***
396 
397  DO ax = 0, la_max
398  fax = f2*real(ax, dp)
399  DO ay = 0, la_max - ax
400  fay = f2*real(ay, dp)
401  az = la_max - ax - ay
402  faz = f2*real(az, dp)
403  coa = coset(ax, ay, az)
404  coamx = coset(ax - 1, ay, az)
405  coamy = coset(ax, ay - 1, az)
406  coamz = coset(ax, ay, az - 1)
407 
408 ! *** Increase the angular momentum component z of b ***
409 
410  f3 = f2*real(lb - 1, dp)
411  cob = coset(0, 0, lb)
412  cobmz = coset(0, 0, lb - 1)
413  cobm2z = coset(0, 0, lb - 2)
414  s(coa, cob, 1) = rbp(3)*s(coa, cobmz, 1) + &
415  faz*s(coamz, cobmz, 1) + &
416  f3*s(coa, cobm2z, 1)
417 
418 ! *** Increase the angular momentum component y of b ***
419 
420  bz = lb - 1
421  cob = coset(0, 1, bz)
422  cobmy = coset(0, 0, bz)
423  s(coa, cob, 1) = rbp(2)*s(coa, cobmy, 1) + &
424  fay*s(coamy, cobmy, 1)
425  DO by = 2, lb
426  bz = lb - by
427  f3 = f2*real(by - 1, dp)
428  cob = coset(0, by, bz)
429  cobmy = coset(0, by - 1, bz)
430  cobm2y = coset(0, by - 2, bz)
431  s(coa, cob, 1) = rbp(2)*s(coa, cobmy, 1) + &
432  fay*s(coamy, cobmy, 1) + &
433  f3*s(coa, cobm2y, 1)
434  END DO
435 
436 ! *** Increase the angular momentum component x of b ***
437 
438  DO by = 0, lb - 1
439  bz = lb - 1 - by
440  cob = coset(1, by, bz)
441  cobmx = coset(0, by, bz)
442  s(coa, cob, 1) = rbp(1)*s(coa, cobmx, 1) + &
443  fax*s(coamx, cobmx, 1)
444  END DO
445  DO bx = 2, lb
446  f3 = f2*real(bx - 1, dp)
447  DO by = 0, lb - bx
448  bz = lb - bx - by
449  cob = coset(bx, by, bz)
450  cobmx = coset(bx - 1, by, bz)
451  cobm2x = coset(bx - 2, by, bz)
452  s(coa, cob, 1) = rbp(1)*s(coa, cobmx, 1) + &
453  fax*s(coamx, cobmx, 1) + &
454  f3*s(coa, cobm2x, 1)
455  END DO
456  END DO
457 
458  END DO
459  END DO
460 
461  END DO
462 
463  END IF
464 
465  ELSE
466 
467  IF (lb_max > 0) THEN
468 
469 ! *** Vertical recurrence steps: [s|s] -> [s|b] ***
470 
471  rbp(:) = (f1 - 1.0_dp)*rab(:)
472 
473 ! *** [s|p] = (Pi - Bi)*[s|s] ***
474 
475  s(1, 2, 1) = rbp(1)*s(1, 1, 1) ! [s|px]
476  s(1, 3, 1) = rbp(2)*s(1, 1, 1) ! [s|py]
477  s(1, 4, 1) = rbp(3)*s(1, 1, 1) ! [s|pz]
478 
479  IF (lb_max > 1) THEN
480 
481 ! *** [s|d] ***
482 
483  f3 = f2*s(1, 1, 1)
484 
485  s(1, 5, 1) = rbp(1)*s(1, 2, 1) + f3 ! [s|dx2]
486  s(1, 6, 1) = rbp(1)*s(1, 3, 1) ! [s|dxy]
487  s(1, 7, 1) = rbp(1)*s(1, 4, 1) ! [s|dxz]
488  s(1, 8, 1) = rbp(2)*s(1, 3, 1) + f3 ! [s|dy2]
489  s(1, 9, 1) = rbp(2)*s(1, 4, 1) ! [s|dyz]
490  s(1, 10, 1) = rbp(3)*s(1, 4, 1) + f3 ! [s|dz2]
491 
492 ! *** [s|b] = (Pi - Bi)*[s|b-1i] + f2*Ni(b-1i)*[s|b-2i] ***
493 
494  DO lb = 3, lb_max
495 
496 ! *** Increase the angular momentum component z of b ***
497 
498  s(1, coset(0, 0, lb), 1) = &
499  rbp(3)*s(1, coset(0, 0, lb - 1), 1) + &
500  f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2), 1)
501 
502 ! *** Increase the angular momentum component y of b ***
503 
504  bz = lb - 1
505  s(1, coset(0, 1, bz), 1) = rbp(2)*s(1, coset(0, 0, bz), 1)
506  DO by = 2, lb
507  bz = lb - by
508  s(1, coset(0, by, bz), 1) = &
509  rbp(2)*s(1, coset(0, by - 1, bz), 1) + &
510  f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz), 1)
511  END DO
512 
513 ! *** Increase the angular momentum component x of b ***
514 
515  DO by = 0, lb - 1
516  bz = lb - 1 - by
517  s(1, coset(1, by, bz), 1) = rbp(1)*s(1, coset(0, by, bz), 1)
518  END DO
519  DO bx = 2, lb
520  f3 = f2*real(bx - 1, dp)
521  DO by = 0, lb - bx
522  bz = lb - bx - by
523  s(1, coset(bx, by, bz), 1) = &
524  rbp(1)*s(1, coset(bx - 1, by, bz), 1) + &
525  f3*s(1, coset(bx - 2, by, bz), 1)
526  END DO
527  END DO
528 
529  END DO
530 
531  END IF
532 
533  END IF
534 
535  END IF
536 
537 ! *** Store the primitive overlap integrals ***
538 
539  DO j = 1, ncoset(lb_max_set)
540  DO i = 1, ncoset(la_max_set)
541  sab(na + i, nb + j) = s(i, j, 1)
542  END DO
543  END DO
544 
545 ! *** Calculate the requested derivatives with respect ***
546 ! *** to the nuclear coordinates of the atomic center a ***
547 
548  IF (PRESENT(sdab) .OR. return_derivatives) THEN
549  la_start = 0
550  lb_start = 0
551  ELSE
552  la_start = la_min_set
553  lb_start = lb_min_set
554  END IF
555 
556  DO da = 0, da_max - 1
557  ftz = 2.0_dp*zeta(ipgf)
558  DO dax = 0, da
559  DO day = 0, da - dax
560  daz = da - dax - day
561  cda = coset(dax, day, daz)
562  cdax = coset(dax + 1, day, daz)
563  cday = coset(dax, day + 1, daz)
564  cdaz = coset(dax, day, daz + 1)
565 
566 ! *** [da/dAi|b] = 2*zeta*[a+1i|b] - Ni(a)[a-1i|b] ***
567 
568  DO la = la_start, la_max - da - 1
569  DO ax = 0, la
570  fax = real(ax, dp)
571  DO ay = 0, la - ax
572  fay = real(ay, dp)
573  az = la - ax - ay
574  faz = real(az, dp)
575  coa = coset(ax, ay, az)
576  coamx = coset(ax - 1, ay, az)
577  coamy = coset(ax, ay - 1, az)
578  coamz = coset(ax, ay, az - 1)
579  coapx = coset(ax + 1, ay, az)
580  coapy = coset(ax, ay + 1, az)
581  coapz = coset(ax, ay, az + 1)
582  DO lb = lb_start, lb_max_set
583  DO bx = 0, lb
584  DO by = 0, lb - bx
585  bz = lb - bx - by
586  cob = coset(bx, by, bz)
587  s(coa, cob, cdax) = ftz*s(coapx, cob, cda) - &
588  fax*s(coamx, cob, cda)
589  s(coa, cob, cday) = ftz*s(coapy, cob, cda) - &
590  fay*s(coamy, cob, cda)
591  s(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - &
592  faz*s(coamz, cob, cda)
593  END DO
594  END DO
595  END DO
596  END DO
597  END DO
598  END DO
599 
600  END DO
601  END DO
602  END DO
603 
604 ! *** Return all the calculated derivatives of the ***
605 ! *** primitive overlap integrals, if requested ***
606 
607  IF (return_derivatives) THEN
608  DO k = 2, ncoset(da_max_set)
609  jstart = (k - 1)*SIZE(sab, 1)
610  DO j = 1, ncoset(lb_max_set)
611  jk = jstart + j
612  DO i = 1, ncoset(la_max_set)
613  sab(na + i, nb + jk) = s(i, j, k)
614  END DO
615  END DO
616  END DO
617  END IF
618 
619 ! *** Calculate the force contribution for the atomic center a ***
620 
621  IF (calculate_force_a) THEN
622  DO k = 1, 3
623  DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
624  DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
625  force_a(k) = force_a(k) + pab(na + i, nb + j)*s(i, j, k + 1)
626  END DO
627  END DO
628  END DO
629  END IF
630 
631 ! *** Store the first derivatives of the primitive overlap integrals ***
632 ! *** which are used as auxiliary integrals for the calculation of ***
633 ! *** the kinetic energy integrals if requested ***
634 
635  IF (PRESENT(sdab)) THEN
636  sdab(nda + 1, nb + 1, 1) = s(1, 1, 1)
637  DO k = 2, 4
638  DO j = 1, ncoset(lb_max_set)
639  DO i = 1, ncoset(la_max - 1)
640  sdab(nda + i, nb + j, k) = s(i, j, k)
641  END DO
642  END DO
643  END DO
644  END IF
645 
646  nb = nb + ncoset(lb_max_set)
647 
648  END DO
649 
650  na = na + ncoset(la_max_set)
651  nda = nda + ncoset(la_max - 1)
652 
653  END DO
654 
655  END SUBROUTINE overlap
656 
657 ! **************************************************************************************************
658 !> \brief Calculation of the two-center overlap integrals [a|b] over
659 !> Cartesian Gaussian-type functions. First and second derivatives
660 !> \param la_max Max L on center A
661 !> \param la_min Min L on center A
662 !> \param npgfa Number of primitives on center A
663 !> \param rpgfa Range of functions on A, used for screening
664 !> \param zeta Exponents on center A
665 !> \param lb_max Max L on center B
666 !> \param lb_min Min L on center B
667 !> \param npgfb Number of primitives on center B
668 !> \param rpgfb Range of functions on B, used for screening
669 !> \param zetb Exponents on center B
670 !> \param rab Distance vector A-B
671 !> \param sab Final overlap integrals
672 !> \param dab First derivative overlap integrals
673 !> \param ddab Second derivative overlap integrals
674 !> \date 01.07.2014
675 !> \author JGH
676 ! **************************************************************************************************
677  SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
678  lb_max, lb_min, npgfb, rpgfb, zetb, &
679  rab, sab, dab, ddab)
680  INTEGER, INTENT(IN) :: la_max, la_min, npgfa
681  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
682  INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
683  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
684  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
685  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
686  OPTIONAL :: sab
687  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
688  OPTIONAL :: dab, ddab
689 
690  INTEGER :: ax, ay, az, bx, by, bz, coa, cob, ia, &
691  ib, ipgf, jpgf, la, lb, ldrr, lma, &
692  lmb, ma, mb, na, nb, ofa, ofb
693  REAL(kind=dp) :: a, ambm, ambp, apbm, apbp, b, dumx, &
694  dumy, dumz, f0, rab2, tab, xhi, zet
695  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
696  REAL(kind=dp), DIMENSION(3) :: rap, rbp
697 
698  ! Distance of the centers a and b
699 
700  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
701  tab = sqrt(rab2)
702 
703  ! Maximum l for auxiliary integrals
704  IF (PRESENT(sab)) THEN
705  lma = la_max
706  lmb = lb_max
707  END IF
708  IF (PRESENT(dab)) THEN
709  lma = la_max + 1
710  lmb = lb_max
711  END IF
712  IF (PRESENT(ddab)) THEN
713  lma = la_max + 1
714  lmb = lb_max + 1
715  END IF
716  ldrr = max(lma, lmb) + 1
717 
718  ! Allocate space for auxiliary integrals
719  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
720 
721  ! Number of integrals, check size of arrays
722  ofa = ncoset(la_min - 1)
723  ofb = ncoset(lb_min - 1)
724  na = ncoset(la_max) - ofa
725  nb = ncoset(lb_max) - ofb
726  IF (PRESENT(sab)) THEN
727  cpassert((SIZE(sab, 1) >= na*npgfa))
728  cpassert((SIZE(sab, 2) >= nb*npgfb))
729  END IF
730  IF (PRESENT(dab)) THEN
731  cpassert((SIZE(dab, 1) >= na*npgfa))
732  cpassert((SIZE(dab, 2) >= nb*npgfb))
733  cpassert((SIZE(dab, 3) >= 3))
734  END IF
735  IF (PRESENT(ddab)) THEN
736  cpassert((SIZE(ddab, 1) >= na*npgfa))
737  cpassert((SIZE(ddab, 2) >= nb*npgfb))
738  cpassert((SIZE(ddab, 3) >= 6))
739  END IF
740 
741  ! Loops over all pairs of primitive Gaussian-type functions
742  ma = 0
743  DO ipgf = 1, npgfa
744  mb = 0
745  DO jpgf = 1, npgfb
746  ! Distance Screening
747  IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
748  IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
749  IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
750  IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
751  mb = mb + nb
752  cycle
753  END IF
754 
755  ! Calculate some prefactors
756  a = zeta(ipgf)
757  b = zetb(jpgf)
758  zet = a + b
759  xhi = a*b/zet
760  rap = b*rab/zet
761  rbp = -a*rab/zet
762 
763  ! [s|s] integral
764  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
765 
766  ! Calculate the recurrence relation
767  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
768 
769  DO lb = lb_min, lb_max
770  DO bx = 0, lb
771  DO by = 0, lb - bx
772  bz = lb - bx - by
773  cob = coset(bx, by, bz) - ofb
774  ib = mb + cob
775  DO la = la_min, la_max
776  DO ax = 0, la
777  DO ay = 0, la - ax
778  az = la - ax - ay
779  coa = coset(ax, ay, az) - ofa
780  ia = ma + coa
781  ! integrals
782  IF (PRESENT(sab)) THEN
783  sab(ia, ib) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
784  END IF
785  ! first derivatives
786  IF (PRESENT(dab)) THEN
787  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
788  ! dx
789  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
790  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
791  dab(ia, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
792  ! dy
793  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
794  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
795  dab(ia, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
796  ! dz
797  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
798  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
799  dab(ia, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
800  END IF
801  ! 2nd derivatives
802  IF (PRESENT(ddab)) THEN
803  ! (dda|b) = -4*a*b*(a+1|b+1) + 2*a*N(b)*(a+1|b-1)
804  ! + 2*b*N(a)*(a-1|b+1) - N(a)*N(b)*(a-1|b-1)
805  ! dx dx
806  apbp = f0*rr(ax + 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
807  IF (bx > 0) THEN
808  apbm = f0*rr(ax + 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
809  ELSE
810  apbm = 0.0_dp
811  END IF
812  IF (ax > 0) THEN
813  ambp = f0*rr(ax - 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
814  ELSE
815  ambp = 0.0_dp
816  END IF
817  IF (ax > 0 .AND. bx > 0) THEN
818  ambm = f0*rr(ax - 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
819  ELSE
820  ambm = 0.0_dp
821  END IF
822  ddab(ia, ib, 1) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bx, dp)*apbm &
823  + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bx, dp)*ambm
824  ! dx dy
825  apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
826  IF (by > 0) THEN
827  apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
828  ELSE
829  apbm = 0.0_dp
830  END IF
831  IF (ax > 0) THEN
832  ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
833  ELSE
834  ambp = 0.0_dp
835  END IF
836  IF (ax > 0 .AND. by > 0) THEN
837  ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
838  ELSE
839  ambm = 0.0_dp
840  END IF
841  ddab(ia, ib, 2) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(by, dp)*apbm &
842  + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(by, dp)*ambm
843  ! dx dz
844  apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
845  IF (bz > 0) THEN
846  apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
847  ELSE
848  apbm = 0.0_dp
849  END IF
850  IF (ax > 0) THEN
851  ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
852  ELSE
853  ambp = 0.0_dp
854  END IF
855  IF (ax > 0 .AND. bz > 0) THEN
856  ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
857  ELSE
858  ambm = 0.0_dp
859  END IF
860  ddab(ia, ib, 3) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
861  + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bz, dp)*ambm
862  ! dy dy
863  apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by + 1, 2)*rr(az, bz, 3)
864  IF (by > 0) THEN
865  apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by - 1, 2)*rr(az, bz, 3)
866  ELSE
867  apbm = 0.0_dp
868  END IF
869  IF (ay > 0) THEN
870  ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by + 1, 2)*rr(az, bz, 3)
871  ELSE
872  ambp = 0.0_dp
873  END IF
874  IF (ay > 0 .AND. by > 0) THEN
875  ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by - 1, 2)*rr(az, bz, 3)
876  ELSE
877  ambm = 0.0_dp
878  END IF
879  ddab(ia, ib, 4) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(by, dp)*apbm &
880  + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(by, dp)*ambm
881  ! dy dz
882  apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz + 1, 3)
883  IF (bz > 0) THEN
884  apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz - 1, 3)
885  ELSE
886  apbm = 0.0_dp
887  END IF
888  IF (ay > 0) THEN
889  ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz + 1, 3)
890  ELSE
891  ambp = 0.0_dp
892  END IF
893  IF (ay > 0 .AND. bz > 0) THEN
894  ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz - 1, 3)
895  ELSE
896  ambm = 0.0_dp
897  END IF
898  ddab(ia, ib, 5) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
899  + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(bz, dp)*ambm
900  ! dz dz
901  apbp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz + 1, 3)
902  IF (bz > 0) THEN
903  apbm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz - 1, 3)
904  ELSE
905  apbm = 0.0_dp
906  END IF
907  IF (az > 0) THEN
908  ambp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz + 1, 3)
909  ELSE
910  ambp = 0.0_dp
911  END IF
912  IF (az > 0 .AND. bz > 0) THEN
913  ambm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz - 1, 3)
914  ELSE
915  ambm = 0.0_dp
916  END IF
917  ddab(ia, ib, 6) = -4.0_dp*a*b*apbm + 2.0_dp*a*real(bz, dp)*apbm &
918  + 2.0_dp*b*real(az, dp)*ambp - real(az, dp)*real(bz, dp)*ambm
919  END IF
920  !
921  END DO
922  END DO
923  END DO !la
924  END DO
925  END DO
926  END DO !lb
927 
928  mb = mb + nb
929  END DO
930  ma = ma + na
931  END DO
932 
933  DEALLOCATE (rr)
934 
935  END SUBROUTINE overlap_ab
936 
937 ! **************************************************************************************************
938 !> \brief Calculation of the two-center overlap integrals [aa|b] over
939 !> Cartesian Gaussian-type functions.
940 !> \param la1_max Max L on center A (basis 1)
941 !> \param la1_min Min L on center A (basis 1)
942 !> \param npgfa1 Number of primitives on center A (basis 1)
943 !> \param rpgfa1 Range of functions on A, used for screening (basis 1)
944 !> \param zeta1 Exponents on center A (basis 1)
945 !> \param la2_max Max L on center A (basis 2)
946 !> \param la2_min Min L on center A (basis 2)
947 !> \param npgfa2 Number of primitives on center A (basis 2)
948 !> \param rpgfa2 Range of functions on A, used for screening (basis 2)
949 !> \param zeta2 Exponents on center A (basis 2)
950 !> \param lb_max Max L on center B
951 !> \param lb_min Min L on center B
952 !> \param npgfb Number of primitives on center B
953 !> \param rpgfb Range of functions on B, used for screening
954 !> \param zetb Exponents on center B
955 !> \param rab Distance vector A-B
956 !> \param saab Final overlap integrals
957 !> \param daab First derivative overlap integrals
958 !> \param saba Final overlap integrals; different order
959 !> \param daba First derivative overlap integrals; different order
960 !> \date 01.07.2014
961 !> \author JGH
962 ! **************************************************************************************************
963  SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
964  la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
965  lb_max, lb_min, npgfb, rpgfb, zetb, &
966  rab, saab, daab, saba, daba)
967  INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
968  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
969  INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
970  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
971  INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
972  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
973  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
974  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
975  OPTIONAL :: saab
976  REAL(kind=dp), DIMENSION(:, :, :, :), &
977  INTENT(INOUT), OPTIONAL :: daab
978  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
979  OPTIONAL :: saba
980  REAL(kind=dp), DIMENSION(:, :, :, :), &
981  INTENT(INOUT), OPTIONAL :: daba
982 
983  INTEGER :: ax, ax1, ax2, ay, ay1, ay2, az, az1, az2, bx, by, bz, coa1, coa2, cob, i1pgf, &
984  i2pgf, ia1, ia2, ib, jpgf, la1, la2, lb, ldrr, lma, lmb, ma1, ma2, mb, na1, na2, nb, &
985  ofa1, ofa2, ofb
986  REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
987  tab, xhi, zet
988  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
989  REAL(kind=dp), DIMENSION(3) :: rap, rbp
990 
991  ! Distance of the centers a and b
992 
993  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
994  tab = sqrt(rab2)
995 
996  ! Maximum l for auxiliary integrals
997  IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
998  lma = la1_max + la2_max
999  lmb = lb_max
1000  END IF
1001  IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
1002  lma = la1_max + la2_max + 1
1003  lmb = lb_max
1004  END IF
1005  ldrr = max(lma, lmb) + 1
1006 
1007  ! Allocate space for auxiliary integrals
1008  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1009 
1010  ! Number of integrals, check size of arrays
1011  ofa1 = ncoset(la1_min - 1)
1012  ofa2 = ncoset(la2_min - 1)
1013  ofb = ncoset(lb_min - 1)
1014  na1 = ncoset(la1_max) - ofa1
1015  na2 = ncoset(la2_max) - ofa2
1016  nb = ncoset(lb_max) - ofb
1017  IF (PRESENT(saab)) THEN
1018  cpassert((SIZE(saab, 1) >= na1*npgfa1))
1019  cpassert((SIZE(saab, 2) >= na2*npgfa2))
1020  cpassert((SIZE(saab, 3) >= nb*npgfb))
1021  END IF
1022  IF (PRESENT(daab)) THEN
1023  cpassert((SIZE(daab, 1) >= na1*npgfa1))
1024  cpassert((SIZE(daab, 2) >= na2*npgfa2))
1025  cpassert((SIZE(daab, 3) >= nb*npgfb))
1026  cpassert((SIZE(daab, 4) >= 3))
1027  END IF
1028  IF (PRESENT(saba)) THEN
1029  cpassert((SIZE(saba, 1) >= na1*npgfa1))
1030  cpassert((SIZE(saba, 2) >= nb*npgfb))
1031  cpassert((SIZE(saba, 3) >= na2*npgfa2))
1032  END IF
1033  IF (PRESENT(daba)) THEN
1034  cpassert((SIZE(daba, 1) >= na1*npgfa1))
1035  cpassert((SIZE(daba, 2) >= nb*npgfb))
1036  cpassert((SIZE(daba, 3) >= na2*npgfa2))
1037  cpassert((SIZE(daba, 4) >= 3))
1038  END IF
1039 
1040  ! Loops over all primitive Gaussian-type functions
1041  ma1 = 0
1042  DO i1pgf = 1, npgfa1
1043  ma2 = 0
1044  DO i2pgf = 1, npgfa2
1045  rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
1046  mb = 0
1047  DO jpgf = 1, npgfb
1048  ! Distance Screening
1049  IF (rpgfa + rpgfb(jpgf) < tab) THEN
1050  IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
1051  IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
1052  IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
1053  IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
1054  mb = mb + nb
1055  cycle
1056  END IF
1057 
1058  ! Calculate some prefactors
1059  a = zeta1(i1pgf) + zeta2(i2pgf)
1060  b = zetb(jpgf)
1061  zet = a + b
1062  xhi = a*b/zet
1063  rap = b*rab/zet
1064  rbp = -a*rab/zet
1065 
1066  ! [ss|s] integral
1067  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1068 
1069  ! Calculate the recurrence relation
1070  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1071 
1072  DO lb = lb_min, lb_max
1073  DO bx = 0, lb
1074  DO by = 0, lb - bx
1075  bz = lb - bx - by
1076  cob = coset(bx, by, bz) - ofb
1077  ib = mb + cob
1078  DO la2 = la2_min, la2_max
1079  DO ax2 = 0, la2
1080  DO ay2 = 0, la2 - ax2
1081  az2 = la2 - ax2 - ay2
1082  coa2 = coset(ax2, ay2, az2) - ofa2
1083  ia2 = ma2 + coa2
1084  DO la1 = la1_min, la1_max
1085  DO ax1 = 0, la1
1086  DO ay1 = 0, la1 - ax1
1087  az1 = la1 - ax1 - ay1
1088  coa1 = coset(ax1, ay1, az1) - ofa1
1089  ia1 = ma1 + coa1
1090  ! integrals
1091  IF (PRESENT(saab)) THEN
1092  saab(ia1, ia2, ib) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
1093  END IF
1094  IF (PRESENT(saba)) THEN
1095  saba(ia1, ib, ia2) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
1096  END IF
1097  ! first derivatives
1098  IF (PRESENT(daab)) THEN
1099  ax = ax1 + ax2
1100  ay = ay1 + ay2
1101  az = az1 + az2
1102  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1103  ! dx
1104  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1105  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1106  daab(ia1, ia2, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1107  ! dy
1108  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1109  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1110  daab(ia1, ia2, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1111  ! dz
1112  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1113  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1114  daab(ia1, ia2, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1115  END IF
1116  IF (PRESENT(daba)) THEN
1117  ax = ax1 + ax2
1118  ay = ay1 + ay2
1119  az = az1 + az2
1120  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1121  ! dx
1122  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1123  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1124  daba(ia1, ib, ia2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1125  ! dy
1126  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1127  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1128  daba(ia1, ib, ia2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1129  ! dz
1130  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1131  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1132  daba(ia1, ib, ia2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1133  END IF
1134  !
1135  END DO
1136  END DO
1137  END DO !la1
1138  END DO
1139  END DO
1140  END DO !la2
1141  END DO
1142  END DO
1143  END DO !lb
1144 
1145  mb = mb + nb
1146  END DO
1147  ma2 = ma2 + na2
1148  END DO
1149  ma1 = ma1 + na1
1150  END DO
1151 
1152  DEALLOCATE (rr)
1153 
1154  END SUBROUTINE overlap_aab
1155 
1156 ! **************************************************************************************************
1157 !> \brief Calculation of the two-center overlap integrals [a|bb] over
1158 !> Cartesian Gaussian-type functions.
1159 !> \param la_max Max L on center A
1160 !> \param la_min Min L on center A
1161 !> \param npgfa Number of primitives on center A
1162 !> \param rpgfa Range of functions on A, used for screening
1163 !> \param zeta Exponents on center A
1164 !> \param lb1_max Max L on center B (basis 1)
1165 !> \param lb1_min Min L on center B (basis 1)
1166 !> \param npgfb1 Number of primitives on center B (basis 1)
1167 !> \param rpgfb1 Range of functions on B, used for screening (basis 1)
1168 !> \param zetb1 Exponents on center B (basis 1)
1169 !> \param lb2_max Max L on center B (basis 2)
1170 !> \param lb2_min Min L on center B (basis 2)
1171 !> \param npgfb2 Number of primitives on center B (basis 2)
1172 !> \param rpgfb2 Range of functions on B, used for screening (basis 2)
1173 !> \param zetb2 Exponents on center B (basis 2)
1174 !> \param rab Distance vector A-B
1175 !> \param sabb Final overlap integrals
1176 !> \param dabb First derivative overlap integrals
1177 !> \date 01.07.2014
1178 !> \author JGH
1179 ! **************************************************************************************************
1180  SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
1181  lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1182  lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1183  rab, sabb, dabb)
1184  INTEGER, INTENT(IN) :: la_max, la_min, npgfa
1185  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
1186  INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1187  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1188  INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1189  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1190  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1191  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
1192  OPTIONAL :: sabb
1193  REAL(kind=dp), DIMENSION(:, :, :, :), &
1194  INTENT(INOUT), OPTIONAL :: dabb
1195 
1196  INTEGER :: ax, ay, az, bx, bx1, bx2, by, by1, by2, bz, bz1, bz2, coa, cob1, cob2, ia, ib1, &
1197  ib2, ipgf, j1pgf, j2pgf, la, lb1, lb2, ldrr, lma, lmb, ma, mb1, mb2, na, nb1, nb2, ofa, &
1198  ofb1, ofb2
1199  REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfb, &
1200  tab, xhi, zet
1201  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1202  REAL(kind=dp), DIMENSION(3) :: rap, rbp
1203 
1204  ! Distance of the centers a and b
1205 
1206  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1207  tab = sqrt(rab2)
1208 
1209  ! Maximum l for auxiliary integrals
1210  IF (PRESENT(sabb)) THEN
1211  lma = la_max
1212  lmb = lb1_max + lb2_max
1213  END IF
1214  IF (PRESENT(dabb)) THEN
1215  lma = la_max + 1
1216  lmb = lb1_max + lb2_max
1217  END IF
1218  ldrr = max(lma, lmb) + 1
1219 
1220  ! Allocate space for auxiliary integrals
1221  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1222 
1223  ! Number of integrals, check size of arrays
1224  ofa = ncoset(la_min - 1)
1225  ofb1 = ncoset(lb1_min - 1)
1226  ofb2 = ncoset(lb2_min - 1)
1227  na = ncoset(la_max) - ofa
1228  nb1 = ncoset(lb1_max) - ofb1
1229  nb2 = ncoset(lb2_max) - ofb2
1230  IF (PRESENT(sabb)) THEN
1231  cpassert((SIZE(sabb, 1) >= na*npgfa))
1232  cpassert((SIZE(sabb, 2) >= nb1*npgfb1))
1233  cpassert((SIZE(sabb, 3) >= nb2*npgfb2))
1234  END IF
1235  IF (PRESENT(dabb)) THEN
1236  cpassert((SIZE(dabb, 1) >= na*npgfa))
1237  cpassert((SIZE(dabb, 2) >= nb1*npgfb1))
1238  cpassert((SIZE(dabb, 3) >= nb2*npgfb2))
1239  cpassert((SIZE(dabb, 4) >= 3))
1240  END IF
1241 
1242  ! Loops over all pairs of primitive Gaussian-type functions
1243  ma = 0
1244  DO ipgf = 1, npgfa
1245  mb1 = 0
1246  DO j1pgf = 1, npgfb1
1247  mb2 = 0
1248  DO j2pgf = 1, npgfb2
1249  ! Distance Screening
1250  rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
1251  IF (rpgfa(ipgf) + rpgfb < tab) THEN
1252  IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
1253  IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
1254  mb2 = mb2 + nb2
1255  cycle
1256  END IF
1257 
1258  ! Calculate some prefactors
1259  a = zeta(ipgf)
1260  b = zetb1(j1pgf) + zetb2(j2pgf)
1261  zet = a + b
1262  xhi = a*b/zet
1263  rap = b*rab/zet
1264  rbp = -a*rab/zet
1265 
1266  ! [s|s] integral
1267  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1268 
1269  ! Calculate the recurrence relation
1270  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1271 
1272  DO lb2 = lb2_min, lb2_max
1273  DO bx2 = 0, lb2
1274  DO by2 = 0, lb2 - bx2
1275  bz2 = lb2 - bx2 - by2
1276  cob2 = coset(bx2, by2, bz2) - ofb2
1277  ib2 = mb2 + cob2
1278  DO lb1 = lb1_min, lb1_max
1279  DO bx1 = 0, lb1
1280  DO by1 = 0, lb1 - bx1
1281  bz1 = lb1 - bx1 - by1
1282  cob1 = coset(bx1, by1, bz1) - ofb1
1283  ib1 = mb1 + cob1
1284  DO la = la_min, la_max
1285  DO ax = 0, la
1286  DO ay = 0, la - ax
1287  az = la - ax - ay
1288  coa = coset(ax, ay, az) - ofa
1289  ia = ma + coa
1290  ! integrals
1291  IF (PRESENT(sabb)) THEN
1292  sabb(ia, ib1, ib2) = f0*rr(ax, bx1 + bx2, 1)*rr(ay, by1 + by2, 2)*rr(az, bz1 + bz2, 3)
1293  END IF
1294  ! first derivatives
1295  IF (PRESENT(dabb)) THEN
1296  bx = bx1 + bx2
1297  by = by1 + by2
1298  bz = bz1 + bz2
1299  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1300  ! dx
1301  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1302  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1303  dabb(ia, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1304  ! dy
1305  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1306  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1307  dabb(ia, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1308  ! dz
1309  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1310  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1311  dabb(ia, ib1, ib2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1312  END IF
1313  !
1314  END DO
1315  END DO
1316  END DO !la
1317  END DO
1318  END DO
1319  END DO !lb1
1320  END DO
1321  END DO
1322  END DO !lb2
1323 
1324  mb2 = mb2 + nb2
1325  END DO
1326  mb1 = mb1 + nb1
1327  END DO
1328  ma = ma + na
1329  END DO
1330 
1331  DEALLOCATE (rr)
1332 
1333  END SUBROUTINE overlap_abb
1334 
1335 ! **************************************************************************************************
1336 !> \brief Calculation of the two-center overlap integrals [aaa|b] over
1337 !> Cartesian Gaussian-type functions.
1338 !> \param la1_max Max L on center A (basis 1)
1339 !> \param la1_min Min L on center A (basis 1)
1340 !> \param npgfa1 Number of primitives on center A (basis 1)
1341 !> \param rpgfa1 Range of functions on A, used for screening (basis 1)
1342 !> \param zeta1 Exponents on center A (basis 1)
1343 !> \param la2_max Max L on center A (basis 2)
1344 !> \param la2_min Min L on center A (basis 2)
1345 !> \param npgfa2 Number of primitives on center A (basis 2)
1346 !> \param rpgfa2 Range of functions on A, used for screening (basis 2)
1347 !> \param zeta2 Exponents on center A (basis 2)
1348 !> \param la3_max Max L on center A (basis 3)
1349 !> \param la3_min Min L on center A (basis 3)
1350 !> \param npgfa3 Number of primitives on center A (basis 3)
1351 !> \param rpgfa3 Range of functions on A, used for screening (basis 3)
1352 !> \param zeta3 Exponents on center A (basis 3)
1353 !> \param lb_max Max L on center B
1354 !> \param lb_min Min L on center B
1355 !> \param npgfb Number of primitives on center B
1356 !> \param rpgfb Range of functions on B, used for screening
1357 !> \param zetb Exponents on center B
1358 !> \param rab Distance vector A-B
1359 !> \param saaab Final overlap integrals
1360 !> \param daaab First derivative overlap integrals
1361 !> \date 01.07.2014
1362 !> \author JGH
1363 ! **************************************************************************************************
1364  SUBROUTINE overlap_aaab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
1365  la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
1366  la3_max, la3_min, npgfa3, rpgfa3, zeta3, &
1367  lb_max, lb_min, npgfb, rpgfb, zetb, &
1368  rab, saaab, daaab)
1369  INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
1370  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
1371  INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
1372  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
1373  INTEGER, INTENT(IN) :: la3_max, la3_min, npgfa3
1374  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa3, zeta3
1375  INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
1376  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
1377  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1378  REAL(kind=dp), DIMENSION(:, :, :, :), &
1379  INTENT(INOUT), OPTIONAL :: saaab
1380  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1381  INTENT(INOUT), OPTIONAL :: daaab
1382 
1383  INTEGER :: ax, ax1, ax2, ax3, ay, ay1, ay2, ay3, az, az1, az2, az3, bx, by, bz, coa1, coa2, &
1384  coa3, cob, i1pgf, i2pgf, i3pgf, ia1, ia2, ia3, ib, jpgf, la1, la2, la3, lb, ldrr, lma, &
1385  lmb, ma1, ma2, ma3, mb, na1, na2, na3, nb, ofa1, ofa2, ofa3, ofb
1386  REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
1387  tab, xhi, zet
1388  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1389  REAL(kind=dp), DIMENSION(3) :: rap, rbp
1390 
1391 ! Distance of the centers a and b
1392 
1393  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1394  tab = sqrt(rab2)
1395 
1396  ! Maximum l for auxiliary integrals
1397  IF (PRESENT(saaab)) THEN
1398  lma = la1_max + la2_max + la3_max
1399  lmb = lb_max
1400  END IF
1401  IF (PRESENT(daaab)) THEN
1402  lma = la1_max + la2_max + la3_max + 1
1403  lmb = lb_max
1404  END IF
1405  ldrr = max(lma, lmb) + 1
1406 
1407  ! Allocate space for auxiliary integrals
1408  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1409 
1410  ! Number of integrals, check size of arrays
1411  ofa1 = ncoset(la1_min - 1)
1412  ofa2 = ncoset(la2_min - 1)
1413  ofa3 = ncoset(la3_min - 1)
1414  ofb = ncoset(lb_min - 1)
1415  na1 = ncoset(la1_max) - ofa1
1416  na2 = ncoset(la2_max) - ofa2
1417  na3 = ncoset(la3_max) - ofa3
1418  nb = ncoset(lb_max) - ofb
1419  IF (PRESENT(saaab)) THEN
1420  cpassert((SIZE(saaab, 1) >= na1*npgfa1))
1421  cpassert((SIZE(saaab, 2) >= na2*npgfa2))
1422  cpassert((SIZE(saaab, 3) >= na3*npgfa3))
1423  cpassert((SIZE(saaab, 4) >= nb*npgfb))
1424  END IF
1425  IF (PRESENT(daaab)) THEN
1426  cpassert((SIZE(daaab, 1) >= na1*npgfa1))
1427  cpassert((SIZE(daaab, 2) >= na2*npgfa2))
1428  cpassert((SIZE(daaab, 3) >= na3*npgfa3))
1429  cpassert((SIZE(daaab, 4) >= nb*npgfb))
1430  cpassert((SIZE(daaab, 5) >= 3))
1431  END IF
1432 
1433  ! Loops over all primitive Gaussian-type functions
1434  ma1 = 0
1435  DO i1pgf = 1, npgfa1
1436  ma2 = 0
1437  DO i2pgf = 1, npgfa2
1438  ma3 = 0
1439  DO i3pgf = 1, npgfa3
1440  rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf), rpgfa3(i3pgf))
1441  mb = 0
1442  DO jpgf = 1, npgfb
1443  ! Distance Screening
1444  IF (rpgfa + rpgfb(jpgf) < tab) THEN
1445  IF (PRESENT(saaab)) saaab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, ma3 + 1:ma3 + na3, mb + 1:mb + nb) = 0.0_dp
1446  IF (PRESENT(daaab)) daaab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, ma3 + 1:ma3 + na3, mb + 1:mb + nb, 1:3) = 0.0_dp
1447  mb = mb + nb
1448  cycle
1449  END IF
1450 
1451  ! Calculate some prefactors
1452  a = zeta1(i1pgf) + zeta2(i2pgf) + zeta3(i3pgf)
1453  b = zetb(jpgf)
1454  zet = a + b
1455  xhi = a*b/zet
1456  rap = b*rab/zet
1457  rbp = -a*rab/zet
1458 
1459  ! [sss|s] integral
1460  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1461 
1462  ! Calculate the recurrence relation
1463  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1464 
1465  DO lb = lb_min, lb_max
1466  DO bx = 0, lb
1467  DO by = 0, lb - bx
1468  bz = lb - bx - by
1469  cob = coset(bx, by, bz) - ofb
1470  ib = mb + cob
1471  DO la3 = la3_min, la3_max
1472  DO ax3 = 0, la3
1473  DO ay3 = 0, la3 - ax3
1474  az3 = la3 - ax3 - ay3
1475  coa3 = coset(ax3, ay3, az3) - ofa3
1476  ia3 = ma3 + coa3
1477  DO la2 = la2_min, la2_max
1478  DO ax2 = 0, la2
1479  DO ay2 = 0, la2 - ax2
1480  az2 = la2 - ax2 - ay2
1481  coa2 = coset(ax2, ay2, az2) - ofa2
1482  ia2 = ma2 + coa2
1483  DO la1 = la1_min, la1_max
1484  DO ax1 = 0, la1
1485  DO ay1 = 0, la1 - ax1
1486  az1 = la1 - ax1 - ay1
1487  coa1 = coset(ax1, ay1, az1) - ofa1
1488  ia1 = ma1 + coa1
1489  ! integrals
1490  IF (PRESENT(saaab)) THEN
1491  saaab(ia1, ia2, ia3, ib) = f0*rr(ax1 + ax2 + ax3, bx, 1)* &
1492  rr(ay1 + ay2 + ay3, by, 2)*rr(az1 + az2 + az3, bz, 3)
1493  END IF
1494  ! first derivatives
1495  IF (PRESENT(daaab)) THEN
1496  ax = ax1 + ax2 + ax3
1497  ay = ay1 + ay2 + ay3
1498  az = az1 + az2 + az3
1499  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1500  ! dx
1501  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1502  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1503  daaab(ia1, ia2, ia3, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1504  ! dy
1505  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1506  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1507  daaab(ia1, ia2, ia3, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1508  ! dz
1509  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1510  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1511  daaab(ia1, ia2, ia3, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1512  END IF
1513  !
1514  END DO
1515  END DO
1516  END DO !la1
1517  END DO
1518  END DO
1519  END DO !la2
1520  END DO
1521  END DO
1522  END DO !la3
1523  END DO
1524  END DO
1525  END DO !lb
1526 
1527  mb = mb + nb
1528  END DO
1529  ma3 = ma3 + na3
1530  END DO
1531  ma2 = ma2 + na2
1532  END DO
1533  ma1 = ma1 + na1
1534  END DO
1535 
1536  DEALLOCATE (rr)
1537 
1538  END SUBROUTINE overlap_aaab
1539 ! **************************************************************************************************
1540 !> \brief Calculation of the two-center overlap integrals [aa|bb] over
1541 !> Cartesian Gaussian-type functions.
1542 !> \param la1_max Max L on center A (basis 1)
1543 !> \param la1_min Min L on center A (basis 1)
1544 !> \param npgfa1 Number of primitives on center A (basis 1)
1545 !> \param rpgfa1 Range of functions on A, used for screening (basis 1)
1546 !> \param zeta1 Exponents on center A (basis 1)
1547 !> \param la2_max Max L on center A (basis 2)
1548 !> \param la2_min Min L on center A (basis 2)
1549 !> \param npgfa2 Number of primitives on center A (basis 2)
1550 !> \param rpgfa2 Range of functions on A, used for screening (basis 2)
1551 !> \param zeta2 Exponents on center A (basis 2)
1552 !> \param lb1_max Max L on center B (basis 3)
1553 !> \param lb1_min Min L on center B (basis 3)
1554 !> \param npgfb1 Number of primitives on center B (basis 3)
1555 !> \param rpgfb1 Range of functions on B, used for screening (basis 3)
1556 !> \param zetb1 Exponents on center B (basis 3)
1557 !> \param lb2_max Max L on center B (basis 4)
1558 !> \param lb2_min Min L on center B (basis 4)
1559 !> \param npgfb2 Number of primitives on center B (basis 4)
1560 !> \param rpgfb2 Range of functions on B, used for screening (basis 4)
1561 !> \param zetb2 Exponents on center B (basis 4)
1562 !> \param rab Distance vector A-B
1563 !> \param saabb Final overlap integrals
1564 !> \param daabb First derivative overlap integrals
1565 !> \date 01.07.2014
1566 !> \author JGH
1567 ! **************************************************************************************************
1568  SUBROUTINE overlap_aabb(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
1569  la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
1570  lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1571  lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1572  rab, saabb, daabb)
1573  INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
1574  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
1575  INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
1576  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
1577  INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1578  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1579  INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1580  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1581  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1582  REAL(kind=dp), DIMENSION(:, :, :, :), &
1583  INTENT(INOUT), OPTIONAL :: saabb
1584  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1585  INTENT(INOUT), OPTIONAL :: daabb
1586 
1587  INTEGER :: ax, ax1, ax2, ay, ay1, ay2, az, az1, az2, bx, bx1, bx2, by, by1, by2, bz, bz1, &
1588  bz2, coa1, coa2, cob1, cob2, i1pgf, i2pgf, ia1, ia2, ib1, ib2, j1pgf, j2pgf, la1, la2, &
1589  lb1, lb2, ldrr, lma, lmb, ma1, ma2, mb1, mb2, na1, na2, nb1, nb2, ofa1, ofa2, ofb1, ofb2
1590  REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
1591  rpgfb, tab, xhi, zet
1592  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1593  REAL(kind=dp), DIMENSION(3) :: rap, rbp
1594 
1595 ! Distance of the centers a and b
1596 
1597  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1598  tab = sqrt(rab2)
1599 
1600  ! Maximum l for auxiliary integrals
1601  IF (PRESENT(saabb)) THEN
1602  lma = la1_max + la2_max
1603  lmb = lb1_max + lb2_max
1604  END IF
1605  IF (PRESENT(daabb)) THEN
1606  lma = la1_max + la2_max + 1
1607  lmb = lb1_max + lb2_max
1608  END IF
1609  ldrr = max(lma, lmb) + 1
1610 
1611  ! Allocate space for auxiliary integrals
1612  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1613 
1614  ! Number of integrals, check size of arrays
1615  ofa1 = ncoset(la1_min - 1)
1616  ofa2 = ncoset(la2_min - 1)
1617  ofb1 = ncoset(lb1_min - 1)
1618  ofb2 = ncoset(lb2_min - 1)
1619  na1 = ncoset(la1_max) - ofa1
1620  na2 = ncoset(la2_max) - ofa2
1621  nb1 = ncoset(lb1_max) - ofb1
1622  nb2 = ncoset(lb2_max) - ofb2
1623  IF (PRESENT(saabb)) THEN
1624  cpassert((SIZE(saabb, 1) >= na1*npgfa1))
1625  cpassert((SIZE(saabb, 2) >= na2*npgfa2))
1626  cpassert((SIZE(saabb, 3) >= nb1*npgfb1))
1627  cpassert((SIZE(saabb, 4) >= nb2*npgfb2))
1628  END IF
1629  IF (PRESENT(daabb)) THEN
1630  cpassert((SIZE(daabb, 1) >= na1*npgfa1))
1631  cpassert((SIZE(daabb, 2) >= na2*npgfa2))
1632  cpassert((SIZE(daabb, 3) >= nb1*npgfb1))
1633  cpassert((SIZE(daabb, 4) >= nb2*npgfb2))
1634  cpassert((SIZE(daabb, 5) >= 3))
1635  END IF
1636 
1637  ! Loops over all primitive Gaussian-type functions
1638  ma1 = 0
1639  DO i1pgf = 1, npgfa1
1640  ma2 = 0
1641  DO i2pgf = 1, npgfa2
1642  mb1 = 0
1643  rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
1644  DO j1pgf = 1, npgfb1
1645  mb2 = 0
1646  DO j2pgf = 1, npgfb2
1647  rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
1648  ! Distance Screening
1649  IF (rpgfa + rpgfb < tab) THEN
1650  IF (PRESENT(saabb)) saabb(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
1651  IF (PRESENT(daabb)) daabb(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
1652  mb2 = mb2 + nb2
1653  cycle
1654  END IF
1655 
1656  ! Calculate some prefactors
1657  a = zeta1(i1pgf) + zeta2(i2pgf)
1658  b = zetb1(j1pgf) + zetb2(j2pgf)
1659  zet = a + b
1660  xhi = a*b/zet
1661  rap = b*rab/zet
1662  rbp = -a*rab/zet
1663 
1664  ! [ss|ss] integral
1665  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1666 
1667  ! Calculate the recurrence relation
1668  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1669 
1670  DO lb2 = lb2_min, lb2_max
1671  DO bx2 = 0, lb2
1672  DO by2 = 0, lb2 - bx2
1673  bz2 = lb2 - bx2 - by2
1674  cob2 = coset(bx2, by2, bz2) - ofb2
1675  ib2 = mb2 + cob2
1676  DO lb1 = lb1_min, lb1_max
1677  DO bx1 = 0, lb1
1678  DO by1 = 0, lb1 - bx1
1679  bz1 = lb1 - bx1 - by1
1680  cob1 = coset(bx1, by1, bz1) - ofb1
1681  ib1 = mb1 + cob1
1682  DO la2 = la2_min, la2_max
1683  DO ax2 = 0, la2
1684  DO ay2 = 0, la2 - ax2
1685  az2 = la2 - ax2 - ay2
1686  coa2 = coset(ax2, ay2, az2) - ofa2
1687  ia2 = ma2 + coa2
1688  DO la1 = la1_min, la1_max
1689  DO ax1 = 0, la1
1690  DO ay1 = 0, la1 - ax1
1691  az1 = la1 - ax1 - ay1
1692  coa1 = coset(ax1, ay1, az1) - ofa1
1693  ia1 = ma1 + coa1
1694  ! integrals
1695  IF (PRESENT(saabb)) THEN
1696  saabb(ia1, ia2, ib1, ib2) = f0*rr(ax1 + ax2, bx1 + bx2, 1)* &
1697  rr(ay1 + ay2, by1 + by2, 2)*rr(az1 + az2, bz1 + bz2, 3)
1698  END IF
1699  ! first derivatives
1700  IF (PRESENT(daabb)) THEN
1701  ax = ax1 + ax2
1702  ay = ay1 + ay2
1703  az = az1 + az2
1704  bx = bx1 + bx2
1705  by = by1 + by2
1706  bz = bz1 + bz2
1707  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1708  ! dx
1709  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1710  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1711  daabb(ia1, ia2, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1712  ! dy
1713  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1714  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1715  daabb(ia1, ia2, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1716  ! dz
1717  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1718  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1719  daabb(ia1, ia2, ib1, ib2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1720  END IF
1721  !
1722  END DO
1723  END DO
1724  END DO !la1
1725  END DO
1726  END DO
1727  END DO !la2
1728  END DO
1729  END DO
1730  END DO !lb1
1731  END DO
1732  END DO
1733  END DO !lb2
1734 
1735  mb2 = mb2 + nb2
1736  END DO
1737  mb1 = mb1 + nb1
1738  END DO
1739  ma2 = ma2 + na2
1740  END DO
1741  ma1 = ma1 + na1
1742  END DO
1743 
1744  DEALLOCATE (rr)
1745 
1746  END SUBROUTINE overlap_aabb
1747 ! **************************************************************************************************
1748 !> \brief Calculation of the two-center overlap integrals [a|bbb] over
1749 !> Cartesian Gaussian-type functions.
1750 !> \param la_max Max L on center A
1751 !> \param la_min Min L on center A
1752 !> \param npgfa Number of primitives on center A
1753 !> \param rpgfa Range of functions on A, used for screening
1754 !> \param zeta Exponents on center A
1755 !> \param lb1_max Max L on center B (basis 1)
1756 !> \param lb1_min Min L on center B (basis 1)
1757 !> \param npgfb1 Number of primitives on center B (basis 1)
1758 !> \param rpgfb1 Range of functions on B, used for screening (basis 1)
1759 !> \param zetb1 Exponents on center B (basis 1)
1760 !> \param lb2_max Max L on center B (basis 2)
1761 !> \param lb2_min Min L on center B (basis 2)
1762 !> \param npgfb2 Number of primitives on center B (basis 2)
1763 !> \param rpgfb2 Range of functions on B, used for screening (basis 2)
1764 !> \param zetb2 Exponents on center B (basis 2)
1765 !> \param lb3_max Max L on center B (basis 3)
1766 !> \param lb3_min Min L on center B (basis 3)
1767 !> \param npgfb3 Number of primitives on center B (basis 3)
1768 !> \param rpgfb3 Range of functions on B, used for screening (basis 3)
1769 !> \param zetb3 Exponents on center B (basis 3)
1770 !> \param rab Distance vector A-B
1771 !> \param sabbb Final overlap integrals
1772 !> \param dabbb First derivative overlap integrals
1773 !> \date 01.07.2014
1774 !> \author JGH
1775 ! **************************************************************************************************
1776  SUBROUTINE overlap_abbb(la_max, la_min, npgfa, rpgfa, zeta, &
1777  lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
1778  lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
1779  lb3_max, lb3_min, npgfb3, rpgfb3, zetb3, &
1780  rab, sabbb, dabbb)
1781  INTEGER, INTENT(IN) :: la_max, la_min, npgfa
1782  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
1783  INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
1784  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
1785  INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
1786  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
1787  INTEGER, INTENT(IN) :: lb3_max, lb3_min, npgfb3
1788  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb3, zetb3
1789  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1790  REAL(kind=dp), DIMENSION(:, :, :, :), &
1791  INTENT(INOUT), OPTIONAL :: sabbb
1792  REAL(kind=dp), DIMENSION(:, :, :, :, :), &
1793  INTENT(INOUT), OPTIONAL :: dabbb
1794 
1795  INTEGER :: ax, ay, az, bx, bx1, bx2, bx3, by, by1, by2, by3, bz, bz1, bz2, bz3, coa, cob1, &
1796  cob2, cob3, ia, ib1, ib2, ib3, ipgf, j1pgf, j2pgf, j3pgf, la, lb1, lb2, lb3, ldrr, lma, &
1797  lmb, ma, mb1, mb2, mb3, na, nb1, nb2, nb3, ofa, ofb1, ofb2, ofb3
1798  REAL(kind=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfb, &
1799  tab, xhi, zet
1800  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
1801  REAL(kind=dp), DIMENSION(3) :: rap, rbp
1802 
1803 ! Distance of the centers a and b
1804 
1805  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1806  tab = sqrt(rab2)
1807 
1808  ! Maximum l for auxiliary integrals
1809  IF (PRESENT(sabbb)) THEN
1810  lma = la_max
1811  lmb = lb1_max + lb2_max + lb3_max
1812  END IF
1813  IF (PRESENT(dabbb)) THEN
1814  lma = la_max + 1
1815  lmb = lb1_max + lb2_max + lb3_max
1816  END IF
1817  ldrr = max(lma, lmb) + 1
1818 
1819  ! Allocate space for auxiliary integrals
1820  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
1821 
1822  ! Number of integrals, check size of arrays
1823  ofa = ncoset(la_min - 1)
1824  ofb1 = ncoset(lb1_min - 1)
1825  ofb2 = ncoset(lb2_min - 1)
1826  ofb3 = ncoset(lb3_min - 1)
1827  na = ncoset(la_max) - ofa
1828  nb1 = ncoset(lb1_max) - ofb1
1829  nb2 = ncoset(lb2_max) - ofb2
1830  nb3 = ncoset(lb3_max) - ofb3
1831  IF (PRESENT(sabbb)) THEN
1832  cpassert((SIZE(sabbb, 1) >= na*npgfa))
1833  cpassert((SIZE(sabbb, 2) >= nb1*npgfb1))
1834  cpassert((SIZE(sabbb, 3) >= nb2*npgfb2))
1835  cpassert((SIZE(sabbb, 4) >= nb3*npgfb3))
1836  END IF
1837  IF (PRESENT(dabbb)) THEN
1838  cpassert((SIZE(dabbb, 1) >= na*npgfa))
1839  cpassert((SIZE(dabbb, 2) >= nb1*npgfb1))
1840  cpassert((SIZE(dabbb, 3) >= nb2*npgfb2))
1841  cpassert((SIZE(dabbb, 4) >= nb3*npgfb3))
1842  cpassert((SIZE(dabbb, 5) >= 3))
1843  END IF
1844 
1845  ! Loops over all pairs of primitive Gaussian-type functions
1846  ma = 0
1847  DO ipgf = 1, npgfa
1848  mb1 = 0
1849  DO j1pgf = 1, npgfb1
1850  mb2 = 0
1851  DO j2pgf = 1, npgfb2
1852  mb3 = 0
1853  DO j3pgf = 1, npgfb3
1854  ! Distance Screening
1855  rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf), rpgfb3(j3pgf))
1856  IF (rpgfa(ipgf) + rpgfb < tab) THEN
1857  IF (PRESENT(sabbb)) sabbb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, mb3 + 1:mb3 + nb3) = 0.0_dp
1858  IF (PRESENT(dabbb)) dabbb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, mb3 + 1:mb3 + nb3, 1:3) = 0.0_dp
1859  mb3 = mb3 + nb3
1860  cycle
1861  END IF
1862 
1863  ! Calculate some prefactors
1864  a = zeta(ipgf)
1865  b = zetb1(j1pgf) + zetb2(j2pgf) + zetb3(j3pgf)
1866  zet = a + b
1867  xhi = a*b/zet
1868  rap = b*rab/zet
1869  rbp = -a*rab/zet
1870 
1871  ! [s|s] integral
1872  f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
1873 
1874  ! Calculate the recurrence relation
1875  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
1876 
1877  DO lb3 = lb3_min, lb3_max
1878  DO bx3 = 0, lb3
1879  DO by3 = 0, lb3 - bx3
1880  bz3 = lb3 - bx3 - by3
1881  cob3 = coset(bx3, by3, bz3) - ofb3
1882  ib3 = mb3 + cob3
1883  DO lb2 = lb2_min, lb2_max
1884  DO bx2 = 0, lb2
1885  DO by2 = 0, lb2 - bx2
1886  bz2 = lb2 - bx2 - by2
1887  cob2 = coset(bx2, by2, bz2) - ofb2
1888  ib2 = mb2 + cob2
1889  DO lb1 = lb1_min, lb1_max
1890  DO bx1 = 0, lb1
1891  DO by1 = 0, lb1 - bx1
1892  bz1 = lb1 - bx1 - by1
1893  cob1 = coset(bx1, by1, bz1) - ofb1
1894  ib1 = mb1 + cob1
1895  DO la = la_min, la_max
1896  DO ax = 0, la
1897  DO ay = 0, la - ax
1898  az = la - ax - ay
1899  coa = coset(ax, ay, az) - ofa
1900  ia = ma + coa
1901  ! integrals
1902  IF (PRESENT(sabbb)) THEN
1903  sabbb(ia, ib1, ib2, ib3) = f0*rr(ax, bx1 + bx2 + bx3, 1)* &
1904  rr(ay, by1 + by2 + by3, 2)*rr(az, bz1 + bz2 + bz3, 3)
1905  END IF
1906  ! first derivatives
1907  IF (PRESENT(dabbb)) THEN
1908  bx = bx1 + bx2 + bx3
1909  by = by1 + by2 + by3
1910  bz = bz1 + bz2 + bz3
1911  ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
1912  ! dx
1913  dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
1914  IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
1915  dabbb(ia, ib1, ib2, ib3, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
1916  ! dy
1917  dumy = 2.0_dp*a*rr(ay + 1, by, 2)
1918  IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
1919  dabbb(ia, ib1, ib2, ib3, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
1920  ! dz
1921  dumz = 2.0_dp*a*rr(az + 1, bz, 3)
1922  IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
1923  dabbb(ia, ib1, ib2, ib3, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
1924  END IF
1925  !
1926  END DO
1927  END DO
1928  END DO !la
1929  END DO
1930  END DO
1931  END DO !lb1
1932  END DO
1933  END DO
1934  END DO !lb2
1935  END DO
1936  END DO
1937  END DO !lb3
1938 
1939  mb3 = mb3 + nb3
1940  END DO
1941  mb2 = mb2 + nb2
1942  END DO
1943  mb1 = mb1 + nb1
1944  END DO
1945  ma = ma + na
1946  END DO
1947 
1948  DEALLOCATE (rr)
1949 
1950  END SUBROUTINE overlap_abbb
1951 
1952 ! **************************************************************************************************
1953 !> \brief Calculation of the two-center overlap integrals [a|b] over
1954 !> Spherical Gaussian-type functions.
1955 !> \param la Max L on center A
1956 !> \param zeta Exponents on center A
1957 !> \param lb Max L on center B
1958 !> \param zetb Exponents on center B
1959 !> \param rab Distance vector A-B
1960 !> \param sab Final overlap integrals
1961 !> \date 01.03.2016
1962 !> \author JGH
1963 ! **************************************************************************************************
1964  SUBROUTINE overlap_ab_s(la, zeta, lb, zetb, rab, sab)
1965  INTEGER, INTENT(IN) :: la
1966  REAL(kind=dp), INTENT(IN) :: zeta
1967  INTEGER, INTENT(IN) :: lb
1968  REAL(kind=dp), INTENT(IN) :: zetb
1969  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1970  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
1971 
1972  REAL(kind=dp), PARAMETER :: huge4 = huge(1._dp)/4._dp
1973 
1974  INTEGER :: nca, ncb, nsa, nsb
1975  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
1976  REAL(kind=dp), DIMENSION(1) :: rpgf, za, zb
1977  REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
1978 
1979  rpgf(1) = huge4
1980  za(1) = zeta
1981  zb(1) = zetb
1982 
1983  nca = nco(la)
1984  ncb = nco(lb)
1985  ALLOCATE (cab(nca, ncb))
1986  nsa = nso(la)
1987  nsb = nso(lb)
1988 
1989  CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
1990 
1991  c2sa => orbtramat(la)%c2s
1992  c2sb => orbtramat(lb)%c2s
1993  sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
1994  matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
1995 
1996  DEALLOCATE (cab)
1997 
1998  END SUBROUTINE overlap_ab_s
1999 
2000 ! **************************************************************************************************
2001 !> \brief Calculation of the overlap integrals [a|b] over
2002 !> cubic periodic Spherical Gaussian-type functions.
2003 !> \param la Max L on center A
2004 !> \param zeta Exponents on center A
2005 !> \param lb Max L on center B
2006 !> \param zetb Exponents on center B
2007 !> \param alat Lattice constant
2008 !> \param sab Final overlap integrals
2009 !> \date 01.03.2016
2010 !> \author JGH
2011 ! **************************************************************************************************
2012  SUBROUTINE overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
2013  INTEGER, INTENT(IN) :: la
2014  REAL(kind=dp), INTENT(IN) :: zeta
2015  INTEGER, INTENT(IN) :: lb
2016  REAL(kind=dp), INTENT(IN) :: zetb, alat
2017  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
2018 
2019  COMPLEX(KIND=dp) :: zfg
2020  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fun, gun
2021  INTEGER :: ax, ay, az, bx, by, bz, i, ia, ib, l, &
2022  l1, l2, na, nb, nca, ncb, nmax, nsa, &
2023  nsb
2024  REAL(kind=dp) :: oa, ob, ovol, zm
2025  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
2026  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
2027  REAL(kind=dp), DIMENSION(0:3, 0:3) :: fgsum
2028  REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
2029 
2030  nca = nco(la)
2031  ncb = nco(lb)
2032  ALLOCATE (cab(nca, ncb))
2033  cab = 0.0_dp
2034  nsa = nso(la)
2035  nsb = nso(lb)
2036 
2037  zm = min(zeta, zetb)
2038  nmax = nint(1.81_dp*alat*sqrt(zm) + 1.0_dp)
2039  ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
2040  fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
2041 
2042  oa = 1._dp/zeta
2043  ob = 1._dp/zetb
2044  DO i = -nmax, nmax
2045  gval(i) = twopi/alat*real(i, kind=dp)
2046  fexp(i) = sqrt(oa*pi)*exp(-0.25_dp*oa*gval(i)**2)
2047  gexp(i) = sqrt(ob*pi)*exp(-0.25_dp*ob*gval(i)**2)
2048  END DO
2049  DO l = 0, la
2050  IF (l == 0) THEN
2051  fun(:, l) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2052  ELSEIF (l == 1) THEN
2053  fun(:, l) = cmplx(0.0_dp, 0.5_dp*oa*gval(:), kind=dp)
2054  ELSEIF (l == 2) THEN
2055  fun(:, l) = cmplx(-(0.5_dp*oa*gval(:))**2, 0.0_dp, kind=dp)
2056  fun(:, l) = fun(:, l) + cmplx(0.5_dp*oa, 0.0_dp, kind=dp)
2057  ELSEIF (l == 3) THEN
2058  fun(:, l) = cmplx(0.0_dp, -(0.5_dp*oa*gval(:))**3, kind=dp)
2059  fun(:, l) = fun(:, l) + cmplx(0.0_dp, 0.75_dp*oa*oa*gval(:), kind=dp)
2060  ELSE
2061  cpabort("l value too high")
2062  END IF
2063  END DO
2064  DO l = 0, lb
2065  IF (l == 0) THEN
2066  gun(:, l) = cmplx(1.0_dp, 0.0_dp, kind=dp)
2067  ELSEIF (l == 1) THEN
2068  gun(:, l) = cmplx(0.0_dp, 0.5_dp*ob*gval(:), kind=dp)
2069  ELSEIF (l == 2) THEN
2070  gun(:, l) = cmplx(-(0.5_dp*ob*gval(:))**2, 0.0_dp, kind=dp)
2071  gun(:, l) = gun(:, l) + cmplx(0.5_dp*ob, 0.0_dp, kind=dp)
2072  ELSEIF (l == 3) THEN
2073  gun(:, l) = cmplx(0.0_dp, -(0.5_dp*ob*gval(:))**3, kind=dp)
2074  gun(:, l) = gun(:, l) + cmplx(0.0_dp, 0.75_dp*ob*ob*gval(:), kind=dp)
2075  ELSE
2076  cpabort("l value too high")
2077  END IF
2078  END DO
2079 
2080  fgsum = 0.0_dp
2081  DO l1 = 0, la
2082  DO l2 = 0, lb
2083  zfg = sum(conjg(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
2084  fgsum(l1, l2) = real(zfg, kind=dp)
2085  END DO
2086  END DO
2087 
2088  na = ncoset(la - 1)
2089  nb = ncoset(lb - 1)
2090  DO ax = 0, la
2091  DO ay = 0, la - ax
2092  az = la - ax - ay
2093  ia = coset(ax, ay, az) - na
2094  DO bx = 0, lb
2095  DO by = 0, lb - bx
2096  bz = lb - bx - by
2097  ib = coset(bx, by, bz) - nb
2098  cab(ia, ib) = fgsum(ax, bx)*fgsum(ay, by)*fgsum(az, bz)
2099  END DO
2100  END DO
2101  END DO
2102  END DO
2103 
2104  c2sa => orbtramat(la)%c2s
2105  c2sb => orbtramat(lb)%c2s
2106  sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
2107  matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
2108  ovol = 1._dp/(alat**3)
2109  sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
2110 
2111  DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
2112 
2113  END SUBROUTINE overlap_ab_sp
2114 
2115 END MODULE ai_overlap
ai_os_rr
Definition: ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition: ai_os_rr.F:39
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition: ai_overlap.F:680
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition: ai_overlap.F:73
ai_overlap::overlap_ab_s
subroutine, public overlap_ab_s(la, zeta, lb, zetb, rab, sab)
Calculation of the two-center overlap integrals [a|b] over Spherical Gaussian-type functions.
Definition: ai_overlap.F:1965
ai_overlap::overlap_ab_sp
subroutine, public overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
Calculation of the overlap integrals [a|b] over cubic periodic Spherical Gaussian-type functions.
Definition: ai_overlap.F:2013
ai_overlap::overlap_aab
subroutine, public overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, la2_max, la2_min, npgfa2, rpgfa2, zeta2, lb_max, lb_min, npgfb, rpgfb, zetb, rab, saab, daab, saba, daba)
Calculation of the two-center overlap integrals [aa|b] over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:967
ai_overlap::overlap_abb
subroutine, public overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, rab, sabb, dabb)
Calculation of the two-center overlap integrals [a|bb] over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:1184
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
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
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
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition: orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition: orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition: orbital_transformation_matrices.F:47