Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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 9 : SUBROUTINE overlap_aabb(la_max_set1, la_min_set1, npgfa1, rpgfa1, zeta1, &
81 9 : la_max_set2, la_min_set2, npgfa2, rpgfa2, zeta2, &
82 9 : lb_max_set1, lb_min_set1, npgfb1, rpgfb1, zetb1, &
83 18 : lb_max_set2, lb_min_set2, npgfb2, rpgfb2, zetb2, &
84 9 : 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 9 : CALL timeset(routineN, handle)
114 :
115 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
116 :
117 9 : ncoa1 = 0
118 9 : ncoa2 = 0
119 9 : ncob1 = 0
120 9 : ncob2 = 0
121 :
122 72 : DO ipgf = 1, npgfa1
123 :
124 63 : ncoa2 = 0
125 :
126 63 : IF (asets_equal) THEN
127 252 : jpgf_start = ipgf
128 252 : DO i = 1, jpgf_start - 1
129 252 : ncoa2 = ncoa2 + ncoset(la_max_set2)
130 : END DO
131 : ELSE
132 : jpgf_start = 1
133 : END IF
134 :
135 315 : DO jpgf = jpgf_start, npgfa2
136 :
137 252 : ncob1 = 0
138 252 : zeta = zeta1(ipgf) + zeta2(jpgf)
139 252 : la_max = la_max_set1 + la_max_set2
140 252 : la_min = la_min_set1 + la_min_set2
141 :
142 2016 : DO kpgf = 1, npgfb1
143 :
144 1764 : ncob2 = 0
145 :
146 1764 : IF (bsets_equal) THEN
147 7056 : lpgf_start = kpgf
148 7056 : DO i = 1, lpgf_start - 1
149 7056 : ncob2 = ncob2 + ncoset(lb_max_set2)
150 : END DO
151 : ELSE
152 : lpgf_start = 1
153 : END IF
154 :
155 8820 : DO lpgf = lpgf_start, npgfb2
156 :
157 : ! *** Screening ***
158 : IF ((rpgfa1(ipgf) + rpgfb1(kpgf) < dab) .OR. &
159 : (rpgfa2(jpgf) + rpgfb1(kpgf) < dab) .OR. &
160 7056 : (rpgfa1(ipgf) + rpgfb2(lpgf) < dab) .OR. &
161 : (rpgfa2(jpgf) + rpgfb2(lpgf) < dab)) THEN
162 735 : DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
163 3087 : DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
164 12348 : DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
165 49392 : DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
166 37632 : saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = 0._dp
167 37632 : IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = 0._dp
168 37632 : IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = 0._dp
169 47040 : IF (asets_equal .AND. bsets_equal) THEN
170 37632 : saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = 0._dp
171 37632 : saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = 0._dp
172 37632 : 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 147 : ncob2 = ncob2 + ncoset(lb_max_set2)
179 147 : CYCLE
180 : END IF
181 :
182 6909 : zetb = zetb1(kpgf) + zetb2(lpgf)
183 6909 : lb_max = lb_max_set1 + lb_max_set2
184 6909 : lb_min = lb_min_set1 + lb_min_set2
185 :
186 : ! *** Calculate some prefactors ***
187 :
188 6909 : zetp = 1.0_dp/(zeta + zetb)
189 :
190 6909 : f0 = SQRT((pi*zetp)**3)
191 6909 : f1 = zetb*zetp
192 6909 : f2 = 0.5_dp*zetp
193 :
194 : ! *** Calculate the basic two-center overlap integral [s|s] ***
195 :
196 6909 : s(1, 1) = f0*EXP(-zeta*f1*dab*dab)
197 :
198 : ! *** Recurrence steps: [s|s] -> [a|b] ***
199 :
200 6909 : IF (la_max > 0) THEN
201 :
202 : ! *** Vertical recurrence steps: [s|s] -> [a|s] ***
203 :
204 27636 : rap(:) = f1*rab(:)
205 :
206 : ! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) ***
207 :
208 6909 : s(2, 1) = rap(1)*s(1, 1) ! [px|s]
209 6909 : s(3, 1) = rap(2)*s(1, 1) ! [py|s]
210 6909 : s(4, 1) = rap(3)*s(1, 1) ! [pz|s]
211 :
212 6909 : IF (la_max > 1) THEN
213 :
214 : ! *** [d|s] ***
215 :
216 6909 : f3 = f2*s(1, 1)
217 :
218 6909 : s(5, 1) = rap(1)*s(2, 1) + f3 ! [dx2|s]
219 6909 : s(6, 1) = rap(1)*s(3, 1) ! [dxy|s]
220 6909 : s(7, 1) = rap(1)*s(4, 1) ! [dxz|s]
221 6909 : s(8, 1) = rap(2)*s(3, 1) + f3 ! [dy2|s]
222 6909 : s(9, 1) = rap(2)*s(4, 1) ! [dyz|s]
223 6909 : s(10, 1) = rap(3)*s(4, 1) + f3 ! [dz2|s]
224 :
225 6909 : IF (la_max > 2) THEN
226 :
227 : ! *** [f|s] ***
228 :
229 0 : f3 = 2.0_dp*f2
230 :
231 0 : s(11, 1) = rap(1)*s(5, 1) + f3*s(2, 1) ! [fx3 |s]
232 0 : s(12, 1) = rap(1)*s(6, 1) + f2*s(3, 1) ! [fx2y|s]
233 0 : s(13, 1) = rap(1)*s(7, 1) + f2*s(4, 1) ! [fx2z|s]
234 0 : s(14, 1) = rap(1)*s(8, 1) ! [fxy2|s]
235 0 : s(15, 1) = rap(1)*s(9, 1) ! [fxyz|s]
236 0 : s(16, 1) = rap(1)*s(10, 1) ! [fxz2|s]
237 0 : s(17, 1) = rap(2)*s(8, 1) + f3*s(3, 1) ! [fy3 |s]
238 0 : s(18, 1) = rap(2)*s(9, 1) + f2*s(4, 1) ! [fy2z|s]
239 0 : s(19, 1) = rap(2)*s(10, 1) ! [fyz2|s]
240 0 : s(20, 1) = rap(3)*s(10, 1) + f3*s(4, 1) ! [fz3 |s]
241 :
242 0 : IF (la_max > 3) THEN
243 :
244 : ! *** [g|s] ***
245 :
246 0 : f4 = 3.0_dp*f2
247 :
248 0 : s(21, 1) = rap(1)*s(11, 1) + f4*s(5, 1) ! [gx4 |s]
249 0 : s(22, 1) = rap(1)*s(12, 1) + f3*s(6, 1) ! [gx3y |s]
250 0 : s(23, 1) = rap(1)*s(13, 1) + f3*s(7, 1) ! [gx3z |s]
251 0 : s(24, 1) = rap(1)*s(14, 1) + f2*s(8, 1) ! [gx2y2|s]
252 0 : s(25, 1) = rap(1)*s(15, 1) + f2*s(9, 1) ! [gx2yz|s]
253 0 : s(26, 1) = rap(1)*s(16, 1) + f2*s(10, 1) ! [gx2z2|s]
254 0 : s(27, 1) = rap(1)*s(17, 1) ! [gxy3 |s]
255 0 : s(28, 1) = rap(1)*s(18, 1) ! [gxy2z|s]
256 0 : s(29, 1) = rap(1)*s(19, 1) ! [gxyz2|s]
257 0 : s(30, 1) = rap(1)*s(20, 1) ! [gxz3 |s]
258 0 : s(31, 1) = rap(2)*s(17, 1) + f4*s(8, 1) ! [gy4 |s]
259 0 : s(32, 1) = rap(2)*s(18, 1) + f3*s(9, 1) ! [gy3z |s]
260 0 : s(33, 1) = rap(2)*s(19, 1) + f2*s(10, 1) ! [gy2z2|s]
261 0 : s(34, 1) = rap(2)*s(20, 1) ! [gyz3 |s]
262 0 : 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 0 : 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 0 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1)
273 :
274 : ! *** Increase the angular momentum component y of a ***
275 :
276 0 : az = la - 1
277 0 : s(coset(0, 1, az), 1) = rap(2)*s(coset(0, 0, az), 1)
278 0 : DO ay = 2, la
279 0 : az = la - ay
280 : s(coset(0, ay, az), 1) = &
281 : rap(2)*s(coset(0, ay - 1, az), 1) + &
282 0 : 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 0 : DO ay = 0, la - 1
288 0 : az = la - 1 - ay
289 0 : s(coset(1, ay, az), 1) = rap(1)*s(coset(0, ay, az), 1)
290 : END DO
291 0 : DO ax = 2, la
292 0 : f3 = f2*REAL(ax - 1, dp)
293 0 : DO ay = 0, la - ax
294 0 : az = la - ax - ay
295 : s(coset(ax, ay, az), 1) = &
296 : rap(1)*s(coset(ax - 1, ay, az), 1) + &
297 0 : 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 6909 : IF (lb_max > 0) THEN
312 :
313 69090 : DO j = 2, ncoset(lb_max)
314 131271 : DO i = 1, ncoset(la_min)
315 124362 : s(i, j) = 0.0_dp
316 : END DO
317 : END DO
318 :
319 : ! *** Horizontal recurrence steps ***
320 :
321 27636 : rbp(:) = rap(:) - rab(:)
322 :
323 : ! *** [a|p] = [a+1i|s] - (Bi - Ai)*[a|s] ***
324 :
325 6909 : IF (lb_max == 1) THEN
326 : la_start = la_min
327 : ELSE
328 6909 : la_start = MAX(0, la_min - 1)
329 : END IF
330 :
331 20727 : DO la = la_start, la_max - 1
332 41454 : DO ax = 0, la
333 62181 : DO ay = 0, la - ax
334 27636 : az = la - ax - ay
335 27636 : coa = coset(ax, ay, az)
336 27636 : coapx = coset(ax + 1, ay, az)
337 27636 : coapy = coset(ax, ay + 1, az)
338 27636 : coapz = coset(ax, ay, az + 1)
339 27636 : s(coa, 2) = s(coapx, 1) - rab(1)*s(coa, 1)
340 27636 : s(coa, 3) = s(coapy, 1) - rab(2)*s(coa, 1)
341 48363 : 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 27636 : DO ax = 0, la_max
351 20727 : fax = f2*REAL(ax, dp)
352 69090 : DO ay = 0, la_max - ax
353 41454 : fay = f2*REAL(ay, dp)
354 41454 : az = la_max - ax - ay
355 41454 : faz = f2*REAL(az, dp)
356 41454 : coa = coset(ax, ay, az)
357 41454 : coamx = coset(ax - 1, ay, az)
358 41454 : coamy = coset(ax, ay - 1, az)
359 41454 : coamz = coset(ax, ay, az - 1)
360 41454 : s(coa, 2) = rbp(1)*s(coa, 1) + fax*s(coamx, 1)
361 41454 : s(coa, 3) = rbp(2)*s(coa, 1) + fay*s(coamy, 1)
362 62181 : 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 13818 : 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 6909 : IF (lb == lb_max) THEN
375 : la_start = la_min
376 : ELSE
377 0 : la_start = MAX(0, la_min - 1)
378 : END IF
379 :
380 20727 : DO la = la_start, la_max - 1
381 41454 : DO ax = 0, la
382 62181 : DO ay = 0, la - ax
383 27636 : az = la - ax - ay
384 27636 : coa = coset(ax, ay, az)
385 27636 : coapx = coset(ax + 1, ay, az)
386 27636 : coapy = coset(ax, ay + 1, az)
387 27636 : coapz = coset(ax, ay, az + 1)
388 :
389 : ! *** Shift of angular momentum component z from a to b ***
390 :
391 27636 : cob = coset(0, 0, lb)
392 27636 : cobmz = coset(0, 0, lb - 1)
393 27636 : 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 82908 : DO by = 1, lb
398 55272 : bz = lb - by
399 55272 : cob = coset(0, by, bz)
400 55272 : cobmy = coset(0, by - 1, bz)
401 82908 : 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 103635 : DO bx = 1, lb
407 165816 : DO by = 0, lb - bx
408 82908 : bz = lb - bx - by
409 82908 : cob = coset(bx, by, bz)
410 82908 : cobmx = coset(bx - 1, by, bz)
411 138180 : 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 34545 : DO ax = 0, la_max
425 20727 : fax = f2*REAL(ax, dp)
426 69090 : DO ay = 0, la_max - ax
427 41454 : fay = f2*REAL(ay, dp)
428 41454 : az = la_max - ax - ay
429 41454 : faz = f2*REAL(az, dp)
430 41454 : coa = coset(ax, ay, az)
431 41454 : coamx = coset(ax - 1, ay, az)
432 41454 : coamy = coset(ax, ay - 1, az)
433 41454 : coamz = coset(ax, ay, az - 1)
434 :
435 : ! *** Increase the angular momentum component z of b ***
436 :
437 41454 : f3 = f2*REAL(lb - 1, dp)
438 41454 : cob = coset(0, 0, lb)
439 41454 : cobmz = coset(0, 0, lb - 1)
440 41454 : cobm2z = coset(0, 0, lb - 2)
441 : s(coa, cob) = rbp(3)*s(coa, cobmz) + &
442 : faz*s(coamz, cobmz) + &
443 41454 : f3*s(coa, cobm2z)
444 :
445 : ! *** Increase the angular momentum component y of b ***
446 :
447 41454 : bz = lb - 1
448 41454 : cob = coset(0, 1, bz)
449 41454 : cobmy = coset(0, 0, bz)
450 : s(coa, cob) = rbp(2)*s(coa, cobmy) + &
451 41454 : fay*s(coamy, cobmy)
452 82908 : DO by = 2, lb
453 41454 : bz = lb - by
454 41454 : f3 = f2*REAL(by - 1, dp)
455 41454 : cob = coset(0, by, bz)
456 41454 : cobmy = coset(0, by - 1, bz)
457 41454 : cobm2y = coset(0, by - 2, bz)
458 : s(coa, cob) = rbp(2)*s(coa, cobmy) + &
459 : fay*s(coamy, cobmy) + &
460 82908 : f3*s(coa, cobm2y)
461 : END DO
462 :
463 : ! *** Increase the angular momentum component x of b ***
464 :
465 124362 : DO by = 0, lb - 1
466 82908 : bz = lb - 1 - by
467 82908 : cob = coset(1, by, bz)
468 82908 : cobmx = coset(0, by, bz)
469 : s(coa, cob) = rbp(1)*s(coa, cobmx) + &
470 124362 : fax*s(coamx, cobmx)
471 : END DO
472 103635 : DO bx = 2, lb
473 41454 : f3 = f2*REAL(bx - 1, dp)
474 124362 : DO by = 0, lb - bx
475 41454 : bz = lb - bx - by
476 41454 : cob = coset(bx, by, bz)
477 41454 : cobmx = coset(bx - 1, by, bz)
478 41454 : cobm2x = coset(bx - 2, by, bz)
479 : s(coa, cob) = rbp(1)*s(coa, cobmx) + &
480 : fax*s(coamx, cobmx) + &
481 82908 : 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 0 : IF (lb_max > 0) THEN
495 :
496 : ! *** Vertical recurrence steps: [s|s] -> [s|b] ***
497 :
498 0 : rbp(:) = (f1 - 1.0_dp)*rab(:)
499 :
500 : ! *** [s|p] = (Pi - Bi)*[s|s] ***
501 :
502 0 : s(1, 2) = rbp(1)*s(1, 1) ! [s|px]
503 0 : s(1, 3) = rbp(2)*s(1, 1) ! [s|py]
504 0 : s(1, 4) = rbp(3)*s(1, 1) ! [s|pz]
505 :
506 0 : IF (lb_max > 1) THEN
507 :
508 : ! *** [s|d] ***
509 :
510 0 : f3 = f2*s(1, 1)
511 :
512 0 : s(1, 5) = rbp(1)*s(1, 2) + f3 ! [s|dx2]
513 0 : s(1, 6) = rbp(1)*s(1, 3) ! [s|dxy]
514 0 : s(1, 7) = rbp(1)*s(1, 4) ! [s|dxz]
515 0 : s(1, 8) = rbp(2)*s(1, 3) + f3 ! [s|dy2]
516 0 : s(1, 9) = rbp(2)*s(1, 4) ! [s|dyz]
517 0 : 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 0 : 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 0 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2))
528 :
529 : ! *** Increase the angular momentum component y of b ***
530 :
531 0 : bz = lb - 1
532 0 : s(1, coset(0, 1, bz)) = rbp(2)*s(1, coset(0, 0, bz))
533 0 : DO by = 2, lb
534 0 : bz = lb - by
535 : s(1, coset(0, by, bz)) = &
536 : rbp(2)*s(1, coset(0, by - 1, bz)) + &
537 0 : 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 0 : DO by = 0, lb - 1
543 0 : bz = lb - 1 - by
544 0 : s(1, coset(1, by, bz)) = rbp(1)*s(1, coset(0, by, bz))
545 : END DO
546 0 : DO bx = 2, lb
547 0 : f3 = f2*REAL(bx - 1, dp)
548 0 : DO by = 0, lb - bx
549 0 : bz = lb - bx - by
550 : s(1, coset(bx, by, bz)) = &
551 : rbp(1)*s(1, coset(bx - 1, by, bz)) + &
552 0 : 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 34545 : DO jb = ncoset(lb_min_set2 - 1) + 1, ncoset(lb_max_set2)
566 110544 : njb(1:3) = indco(1:3, jb)
567 145089 : DO ib = ncoset(lb_min_set1 - 1) + 1, ncoset(lb_max_set1)
568 442176 : nib(1:3) = indco(1:3, ib)
569 442176 : nbb(1:3) = nib + njb
570 580356 : DO ja = ncoset(la_min_set2 - 1) + 1, ncoset(la_max_set2)
571 1768704 : nja(1:3) = indco(1:3, ja)
572 2321424 : DO ia = ncoset(la_min_set1 - 1) + 1, ncoset(la_max_set1)
573 7074816 : nia(1:3) = indco(1:3, ia)
574 7074816 : naa(1:3) = nia + nja
575 : ! now loop over all elements of s
576 19897920 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
577 70748160 : nb(1:3) = indco(1:3, j)
578 196326144 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max)
579 707481600 : na(1:3) = indco(1:3, i)
580 638944320 : IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN
581 1768704 : saabb(ncoa1 + ia, ncoa2 + ja, ncob1 + ib, ncob2 + jb) = s(i, j)
582 1768704 : IF (asets_equal) saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(i, j)
583 1768704 : IF (bsets_equal) saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(i, j)
584 1768704 : IF (asets_equal .AND. bsets_equal) THEN
585 1768704 : saabb(ncoa2 + ja, ncoa1 + ia, ncob1 + ib, ncob2 + jb) = s(i, j)
586 1768704 : saabb(ncoa1 + ia, ncoa2 + ja, ncob2 + jb, ncob1 + ib) = s(i, j)
587 1768704 : 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 8673 : ncob2 = ncob2 + ncoset(lb_max_set2)
598 :
599 : END DO
600 :
601 2016 : ncob1 = ncob1 + ncoset(lb_max_set1)
602 :
603 : END DO
604 :
605 315 : ncoa2 = ncoa2 + ncoset(la_max_set2)
606 :
607 : END DO
608 :
609 72 : ncoa1 = ncoa1 + ncoset(la_max_set1)
610 :
611 : END DO
612 :
613 9 : CALL timestop(handle)
614 :
615 9 : END SUBROUTINE overlap_aabb
616 :
617 : END MODULE ai_overlap_aabb
|
