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_cartesian_operators, ONLY: apply_cartesian_1d,&
20 : apply_cartesian_operator_pair,&
21 : cartesian_ladder,&
22 : cartesian_ladder_map
23 : USE ai_os_rr, ONLY: os_rr_ovlp
24 : USE kinds, ONLY: dp
25 : USE mathconstants, ONLY: pi,&
26 : twopi,&
27 : z_one
28 : USE orbital_pointers, ONLY: coset,&
29 : indco,&
30 : nco,&
31 : ncoset,&
32 : nso
33 : USE orbital_transformation_matrices, ONLY: orbtramat
34 : #include "../base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap'
41 :
42 : ! *** Public subroutines ***
43 : PUBLIC :: overlap, overlap_ab, overlap_aab, overlap_ab_s, overlap_ab_sp, &
44 : overlap_abb
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Purpose: Calculation of the two-center overlap integrals [a|b] over
50 : !> Cartesian Gaussian-type functions.
51 : !> \param la_max_set Max L on center A
52 : !> \param la_min_set Min L on center A
53 : !> \param npgfa Number of primitives on center A
54 : !> \param rpgfa Range of functions on A, used for screening
55 : !> \param zeta Exponents on center A
56 : !> \param lb_max_set Max L on center B
57 : !> \param lb_min_set Min L on center B
58 : !> \param npgfb Number of primitives on center B
59 : !> \param rpgfb Range of functions on B, used for screening
60 : !> \param zetb Exponents on center B
61 : !> \param rab Distance vector A-B
62 : !> \param dab Distance A-B
63 : !> \param sab Final Integrals, basic and derivatives
64 : !> \param da_max_set Some additional derivative information
65 : !> \param return_derivatives Return integral derivatives
66 : !> \param s Work space
67 : !> \param lds Leading dimension of s
68 : !> \date 19.09.2000
69 : !> \author MK
70 : !> \version 1.0
71 : ! **************************************************************************************************
72 2055200 : SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
73 4110400 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
74 2055200 : rab, dab, sab, da_max_set, return_derivatives, s, lds)
75 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
76 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
77 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
78 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
79 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
80 : REAL(KIND=dp), INTENT(IN) :: dab
81 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
82 : INTEGER, INTENT(IN) :: da_max_set
83 : LOGICAL, INTENT(IN) :: return_derivatives
84 : INTEGER, INTENT(IN) :: lds
85 : REAL(KIND=dp), DIMENSION(lds, lds, *), &
86 : INTENT(INOUT) :: s
87 :
88 : INTEGER :: ax, axis, ay, az, bx, by, bz, cda, coa, cob, da, da_max, first, first_b, i, ipgf, &
89 : j, jk, jpgf, jstart, k, la_max, la_start, last, last_b, lb_max, lb_start, ldrr, na, nb, &
90 : parent, power(3)
91 2055200 : INTEGER, ALLOCATABLE :: derivative_map(:, :, :)
92 : REAL(KIND=dp) :: f0, zetp
93 2055200 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
94 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
95 :
96 2055200 : da_max = da_max_set
97 2055200 : la_max = la_max_set + da_max_set
98 :
99 2055200 : lb_max = lb_max_set
100 2055200 : ldrr = MAX(la_max, lb_max) + 1
101 10276000 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
102 2055200 : IF (da_max > 0) THEN
103 2557544 : ALLOCATE (derivative_map(4, ncoset(la_max - 1), 3))
104 2557544 : DO axis = 1, 3
105 2557544 : CALL cartesian_ladder_map(la_max - 1, 0, axis, derivative_map(:, :, axis))
106 : END DO
107 : END IF
108 :
109 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
110 :
111 2055200 : na = 0
112 10527612 : DO ipgf = 1, npgfa
113 :
114 8472412 : nb = 0
115 :
116 22873998 : DO jpgf = 1, npgfb
117 :
118 : ! *** Screening ***
119 :
120 14401586 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
121 51350752 : DO j = nb + 1, nb + ncoset(lb_max_set)
122 258651398 : DO i = na + 1, na + ncoset(la_max_set)
123 249293258 : sab(i, j) = 0.0_dp
124 : END DO
125 : END DO
126 9358140 : IF (return_derivatives) THEN
127 20027010 : DO k = 2, ncoset(da_max_set)
128 10668870 : jstart = (k - 1)*SIZE(sab, 1)
129 62639940 : DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
130 228142932 : DO i = na + 1, na + ncoset(la_max_set)
131 217474062 : sab(i, j) = 0.0_dp
132 : END DO
133 : END DO
134 : END DO
135 : END IF
136 9358140 : nb = nb + ncoset(lb_max_set)
137 9358140 : CYCLE
138 : END IF
139 :
140 : ! *** Calculate some prefactors ***
141 :
142 5043446 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
143 :
144 5043446 : f0 = SQRT((pi*zetp)**3)*EXP(-zeta(ipgf)*zetb(jpgf)*zetp*dab*dab)
145 20173784 : rap(:) = zetb(jpgf)*zetp*rab(:)
146 20173784 : rbp(:) = -zeta(ipgf)*zetp*rab(:)
147 :
148 5043446 : CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
149 :
150 27949398 : DO cob = 1, ncoset(lb_max)
151 22905952 : bx = indco(1, cob)
152 22905952 : by = indco(2, cob)
153 22905952 : bz = indco(3, cob)
154 195509494 : DO coa = 1, ncoset(la_max)
155 167560096 : ax = indco(1, coa)
156 167560096 : ay = indco(2, coa)
157 167560096 : az = indco(3, coa)
158 190466048 : s(coa, cob, 1) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
159 : END DO
160 : END DO
161 :
162 : ! *** Store the primitive overlap integrals ***
163 :
164 27949398 : DO j = 1, ncoset(lb_max_set)
165 151019308 : DO i = 1, ncoset(la_max_set)
166 145975862 : sab(na + i, nb + j) = s(i, j, 1)
167 : END DO
168 : END DO
169 :
170 : ! *** Calculate the requested derivatives with respect ***
171 : ! *** to the nuclear coordinates of the atomic center a ***
172 :
173 5043446 : IF (return_derivatives) THEN
174 : la_start = 0
175 : lb_start = 0
176 : ELSE
177 0 : la_start = la_min_set
178 0 : lb_start = lb_min_set
179 : END IF
180 :
181 5043446 : first_b = ncoset(lb_start - 1) + 1
182 5043446 : last_b = ncoset(lb_max_set)
183 6785847 : DO da = 1, da_max
184 : ! Intermediate shells must cover lowering in the remaining derivative steps.
185 1742401 : first = ncoset(MAX(0, la_start - da_max + da) - 1) + 1
186 1742401 : last = ncoset(la_max - da)
187 12014850 : DO cda = ncoset(da - 1) + 1, ncoset(da)
188 20916012 : power = indco(:, cda)
189 5229003 : axis = 3
190 5229003 : IF (power(3) == 0) axis = 2
191 5229003 : IF (power(axis) == 0) axis = 1
192 5229003 : power(axis) = power(axis) - 1
193 5229003 : parent = coset(power(1), power(2), power(3))
194 : CALL cartesian_ladder(s(:, first_b:last_b, parent), s(first:last, first_b:last_b, cda), &
195 6971404 : derivative_map(:, first:last, axis), .TRUE., zeta(ipgf), 0.0_dp, .FALSE.)
196 : END DO
197 : END DO
198 :
199 : ! *** Return all the calculated derivatives of the ***
200 : ! *** primitive overlap integrals, if requested ***
201 :
202 5043446 : IF (return_derivatives) THEN
203 10272449 : DO k = 2, ncoset(da_max_set)
204 5229003 : jstart = (k - 1)*SIZE(sab, 1)
205 31419164 : DO j = 1, ncoset(lb_max_set)
206 21146715 : jk = jstart + j
207 132083964 : DO i = 1, ncoset(la_max_set)
208 126854961 : sab(na + i, nb + jk) = s(i, j, k)
209 : END DO
210 : END DO
211 : END DO
212 : END IF
213 :
214 13515858 : nb = nb + ncoset(lb_max_set)
215 :
216 : END DO
217 :
218 10527612 : na = na + ncoset(la_max_set)
219 : END DO
220 :
221 2055200 : DEALLOCATE (rr)
222 :
223 2055200 : END SUBROUTINE overlap
224 :
225 : ! **************************************************************************************************
226 : !> \brief Calculation of the two-center overlap integrals [a|b] over
227 : !> Cartesian Gaussian-type functions. First and second derivatives
228 : !> \param la_max Max L on center A
229 : !> \param la_min Min L on center A
230 : !> \param npgfa Number of primitives on center A
231 : !> \param rpgfa Range of functions on A, used for screening
232 : !> \param zeta Exponents on center A
233 : !> \param lb_max Max L on center B
234 : !> \param lb_min Min L on center B
235 : !> \param npgfb Number of primitives on center B
236 : !> \param rpgfb Range of functions on B, used for screening
237 : !> \param zetb Exponents on center B
238 : !> \param rab Distance vector A-B
239 : !> \param sab Final overlap integrals
240 : !> \param dab First derivative overlap integrals
241 : !> \param ddab Second derivative overlap integrals ordered xx, xy, xz, yy, yz, zz
242 : !> \param rr_work Optional caller-provided one-dimensional overlap recurrence workspace
243 : !> \date 01.07.2014
244 : !> \author JGH
245 : ! **************************************************************************************************
246 18866245 : SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
247 18866245 : lb_max, lb_min, npgfb, rpgfb, zetb, &
248 18866245 : rab, sab, dab, ddab, rr_work)
249 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
250 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
251 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
252 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
253 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
254 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
255 : OPTIONAL :: sab
256 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
257 : OPTIONAL :: dab, ddab
258 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
259 : OPTIONAL, TARGET :: rr_work
260 :
261 : INTEGER, PARAMETER :: component(2, 1) = RESHAPE([1, 1], [2, 1]), &
262 : pairs(2, 6) = RESHAPE([1, 1, 1, 2, 1, 3, 2, 2, 2, 3, 3, 3], [2, 6])
263 :
264 : INTEGER :: avec(3), bvec(3), coa, cob, ia, ib, icomp, idir, ipgf, jdir, jpgf, kdir, ldrr, &
265 37732490 : lma, lmb, ma, mb, na, nb, ofa, ofb, map_b(4, lb_max + 1, 1), map_a(4, la_max + 1, 1)
266 : REAL(KIND=dp) :: a, b, f0, factor(3), rab2, sval(3), tab, &
267 : xhi, zet
268 18866245 : REAL(KIND=dp), ALLOCATABLE :: drr_a(:, :, :), drr_ab(:, :, :)
269 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
270 18866245 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rr
271 :
272 : ! Distance of the centers a and b
273 :
274 18866245 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
275 18866245 : tab = SQRT(rab2)
276 :
277 : ! Maximum l for auxiliary integrals
278 18866245 : CPASSERT(PRESENT(sab) .OR. PRESENT(dab) .OR. PRESENT(ddab))
279 18866245 : IF (PRESENT(sab)) THEN
280 18569413 : lma = la_max
281 18569413 : lmb = lb_max
282 : END IF
283 18866245 : IF (PRESENT(dab)) THEN
284 4966681 : lma = la_max + 1
285 4966681 : lmb = lb_max
286 : END IF
287 18866245 : IF (PRESENT(ddab)) THEN
288 13855 : lma = la_max + 1
289 13855 : lmb = lb_max + 1
290 : END IF
291 18866245 : ldrr = MAX(lma, lmb) + 1
292 :
293 : ! Allocate or attach the workspace for auxiliary integrals
294 : NULLIFY (rr)
295 18866245 : IF (PRESENT(rr_work)) THEN
296 157518 : CPASSERT(SIZE(rr_work) >= ldrr*ldrr*3)
297 157518 : rr(0:ldrr - 1, 0:ldrr - 1, 1:3) => rr_work(1:ldrr*ldrr*3)
298 : ELSE
299 93543635 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
300 : END IF
301 :
302 38732969 : IF (PRESENT(dab) .OR. PRESENT(ddab)) ALLOCATE (drr_a(0:ldrr - 1, 0:ldrr - 1, 3))
303 4966681 : IF (PRESENT(ddab)) THEN
304 41565 : ALLOCATE (drr_ab(0:ldrr - 1, 0:ldrr - 1, 3))
305 13855 : CALL cartesian_ladder_map(la_max, 0, 0, map_a(:, :, 1))
306 13855 : CALL cartesian_ladder_map(lb_max, 0, 0, map_b(:, :, 1))
307 : END IF
308 :
309 : ! Number of integrals, check size of arrays
310 18866245 : ofa = ncoset(la_min - 1)
311 18866245 : ofb = ncoset(lb_min - 1)
312 18866245 : na = ncoset(la_max) - ofa
313 18866245 : nb = ncoset(lb_max) - ofb
314 18866245 : IF (PRESENT(sab)) THEN
315 18569413 : CPASSERT((SIZE(sab, 1) >= na*npgfa))
316 18569413 : CPASSERT((SIZE(sab, 2) >= nb*npgfb))
317 : END IF
318 18866245 : IF (PRESENT(dab)) THEN
319 4966681 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
320 4966681 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
321 4966681 : CPASSERT((SIZE(dab, 3) >= 3))
322 : END IF
323 18866245 : IF (PRESENT(ddab)) THEN
324 13855 : CPASSERT((SIZE(ddab, 1) >= na*npgfa))
325 13855 : CPASSERT((SIZE(ddab, 2) >= nb*npgfb))
326 13855 : CPASSERT((SIZE(ddab, 3) >= 6))
327 : END IF
328 :
329 : ! Loops over all pairs of primitive Gaussian-type functions
330 18866245 : ma = 0
331 114179650 : DO ipgf = 1, npgfa
332 95313405 : mb = 0
333 648470381 : DO jpgf = 1, npgfb
334 : ! Distance Screening
335 553156976 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
336 4627628707 : IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
337 3587358114 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
338 430362207 : IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
339 423653793 : mb = mb + nb
340 423653793 : CYCLE
341 : END IF
342 :
343 : ! Calculate some prefactors
344 129503183 : a = zeta(ipgf)
345 129503183 : b = zetb(jpgf)
346 129503183 : zet = a + b
347 129503183 : xhi = a*b/zet
348 518012732 : rap = b*rab/zet
349 518012732 : rbp = -a*rab/zet
350 :
351 : ! [s|s] integral
352 129503183 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
353 :
354 : ! Calculate the recurrence relation
355 129503183 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
356 129503183 : IF (PRESENT(ddab)) THEN
357 223420 : DO idir = 1, 3
358 : CALL apply_cartesian_operator_pair(rr(:, :, idir), drr_ab(0:la_max, 0:lb_max, idir:idir), &
359 : map_a, map_b, component, [.TRUE., .TRUE.], [a, b], &
360 558550 : [.FALSE., .TRUE.], drr_a(0:la_max, 0:lmb, idir))
361 : END DO
362 129447328 : ELSE IF (PRESENT(dab)) THEN
363 32109437 : CALL apply_cartesian_1d(rr, drr_a, la_max, lmb, .FALSE., .TRUE., a, [0.0_dp, 0.0_dp, 0.0_dp])
364 : END IF
365 :
366 465976551 : DO cob = ofb + 1, ncoset(lb_max)
367 1345893472 : bvec = indco(:, cob)
368 336473368 : ib = mb + cob - ofb
369 1534507102 : DO coa = ofa + 1, ncoset(la_max)
370 4274122204 : avec = indco(:, coa)
371 1068530551 : ia = ma + coa - ofa
372 4274122204 : sval = [(rr(avec(idir), bvec(idir), idir), idir=1, 3)]
373 : ! integrals
374 1068530551 : IF (PRESENT(sab)) THEN
375 1063871290 : sab(ia, ib) = f0*sval(1)*sval(2)*sval(3)
376 : END IF
377 : ! first derivatives
378 1068530551 : IF (PRESENT(dab)) THEN
379 867492952 : DO idir = 1, 3
380 650619714 : jdir = MOD(idir, 3) + 1
381 650619714 : kdir = MOD(idir + 1, 3) + 1
382 867492952 : dab(ia, ib, idir) = f0*drr_a(avec(idir), bvec(idir), idir)*sval(jdir)*sval(kdir)
383 : END DO
384 : END IF
385 : ! Translational invariance gives D_B S = -D_A S and D_A D_A S = -D_A D_B S.
386 1405003919 : IF (PRESENT(ddab)) THEN
387 2444281 : DO icomp = 1, 6
388 2095098 : idir = pairs(1, icomp)
389 2095098 : jdir = pairs(2, icomp)
390 2095098 : factor = sval
391 2095098 : IF (idir == jdir) THEN
392 1047549 : factor(idir) = -drr_ab(avec(idir), bvec(idir), idir)
393 : ELSE
394 1047549 : factor(idir) = drr_a(avec(idir), bvec(idir), idir)
395 1047549 : factor(jdir) = drr_a(avec(jdir), bvec(jdir), jdir)
396 : END IF
397 2444281 : ddab(ia, ib, icomp) = f0*factor(1)*factor(2)*factor(3)
398 : END DO
399 : END IF
400 : !
401 : END DO
402 : END DO
403 :
404 224816588 : mb = mb + nb
405 : END DO
406 114179650 : ma = ma + na
407 : END DO
408 :
409 18866245 : IF (.NOT. PRESENT(rr_work)) DEALLOCATE (rr)
410 : NULLIFY (rr)
411 :
412 18866245 : END SUBROUTINE overlap_ab
413 :
414 : ! **************************************************************************************************
415 : !> \brief Calculation of the two-center overlap integrals [aa|b] over
416 : !> Cartesian Gaussian-type functions.
417 : !> \param la1_max Max L on center A (basis 1)
418 : !> \param la1_min Min L on center A (basis 1)
419 : !> \param npgfa1 Number of primitives on center A (basis 1)
420 : !> \param rpgfa1 Range of functions on A, used for screening (basis 1)
421 : !> \param zeta1 Exponents on center A (basis 1)
422 : !> \param la2_max Max L on center A (basis 2)
423 : !> \param la2_min Min L on center A (basis 2)
424 : !> \param npgfa2 Number of primitives on center A (basis 2)
425 : !> \param rpgfa2 Range of functions on A, used for screening (basis 2)
426 : !> \param zeta2 Exponents on center A (basis 2)
427 : !> \param lb_max Max L on center B
428 : !> \param lb_min Min L on center B
429 : !> \param npgfb Number of primitives on center B
430 : !> \param rpgfb Range of functions on B, used for screening
431 : !> \param zetb Exponents on center B
432 : !> \param rab Distance vector A-B
433 : !> \param saab Final overlap integrals
434 : !> \param daab First derivatives with respect to the shared center of both A functions
435 : !> \param saba Final overlap integrals; different order
436 : !> \param daba Same shared-center derivatives as daab, with A-B-A output ordering
437 : !> \date 01.07.2014
438 : !> \author JGH
439 : ! **************************************************************************************************
440 9695 : SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
441 19390 : la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
442 19390 : lb_max, lb_min, npgfb, rpgfb, zetb, &
443 9695 : rab, saab, daab, saba, daba)
444 : INTEGER, INTENT(IN) :: la1_max, la1_min, npgfa1
445 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa1, zeta1
446 : INTEGER, INTENT(IN) :: la2_max, la2_min, npgfa2
447 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa2, zeta2
448 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
449 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
450 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
451 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
452 : OPTIONAL :: saab
453 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
454 : INTENT(INOUT), OPTIONAL :: daab
455 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
456 : OPTIONAL :: saba
457 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
458 : INTENT(INOUT), OPTIONAL :: daba
459 :
460 : INTEGER :: coa1, coa2, cob, i1pgf, i2pgf, ia1, ia2, &
461 : ib, jpgf, ldrr, lma, lmb, ma1, ma2, &
462 : mb, na1, na2, nb, ofa1, ofa2, ofb
463 : INTEGER, DIMENSION(3) :: a2vec, avec, bvec
464 : REAL(KIND=dp) :: a, b, dum(3), f0, rab2, rpgfa, tab, &
465 : value, xhi, zet
466 9695 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drr, rr
467 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
468 :
469 : ! Distance of the centers a and b
470 :
471 9695 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
472 9695 : tab = SQRT(rab2)
473 :
474 : ! Maximum l for auxiliary integrals
475 9695 : CPASSERT(PRESENT(saab) .OR. PRESENT(daab) .OR. PRESENT(saba) .OR. PRESENT(daba))
476 9695 : IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
477 9567 : lma = la1_max + la2_max
478 9567 : lmb = lb_max
479 : END IF
480 9695 : IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
481 3525 : lma = la1_max + la2_max + 1
482 3525 : lmb = lb_max
483 : END IF
484 9695 : ldrr = MAX(lma, lmb) + 1
485 :
486 : ! Allocate space for auxiliary integrals
487 48475 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
488 23795 : IF (PRESENT(daab) .OR. PRESENT(daba)) ALLOCATE (drr(0:lma - 1, 0:lmb, 3))
489 :
490 : ! Number of integrals, check size of arrays
491 9695 : ofa1 = ncoset(la1_min - 1)
492 9695 : ofa2 = ncoset(la2_min - 1)
493 9695 : ofb = ncoset(lb_min - 1)
494 9695 : na1 = ncoset(la1_max) - ofa1
495 9695 : na2 = ncoset(la2_max) - ofa2
496 9695 : nb = ncoset(lb_max) - ofb
497 9695 : IF (PRESENT(saab)) THEN
498 1570 : CPASSERT((SIZE(saab, 1) >= na1*npgfa1))
499 1570 : CPASSERT((SIZE(saab, 2) >= na2*npgfa2))
500 1570 : CPASSERT((SIZE(saab, 3) >= nb*npgfb))
501 : END IF
502 9695 : IF (PRESENT(daab)) THEN
503 128 : CPASSERT((SIZE(daab, 1) >= na1*npgfa1))
504 128 : CPASSERT((SIZE(daab, 2) >= na2*npgfa2))
505 128 : CPASSERT((SIZE(daab, 3) >= nb*npgfb))
506 128 : CPASSERT((SIZE(daab, 4) >= 3))
507 : END IF
508 9695 : IF (PRESENT(saba)) THEN
509 7997 : CPASSERT((SIZE(saba, 1) >= na1*npgfa1))
510 7997 : CPASSERT((SIZE(saba, 2) >= nb*npgfb))
511 7997 : CPASSERT((SIZE(saba, 3) >= na2*npgfa2))
512 : END IF
513 9695 : IF (PRESENT(daba)) THEN
514 3397 : CPASSERT((SIZE(daba, 1) >= na1*npgfa1))
515 3397 : CPASSERT((SIZE(daba, 2) >= nb*npgfb))
516 3397 : CPASSERT((SIZE(daba, 3) >= na2*npgfa2))
517 3397 : CPASSERT((SIZE(daba, 4) >= 3))
518 : END IF
519 :
520 : ! Loops over all primitive Gaussian-type functions
521 9695 : ma1 = 0
522 72326 : DO i1pgf = 1, npgfa1
523 62631 : ma2 = 0
524 302786 : DO i2pgf = 1, npgfa2
525 240155 : rpgfa = MIN(rpgfa1(i1pgf), rpgfa2(i2pgf))
526 240155 : mb = 0
527 1362964 : DO jpgf = 1, npgfb
528 : ! Distance Screening
529 1122809 : IF (rpgfa + rpgfb(jpgf) < tab) THEN
530 251494 : IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
531 251494 : IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
532 12150124 : IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
533 15308095 : IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
534 251494 : mb = mb + nb
535 251494 : CYCLE
536 : END IF
537 :
538 : ! Calculate some prefactors
539 871315 : a = zeta1(i1pgf) + zeta2(i2pgf)
540 871315 : b = zetb(jpgf)
541 871315 : zet = a + b
542 871315 : xhi = a*b/zet
543 3485260 : rap = b*rab/zet
544 3485260 : rbp = -a*rab/zet
545 :
546 : ! [ss|s] integral
547 871315 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
548 :
549 : ! Calculate the recurrence relation
550 871315 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
551 871315 : IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
552 : CALL apply_cartesian_1d(rr, drr, la1_max + la2_max, lmb, .FALSE., .TRUE., a, &
553 348900 : [0.0_dp, 0.0_dp, 0.0_dp])
554 : END IF
555 :
556 3125262 : DO cob = ofb + 1, ncoset(lb_max)
557 9015788 : bvec = indco(:, cob)
558 2253947 : ib = mb + cob - ofb
559 15740535 : DO coa2 = ofa2 + 1, ncoset(la2_max)
560 50461092 : a2vec = indco(:, coa2)
561 12615273 : ia2 = ma2 + coa2 - ofa2
562 51499096 : DO coa1 = ofa1 + 1, ncoset(la1_max)
563 146519504 : avec = indco(:, coa1) + a2vec
564 36629876 : ia1 = ma1 + coa1 - ofa1
565 : ! integrals
566 36629876 : IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
567 36598584 : value = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
568 36598584 : IF (PRESENT(saab)) saab(ia1, ia2, ib) = value
569 36598584 : IF (PRESENT(saba)) saba(ia1, ib, ia2) = value
570 : END IF
571 : ! first derivatives
572 49245149 : IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
573 19863985 : dum(1) = f0*drr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
574 19863985 : dum(2) = f0*rr(avec(1), bvec(1), 1)*drr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
575 19863985 : dum(3) = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*drr(avec(3), bvec(3), 3)
576 19957861 : IF (PRESENT(daab)) daab(ia1, ia2, ib, 1:3) = dum
577 79362064 : IF (PRESENT(daba)) daba(ia1, ib, ia2, 1:3) = dum
578 : END IF
579 : END DO
580 : END DO
581 : END DO
582 :
583 1111470 : mb = mb + nb
584 : END DO
585 302786 : ma2 = ma2 + na2
586 : END DO
587 72326 : ma1 = ma1 + na1
588 : END DO
589 :
590 9695 : DEALLOCATE (rr)
591 :
592 9695 : END SUBROUTINE overlap_aab
593 :
594 : ! **************************************************************************************************
595 : !> \brief Calculation of the two-center overlap integrals [a|bb] over
596 : !> Cartesian Gaussian-type functions.
597 : !> \param la_max Max L on center A
598 : !> \param la_min Min L on center A
599 : !> \param npgfa Number of primitives on center A
600 : !> \param rpgfa Range of functions on A, used for screening
601 : !> \param zeta Exponents on center A
602 : !> \param lb1_max Max L on center B (basis 1)
603 : !> \param lb1_min Min L on center B (basis 1)
604 : !> \param npgfb1 Number of primitives on center B (basis 1)
605 : !> \param rpgfb1 Range of functions on B, used for screening (basis 1)
606 : !> \param zetb1 Exponents on center B (basis 1)
607 : !> \param lb2_max Max L on center B (basis 2)
608 : !> \param lb2_min Min L on center B (basis 2)
609 : !> \param npgfb2 Number of primitives on center B (basis 2)
610 : !> \param rpgfb2 Range of functions on B, used for screening (basis 2)
611 : !> \param zetb2 Exponents on center B (basis 2)
612 : !> \param rab Distance vector A-B
613 : !> \param sabb Final overlap integrals
614 : !> \param dabb First derivatives with respect to center A
615 : !> \date 01.07.2014
616 : !> \author JGH
617 : ! **************************************************************************************************
618 7997 : SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
619 15994 : lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
620 15994 : lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
621 7997 : rab, sabb, dabb)
622 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
623 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
624 : INTEGER, INTENT(IN) :: lb1_max, lb1_min, npgfb1
625 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb1, zetb1
626 : INTEGER, INTENT(IN) :: lb2_max, lb2_min, npgfb2
627 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb2, zetb2
628 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
629 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
630 : OPTIONAL :: sabb
631 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
632 : INTENT(INOUT), OPTIONAL :: dabb
633 :
634 : INTEGER :: coa, cob1, cob2, ia, ib1, ib2, ipgf, &
635 : j1pgf, j2pgf, ldrr, lma, lmb, ma, mb1, &
636 : mb2, na, nb1, nb2, ofa, ofb1, ofb2
637 : INTEGER, DIMENSION(3) :: avec, b2vec, bvec
638 : REAL(KIND=dp) :: a, b, f0, rab2, rpgfb, tab, xhi, zet
639 7997 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drr, rr
640 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
641 :
642 : ! Distance of the centers a and b
643 :
644 7997 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
645 7997 : tab = SQRT(rab2)
646 :
647 : ! Maximum l for auxiliary integrals
648 7997 : CPASSERT(PRESENT(sabb) .OR. PRESENT(dabb))
649 7997 : IF (PRESENT(sabb)) THEN
650 7997 : lma = la_max
651 7997 : lmb = lb1_max + lb2_max
652 : END IF
653 7997 : IF (PRESENT(dabb)) THEN
654 3397 : lma = la_max + 1
655 3397 : lmb = lb1_max + lb2_max
656 : END IF
657 7997 : ldrr = MAX(lma, lmb) + 1
658 :
659 : ! Allocate space for auxiliary integrals
660 39985 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
661 21585 : IF (PRESENT(dabb)) ALLOCATE (drr(0:lma - 1, 0:lmb, 3))
662 :
663 : ! Number of integrals, check size of arrays
664 7997 : ofa = ncoset(la_min - 1)
665 7997 : ofb1 = ncoset(lb1_min - 1)
666 7997 : ofb2 = ncoset(lb2_min - 1)
667 7997 : na = ncoset(la_max) - ofa
668 7997 : nb1 = ncoset(lb1_max) - ofb1
669 7997 : nb2 = ncoset(lb2_max) - ofb2
670 7997 : IF (PRESENT(sabb)) THEN
671 7997 : CPASSERT((SIZE(sabb, 1) >= na*npgfa))
672 7997 : CPASSERT((SIZE(sabb, 2) >= nb1*npgfb1))
673 7997 : CPASSERT((SIZE(sabb, 3) >= nb2*npgfb2))
674 : END IF
675 7997 : IF (PRESENT(dabb)) THEN
676 3397 : CPASSERT((SIZE(dabb, 1) >= na*npgfa))
677 3397 : CPASSERT((SIZE(dabb, 2) >= nb1*npgfb1))
678 3397 : CPASSERT((SIZE(dabb, 3) >= nb2*npgfb2))
679 3397 : CPASSERT((SIZE(dabb, 4) >= 3))
680 : END IF
681 :
682 : ! Loops over all pairs of primitive Gaussian-type functions
683 7997 : ma = 0
684 60026 : DO ipgf = 1, npgfa
685 52029 : mb1 = 0
686 392532 : DO j1pgf = 1, npgfb1
687 340503 : mb2 = 0
688 1393966 : DO j2pgf = 1, npgfb2
689 : ! Distance Screening
690 1053463 : rpgfb = MIN(rpgfb1(j1pgf), rpgfb2(j2pgf))
691 1053463 : IF (rpgfa(ipgf) + rpgfb < tab) THEN
692 11929000 : IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
693 15070032 : IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
694 253218 : mb2 = mb2 + nb2
695 253218 : CYCLE
696 : END IF
697 :
698 : ! Calculate some prefactors
699 800245 : a = zeta(ipgf)
700 800245 : b = zetb1(j1pgf) + zetb2(j2pgf)
701 800245 : zet = a + b
702 800245 : xhi = a*b/zet
703 3200980 : rap = b*rab/zet
704 3200980 : rbp = -a*rab/zet
705 :
706 : ! [s|s] integral
707 800245 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
708 :
709 : ! Calculate the recurrence relation
710 800245 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
711 800245 : IF (PRESENT(dabb)) THEN
712 : CALL apply_cartesian_1d(rr, drr, la_max, lmb, .FALSE., .TRUE., a, &
713 342382 : [0.0_dp, 0.0_dp, 0.0_dp])
714 : END IF
715 :
716 4903571 : DO cob2 = ofb2 + 1, ncoset(lb2_max)
717 16413304 : b2vec = indco(:, cob2)
718 4103326 : ib2 = mb2 + cob2 - ofb2
719 17227305 : DO cob1 = ofb1 + 1, ncoset(lb1_max)
720 49294936 : bvec = indco(:, cob1) + b2vec
721 12323734 : ib1 = mb1 + cob1 - ofb1
722 52406636 : DO coa = ofa + 1, ncoset(la_max)
723 143918304 : avec = indco(:, coa)
724 35979576 : ia = ma + coa - ofa
725 : ! integrals
726 35979576 : IF (PRESENT(sabb)) THEN
727 35979576 : sabb(ia, ib1, ib2) = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
728 : END IF
729 : ! first derivatives
730 48303310 : IF (PRESENT(dabb)) THEN
731 19827139 : dabb(ia, ib1, ib2, 1) = f0*drr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
732 19827139 : dabb(ia, ib1, ib2, 2) = f0*rr(avec(1), bvec(1), 1)*drr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
733 19827139 : dabb(ia, ib1, ib2, 3) = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*drr(avec(3), bvec(3), 3)
734 : END IF
735 : END DO
736 : END DO
737 : END DO
738 :
739 1140748 : mb2 = mb2 + nb2
740 : END DO
741 392532 : mb1 = mb1 + nb1
742 : END DO
743 60026 : ma = ma + na
744 : END DO
745 :
746 7997 : DEALLOCATE (rr)
747 :
748 7997 : END SUBROUTINE overlap_abb
749 :
750 : ! **************************************************************************************************
751 :
752 : ! **************************************************************************************************
753 : !> \brief Calculation of the two-center overlap integrals [a|b] over
754 : !> Spherical Gaussian-type functions.
755 : !> \param la Max L on center A
756 : !> \param zeta Exponents on center A
757 : !> \param lb Max L on center B
758 : !> \param zetb Exponents on center B
759 : !> \param rab Distance vector A-B
760 : !> \param sab Final overlap integrals
761 : !> \date 01.03.2016
762 : !> \author JGH
763 : ! **************************************************************************************************
764 275350 : SUBROUTINE overlap_ab_s(la, zeta, lb, zetb, rab, sab)
765 : INTEGER, INTENT(IN) :: la
766 : REAL(KIND=dp), INTENT(IN) :: zeta
767 : INTEGER, INTENT(IN) :: lb
768 : REAL(KIND=dp), INTENT(IN) :: zetb
769 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
770 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
771 :
772 : REAL(KIND=dp), PARAMETER :: huge4 = HUGE(1._dp)/4._dp
773 :
774 : INTEGER :: nca, ncb, nsa, nsb
775 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
776 : REAL(KIND=dp), DIMENSION(1) :: rpgf, za, zb
777 275350 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
778 :
779 275350 : rpgf(1) = huge4
780 275350 : za(1) = zeta
781 275350 : zb(1) = zetb
782 :
783 275350 : nca = nco(la)
784 275350 : ncb = nco(lb)
785 1101400 : ALLOCATE (cab(nca, ncb))
786 275350 : nsa = nso(la)
787 275350 : nsb = nso(lb)
788 :
789 275350 : CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
790 :
791 275350 : c2sa => orbtramat(la)%c2s
792 275350 : c2sb => orbtramat(lb)%c2s
793 275350 : sab(1:nsa, 1:nsb) = MATMUL(c2sa(1:nsa, 1:nca), &
794 11068283 : MATMUL(cab(1:nca, 1:ncb), TRANSPOSE(c2sb(1:nsb, 1:ncb))))
795 :
796 275350 : DEALLOCATE (cab)
797 :
798 275350 : END SUBROUTINE overlap_ab_s
799 :
800 : ! **************************************************************************************************
801 : !> \brief Calculation of the overlap integrals [a|b] over
802 : !> cubic periodic Spherical Gaussian-type functions.
803 : !> \param la Max L on center A
804 : !> \param zeta Exponents on center A
805 : !> \param lb Max L on center B
806 : !> \param zetb Exponents on center B
807 : !> \param alat Lattice constant
808 : !> \param sab Final overlap integrals
809 : !> \date 01.03.2016
810 : !> \author JGH
811 : ! **************************************************************************************************
812 36066 : SUBROUTINE overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
813 : INTEGER, INTENT(IN) :: la
814 : REAL(KIND=dp), INTENT(IN) :: zeta
815 : INTEGER, INTENT(IN) :: lb
816 : REAL(KIND=dp), INTENT(IN) :: zetb, alat
817 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
818 :
819 : COMPLEX(KIND=dp) :: zfg
820 36066 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fun, gun
821 : INTEGER :: ax, ay, az, bx, by, bz, i, ia, ib, l, &
822 : l1, l2, na, nb, nca, ncb, nmax, nsa, &
823 : nsb
824 : REAL(KIND=dp) :: oa, ob, ovol, zm
825 36066 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
826 36066 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
827 : REAL(KIND=dp), DIMENSION(0:3, 0:3) :: fgsum
828 36066 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
829 :
830 36066 : nca = nco(la)
831 36066 : ncb = nco(lb)
832 144264 : ALLOCATE (cab(nca, ncb))
833 36066 : cab = 0.0_dp
834 36066 : nsa = nso(la)
835 36066 : nsb = nso(lb)
836 :
837 36066 : zm = MIN(zeta, zetb)
838 36066 : nmax = NINT(1.81_dp*alat*SQRT(zm) + 1.0_dp)
839 : ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
840 396726 : fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
841 :
842 36066 : oa = 1._dp/zeta
843 36066 : ob = 1._dp/zetb
844 522120 : DO i = -nmax, nmax
845 486054 : gval(i) = twopi/alat*REAL(i, KIND=dp)
846 486054 : fexp(i) = SQRT(oa*pi)*EXP(-0.25_dp*oa*gval(i)**2)
847 522120 : gexp(i) = SQRT(ob*pi)*EXP(-0.25_dp*ob*gval(i)**2)
848 : END DO
849 91966 : DO l = 0, la
850 91966 : IF (l == 0) THEN
851 522120 : fun(:, l) = z_one
852 19834 : ELSE IF (l == 1) THEN
853 252786 : fun(:, l) = CMPLX(0.0_dp, 0.5_dp*oa*gval(:), KIND=dp)
854 2073 : ELSE IF (l == 2) THEN
855 29772 : fun(:, l) = CMPLX(-(0.5_dp*oa*gval(:))**2, 0.0_dp, KIND=dp)
856 29772 : fun(:, l) = fun(:, l) + CMPLX(0.5_dp*oa, 0.0_dp, KIND=dp)
857 0 : ELSE IF (l == 3) THEN
858 0 : fun(:, l) = CMPLX(0.0_dp, -(0.5_dp*oa*gval(:))**3, KIND=dp)
859 0 : fun(:, l) = fun(:, l) + CMPLX(0.0_dp, 0.75_dp*oa*oa*gval(:), KIND=dp)
860 : ELSE
861 0 : CPABORT("l value too high")
862 : END IF
863 : END DO
864 91966 : DO l = 0, lb
865 91966 : IF (l == 0) THEN
866 522120 : gun(:, l) = z_one
867 19834 : ELSE IF (l == 1) THEN
868 252786 : gun(:, l) = CMPLX(0.0_dp, 0.5_dp*ob*gval(:), KIND=dp)
869 2073 : ELSE IF (l == 2) THEN
870 29772 : gun(:, l) = CMPLX(-(0.5_dp*ob*gval(:))**2, 0.0_dp, KIND=dp)
871 29772 : gun(:, l) = gun(:, l) + CMPLX(0.5_dp*ob, 0.0_dp, KIND=dp)
872 0 : ELSE IF (l == 3) THEN
873 0 : gun(:, l) = CMPLX(0.0_dp, -(0.5_dp*ob*gval(:))**3, KIND=dp)
874 0 : gun(:, l) = gun(:, l) + CMPLX(0.0_dp, 0.75_dp*ob*ob*gval(:), KIND=dp)
875 : ELSE
876 0 : CPABORT("l value too high")
877 : END IF
878 : END DO
879 :
880 36066 : fgsum = 0.0_dp
881 91966 : DO l1 = 0, la
882 179954 : DO l2 = 0, lb
883 1260466 : zfg = SUM(CONJG(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
884 143888 : fgsum(l1, l2) = REAL(zfg, KIND=dp)
885 : END DO
886 : END DO
887 :
888 36066 : na = ncoset(la - 1)
889 36066 : nb = ncoset(lb - 1)
890 91966 : DO ax = 0, la
891 169773 : DO ay = 0, la - ax
892 77807 : az = la - ax - ay
893 77807 : ia = coset(ax, ay, az) - na
894 257848 : DO bx = 0, lb
895 379135 : DO by = 0, lb - bx
896 177187 : bz = lb - bx - by
897 177187 : ib = coset(bx, by, bz) - nb
898 301328 : cab(ia, ib) = fgsum(ax, bx)*fgsum(ay, by)*fgsum(az, bz)
899 : END DO
900 : END DO
901 : END DO
902 : END DO
903 :
904 36066 : c2sa => orbtramat(la)%c2s
905 36066 : c2sb => orbtramat(lb)%c2s
906 36066 : sab(1:nsa, 1:nsb) = MATMUL(c2sa(1:nsa, 1:nca), &
907 2482203 : MATMUL(cab(1:nca, 1:ncb), TRANSPOSE(c2sb(1:nsb, 1:ncb))))
908 36066 : ovol = 1._dp/(alat**3)
909 276218 : sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
910 :
911 36066 : DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
912 :
913 36066 : END SUBROUTINE overlap_ab_sp
914 :
915 311416 : END MODULE ai_overlap
|
