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