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