(git:b279b6b)
ai_overlap_aabb.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 Parameters
14 !> - ax,ay,az : Angular momentum index numbers of orbital a.
15 !> - bx,by,bz : Angular momentum index numbers of orbital b.
16 !> - coset : Cartesian orbital set pointer.
17 !> - dab : Distance between the atomic centers a and b.
18 !> - l{a,b} : Angular momentum quantum number of shell a or b.
19 !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
20 !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
21 !> - rab : Distance vector between the atomic centers a and b.
22 !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
23 !> - sab : Shell set of overlap integrals.
24 !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
25 !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
26 ! **************************************************************************************************
27 MODULE ai_overlap_aabb
28 
29  USE kinds, ONLY: dp
30  USE mathconstants, ONLY: pi
31  USE orbital_pointers, ONLY: coset,&
32  indco,&
33  ncoset
34 #include "../base/base_uses.f90"
35 
36  IMPLICIT NONE
37 
38  PRIVATE
39 
40  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_aabb'
41 
42 ! *** Public subroutines ***
43  PUBLIC :: overlap_aabb
44 
45 CONTAINS
46 
47 ! **************************************************************************************************
48 !> \brief Purpose: Calculation of the two-center overlap integrals [aa|bb]
49 !> over Cartesian Gaussian-type functions.
50 !> \param la_max_set1 ...
51 !> \param la_min_set1 ...
52 !> \param npgfa1 ...
53 !> \param rpgfa1 ...
54 !> \param zeta1 ...
55 !> \param la_max_set2 ...
56 !> \param la_min_set2 ...
57 !> \param npgfa2 ...
58 !> \param rpgfa2 ...
59 !> \param zeta2 ...
60 !> \param lb_max_set1 ...
61 !> \param lb_min_set1 ...
62 !> \param npgfb1 ...
63 !> \param rpgfb1 ...
64 !> \param zetb1 ...
65 !> \param lb_max_set2 ...
66 !> \param lb_min_set2 ...
67 !> \param npgfb2 ...
68 !> \param rpgfb2 ...
69 !> \param zetb2 ...
70 !> \param asets_equal ...
71 !> \param bsets_equal ...
72 !> \param rab ...
73 !> \param dab ...
74 !> \param saabb ...
75 !> \param s ...
76 !> \param lds ...
77 !> \date 06.2014
78 !> \author Dorothea Golze
79 ! **************************************************************************************************
80  SUBROUTINE overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, &
81  la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, &
82  lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, &
83  lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, &
84  asets_equal, bsets_equal, rab, dab, saabb, s, lds)
85 
86  INTEGER, INTENT(IN) :: la_max_set1, la_min_set1, npgfa1
87  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
88  INTEGER, INTENT(IN) :: la_max_set2, la_min_set2, npgfa2
89  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
90  INTEGER, INTENT(IN) :: lb_max_set1, lb_min_set1, npgfb1
91  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
92  INTEGER, INTENT(IN) :: lb_max_set2, lb_min_set2, npgfb2
93  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
94  LOGICAL, INTENT(IN) :: asets_equal, bsets_equal
95  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
96  REAL(kind=dp), INTENT(IN) :: dab
97  REAL(kind=dp), DIMENSION(:, :, :, :), &
98  INTENT(INOUT) :: saabb
99  INTEGER, INTENT(IN) :: lds
100  REAL(kind=dp), DIMENSION(lds, lds), INTENT(INOUT) :: s
101 
102  CHARACTER(len=*), PARAMETER :: routinen = 'overlap_aabb'
103 
104  INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, coamy, coamz, coapx, coapy, coapz, cob, &
105  cobm2x, cobm2y, cobm2z, cobmx, cobmy, cobmz, handle, i, ia, ib, ipgf, j, ja, jb, jpgf, &
106  jpgf_start, kpgf, la, la_max, la_min, la_start, lb, lb_max, lb_min, lpgf, lpgf_start, &
107  ncoa1, ncoa2, ncob1, ncob2
108  INTEGER, DIMENSION(3) :: na, naa, nb, nbb, nia, nib, nja, njb
109  REAL(kind=dp) :: f0, f1, f2, f3, f4, fax, fay, faz, zeta, &
110  zetb, zetp
111  REAL(kind=dp), DIMENSION(3) :: rap, rbp
112 
113  CALL timeset(routinen, handle)
114 
115 ! *** Loop over all pairs of primitive Gaussian-type functions ***
116 
117  ncoa1 = 0
118  ncoa2 = 0
119  ncob1 = 0
120  ncob2 = 0
121 
122  DO ipgf = 1, npgfa1
123 
124  ncoa2 = 0
125 
126  IF (asets_equal) THEN
127  jpgf_start = ipgf
128  DO i = 1, jpgf_start - 1
129  ncoa2 = ncoa2 + ncoset(la_max_set2)
130  END DO
131  ELSE
132  jpgf_start = 1
133  END IF
134 
135  DO jpgf = jpgf_start, npgfa2
136 
137  ncob1 = 0
138  zeta = zeta1(ipgf) + zeta2(jpgf)
139  la_max = la_max_set1 + la_max_set2
140  la_min = la_min_set1 + la_min_set2
141 
142  DO kpgf = 1, npgfb1
143 
144  ncob2 = 0
145 
146  IF (bsets_equal) THEN
147  lpgf_start = kpgf
148  DO i = 1, lpgf_start - 1
149  ncob2 = ncob2 + ncoset(lb_max_set2)
150  END DO
151  ELSE
152  lpgf_start = 1
153  END IF
154 
155  DO lpgf = lpgf_start, npgfb2
156 
157  ! *** Screening ***
158  IF ((rpgfa1(ipgf) + rpgfb1(kpgf) < dab) .OR. &
159  (rpgfa2(jpgf) + rpgfb1(kpgf) < dab) .OR. &
160  (rpgfa1(ipgf) + rpgfb2(lpgf) < dab) .OR. &
161  (rpgfa2(jpgf) + rpgfb2(lpgf) < dab)) THEN
162  DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
163  DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
164  DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
165  DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
166  saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = 0._dp
167  IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = 0._dp
168  IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = 0._dp
169  IF (asets_equal .AND. bsets_equal) THEN
170  saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = 0._dp
171  saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = 0._dp
172  saabb(ncoa2 + ja, ncoa1 + ia, ncob2 + jb, ncob1 + ib) = 0._dp
173  END IF
174  END DO
175  END DO
176  END DO
177  END DO
178  ncob2 = ncob2 + ncoset(lb_max_set2)
179  cycle
180  END IF
181 
182  zetb = zetb1(kpgf) + zetb2(lpgf)
183  lb_max = lb_max_set1 + lb_max_set2
184  lb_min = lb_min_set1 + lb_min_set2
185 
186 ! *** Calculate some prefactors ***
187 
188  zetp = 1.0_dp/(zeta + zetb)
189 
190  f0 = sqrt((pi*zetp)**3)
191  f1 = zetb*zetp
192  f2 = 0.5_dp*zetp
193 
194 ! *** Calculate the basic two-center overlap integral [s|s] ***
195 
196  s(1, 1) = f0*exp(-zeta*f1*dab*dab)
197 
198 ! *** Recurrence steps: [s|s] -> [a|b] ***
199 
200  IF (la_max > 0) THEN
201 
202 ! *** Vertical recurrence steps: [s|s] -> [a|s] ***
203 
204  rap(:) = f1*rab(:)
205 
206 ! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) ***
207 
208  s(2, 1) = rap(1)*s(1, 1) ! [px|s]
209  s(3, 1) = rap(2)*s(1, 1) ! [py|s]
210  s(4, 1) = rap(3)*s(1, 1) ! [pz|s]
211 
212  IF (la_max > 1) THEN
213 
214 ! *** [d|s] ***
215 
216  f3 = f2*s(1, 1)
217 
218  s(5, 1) = rap(1)*s(2, 1) + f3 ! [dx2|s]
219  s(6, 1) = rap(1)*s(3, 1) ! [dxy|s]
220  s(7, 1) = rap(1)*s(4, 1) ! [dxz|s]
221  s(8, 1) = rap(2)*s(3, 1) + f3 ! [dy2|s]
222  s(9, 1) = rap(2)*s(4, 1) ! [dyz|s]
223  s(10, 1) = rap(3)*s(4, 1) + f3 ! [dz2|s]
224 
225  IF (la_max > 2) THEN
226 
227 ! *** [f|s] ***
228 
229  f3 = 2.0_dp*f2
230 
231  s(11, 1) = rap(1)*s(5, 1) + f3*s(2, 1) ! [fx3 |s]
232  s(12, 1) = rap(1)*s(6, 1) + f2*s(3, 1) ! [fx2y|s]
233  s(13, 1) = rap(1)*s(7, 1) + f2*s(4, 1) ! [fx2z|s]
234  s(14, 1) = rap(1)*s(8, 1) ! [fxy2|s]
235  s(15, 1) = rap(1)*s(9, 1) ! [fxyz|s]
236  s(16, 1) = rap(1)*s(10, 1) ! [fxz2|s]
237  s(17, 1) = rap(2)*s(8, 1) + f3*s(3, 1) ! [fy3 |s]
238  s(18, 1) = rap(2)*s(9, 1) + f2*s(4, 1) ! [fy2z|s]
239  s(19, 1) = rap(2)*s(10, 1) ! [fyz2|s]
240  s(20, 1) = rap(3)*s(10, 1) + f3*s(4, 1) ! [fz3 |s]
241 
242  IF (la_max > 3) THEN
243 
244 ! *** [g|s] ***
245 
246  f4 = 3.0_dp*f2
247 
248  s(21, 1) = rap(1)*s(11, 1) + f4*s(5, 1) ! [gx4 |s]
249  s(22, 1) = rap(1)*s(12, 1) + f3*s(6, 1) ! [gx3y |s]
250  s(23, 1) = rap(1)*s(13, 1) + f3*s(7, 1) ! [gx3z |s]
251  s(24, 1) = rap(1)*s(14, 1) + f2*s(8, 1) ! [gx2y2|s]
252  s(25, 1) = rap(1)*s(15, 1) + f2*s(9, 1) ! [gx2yz|s]
253  s(26, 1) = rap(1)*s(16, 1) + f2*s(10, 1) ! [gx2z2|s]
254  s(27, 1) = rap(1)*s(17, 1) ! [gxy3 |s]
255  s(28, 1) = rap(1)*s(18, 1) ! [gxy2z|s]
256  s(29, 1) = rap(1)*s(19, 1) ! [gxyz2|s]
257  s(30, 1) = rap(1)*s(20, 1) ! [gxz3 |s]
258  s(31, 1) = rap(2)*s(17, 1) + f4*s(8, 1) ! [gy4 |s]
259  s(32, 1) = rap(2)*s(18, 1) + f3*s(9, 1) ! [gy3z |s]
260  s(33, 1) = rap(2)*s(19, 1) + f2*s(10, 1) ! [gy2z2|s]
261  s(34, 1) = rap(2)*s(20, 1) ! [gyz3 |s]
262  s(35, 1) = rap(3)*s(20, 1) + f4*s(10, 1) ! [gz4 |s]
263 
264 ! *** [a|s] = (Pi - Ai)*[a-1i|s] + f2*Ni(a-1i)*[a-2i|s] ***
265 
266  DO la = 5, la_max
267 
268 ! *** Increase the angular momentum component z of a ***
269 
270  s(coset(0, 0, la), 1) = &
271  rap(3)*s(coset(0, 0, la - 1), 1) + &
272  f2*real(la - 1, dp)*s(coset(0, 0, la - 2), 1)
273 
274 ! *** Increase the angular momentum component y of a ***
275 
276  az = la - 1
277  s(coset(0, 1, az), 1) = rap(2)*s(coset(0, 0, az), 1)
278  DO ay = 2, la
279  az = la - ay
280  s(coset(0, ay, az), 1) = &
281  rap(2)*s(coset(0, ay - 1, az), 1) + &
282  f2*real(ay - 1, dp)*s(coset(0, ay - 2, az), 1)
283  END DO
284 
285 ! *** Increase the angular momentum component x of a ***
286 
287  DO ay = 0, la - 1
288  az = la - 1 - ay
289  s(coset(1, ay, az), 1) = rap(1)*s(coset(0, ay, az), 1)
290  END DO
291  DO ax = 2, la
292  f3 = f2*real(ax - 1, dp)
293  DO ay = 0, la - ax
294  az = la - ax - ay
295  s(coset(ax, ay, az), 1) = &
296  rap(1)*s(coset(ax - 1, ay, az), 1) + &
297  f3*s(coset(ax - 2, ay, az), 1)
298  END DO
299  END DO
300 
301  END DO
302 
303  END IF
304 
305  END IF
306 
307  END IF
308 
309 ! *** Recurrence steps: [a|s] -> [a|b] ***
310 
311  IF (lb_max > 0) THEN
312 
313  DO j = 2, ncoset(lb_max)
314  DO i = 1, ncoset(la_min)
315  s(i, j) = 0.0_dp
316  END DO
317  END DO
318 
319 ! *** Horizontal recurrence steps ***
320 
321  rbp(:) = rap(:) - rab(:)
322 
323 ! *** [a|p] = [a+1i|s] - (Bi - Ai)*[a|s] ***
324 
325  IF (lb_max == 1) THEN
326  la_start = la_min
327  ELSE
328  la_start = max(0, la_min - 1)
329  END IF
330 
331  DO la = la_start, la_max - 1
332  DO ax = 0, la
333  DO ay = 0, la - ax
334  az = la - ax - ay
335  coa = coset(ax, ay, az)
336  coapx = coset(ax + 1, ay, az)
337  coapy = coset(ax, ay + 1, az)
338  coapz = coset(ax, ay, az + 1)
339  s(coa, 2) = s(coapx, 1) - rab(1)*s(coa, 1)
340  s(coa, 3) = s(coapy, 1) - rab(2)*s(coa, 1)
341  s(coa, 4) = s(coapz, 1) - rab(3)*s(coa, 1)
342  END DO
343  END DO
344  END DO
345 
346 ! *** Vertical recurrence step ***
347 
348 ! *** [a|p] = (Pi - Bi)*[a|s] + f2*Ni(a)*[a-1i|s] ***
349 
350  DO ax = 0, la_max
351  fax = f2*real(ax, dp)
352  DO ay = 0, la_max - ax
353  fay = f2*real(ay, dp)
354  az = la_max - ax - ay
355  faz = f2*real(az, dp)
356  coa = coset(ax, ay, az)
357  coamx = coset(ax - 1, ay, az)
358  coamy = coset(ax, ay - 1, az)
359  coamz = coset(ax, ay, az - 1)
360  s(coa, 2) = rbp(1)*s(coa, 1) + fax*s(coamx, 1)
361  s(coa, 3) = rbp(2)*s(coa, 1) + fay*s(coamy, 1)
362  s(coa, 4) = rbp(3)*s(coa, 1) + faz*s(coamz, 1)
363  END DO
364  END DO
365 
366 ! *** Recurrence steps: [a|p] -> [a|b] ***
367 
368  DO lb = 2, lb_max
369 
370 ! *** Horizontal recurrence steps ***
371 
372 ! *** [a|b] = [a+1i|b-1i] - (Bi - Ai)*[a|b-1i] ***
373 
374  IF (lb == lb_max) THEN
375  la_start = la_min
376  ELSE
377  la_start = max(0, la_min - 1)
378  END IF
379 
380  DO la = la_start, la_max - 1
381  DO ax = 0, la
382  DO ay = 0, la - ax
383  az = la - ax - ay
384  coa = coset(ax, ay, az)
385  coapx = coset(ax + 1, ay, az)
386  coapy = coset(ax, ay + 1, az)
387  coapz = coset(ax, ay, az + 1)
388 
389 ! *** Shift of angular momentum component z from a to b ***
390 
391  cob = coset(0, 0, lb)
392  cobmz = coset(0, 0, lb - 1)
393  s(coa, cob) = s(coapz, cobmz) - rab(3)*s(coa, cobmz)
394 
395 ! *** Shift of angular momentum component y from a to b ***
396 
397  DO by = 1, lb
398  bz = lb - by
399  cob = coset(0, by, bz)
400  cobmy = coset(0, by - 1, bz)
401  s(coa, cob) = s(coapy, cobmy) - rab(2)*s(coa, cobmy)
402  END DO
403 
404 ! *** Shift of angular momentum component x from a to b ***
405 
406  DO bx = 1, lb
407  DO by = 0, lb - bx
408  bz = lb - bx - by
409  cob = coset(bx, by, bz)
410  cobmx = coset(bx - 1, by, bz)
411  s(coa, cob) = s(coapx, cobmx) - rab(1)*s(coa, cobmx)
412  END DO
413  END DO
414 
415  END DO
416  END DO
417  END DO
418 
419 ! *** Vertical recurrence step ***
420 
421 ! *** [a|b] = (Pi - Bi)*[a|b-1i] + f2*Ni(a)*[a-1i|b-1i] + ***
422 ! *** f2*Ni(b-1i)*[a|b-2i] ***
423 
424  DO ax = 0, la_max
425  fax = f2*real(ax, dp)
426  DO ay = 0, la_max - ax
427  fay = f2*real(ay, dp)
428  az = la_max - ax - ay
429  faz = f2*real(az, dp)
430  coa = coset(ax, ay, az)
431  coamx = coset(ax - 1, ay, az)
432  coamy = coset(ax, ay - 1, az)
433  coamz = coset(ax, ay, az - 1)
434 
435 ! *** Increase the angular momentum component z of b ***
436 
437  f3 = f2*real(lb - 1, dp)
438  cob = coset(0, 0, lb)
439  cobmz = coset(0, 0, lb - 1)
440  cobm2z = coset(0, 0, lb - 2)
441  s(coa, cob) = rbp(3)*s(coa, cobmz) + &
442  faz*s(coamz, cobmz) + &
443  f3*s(coa, cobm2z)
444 
445 ! *** Increase the angular momentum component y of b ***
446 
447  bz = lb - 1
448  cob = coset(0, 1, bz)
449  cobmy = coset(0, 0, bz)
450  s(coa, cob) = rbp(2)*s(coa, cobmy) + &
451  fay*s(coamy, cobmy)
452  DO by = 2, lb
453  bz = lb - by
454  f3 = f2*real(by - 1, dp)
455  cob = coset(0, by, bz)
456  cobmy = coset(0, by - 1, bz)
457  cobm2y = coset(0, by - 2, bz)
458  s(coa, cob) = rbp(2)*s(coa, cobmy) + &
459  fay*s(coamy, cobmy) + &
460  f3*s(coa, cobm2y)
461  END DO
462 
463 ! *** Increase the angular momentum component x of b ***
464 
465  DO by = 0, lb - 1
466  bz = lb - 1 - by
467  cob = coset(1, by, bz)
468  cobmx = coset(0, by, bz)
469  s(coa, cob) = rbp(1)*s(coa, cobmx) + &
470  fax*s(coamx, cobmx)
471  END DO
472  DO bx = 2, lb
473  f3 = f2*real(bx - 1, dp)
474  DO by = 0, lb - bx
475  bz = lb - bx - by
476  cob = coset(bx, by, bz)
477  cobmx = coset(bx - 1, by, bz)
478  cobm2x = coset(bx - 2, by, bz)
479  s(coa, cob) = rbp(1)*s(coa, cobmx) + &
480  fax*s(coamx, cobmx) + &
481  f3*s(coa, cobm2x)
482  END DO
483  END DO
484 
485  END DO
486  END DO
487 
488  END DO
489 
490  END IF
491 
492  ELSE
493 
494  IF (lb_max > 0) THEN
495 
496 ! *** Vertical recurrence steps: [s|s] -> [s|b] ***
497 
498  rbp(:) = (f1 - 1.0_dp)*rab(:)
499 
500 ! *** [s|p] = (Pi - Bi)*[s|s] ***
501 
502  s(1, 2) = rbp(1)*s(1, 1) ! [s|px]
503  s(1, 3) = rbp(2)*s(1, 1) ! [s|py]
504  s(1, 4) = rbp(3)*s(1, 1) ! [s|pz]
505 
506  IF (lb_max > 1) THEN
507 
508 ! *** [s|d] ***
509 
510  f3 = f2*s(1, 1)
511 
512  s(1, 5) = rbp(1)*s(1, 2) + f3 ! [s|dx2]
513  s(1, 6) = rbp(1)*s(1, 3) ! [s|dxy]
514  s(1, 7) = rbp(1)*s(1, 4) ! [s|dxz]
515  s(1, 8) = rbp(2)*s(1, 3) + f3 ! [s|dy2]
516  s(1, 9) = rbp(2)*s(1, 4) ! [s|dyz]
517  s(1, 10) = rbp(3)*s(1, 4) + f3 ! [s|dz2]
518 
519 ! *** [s|b] = (Pi - Bi)*[s|b-1i] + f2*Ni(b-1i)*[s|b-2i] ***
520 
521  DO lb = 3, lb_max
522 
523 ! *** Increase the angular momentum component z of b ***
524 
525  s(1, coset(0, 0, lb)) = &
526  rbp(3)*s(1, coset(0, 0, lb - 1)) + &
527  f2*real(lb - 1, dp)*s(1, coset(0, 0, lb - 2))
528 
529 ! *** Increase the angular momentum component y of b ***
530 
531  bz = lb - 1
532  s(1, coset(0, 1, bz)) = rbp(2)*s(1, coset(0, 0, bz))
533  DO by = 2, lb
534  bz = lb - by
535  s(1, coset(0, by, bz)) = &
536  rbp(2)*s(1, coset(0, by - 1, bz)) + &
537  f2*real(by - 1, dp)*s(1, coset(0, by - 2, bz))
538  END DO
539 
540 ! *** Increase the angular momentum component x of b ***
541 
542  DO by = 0, lb - 1
543  bz = lb - 1 - by
544  s(1, coset(1, by, bz)) = rbp(1)*s(1, coset(0, by, bz))
545  END DO
546  DO bx = 2, lb
547  f3 = f2*real(bx - 1, dp)
548  DO by = 0, lb - bx
549  bz = lb - bx - by
550  s(1, coset(bx, by, bz)) = &
551  rbp(1)*s(1, coset(bx - 1, by, bz)) + &
552  f3*s(1, coset(bx - 2, by, bz))
553  END DO
554  END DO
555 
556  END DO
557 
558  END IF
559 
560  END IF
561 
562  END IF
563 
564 ! *** Store the primitive overlap integrals ***
565  DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
566  njb(1:3) = indco(1:3, jb)
567  DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
568  nib(1:3) = indco(1:3, ib)
569  nbb(1:3) = nib + njb
570  DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
571  nja(1:3) = indco(1:3, ja)
572  DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
573  nia(1:3) = indco(1:3, ia)
574  naa(1:3) = nia + nja
575  ! now loop over all elements of s
576  DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
577  nb(1:3) = indco(1:3, j)
578  DO i = ncoset(la_min - 1) + 1, ncoset(la_max)
579  na(1:3) = indco(1:3, i)
580  IF (all(na == naa) .AND. all(nb == nbb)) THEN
581  saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = s(i, j)
582  IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(i, j)
583  IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(i, j)
584  IF (asets_equal .AND. bsets_equal) THEN
585  saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(i, j)
586  saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(i, j)
587  saabb(ncoa2 + ja, ncoa1 + ia, ncob2 + jb, ncob1 + ib) = s(i, j)
588  END IF
589  END IF
590  END DO
591  END DO
592  END DO
593  END DO
594  END DO
595  END DO
596 
597  ncob2 = ncob2 + ncoset(lb_max_set2)
598 
599  END DO
600 
601  ncob1 = ncob1 + ncoset(lb_max_set1)
602 
603  END DO
604 
605  ncoa2 = ncoa2 + ncoset(la_max_set2)
606 
607  END DO
608 
609  ncoa1 = ncoa1 + ncoset(la_max_set1)
610 
611  END DO
612 
613  CALL timestop(handle)
614 
615  END SUBROUTINE overlap_aabb
616 
617 END MODULE ai_overlap_aabb
ai_overlap_aabb
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap_aabb.F:27
ai_overlap_aabb::overlap_aabb
subroutine, public overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, asets_equal, bsets_equal, rab, dab, saabb, s, lds)
Purpose: Calculation of the two-center overlap integrals [aa|bb] over Cartesian Gaussian-type functio...
Definition: ai_overlap_aabb.F:85
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
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43