Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : !> \date 19.09.2000
64 : !> \author MK
65 : !> \version 1.0
66 : ! **************************************************************************************************
67 2012596 : SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
68 4025192 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
69 2012596 : rab, dab, sab, da_max_set, return_derivatives, s, lds)
70 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
71 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
72 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
73 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
74 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
75 : REAL(KIND=dp), INTENT(IN) :: dab
76 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
77 : INTEGER, INTENT(IN) :: da_max_set
78 : LOGICAL, INTENT(IN) :: return_derivatives
79 : INTEGER, INTENT(IN) :: lds
80 : REAL(KIND=dp), DIMENSION(lds, lds, *), &
81 : INTENT(INOUT) :: s
82 :
83 : INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, coa, coamx, coamy, coamz, coapx, &
84 : coapy, coapz, cob, da, da_max, dax, day, daz, i, ipgf, j, jk, jpgf, jstart, k, la, &
85 : la_max, la_start, lb, lb_max, lb_start, ldrr, na, nb
86 : REAL(KIND=dp) :: f0, fax, fay, faz, ftz, zetp
87 2012596 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
88 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
89 :
90 2012596 : da_max = da_max_set
91 2012596 : la_max = la_max_set + da_max_set
92 :
93 2012596 : lb_max = lb_max_set
94 2012596 : ldrr = MAX(la_max, lb_max) + 1
95 10062980 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
96 :
97 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
98 :
99 2012596 : na = 0
100 10268382 : DO ipgf = 1, npgfa
101 :
102 8255786 : nb = 0
103 :
104 22287392 : DO jpgf = 1, npgfb
105 :
106 : ! *** Screening ***
107 :
108 14031606 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 49876200 : DO j = nb + 1, nb + ncoset(lb_max_set)
110 248700078 : DO i = na + 1, na + ncoset(la_max_set)
111 239599280 : sab(i, j) = 0.0_dp
112 : END DO
113 : END DO
114 9100798 : IF (return_derivatives) THEN
115 19679854 : DO k = 2, ncoset(da_max_set)
116 10592208 : jstart = (k - 1)*SIZE(sab, 1)
117 61851544 : DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
118 223903776 : DO i = na + 1, na + ncoset(la_max_set)
119 213311568 : sab(i, j) = 0.0_dp
120 : END DO
121 : END DO
122 : END DO
123 : END IF
124 9100798 : nb = nb + ncoset(lb_max_set)
125 9100798 : CYCLE
126 : END IF
127 :
128 : ! *** Calculate some prefactors ***
129 :
130 4930808 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
131 :
132 4930808 : f0 = SQRT((pi*zetp)**3)*EXP(-zeta(ipgf)*zetb(jpgf)*zetp*dab*dab)
133 19723232 : rap(:) = zetb(jpgf)*zetp*rab(:)
134 19723232 : rbp(:) = -zeta(ipgf)*zetp*rab(:)
135 :
136 4930808 : CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
137 :
138 13992856 : DO lb = 0, lb_max
139 28647660 : DO bx = 0, lb
140 45819432 : DO by = 0, lb - bx
141 22102580 : bz = lb - bx - by
142 22102580 : cob = coset(bx, by, bz)
143 89817255 : DO la = 0, la_max
144 173458949 : DO ax = 0, la
145 311384221 : DO ay = 0, la - ax
146 160027852 : az = la - ax - ay
147 160027852 : coa = coset(ax, ay, az)
148 258324350 : s(coa, cob, 1) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
149 : END DO
150 : END DO
151 : END DO
152 : END DO
153 : END DO
154 : END DO
155 :
156 : ! *** Store the primitive overlap integrals ***
157 :
158 27033388 : DO j = 1, ncoset(lb_max_set)
159 143730734 : DO i = 1, ncoset(la_max_set)
160 138799926 : sab(na + i, nb + j) = s(i, j, 1)
161 : END DO
162 : END DO
163 :
164 : ! *** Calculate the requested derivatives with respect ***
165 : ! *** to the nuclear coordinates of the atomic center a ***
166 :
167 4930808 : IF (return_derivatives) THEN
168 : la_start = 0
169 : lb_start = 0
170 : ELSE
171 16782 : la_start = la_min_set
172 16782 : lb_start = lb_min_set
173 : END IF
174 :
175 6659621 : DO da = 0, da_max - 1
176 1728813 : ftz = 2.0_dp*zeta(ipgf)
177 8388434 : DO dax = 0, da
178 5186439 : DO day = 0, da - dax
179 1728813 : daz = da - dax - day
180 1728813 : cda = coset(dax, day, daz)
181 1728813 : cdax = coset(dax + 1, day, daz)
182 1728813 : cday = coset(dax, day + 1, daz)
183 1728813 : cdaz = coset(dax, day, daz + 1)
184 :
185 : ! *** [da/dAi|b] = 2*zeta*[a+1i|b] - Ni(a)[a-1i|b] ***
186 :
187 6573757 : DO la = la_start, la_max - da - 1
188 9607658 : DO ax = 0, la
189 4762714 : fax = REAL(ax, dp)
190 14548351 : DO ay = 0, la - ax
191 6669506 : fay = REAL(ay, dp)
192 6669506 : az = la - ax - ay
193 6669506 : faz = REAL(az, dp)
194 6669506 : coa = coset(ax, ay, az)
195 6669506 : coamx = coset(ax - 1, ay, az)
196 6669506 : coamy = coset(ax, ay - 1, az)
197 6669506 : coamz = coset(ax, ay, az - 1)
198 6669506 : coapx = coset(ax + 1, ay, az)
199 6669506 : coapy = coset(ax, ay + 1, az)
200 6669506 : coapz = coset(ax, ay, az + 1)
201 24052744 : DO lb = lb_start, lb_max_set
202 40529952 : DO bx = 0, lb
203 67290196 : DO by = 0, lb - bx
204 33429750 : bz = lb - bx - by
205 33429750 : cob = coset(bx, by, bz)
206 : s(coa, cob, cdax) = ftz*s(coapx, cob, cda) - &
207 33429750 : fax*s(coamx, cob, cda)
208 : s(coa, cob, cday) = ftz*s(coapy, cob, cda) - &
209 33429750 : fay*s(coamy, cob, cda)
210 : s(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - &
211 54669672 : faz*s(coamz, cob, cda)
212 : END DO
213 : END DO
214 : END DO
215 : END DO
216 : END DO
217 : END DO
218 :
219 : END DO
220 : END DO
221 : END DO
222 :
223 : ! *** Return all the calculated derivatives of the ***
224 : ! *** primitive overlap integrals, if requested ***
225 :
226 4930808 : IF (return_derivatives) THEN
227 10100465 : DO k = 2, ncoset(da_max_set)
228 5186439 : jstart = (k - 1)*SIZE(sab, 1)
229 30936890 : DO j = 1, ncoset(lb_max_set)
230 20836425 : jk = jstart + j
231 126312114 : DO i = 1, ncoset(la_max_set)
232 121125675 : sab(na + i, nb + jk) = s(i, j, k)
233 : END DO
234 : END DO
235 : END DO
236 : END IF
237 :
238 13186594 : nb = nb + ncoset(lb_max_set)
239 :
240 : END DO
241 :
242 10268382 : na = na + ncoset(la_max_set)
243 : END DO
244 :
245 2012596 : DEALLOCATE (rr)
246 :
247 2012596 : END SUBROUTINE overlap
248 :
249 : ! **************************************************************************************************
250 : !> \brief Calculation of the two-center overlap integrals [a|b] over
251 : !> Cartesian Gaussian-type functions. First and second derivatives
252 : !> \param la_max Max L on center A
253 : !> \param la_min Min L on center A
254 : !> \param npgfa Number of primitives on center A
255 : !> \param rpgfa Range of functions on A, used for screening
256 : !> \param zeta Exponents on center A
257 : !> \param lb_max Max L on center B
258 : !> \param lb_min Min L on center B
259 : !> \param npgfb Number of primitives on center B
260 : !> \param rpgfb Range of functions on B, used for screening
261 : !> \param zetb Exponents on center B
262 : !> \param rab Distance vector A-B
263 : !> \param sab Final overlap integrals
264 : !> \param dab First derivative overlap integrals
265 : !> \param ddab Second derivative overlap integrals
266 : !> \param rr_work Optional caller-provided one-dimensional overlap recurrence workspace
267 : !> \date 01.07.2014
268 : !> \author JGH
269 : ! **************************************************************************************************
270 18651647 : SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
271 18651647 : lb_max, lb_min, npgfb, rpgfb, zetb, &
272 18651647 : rab, sab, dab, ddab, rr_work)
273 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
274 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
275 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
276 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
277 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
278 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
279 : OPTIONAL :: sab
280 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
281 : OPTIONAL :: dab, ddab
282 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
283 : OPTIONAL, TARGET :: rr_work
284 :
285 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, ia, &
286 : ib, ipgf, jpgf, la, lb, ldrr, lma, &
287 : lmb, ma, mb, na, nb, ofa, ofb
288 : REAL(KIND=dp) :: a, ambm, ambp, apbm, apbp, b, dumx, &
289 : dumy, dumz, f0, rab2, tab, xhi, zet
290 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
291 18651647 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rr
292 :
293 : ! Distance of the centers a and b
294 :
295 18651647 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
296 18651647 : tab = SQRT(rab2)
297 :
298 : ! Maximum l for auxiliary integrals
299 18651647 : CPASSERT(PRESENT(sab) .OR. PRESENT(dab) .OR. PRESENT(ddab))
300 18651647 : IF (PRESENT(sab)) THEN
301 18349854 : lma = la_max
302 18349854 : lmb = lb_max
303 : END IF
304 18651647 : IF (PRESENT(dab)) THEN
305 4931598 : lma = la_max + 1
306 4931598 : lmb = lb_max
307 : END IF
308 18651647 : IF (PRESENT(ddab)) THEN
309 13855 : lma = la_max + 1
310 13855 : lmb = lb_max + 1
311 : END IF
312 18651647 : ldrr = MAX(lma, lmb) + 1
313 :
314 : ! Allocate or attach the workspace for auxiliary integrals
315 : NULLIFY (rr)
316 18651647 : IF (PRESENT(rr_work)) THEN
317 149387 : CPASSERT(SIZE(rr_work) >= ldrr*ldrr*3)
318 149387 : rr(0:ldrr - 1, 0:ldrr - 1, 1:3) => rr_work(1:ldrr*ldrr*3)
319 : ELSE
320 92511300 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
321 : END IF
322 :
323 : ! Number of integrals, check size of arrays
324 18651647 : ofa = ncoset(la_min - 1)
325 18651647 : ofb = ncoset(lb_min - 1)
326 18651647 : na = ncoset(la_max) - ofa
327 18651647 : nb = ncoset(lb_max) - ofb
328 18651647 : IF (PRESENT(sab)) THEN
329 18349854 : CPASSERT((SIZE(sab, 1) >= na*npgfa))
330 18349854 : CPASSERT((SIZE(sab, 2) >= nb*npgfb))
331 : END IF
332 18651647 : IF (PRESENT(dab)) THEN
333 4931598 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
334 4931598 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
335 4931598 : CPASSERT((SIZE(dab, 3) >= 3))
336 : END IF
337 18651647 : IF (PRESENT(ddab)) THEN
338 13855 : CPASSERT((SIZE(ddab, 1) >= na*npgfa))
339 13855 : CPASSERT((SIZE(ddab, 2) >= nb*npgfb))
340 13855 : CPASSERT((SIZE(ddab, 3) >= 6))
341 : END IF
342 :
343 : ! Loops over all pairs of primitive Gaussian-type functions
344 18651647 : ma = 0
345 112919241 : DO ipgf = 1, npgfa
346 94267594 : mb = 0
347 641915073 : DO jpgf = 1, npgfb
348 : ! Distance Screening
349 547647479 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
350 4521227905 : IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
351 3527545863 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
352 425926419 : IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
353 419218005 : mb = mb + nb
354 419218005 : CYCLE
355 : END IF
356 :
357 : ! Calculate some prefactors
358 128429474 : a = zeta(ipgf)
359 128429474 : b = zetb(jpgf)
360 128429474 : zet = a + b
361 128429474 : xhi = a*b/zet
362 513717896 : rap = b*rab/zet
363 513717896 : rbp = -a*rab/zet
364 :
365 : ! [s|s] integral
366 128429474 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
367 :
368 : ! Calculate the recurrence relation
369 128429474 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
370 :
371 267067143 : DO lb = lb_min, lb_max
372 491340865 : DO bx = 0, lb
373 695185687 : DO by = 0, lb - bx
374 332274296 : bz = lb - bx - by
375 332274296 : cob = coset(bx, by, bz) - ofb
376 332274296 : ib = mb + cob
377 939590113 : DO la = la_min, la_max
378 1388317254 : DO ax = 0, la
379 2102832922 : DO ay = 0, la - ax
380 1046789964 : az = la - ax - ay
381 1046789964 : coa = coset(ax, ay, az) - ofa
382 1046789964 : ia = ma + coa
383 : ! integrals
384 1046789964 : IF (PRESENT(sab)) THEN
385 1042003065 : sab(ia, ib) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
386 : END IF
387 : ! first derivatives
388 1046789964 : IF (PRESENT(dab)) THEN
389 : ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
390 : ! dx
391 213399147 : dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
392 213399147 : IF (ax > 0) dumx = dumx - REAL(ax, dp)*rr(ax - 1, bx, 1)
393 213399147 : dab(ia, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
394 : ! dy
395 213399147 : dumy = 2.0_dp*a*rr(ay + 1, by, 2)
396 213399147 : IF (ay > 0) dumy = dumy - REAL(ay, dp)*rr(ay - 1, by, 2)
397 213399147 : dab(ia, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
398 : ! dz
399 213399147 : dumz = 2.0_dp*a*rr(az + 1, bz, 3)
400 213399147 : IF (az > 0) dumz = dumz - REAL(az, dp)*rr(az - 1, bz, 3)
401 213399147 : dab(ia, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
402 : END IF
403 : ! 2nd derivatives
404 1719790827 : IF (PRESENT(ddab)) THEN
405 : ! (dda|b) = -4*a*b*(a+1|b+1) + 2*a*N(b)*(a+1|b-1)
406 : ! + 2*b*N(a)*(a-1|b+1) - N(a)*N(b)*(a-1|b-1)
407 : ! dx dx
408 349183 : apbp = f0*rr(ax + 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
409 349183 : IF (bx > 0) THEN
410 96759 : apbm = f0*rr(ax + 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
411 : ELSE
412 : apbm = 0.0_dp
413 : END IF
414 349183 : IF (ax > 0) THEN
415 96605 : ambp = f0*rr(ax - 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
416 : ELSE
417 : ambp = 0.0_dp
418 : END IF
419 349183 : IF (ax > 0 .AND. bx > 0) THEN
420 29337 : ambm = f0*rr(ax - 1, bx - 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
421 : ELSE
422 : ambm = 0.0_dp
423 : END IF
424 : ddab(ia, ib, 1) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(bx, dp)*apbm &
425 349183 : + 2.0_dp*b*REAL(ax, dp)*ambp - REAL(ax, dp)*REAL(bx, dp)*ambm
426 : ! dx dy
427 349183 : apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
428 349183 : IF (by > 0) THEN
429 96759 : apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
430 : ELSE
431 : apbm = 0.0_dp
432 : END IF
433 349183 : IF (ax > 0) THEN
434 96605 : ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
435 : ELSE
436 : ambp = 0.0_dp
437 : END IF
438 349183 : IF (ax > 0 .AND. by > 0) THEN
439 29337 : ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by - 1, 2)*rr(az, bz, 3)
440 : ELSE
441 : ambm = 0.0_dp
442 : END IF
443 : ddab(ia, ib, 2) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(by, dp)*apbm &
444 349183 : + 2.0_dp*b*REAL(ax, dp)*ambp - REAL(ax, dp)*REAL(by, dp)*ambm
445 : ! dx dz
446 349183 : apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
447 349183 : IF (bz > 0) THEN
448 96759 : apbm = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
449 : ELSE
450 : apbm = 0.0_dp
451 : END IF
452 349183 : IF (ax > 0) THEN
453 96605 : ambp = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
454 : ELSE
455 : ambp = 0.0_dp
456 : END IF
457 349183 : IF (ax > 0 .AND. bz > 0) THEN
458 29337 : ambm = f0*rr(ax - 1, bx, 1)*rr(ay, by, 2)*rr(az, bz - 1, 3)
459 : ELSE
460 : ambm = 0.0_dp
461 : END IF
462 : ddab(ia, ib, 3) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(bz, dp)*apbm &
463 349183 : + 2.0_dp*b*REAL(ax, dp)*ambp - REAL(ax, dp)*REAL(bz, dp)*ambm
464 : ! dy dy
465 349183 : apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by + 1, 2)*rr(az, bz, 3)
466 349183 : IF (by > 0) THEN
467 96759 : apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by - 1, 2)*rr(az, bz, 3)
468 : ELSE
469 : apbm = 0.0_dp
470 : END IF
471 349183 : IF (ay > 0) THEN
472 96605 : ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by + 1, 2)*rr(az, bz, 3)
473 : ELSE
474 : ambp = 0.0_dp
475 : END IF
476 349183 : IF (ay > 0 .AND. by > 0) THEN
477 29337 : ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by - 1, 2)*rr(az, bz, 3)
478 : ELSE
479 : ambm = 0.0_dp
480 : END IF
481 : ddab(ia, ib, 4) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(by, dp)*apbm &
482 349183 : + 2.0_dp*b*REAL(ay, dp)*ambp - REAL(ay, dp)*REAL(by, dp)*ambm
483 : ! dy dz
484 349183 : apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz + 1, 3)
485 349183 : IF (bz > 0) THEN
486 96759 : apbm = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz - 1, 3)
487 : ELSE
488 : apbm = 0.0_dp
489 : END IF
490 349183 : IF (ay > 0) THEN
491 96605 : ambp = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz + 1, 3)
492 : ELSE
493 : ambp = 0.0_dp
494 : END IF
495 349183 : IF (ay > 0 .AND. bz > 0) THEN
496 29337 : ambm = f0*rr(ax, bx, 1)*rr(ay - 1, by, 2)*rr(az, bz - 1, 3)
497 : ELSE
498 : ambm = 0.0_dp
499 : END IF
500 : ddab(ia, ib, 5) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(bz, dp)*apbm &
501 349183 : + 2.0_dp*b*REAL(ay, dp)*ambp - REAL(ay, dp)*REAL(bz, dp)*ambm
502 : ! dz dz
503 349183 : apbp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz + 1, 3)
504 349183 : IF (bz > 0) THEN
505 96759 : apbm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz - 1, 3)
506 : ELSE
507 : apbm = 0.0_dp
508 : END IF
509 349183 : IF (az > 0) THEN
510 96605 : ambp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz + 1, 3)
511 : ELSE
512 : ambp = 0.0_dp
513 : END IF
514 349183 : IF (az > 0 .AND. bz > 0) THEN
515 29337 : ambm = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az - 1, bz - 1, 3)
516 : ELSE
517 : ambm = 0.0_dp
518 : END IF
519 : ddab(ia, ib, 6) = -4.0_dp*a*b*apbp + 2.0_dp*a*REAL(bz, dp)*apbm &
520 349183 : + 2.0_dp*b*REAL(az, dp)*ambp - REAL(az, dp)*REAL(bz, dp)*ambm
521 : END IF
522 : !
523 : END DO
524 : END DO
525 : END DO !la
526 : END DO
527 : END DO
528 : END DO !lb
529 :
530 222697068 : mb = mb + nb
531 : END DO
532 112919241 : ma = ma + na
533 : END DO
534 :
535 18651647 : IF (.NOT. PRESENT(rr_work)) DEALLOCATE (rr)
536 : NULLIFY (rr)
537 :
538 18651647 : END SUBROUTINE overlap_ab
539 :
540 : ! **************************************************************************************************
541 : !> \brief Calculation of the two-center overlap integrals [aa|b] over
542 : !> Cartesian Gaussian-type functions.
543 : !> \param la1_max Max L on center A (basis 1)
544 : !> \param la1_min Min L on center A (basis 1)
545 : !> \param npgfa1 Number of primitives on center A (basis 1)
546 : !> \param rpgfa1 Range of functions on A, used for screening (basis 1)
547 : !> \param zeta1 Exponents on center A (basis 1)
548 : !> \param la2_max Max L on center A (basis 2)
549 : !> \param la2_min Min L on center A (basis 2)
550 : !> \param npgfa2 Number of primitives on center A (basis 2)
551 : !> \param rpgfa2 Range of functions on A, used for screening (basis 2)
552 : !> \param zeta2 Exponents on center A (basis 2)
553 : !> \param lb_max Max L on center B
554 : !> \param lb_min Min L on center B
555 : !> \param npgfb Number of primitives on center B
556 : !> \param rpgfb Range of functions on B, used for screening
557 : !> \param zetb Exponents on center B
558 : !> \param rab Distance vector A-B
559 : !> \param saab Final overlap integrals
560 : !> \param daab First derivative overlap integrals
561 : !> \param saba Final overlap integrals; different order
562 : !> \param daba First derivative overlap integrals; different order
563 : !> \date 01.07.2014
564 : !> \author JGH
565 : ! **************************************************************************************************
566 11331 : SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
567 22662 : la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
568 22662 : lb_max, lb_min, npgfb, rpgfb, zetb, &
569 11331 : rab, saab, daab, saba, daba)
570 : INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
571 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
572 : INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
573 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
574 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
575 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
576 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
577 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
578 : OPTIONAL :: saab
579 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
580 : INTENT(INOUT), OPTIONAL :: daab
581 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
582 : OPTIONAL :: saba
583 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
584 : INTENT(INOUT), OPTIONAL :: daba
585 :
586 : INTEGER :: ax, ax1, ax2, ay, ay1, ay2, az, az1, az2, bx, by, bz, coa1, coa2, cob, i1pgf, &
587 : i2pgf, ia1, ia2, ib, jpgf, la1, la2, lb, ldrr, lma, lmb, ma1, ma2, mb, na1, na2, nb, &
588 : ofa1, ofa2, ofb
589 : REAL(KIND=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfa, &
590 : tab, xhi, zet
591 11331 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
592 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
593 :
594 : ! Distance of the centers a and b
595 :
596 11331 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
597 11331 : tab = SQRT(rab2)
598 :
599 : ! Maximum l for auxiliary integrals
600 11331 : CPASSERT(PRESENT(saab) .OR. PRESENT(daab) .OR. PRESENT(saba) .OR. PRESENT(daba))
601 11331 : IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
602 11203 : lma = la1_max + la2_max
603 11203 : lmb = lb_max
604 : END IF
605 11331 : IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
606 3525 : lma = la1_max + la2_max + 1
607 3525 : lmb = lb_max
608 : END IF
609 11331 : ldrr = MAX(lma, lmb) + 1
610 :
611 : ! Allocate space for auxiliary integrals
612 56655 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
613 :
614 : ! Number of integrals, check size of arrays
615 11331 : ofa1 = ncoset(la1_min - 1)
616 11331 : ofa2 = ncoset(la2_min - 1)
617 11331 : ofb = ncoset(lb_min - 1)
618 11331 : na1 = ncoset(la1_max) - ofa1
619 11331 : na2 = ncoset(la2_max) - ofa2
620 11331 : nb = ncoset(lb_max) - ofb
621 11331 : IF (PRESENT(saab)) THEN
622 3206 : CPASSERT((SIZE(saab, 1) >= na1*npgfa1))
623 3206 : CPASSERT((SIZE(saab, 2) >= na2*npgfa2))
624 3206 : CPASSERT((SIZE(saab, 3) >= nb*npgfb))
625 : END IF
626 11331 : IF (PRESENT(daab)) THEN
627 128 : CPASSERT((SIZE(daab, 1) >= na1*npgfa1))
628 128 : CPASSERT((SIZE(daab, 2) >= na2*npgfa2))
629 128 : CPASSERT((SIZE(daab, 3) >= nb*npgfb))
630 128 : CPASSERT((SIZE(daab, 4) >= 3))
631 : END IF
632 11331 : IF (PRESENT(saba)) THEN
633 7997 : CPASSERT((SIZE(saba, 1) >= na1*npgfa1))
634 7997 : CPASSERT((SIZE(saba, 2) >= nb*npgfb))
635 7997 : CPASSERT((SIZE(saba, 3) >= na2*npgfa2))
636 : END IF
637 11331 : IF (PRESENT(daba)) THEN
638 3397 : CPASSERT((SIZE(daba, 1) >= na1*npgfa1))
639 3397 : CPASSERT((SIZE(daba, 2) >= nb*npgfb))
640 3397 : CPASSERT((SIZE(daba, 3) >= na2*npgfa2))
641 3397 : CPASSERT((SIZE(daba, 4) >= 3))
642 : END IF
643 :
644 : ! Loops over all primitive Gaussian-type functions
645 11331 : ma1 = 0
646 82142 : DO i1pgf = 1, npgfa1
647 70811 : ma2 = 0
648 351866 : DO i2pgf = 1, npgfa2
649 281055 : rpgfa = MIN(rpgfa1(i1pgf), rpgfa2(i2pgf))
650 281055 : mb = 0
651 1444764 : DO jpgf = 1, npgfb
652 : ! Distance Screening
653 1163709 : IF (rpgfa + rpgfb(jpgf) < tab) THEN
654 251494 : IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
655 251494 : IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
656 12150124 : IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
657 15308095 : IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
658 251494 : mb = mb + nb
659 251494 : CYCLE
660 : END IF
661 :
662 : ! Calculate some prefactors
663 912215 : a = zeta1(i1pgf) + zeta2(i2pgf)
664 912215 : b = zetb(jpgf)
665 912215 : zet = a + b
666 912215 : xhi = a*b/zet
667 3648860 : rap = b*rab/zet
668 3648860 : rbp = -a*rab/zet
669 :
670 : ! [ss|s] integral
671 912215 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
672 :
673 : ! Calculate the recurrence relation
674 912215 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
675 :
676 1874255 : DO lb = lb_min, lb_max
677 3399677 : DO bx = 0, lb
678 4817409 : DO by = 0, lb - bx
679 2329947 : bz = lb - bx - by
680 2329947 : cob = coset(bx, by, bz) - ofb
681 2329947 : ib = mb + cob
682 6743495 : DO la2 = la2_min, la2_max
683 11694698 : DO ax2 = 0, la2
684 22200024 : DO ay2 = 0, la2 - ax2
685 12835273 : az2 = la2 - ax2 - ay2
686 12835273 : coa2 = coset(ax2, ay2, az2) - ofa2
687 12835273 : ia2 = ma2 + coa2
688 36018542 : DO la1 = la1_min, la1_max
689 55660182 : DO ax1 = 0, la1
690 80322185 : DO ay1 = 0, la1 - ax1
691 37497276 : az1 = la1 - ax1 - ay1
692 37497276 : coa1 = coset(ax1, ay1, az1) - ofa1
693 37497276 : ia1 = ma1 + coa1
694 : ! integrals
695 37497276 : IF (PRESENT(saab)) THEN
696 1475300 : saab(ia1, ia2, ib) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
697 : END IF
698 37497276 : IF (PRESENT(saba)) THEN
699 35990684 : saba(ia1, ib, ia2) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
700 : END IF
701 : ! first derivatives
702 63615541 : IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
703 19863985 : ax = ax1 + ax2
704 19863985 : ay = ay1 + ay2
705 19863985 : az = az1 + az2
706 : ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
707 : ! dx
708 19863985 : dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
709 19863985 : IF (ax > 0) dumx = dumx - REAL(ax, dp)*rr(ax - 1, bx, 1)
710 19863985 : dumx = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
711 : ! dy
712 19863985 : dumy = 2.0_dp*a*rr(ay + 1, by, 2)
713 19863985 : IF (ay > 0) dumy = dumy - REAL(ay, dp)*rr(ay - 1, by, 2)
714 19863985 : dumy = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
715 : ! dz
716 19863985 : dumz = 2.0_dp*a*rr(az + 1, bz, 3)
717 19863985 : IF (az > 0) dumz = dumz - REAL(az, dp)*rr(az - 1, bz, 3)
718 19863985 : dumz = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
719 19863985 : IF (PRESENT(daab)) THEN
720 31292 : daab(ia1, ia2, ib, 1) = dumx
721 31292 : daab(ia1, ia2, ib, 2) = dumy
722 31292 : daab(ia1, ia2, ib, 3) = dumz
723 : END IF
724 19863985 : IF (PRESENT(daba)) THEN
725 19832693 : daba(ia1, ib, ia2, 1) = dumx
726 19832693 : daba(ia1, ib, ia2, 2) = dumy
727 19832693 : daba(ia1, ib, ia2, 3) = dumz
728 : END IF
729 : END IF
730 : !
731 : END DO
732 : END DO
733 : END DO !la1
734 : END DO
735 : END DO
736 : END DO !la2
737 : END DO
738 : END DO
739 : END DO !lb
740 :
741 1193270 : mb = mb + nb
742 : END DO
743 351866 : ma2 = ma2 + na2
744 : END DO
745 82142 : ma1 = ma1 + na1
746 : END DO
747 :
748 11331 : DEALLOCATE (rr)
749 :
750 11331 : END SUBROUTINE overlap_aab
751 :
752 : ! **************************************************************************************************
753 : !> \brief Calculation of the two-center overlap integrals [a|bb] over
754 : !> Cartesian Gaussian-type functions.
755 : !> \param la_max Max L on center A
756 : !> \param la_min Min L on center A
757 : !> \param npgfa Number of primitives on center A
758 : !> \param rpgfa Range of functions on A, used for screening
759 : !> \param zeta Exponents on center A
760 : !> \param lb1_max Max L on center B (basis 1)
761 : !> \param lb1_min Min L on center B (basis 1)
762 : !> \param npgfb1 Number of primitives on center B (basis 1)
763 : !> \param rpgfb1 Range of functions on B, used for screening (basis 1)
764 : !> \param zetb1 Exponents on center B (basis 1)
765 : !> \param lb2_max Max L on center B (basis 2)
766 : !> \param lb2_min Min L on center B (basis 2)
767 : !> \param npgfb2 Number of primitives on center B (basis 2)
768 : !> \param rpgfb2 Range of functions on B, used for screening (basis 2)
769 : !> \param zetb2 Exponents on center B (basis 2)
770 : !> \param rab Distance vector A-B
771 : !> \param sabb Final overlap integrals
772 : !> \param dabb First derivative overlap integrals
773 : !> \date 01.07.2014
774 : !> \author JGH
775 : ! **************************************************************************************************
776 7997 : SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
777 15994 : lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
778 15994 : lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
779 7997 : rab, sabb, dabb)
780 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
781 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
782 : INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
783 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
784 : INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
785 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
786 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
787 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
788 : OPTIONAL :: sabb
789 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
790 : INTENT(INOUT), OPTIONAL :: dabb
791 :
792 : INTEGER :: ax, ay, az, bx, bx1, bx2, by, by1, by2, bz, bz1, bz2, coa, cob1, cob2, ia, ib1, &
793 : ib2, ipgf, j1pgf, j2pgf, la, lb1, lb2, ldrr, lma, lmb, ma, mb1, mb2, na, nb1, nb2, ofa, &
794 : ofb1, ofb2
795 : REAL(KIND=dp) :: a, b, dumx, dumy, dumz, f0, rab2, rpgfb, &
796 : tab, xhi, zet
797 7997 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
798 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
799 :
800 : ! Distance of the centers a and b
801 :
802 7997 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
803 7997 : tab = SQRT(rab2)
804 :
805 : ! Maximum l for auxiliary integrals
806 7997 : CPASSERT(PRESENT(sabb) .OR. PRESENT(dabb))
807 7997 : IF (PRESENT(sabb)) THEN
808 7997 : lma = la_max
809 7997 : lmb = lb1_max + lb2_max
810 : END IF
811 7997 : IF (PRESENT(dabb)) THEN
812 3397 : lma = la_max + 1
813 3397 : lmb = lb1_max + lb2_max
814 : END IF
815 7997 : ldrr = MAX(lma, lmb) + 1
816 :
817 : ! Allocate space for auxiliary integrals
818 39985 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
819 :
820 : ! Number of integrals, check size of arrays
821 7997 : ofa = ncoset(la_min - 1)
822 7997 : ofb1 = ncoset(lb1_min - 1)
823 7997 : ofb2 = ncoset(lb2_min - 1)
824 7997 : na = ncoset(la_max) - ofa
825 7997 : nb1 = ncoset(lb1_max) - ofb1
826 7997 : nb2 = ncoset(lb2_max) - ofb2
827 7997 : IF (PRESENT(sabb)) THEN
828 7997 : CPASSERT((SIZE(sabb, 1) >= na*npgfa))
829 7997 : CPASSERT((SIZE(sabb, 2) >= nb1*npgfb1))
830 7997 : CPASSERT((SIZE(sabb, 3) >= nb2*npgfb2))
831 : END IF
832 7997 : IF (PRESENT(dabb)) THEN
833 3397 : CPASSERT((SIZE(dabb, 1) >= na*npgfa))
834 3397 : CPASSERT((SIZE(dabb, 2) >= nb1*npgfb1))
835 3397 : CPASSERT((SIZE(dabb, 3) >= nb2*npgfb2))
836 3397 : CPASSERT((SIZE(dabb, 4) >= 3))
837 : END IF
838 :
839 : ! Loops over all pairs of primitive Gaussian-type functions
840 7997 : ma = 0
841 60026 : DO ipgf = 1, npgfa
842 52029 : mb1 = 0
843 392532 : DO j1pgf = 1, npgfb1
844 340503 : mb2 = 0
845 1393966 : DO j2pgf = 1, npgfb2
846 : ! Distance Screening
847 1053463 : rpgfb = MIN(rpgfb1(j1pgf), rpgfb2(j2pgf))
848 1053463 : IF (rpgfa(ipgf) + rpgfb < tab) THEN
849 11929000 : IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
850 15070032 : IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
851 253218 : mb2 = mb2 + nb2
852 253218 : CYCLE
853 : END IF
854 :
855 : ! Calculate some prefactors
856 800245 : a = zeta(ipgf)
857 800245 : b = zetb1(j1pgf) + zetb2(j2pgf)
858 800245 : zet = a + b
859 800245 : xhi = a*b/zet
860 3200980 : rap = b*rab/zet
861 3200980 : rbp = -a*rab/zet
862 :
863 : ! [s|s] integral
864 800245 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
865 :
866 : ! Calculate the recurrence relation
867 800245 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
868 :
869 1706310 : DO lb2 = lb2_min, lb2_max
870 3765333 : DO bx2 = 0, lb2
871 7068414 : DO by2 = 0, lb2 - bx2
872 4103326 : bz2 = lb2 - bx2 - by2
873 4103326 : cob2 = coset(bx2, by2, bz2) - ofb2
874 4103326 : ib2 = mb2 + cob2
875 11070029 : DO lb1 = lb1_min, lb1_max
876 17121130 : DO bx1 = 0, lb1
877 25341538 : DO by1 = 0, lb1 - bx1
878 12323734 : bz1 = lb1 - bx1 - by1
879 12323734 : cob1 = coset(bx1, by1, bz1) - ofb1
880 12323734 : ib1 = mb1 + cob1
881 36559042 : DO la = la_min, la_max
882 53606848 : DO ax = 0, la
883 77262690 : DO ay = 0, la - ax
884 35979576 : az = la - ax - ay
885 35979576 : coa = coset(ax, ay, az) - ofa
886 35979576 : ia = ma + coa
887 : ! integrals
888 35979576 : IF (PRESENT(sabb)) THEN
889 35979576 : sabb(ia, ib1, ib2) = f0*rr(ax, bx1 + bx2, 1)*rr(ay, by1 + by2, 2)*rr(az, bz1 + bz2, 3)
890 : END IF
891 : ! first derivatives
892 61137506 : IF (PRESENT(dabb)) THEN
893 19827139 : bx = bx1 + bx2
894 19827139 : by = by1 + by2
895 19827139 : bz = bz1 + bz2
896 : ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
897 : ! dx
898 19827139 : dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
899 19827139 : IF (ax > 0) dumx = dumx - REAL(ax, dp)*rr(ax - 1, bx, 1)
900 19827139 : dabb(ia, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
901 : ! dy
902 19827139 : dumy = 2.0_dp*a*rr(ay + 1, by, 2)
903 19827139 : IF (ay > 0) dumy = dumy - REAL(ay, dp)*rr(ay - 1, by, 2)
904 19827139 : dabb(ia, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
905 : ! dz
906 19827139 : dumz = 2.0_dp*a*rr(az + 1, bz, 3)
907 19827139 : IF (az > 0) dumz = dumz - REAL(az, dp)*rr(az - 1, bz, 3)
908 19827139 : dabb(ia, ib1, ib2, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
909 : END IF
910 : !
911 : END DO
912 : END DO
913 : END DO !la
914 : END DO
915 : END DO
916 : END DO !lb1
917 : END DO
918 : END DO
919 : END DO !lb2
920 :
921 1140748 : mb2 = mb2 + nb2
922 : END DO
923 392532 : mb1 = mb1 + nb1
924 : END DO
925 60026 : ma = ma + na
926 : END DO
927 :
928 7997 : DEALLOCATE (rr)
929 :
930 7997 : END SUBROUTINE overlap_abb
931 :
932 : ! **************************************************************************************************
933 :
934 : ! **************************************************************************************************
935 : !> \brief Calculation of the two-center overlap integrals [a|b] over
936 : !> Spherical Gaussian-type functions.
937 : !> \param la Max L on center A
938 : !> \param zeta Exponents on center A
939 : !> \param lb Max L on center B
940 : !> \param zetb Exponents on center B
941 : !> \param rab Distance vector A-B
942 : !> \param sab Final overlap integrals
943 : !> \date 01.03.2016
944 : !> \author JGH
945 : ! **************************************************************************************************
946 275346 : SUBROUTINE overlap_ab_s(la, zeta, lb, zetb, rab, sab)
947 : INTEGER, INTENT(IN) :: la
948 : REAL(KIND=dp), INTENT(IN) :: zeta
949 : INTEGER, INTENT(IN) :: lb
950 : REAL(KIND=dp), INTENT(IN) :: zetb
951 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
952 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
953 :
954 : REAL(KIND=dp), PARAMETER :: huge4 = HUGE(1._dp)/4._dp
955 :
956 : INTEGER :: nca, ncb, nsa, nsb
957 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
958 : REAL(KIND=dp), DIMENSION(1) :: rpgf, za, zb
959 275346 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
960 :
961 275346 : rpgf(1) = huge4
962 275346 : za(1) = zeta
963 275346 : zb(1) = zetb
964 :
965 275346 : nca = nco(la)
966 275346 : ncb = nco(lb)
967 1101384 : ALLOCATE (cab(nca, ncb))
968 275346 : nsa = nso(la)
969 275346 : nsb = nso(lb)
970 :
971 275346 : CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
972 :
973 275346 : c2sa => orbtramat(la)%c2s
974 275346 : c2sb => orbtramat(lb)%c2s
975 275346 : sab(1:nsa, 1:nsb) = MATMUL(c2sa(1:nsa, 1:nca), &
976 11068223 : MATMUL(cab(1:nca, 1:ncb), TRANSPOSE(c2sb(1:nsb, 1:ncb))))
977 :
978 275346 : DEALLOCATE (cab)
979 :
980 275346 : END SUBROUTINE overlap_ab_s
981 :
982 : ! **************************************************************************************************
983 : !> \brief Calculation of the overlap integrals [a|b] over
984 : !> cubic periodic Spherical Gaussian-type functions.
985 : !> \param la Max L on center A
986 : !> \param zeta Exponents on center A
987 : !> \param lb Max L on center B
988 : !> \param zetb Exponents on center B
989 : !> \param alat Lattice constant
990 : !> \param sab Final overlap integrals
991 : !> \date 01.03.2016
992 : !> \author JGH
993 : ! **************************************************************************************************
994 36030 : SUBROUTINE overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
995 : INTEGER, INTENT(IN) :: la
996 : REAL(KIND=dp), INTENT(IN) :: zeta
997 : INTEGER, INTENT(IN) :: lb
998 : REAL(KIND=dp), INTENT(IN) :: zetb, alat
999 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
1000 :
1001 : COMPLEX(KIND=dp) :: zfg
1002 36030 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fun, gun
1003 : INTEGER :: ax, ay, az, bx, by, bz, i, ia, ib, l, &
1004 : l1, l2, na, nb, nca, ncb, nmax, nsa, &
1005 : nsb
1006 : REAL(KIND=dp) :: oa, ob, ovol, zm
1007 36030 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
1008 36030 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
1009 : REAL(KIND=dp), DIMENSION(0:3, 0:3) :: fgsum
1010 36030 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
1011 :
1012 36030 : nca = nco(la)
1013 36030 : ncb = nco(lb)
1014 144120 : ALLOCATE (cab(nca, ncb))
1015 36030 : cab = 0.0_dp
1016 36030 : nsa = nso(la)
1017 36030 : nsb = nso(lb)
1018 :
1019 36030 : zm = MIN(zeta, zetb)
1020 36030 : nmax = NINT(1.81_dp*alat*SQRT(zm) + 1.0_dp)
1021 : ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
1022 396330 : fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
1023 :
1024 36030 : oa = 1._dp/zeta
1025 36030 : ob = 1._dp/zetb
1026 521510 : DO i = -nmax, nmax
1027 485480 : gval(i) = twopi/alat*REAL(i, KIND=dp)
1028 485480 : fexp(i) = SQRT(oa*pi)*EXP(-0.25_dp*oa*gval(i)**2)
1029 521510 : gexp(i) = SQRT(ob*pi)*EXP(-0.25_dp*ob*gval(i)**2)
1030 : END DO
1031 91894 : DO l = 0, la
1032 91894 : IF (l == 0) THEN
1033 521510 : fun(:, l) = z_one
1034 19834 : ELSE IF (l == 1) THEN
1035 252786 : fun(:, l) = CMPLX(0.0_dp, 0.5_dp*oa*gval(:), KIND=dp)
1036 2073 : ELSE IF (l == 2) THEN
1037 29772 : fun(:, l) = CMPLX(-(0.5_dp*oa*gval(:))**2, 0.0_dp, KIND=dp)
1038 29772 : fun(:, l) = fun(:, l) + CMPLX(0.5_dp*oa, 0.0_dp, KIND=dp)
1039 0 : ELSE IF (l == 3) THEN
1040 0 : fun(:, l) = CMPLX(0.0_dp, -(0.5_dp*oa*gval(:))**3, KIND=dp)
1041 0 : fun(:, l) = fun(:, l) + CMPLX(0.0_dp, 0.75_dp*oa*oa*gval(:), KIND=dp)
1042 : ELSE
1043 0 : CPABORT("l value too high")
1044 : END IF
1045 : END DO
1046 91894 : DO l = 0, lb
1047 91894 : IF (l == 0) THEN
1048 521510 : gun(:, l) = z_one
1049 19834 : ELSE IF (l == 1) THEN
1050 252786 : gun(:, l) = CMPLX(0.0_dp, 0.5_dp*ob*gval(:), KIND=dp)
1051 2073 : ELSE IF (l == 2) THEN
1052 29772 : gun(:, l) = CMPLX(-(0.5_dp*ob*gval(:))**2, 0.0_dp, KIND=dp)
1053 29772 : gun(:, l) = gun(:, l) + CMPLX(0.5_dp*ob, 0.0_dp, KIND=dp)
1054 0 : ELSE IF (l == 3) THEN
1055 0 : gun(:, l) = CMPLX(0.0_dp, -(0.5_dp*ob*gval(:))**3, KIND=dp)
1056 0 : gun(:, l) = gun(:, l) + CMPLX(0.0_dp, 0.75_dp*ob*ob*gval(:), KIND=dp)
1057 : ELSE
1058 0 : CPABORT("l value too high")
1059 : END IF
1060 : END DO
1061 :
1062 36030 : fgsum = 0.0_dp
1063 91894 : DO l1 = 0, la
1064 179846 : DO l2 = 0, lb
1065 1259856 : zfg = SUM(CONJG(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
1066 143816 : fgsum(l1, l2) = REAL(zfg, KIND=dp)
1067 : END DO
1068 : END DO
1069 :
1070 36030 : na = ncoset(la - 1)
1071 36030 : nb = ncoset(lb - 1)
1072 91894 : DO ax = 0, la
1073 169665 : DO ay = 0, la - ax
1074 77771 : az = la - ax - ay
1075 77771 : ia = coset(ax, ay, az) - na
1076 257740 : DO bx = 0, lb
1077 379027 : DO by = 0, lb - bx
1078 177151 : bz = lb - bx - by
1079 177151 : ib = coset(bx, by, bz) - nb
1080 301256 : cab(ia, ib) = fgsum(ax, bx)*fgsum(ay, by)*fgsum(az, bz)
1081 : END DO
1082 : END DO
1083 : END DO
1084 : END DO
1085 :
1086 36030 : c2sa => orbtramat(la)%c2s
1087 36030 : c2sb => orbtramat(lb)%c2s
1088 36030 : sab(1:nsa, 1:nsb) = MATMUL(c2sa(1:nsa, 1:nca), &
1089 2481735 : MATMUL(cab(1:nca, 1:ncb), TRANSPOSE(c2sb(1:nsb, 1:ncb))))
1090 36030 : ovol = 1._dp/(alat**3)
1091 276110 : sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
1092 :
1093 36030 : DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
1094 :
1095 36030 : END SUBROUTINE overlap_ab_sp
1096 :
1097 311376 : END MODULE ai_overlap
|
