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