(git:d2a9ebd)
Loading...
Searching...
No Matches
ai_overlap.F
Go to the documentation of this file.
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! **************************************************************************************************
18MODULE 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
41CONTAINS
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 SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
68 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
69 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
88 REAL(kind=dp), DIMENSION(3) :: rap, rbp
89
90 da_max = da_max_set
91 la_max = la_max_set + da_max_set
92
93 lb_max = lb_max_set
94 ldrr = max(la_max, lb_max) + 1
95 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
96
97! *** Loop over all pairs of primitive Gaussian-type functions ***
98
99 na = 0
100 DO ipgf = 1, npgfa
101
102 nb = 0
103
104 DO jpgf = 1, npgfb
105
106! *** Screening ***
107
108 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 DO j = nb + 1, nb + ncoset(lb_max_set)
110 DO i = na + 1, na + ncoset(la_max_set)
111 sab(i, j) = 0.0_dp
112 END DO
113 END DO
114 IF (return_derivatives) THEN
115 DO k = 2, ncoset(da_max_set)
116 jstart = (k - 1)*SIZE(sab, 1)
117 DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
118 DO i = na + 1, na + ncoset(la_max_set)
119 sab(i, j) = 0.0_dp
120 END DO
121 END DO
122 END DO
123 END IF
124 nb = nb + ncoset(lb_max_set)
125 cycle
126 END IF
127
128! *** Calculate some prefactors ***
129
130 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
131
132 f0 = sqrt((pi*zetp)**3)*exp(-zeta(ipgf)*zetb(jpgf)*zetp*dab*dab)
133 rap(:) = zetb(jpgf)*zetp*rab(:)
134 rbp(:) = -zeta(ipgf)*zetp*rab(:)
135
136 CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
137
138 DO lb = 0, lb_max
139 DO bx = 0, lb
140 DO by = 0, lb - bx
141 bz = lb - bx - by
142 cob = coset(bx, by, bz)
143 DO la = 0, la_max
144 DO ax = 0, la
145 DO ay = 0, la - ax
146 az = la - ax - ay
147 coa = coset(ax, ay, az)
148 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 DO j = 1, ncoset(lb_max_set)
159 DO i = 1, ncoset(la_max_set)
160 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 IF (return_derivatives) THEN
168 la_start = 0
169 lb_start = 0
170 ELSE
171 la_start = la_min_set
172 lb_start = lb_min_set
173 END IF
174
175 DO da = 0, da_max - 1
176 ftz = 2.0_dp*zeta(ipgf)
177 DO dax = 0, da
178 DO day = 0, da - dax
179 daz = da - dax - day
180 cda = coset(dax, day, daz)
181 cdax = coset(dax + 1, day, daz)
182 cday = coset(dax, day + 1, daz)
183 cdaz = coset(dax, day, daz + 1)
184
185! *** [da/dAi|b] = 2*zeta*[a+1i|b] - Ni(a)[a-1i|b] ***
186
187 DO la = la_start, la_max - da - 1
188 DO ax = 0, la
189 fax = real(ax, dp)
190 DO ay = 0, la - ax
191 fay = real(ay, dp)
192 az = la - ax - ay
193 faz = real(az, dp)
194 coa = coset(ax, ay, az)
195 coamx = coset(ax - 1, ay, az)
196 coamy = coset(ax, ay - 1, az)
197 coamz = coset(ax, ay, az - 1)
198 coapx = coset(ax + 1, ay, az)
199 coapy = coset(ax, ay + 1, az)
200 coapz = coset(ax, ay, az + 1)
201 DO lb = lb_start, lb_max_set
202 DO bx = 0, lb
203 DO by = 0, lb - bx
204 bz = lb - bx - by
205 cob = coset(bx, by, bz)
206 s(coa, cob, cdax) = ftz*s(coapx, cob, cda) - &
207 fax*s(coamx, cob, cda)
208 s(coa, cob, cday) = ftz*s(coapy, cob, cda) - &
209 fay*s(coamy, cob, cda)
210 s(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - &
211 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 IF (return_derivatives) THEN
227 DO k = 2, ncoset(da_max_set)
228 jstart = (k - 1)*SIZE(sab, 1)
229 DO j = 1, ncoset(lb_max_set)
230 jk = jstart + j
231 DO i = 1, ncoset(la_max_set)
232 sab(na + i, nb + jk) = s(i, j, k)
233 END DO
234 END DO
235 END DO
236 END IF
237
238 nb = nb + ncoset(lb_max_set)
239
240 END DO
241
242 na = na + ncoset(la_max_set)
243 END DO
244
245 DEALLOCATE (rr)
246
247 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 SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
271 lb_max, lb_min, npgfb, rpgfb, zetb, &
272 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 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rr
292
293 ! Distance of the centers a and b
294
295 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
296 tab = sqrt(rab2)
297
298 ! Maximum l for auxiliary integrals
299 cpassert(PRESENT(sab) .OR. PRESENT(dab) .OR. PRESENT(ddab))
300 IF (PRESENT(sab)) THEN
301 lma = la_max
302 lmb = lb_max
303 END IF
304 IF (PRESENT(dab)) THEN
305 lma = la_max + 1
306 lmb = lb_max
307 END IF
308 IF (PRESENT(ddab)) THEN
309 lma = la_max + 1
310 lmb = lb_max + 1
311 END IF
312 ldrr = max(lma, lmb) + 1
313
314 ! Allocate or attach the workspace for auxiliary integrals
315 NULLIFY (rr)
316 IF (PRESENT(rr_work)) THEN
317 cpassert(SIZE(rr_work) >= ldrr*ldrr*3)
318 rr(0:ldrr - 1, 0:ldrr - 1, 1:3) => rr_work(1:ldrr*ldrr*3)
319 ELSE
320 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
321 END IF
322
323 ! Number of integrals, check size of arrays
324 ofa = ncoset(la_min - 1)
325 ofb = ncoset(lb_min - 1)
326 na = ncoset(la_max) - ofa
327 nb = ncoset(lb_max) - ofb
328 IF (PRESENT(sab)) THEN
329 cpassert((SIZE(sab, 1) >= na*npgfa))
330 cpassert((SIZE(sab, 2) >= nb*npgfb))
331 END IF
332 IF (PRESENT(dab)) THEN
333 cpassert((SIZE(dab, 1) >= na*npgfa))
334 cpassert((SIZE(dab, 2) >= nb*npgfb))
335 cpassert((SIZE(dab, 3) >= 3))
336 END IF
337 IF (PRESENT(ddab)) THEN
338 cpassert((SIZE(ddab, 1) >= na*npgfa))
339 cpassert((SIZE(ddab, 2) >= nb*npgfb))
340 cpassert((SIZE(ddab, 3) >= 6))
341 END IF
342
343 ! Loops over all pairs of primitive Gaussian-type functions
344 ma = 0
345 DO ipgf = 1, npgfa
346 mb = 0
347 DO jpgf = 1, npgfb
348 ! Distance Screening
349 IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
350 IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
351 IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
352 IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
353 mb = mb + nb
354 cycle
355 END IF
356
357 ! Calculate some prefactors
358 a = zeta(ipgf)
359 b = zetb(jpgf)
360 zet = a + b
361 xhi = a*b/zet
362 rap = b*rab/zet
363 rbp = -a*rab/zet
364
365 ! [s|s] integral
366 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
367
368 ! Calculate the recurrence relation
369 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
370
371 DO lb = lb_min, lb_max
372 DO bx = 0, lb
373 DO by = 0, lb - bx
374 bz = lb - bx - by
375 cob = coset(bx, by, bz) - ofb
376 ib = mb + cob
377 DO la = la_min, la_max
378 DO ax = 0, la
379 DO ay = 0, la - ax
380 az = la - ax - ay
381 coa = coset(ax, ay, az) - ofa
382 ia = ma + coa
383 ! integrals
384 IF (PRESENT(sab)) THEN
385 sab(ia, ib) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3)
386 END IF
387 ! first derivatives
388 IF (PRESENT(dab)) THEN
389 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
390 ! dx
391 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
392 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
393 dab(ia, ib, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
394 ! dy
395 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
396 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
397 dab(ia, ib, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
398 ! dz
399 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
400 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
401 dab(ia, ib, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
402 END IF
403 ! 2nd derivatives
404 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 apbp = f0*rr(ax + 1, bx + 1, 1)*rr(ay, by, 2)*rr(az, bz, 3)
409 IF (bx > 0) THEN
410 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 IF (ax > 0) THEN
415 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 IF (ax > 0 .AND. bx > 0) THEN
420 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 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bx, dp)*ambm
426 ! dx dy
427 apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by + 1, 2)*rr(az, bz, 3)
428 IF (by > 0) THEN
429 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 IF (ax > 0) THEN
434 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 IF (ax > 0 .AND. by > 0) THEN
439 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 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(by, dp)*ambm
445 ! dx dz
446 apbp = f0*rr(ax + 1, bx, 1)*rr(ay, by, 2)*rr(az, bz + 1, 3)
447 IF (bz > 0) THEN
448 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 IF (ax > 0) THEN
453 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 IF (ax > 0 .AND. bz > 0) THEN
458 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 + 2.0_dp*b*real(ax, dp)*ambp - real(ax, dp)*real(bz, dp)*ambm
464 ! dy dy
465 apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by + 1, 2)*rr(az, bz, 3)
466 IF (by > 0) THEN
467 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 IF (ay > 0) THEN
472 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 IF (ay > 0 .AND. by > 0) THEN
477 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 + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(by, dp)*ambm
483 ! dy dz
484 apbp = f0*rr(ax, bx, 1)*rr(ay + 1, by, 2)*rr(az, bz + 1, 3)
485 IF (bz > 0) THEN
486 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 IF (ay > 0) THEN
491 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 IF (ay > 0 .AND. bz > 0) THEN
496 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 + 2.0_dp*b*real(ay, dp)*ambp - real(ay, dp)*real(bz, dp)*ambm
502 ! dz dz
503 apbp = f0*rr(ax, bx, 1)*rr(ay, by, 2)*rr(az + 1, bz + 1, 3)
504 IF (bz > 0) THEN
505 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 IF (az > 0) THEN
510 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 IF (az > 0 .AND. bz > 0) THEN
515 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 + 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 mb = mb + nb
531 END DO
532 ma = ma + na
533 END DO
534
535 IF (.NOT. PRESENT(rr_work)) DEALLOCATE (rr)
536 NULLIFY (rr)
537
538 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 SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
567 la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
568 lb_max, lb_min, npgfb, rpgfb, zetb, &
569 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
597 tab = sqrt(rab2)
598
599 ! Maximum l for auxiliary integrals
600 cpassert(PRESENT(saab) .OR. PRESENT(daab) .OR. PRESENT(saba) .OR. PRESENT(daba))
601 IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
602 lma = la1_max + la2_max
603 lmb = lb_max
604 END IF
605 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
606 lma = la1_max + la2_max + 1
607 lmb = lb_max
608 END IF
609 ldrr = max(lma, lmb) + 1
610
611 ! Allocate space for auxiliary integrals
612 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
613
614 ! Number of integrals, check size of arrays
615 ofa1 = ncoset(la1_min - 1)
616 ofa2 = ncoset(la2_min - 1)
617 ofb = ncoset(lb_min - 1)
618 na1 = ncoset(la1_max) - ofa1
619 na2 = ncoset(la2_max) - ofa2
620 nb = ncoset(lb_max) - ofb
621 IF (PRESENT(saab)) THEN
622 cpassert((SIZE(saab, 1) >= na1*npgfa1))
623 cpassert((SIZE(saab, 2) >= na2*npgfa2))
624 cpassert((SIZE(saab, 3) >= nb*npgfb))
625 END IF
626 IF (PRESENT(daab)) THEN
627 cpassert((SIZE(daab, 1) >= na1*npgfa1))
628 cpassert((SIZE(daab, 2) >= na2*npgfa2))
629 cpassert((SIZE(daab, 3) >= nb*npgfb))
630 cpassert((SIZE(daab, 4) >= 3))
631 END IF
632 IF (PRESENT(saba)) THEN
633 cpassert((SIZE(saba, 1) >= na1*npgfa1))
634 cpassert((SIZE(saba, 2) >= nb*npgfb))
635 cpassert((SIZE(saba, 3) >= na2*npgfa2))
636 END IF
637 IF (PRESENT(daba)) THEN
638 cpassert((SIZE(daba, 1) >= na1*npgfa1))
639 cpassert((SIZE(daba, 2) >= nb*npgfb))
640 cpassert((SIZE(daba, 3) >= na2*npgfa2))
641 cpassert((SIZE(daba, 4) >= 3))
642 END IF
643
644 ! Loops over all primitive Gaussian-type functions
645 ma1 = 0
646 DO i1pgf = 1, npgfa1
647 ma2 = 0
648 DO i2pgf = 1, npgfa2
649 rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
650 mb = 0
651 DO jpgf = 1, npgfb
652 ! Distance Screening
653 IF (rpgfa + rpgfb(jpgf) < tab) THEN
654 IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
655 IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
656 IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
657 IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
658 mb = mb + nb
659 cycle
660 END IF
661
662 ! Calculate some prefactors
663 a = zeta1(i1pgf) + zeta2(i2pgf)
664 b = zetb(jpgf)
665 zet = a + b
666 xhi = a*b/zet
667 rap = b*rab/zet
668 rbp = -a*rab/zet
669
670 ! [ss|s] integral
671 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
672
673 ! Calculate the recurrence relation
674 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
675
676 DO lb = lb_min, lb_max
677 DO bx = 0, lb
678 DO by = 0, lb - bx
679 bz = lb - bx - by
680 cob = coset(bx, by, bz) - ofb
681 ib = mb + cob
682 DO la2 = la2_min, la2_max
683 DO ax2 = 0, la2
684 DO ay2 = 0, la2 - ax2
685 az2 = la2 - ax2 - ay2
686 coa2 = coset(ax2, ay2, az2) - ofa2
687 ia2 = ma2 + coa2
688 DO la1 = la1_min, la1_max
689 DO ax1 = 0, la1
690 DO ay1 = 0, la1 - ax1
691 az1 = la1 - ax1 - ay1
692 coa1 = coset(ax1, ay1, az1) - ofa1
693 ia1 = ma1 + coa1
694 ! integrals
695 IF (PRESENT(saab)) THEN
696 saab(ia1, ia2, ib) = f0*rr(ax1 + ax2, bx, 1)*rr(ay1 + ay2, by, 2)*rr(az1 + az2, bz, 3)
697 END IF
698 IF (PRESENT(saba)) THEN
699 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 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
703 ax = ax1 + ax2
704 ay = ay1 + ay2
705 az = az1 + az2
706 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
707 ! dx
708 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
709 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
710 dumx = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
711 ! dy
712 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
713 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
714 dumy = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
715 ! dz
716 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
717 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
718 dumz = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
719 IF (PRESENT(daab)) THEN
720 daab(ia1, ia2, ib, 1) = dumx
721 daab(ia1, ia2, ib, 2) = dumy
722 daab(ia1, ia2, ib, 3) = dumz
723 END IF
724 IF (PRESENT(daba)) THEN
725 daba(ia1, ib, ia2, 1) = dumx
726 daba(ia1, ib, ia2, 2) = dumy
727 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 mb = mb + nb
742 END DO
743 ma2 = ma2 + na2
744 END DO
745 ma1 = ma1 + na1
746 END DO
747
748 DEALLOCATE (rr)
749
750 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 SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
777 lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
778 lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
779 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
803 tab = sqrt(rab2)
804
805 ! Maximum l for auxiliary integrals
806 cpassert(PRESENT(sabb) .OR. PRESENT(dabb))
807 IF (PRESENT(sabb)) THEN
808 lma = la_max
809 lmb = lb1_max + lb2_max
810 END IF
811 IF (PRESENT(dabb)) THEN
812 lma = la_max + 1
813 lmb = lb1_max + lb2_max
814 END IF
815 ldrr = max(lma, lmb) + 1
816
817 ! Allocate space for auxiliary integrals
818 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
819
820 ! Number of integrals, check size of arrays
821 ofa = ncoset(la_min - 1)
822 ofb1 = ncoset(lb1_min - 1)
823 ofb2 = ncoset(lb2_min - 1)
824 na = ncoset(la_max) - ofa
825 nb1 = ncoset(lb1_max) - ofb1
826 nb2 = ncoset(lb2_max) - ofb2
827 IF (PRESENT(sabb)) THEN
828 cpassert((SIZE(sabb, 1) >= na*npgfa))
829 cpassert((SIZE(sabb, 2) >= nb1*npgfb1))
830 cpassert((SIZE(sabb, 3) >= nb2*npgfb2))
831 END IF
832 IF (PRESENT(dabb)) THEN
833 cpassert((SIZE(dabb, 1) >= na*npgfa))
834 cpassert((SIZE(dabb, 2) >= nb1*npgfb1))
835 cpassert((SIZE(dabb, 3) >= nb2*npgfb2))
836 cpassert((SIZE(dabb, 4) >= 3))
837 END IF
838
839 ! Loops over all pairs of primitive Gaussian-type functions
840 ma = 0
841 DO ipgf = 1, npgfa
842 mb1 = 0
843 DO j1pgf = 1, npgfb1
844 mb2 = 0
845 DO j2pgf = 1, npgfb2
846 ! Distance Screening
847 rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
848 IF (rpgfa(ipgf) + rpgfb < tab) THEN
849 IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
850 IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
851 mb2 = mb2 + nb2
852 cycle
853 END IF
854
855 ! Calculate some prefactors
856 a = zeta(ipgf)
857 b = zetb1(j1pgf) + zetb2(j2pgf)
858 zet = a + b
859 xhi = a*b/zet
860 rap = b*rab/zet
861 rbp = -a*rab/zet
862
863 ! [s|s] integral
864 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
865
866 ! Calculate the recurrence relation
867 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
868
869 DO lb2 = lb2_min, lb2_max
870 DO bx2 = 0, lb2
871 DO by2 = 0, lb2 - bx2
872 bz2 = lb2 - bx2 - by2
873 cob2 = coset(bx2, by2, bz2) - ofb2
874 ib2 = mb2 + cob2
875 DO lb1 = lb1_min, lb1_max
876 DO bx1 = 0, lb1
877 DO by1 = 0, lb1 - bx1
878 bz1 = lb1 - bx1 - by1
879 cob1 = coset(bx1, by1, bz1) - ofb1
880 ib1 = mb1 + cob1
881 DO la = la_min, la_max
882 DO ax = 0, la
883 DO ay = 0, la - ax
884 az = la - ax - ay
885 coa = coset(ax, ay, az) - ofa
886 ia = ma + coa
887 ! integrals
888 IF (PRESENT(sabb)) THEN
889 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 IF (PRESENT(dabb)) THEN
893 bx = bx1 + bx2
894 by = by1 + by2
895 bz = bz1 + bz2
896 ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
897 ! dx
898 dumx = 2.0_dp*a*rr(ax + 1, bx, 1)
899 IF (ax > 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
900 dabb(ia, ib1, ib2, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
901 ! dy
902 dumy = 2.0_dp*a*rr(ay + 1, by, 2)
903 IF (ay > 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
904 dabb(ia, ib1, ib2, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
905 ! dz
906 dumz = 2.0_dp*a*rr(az + 1, bz, 3)
907 IF (az > 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
908 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 mb2 = mb2 + nb2
922 END DO
923 mb1 = mb1 + nb1
924 END DO
925 ma = ma + na
926 END DO
927
928 DEALLOCATE (rr)
929
930 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
960
961 rpgf(1) = huge4
962 za(1) = zeta
963 zb(1) = zetb
964
965 nca = nco(la)
966 ncb = nco(lb)
967 ALLOCATE (cab(nca, ncb))
968 nsa = nso(la)
969 nsb = nso(lb)
970
971 CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
972
973 c2sa => orbtramat(la)%c2s
974 c2sb => orbtramat(lb)%c2s
975 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
976 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
977
978 DEALLOCATE (cab)
979
980 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
1008 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
1009 REAL(kind=dp), DIMENSION(0:3, 0:3) :: fgsum
1010 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
1011
1012 nca = nco(la)
1013 ncb = nco(lb)
1014 ALLOCATE (cab(nca, ncb))
1015 cab = 0.0_dp
1016 nsa = nso(la)
1017 nsb = nso(lb)
1018
1019 zm = min(zeta, zetb)
1020 nmax = nint(1.81_dp*alat*sqrt(zm) + 1.0_dp)
1021 ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
1022 fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
1023
1024 oa = 1._dp/zeta
1025 ob = 1._dp/zetb
1026 DO i = -nmax, nmax
1027 gval(i) = twopi/alat*real(i, kind=dp)
1028 fexp(i) = sqrt(oa*pi)*exp(-0.25_dp*oa*gval(i)**2)
1029 gexp(i) = sqrt(ob*pi)*exp(-0.25_dp*ob*gval(i)**2)
1030 END DO
1031 DO l = 0, la
1032 IF (l == 0) THEN
1033 fun(:, l) = z_one
1034 ELSE IF (l == 1) THEN
1035 fun(:, l) = cmplx(0.0_dp, 0.5_dp*oa*gval(:), kind=dp)
1036 ELSE IF (l == 2) THEN
1037 fun(:, l) = cmplx(-(0.5_dp*oa*gval(:))**2, 0.0_dp, kind=dp)
1038 fun(:, l) = fun(:, l) + cmplx(0.5_dp*oa, 0.0_dp, kind=dp)
1039 ELSE IF (l == 3) THEN
1040 fun(:, l) = cmplx(0.0_dp, -(0.5_dp*oa*gval(:))**3, kind=dp)
1041 fun(:, l) = fun(:, l) + cmplx(0.0_dp, 0.75_dp*oa*oa*gval(:), kind=dp)
1042 ELSE
1043 cpabort("l value too high")
1044 END IF
1045 END DO
1046 DO l = 0, lb
1047 IF (l == 0) THEN
1048 gun(:, l) = z_one
1049 ELSE IF (l == 1) THEN
1050 gun(:, l) = cmplx(0.0_dp, 0.5_dp*ob*gval(:), kind=dp)
1051 ELSE IF (l == 2) THEN
1052 gun(:, l) = cmplx(-(0.5_dp*ob*gval(:))**2, 0.0_dp, kind=dp)
1053 gun(:, l) = gun(:, l) + cmplx(0.5_dp*ob, 0.0_dp, kind=dp)
1054 ELSE IF (l == 3) THEN
1055 gun(:, l) = cmplx(0.0_dp, -(0.5_dp*ob*gval(:))**3, kind=dp)
1056 gun(:, l) = gun(:, l) + cmplx(0.0_dp, 0.75_dp*ob*ob*gval(:), kind=dp)
1057 ELSE
1058 cpabort("l value too high")
1059 END IF
1060 END DO
1061
1062 fgsum = 0.0_dp
1063 DO l1 = 0, la
1064 DO l2 = 0, lb
1065 zfg = sum(conjg(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
1066 fgsum(l1, l2) = real(zfg, kind=dp)
1067 END DO
1068 END DO
1069
1070 na = ncoset(la - 1)
1071 nb = ncoset(lb - 1)
1072 DO ax = 0, la
1073 DO ay = 0, la - ax
1074 az = la - ax - ay
1075 ia = coset(ax, ay, az) - na
1076 DO bx = 0, lb
1077 DO by = 0, lb - bx
1078 bz = lb - bx - by
1079 ib = coset(bx, by, bz) - nb
1080 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 c2sa => orbtramat(la)%c2s
1087 c2sb => orbtramat(lb)%c2s
1088 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
1089 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
1090 ovol = 1._dp/(alat**3)
1091 sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
1092
1093 DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
1094
1095 END SUBROUTINE overlap_ab_sp
1096
1097END MODULE ai_overlap
ai_os_rr
Definition ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition ai_os_rr.F:39
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition ai_overlap.F:70
ai_overlap::overlap_ab_s
subroutine, public overlap_ab_s(la, zeta, lb, zetb, rab, sab)
Calculation of the two-center overlap integrals [a|b] over Spherical Gaussian-type functions.
Definition ai_overlap.F:947
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab, rr_work)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:273
ai_overlap::overlap_ab_sp
subroutine, public overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
Calculation of the overlap integrals [a|b] over cubic periodic Spherical Gaussian-type functions.
Definition ai_overlap.F:995
ai_overlap::overlap_aab
subroutine, public overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, la2_max, la2_min, npgfa2, rpgfa2, zeta2, lb_max, lb_min, npgfb, rpgfb, zetb, rab, saab, daab, saba, daba)
Calculation of the two-center overlap integrals [aa|b] over Cartesian Gaussian-type functions.
Definition ai_overlap.F:570
ai_overlap::overlap_abb
subroutine, public overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, rab, sabb, dabb)
Calculation of the two-center overlap integrals [a|bb] over Cartesian Gaussian-type functions.
Definition ai_overlap.F:780
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47