(git:9ccafda)
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_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
45CONTAINS
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 SUBROUTINE overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
72 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
73 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 INTEGER, ALLOCATABLE :: derivative_map(:, :, :)
91 REAL(kind=dp) :: f0, zetp
92 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
93 REAL(kind=dp), DIMENSION(3) :: rap, rbp
94
95 da_max = da_max_set
96 la_max = la_max_set + da_max_set
97
98 lb_max = lb_max_set
99 ldrr = max(la_max, lb_max) + 1
100 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
101 IF (da_max > 0) THEN
102 ALLOCATE (derivative_map(4, ncoset(la_max - 1), 3))
103 DO axis = 1, 3
104 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 na = 0
111 DO ipgf = 1, npgfa
112
113 nb = 0
114
115 DO jpgf = 1, npgfb
116
117! *** Screening ***
118
119 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
120 DO j = nb + 1, nb + ncoset(lb_max_set)
121 DO i = na + 1, na + ncoset(la_max_set)
122 sab(i, j) = 0.0_dp
123 END DO
124 END DO
125 IF (return_derivatives) THEN
126 DO k = 2, ncoset(da_max_set)
127 jstart = (k - 1)*SIZE(sab, 1)
128 DO j = jstart + nb + 1, jstart + nb + ncoset(lb_max_set)
129 DO i = na + 1, na + ncoset(la_max_set)
130 sab(i, j) = 0.0_dp
131 END DO
132 END DO
133 END DO
134 END IF
135 nb = nb + ncoset(lb_max_set)
136 cycle
137 END IF
138
139! *** Calculate some prefactors ***
140
141 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
142
143 f0 = sqrt((pi*zetp)**3)*exp(-zeta(ipgf)*zetb(jpgf)*zetp*dab*dab)
144 rap(:) = zetb(jpgf)*zetp*rab(:)
145 rbp(:) = -zeta(ipgf)*zetp*rab(:)
146
147 CALL os_rr_ovlp(rap, la_max, rbp, lb_max, 1.0_dp/zetp, ldrr, rr)
148
149 DO cob = 1, ncoset(lb_max)
150 bx = indco(1, cob)
151 by = indco(2, cob)
152 bz = indco(3, cob)
153 DO coa = 1, ncoset(la_max)
154 ax = indco(1, coa)
155 ay = indco(2, coa)
156 az = indco(3, coa)
157 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 DO j = 1, ncoset(lb_max_set)
164 DO i = 1, ncoset(la_max_set)
165 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 IF (return_derivatives) THEN
173 la_start = 0
174 lb_start = 0
175 ELSE
176 la_start = la_min_set
177 lb_start = lb_min_set
178 END IF
179
180 first = ncoset(la_start - 1) + 1
181 first_b = ncoset(lb_start - 1) + 1
182 last_b = ncoset(lb_max_set)
183 DO da = 0, da_max - 1
184 last = ncoset(la_max - da - 1)
185 DO dax = 0, da
186 DO day = 0, da - dax
187 daz = da - dax - day
188 cda = coset(dax, day, daz)
189 cdax = coset(dax + 1, day, daz)
190 cday = coset(dax, day + 1, daz)
191 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 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 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 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 IF (return_derivatives) THEN
206 DO k = 2, ncoset(da_max_set)
207 jstart = (k - 1)*SIZE(sab, 1)
208 DO j = 1, ncoset(lb_max_set)
209 jk = jstart + j
210 DO i = 1, ncoset(la_max_set)
211 sab(na + i, nb + jk) = s(i, j, k)
212 END DO
213 END DO
214 END DO
215 END IF
216
217 nb = nb + ncoset(lb_max_set)
218
219 END DO
220
221 na = na + ncoset(la_max_set)
222 END DO
223
224 DEALLOCATE (rr)
225
226 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 SUBROUTINE overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, &
250 lb_max, lb_min, npgfb, rpgfb, zetb, &
251 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 REAL(kind=dp), ALLOCATABLE :: drr_a(:, :, :), drr_ab(:, :, :)
273 REAL(kind=dp), DIMENSION(3) :: rap, rbp
274 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rr
275
276 ! Distance of the centers a and b
277
278 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
279 tab = sqrt(rab2)
280
281 ! Maximum l for auxiliary integrals
282 cpassert(PRESENT(sab) .OR. PRESENT(dab) .OR. PRESENT(ddab))
283 IF (PRESENT(sab)) THEN
284 lma = la_max
285 lmb = lb_max
286 END IF
287 IF (PRESENT(dab)) THEN
288 lma = la_max + 1
289 lmb = lb_max
290 END IF
291 IF (PRESENT(ddab)) THEN
292 lma = la_max + 1
293 lmb = lb_max + 1
294 END IF
295 ldrr = max(lma, lmb) + 1
296
297 ! Allocate or attach the workspace for auxiliary integrals
298 NULLIFY (rr)
299 IF (PRESENT(rr_work)) THEN
300 cpassert(SIZE(rr_work) >= ldrr*ldrr*3)
301 rr(0:ldrr - 1, 0:ldrr - 1, 1:3) => rr_work(1:ldrr*ldrr*3)
302 ELSE
303 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
304 END IF
305
306 IF (PRESENT(dab) .OR. PRESENT(ddab)) ALLOCATE (drr_a(0:ldrr - 1, 0:ldrr - 1, 3))
307 IF (PRESENT(ddab)) ALLOCATE (drr_ab(0:ldrr - 1, 0:ldrr - 1, 3))
308
309 ! Number of integrals, check size of arrays
310 ofa = ncoset(la_min - 1)
311 ofb = ncoset(lb_min - 1)
312 na = ncoset(la_max) - ofa
313 nb = ncoset(lb_max) - ofb
314 IF (PRESENT(sab)) THEN
315 cpassert((SIZE(sab, 1) >= na*npgfa))
316 cpassert((SIZE(sab, 2) >= nb*npgfb))
317 END IF
318 IF (PRESENT(dab)) THEN
319 cpassert((SIZE(dab, 1) >= na*npgfa))
320 cpassert((SIZE(dab, 2) >= nb*npgfb))
321 cpassert((SIZE(dab, 3) >= 3))
322 END IF
323 IF (PRESENT(ddab)) THEN
324 cpassert((SIZE(ddab, 1) >= na*npgfa))
325 cpassert((SIZE(ddab, 2) >= nb*npgfb))
326 cpassert((SIZE(ddab, 3) >= 6))
327 END IF
328
329 ! Loops over all pairs of primitive Gaussian-type functions
330 ma = 0
331 DO ipgf = 1, npgfa
332 mb = 0
333 DO jpgf = 1, npgfb
334 ! Distance Screening
335 IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
336 IF (PRESENT(sab)) sab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
337 IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
338 IF (PRESENT(ddab)) ddab(ma + 1:ma + na, mb + 1:mb + nb, 1:6) = 0.0_dp
339 mb = mb + nb
340 cycle
341 END IF
342
343 ! Calculate some prefactors
344 a = zeta(ipgf)
345 b = zetb(jpgf)
346 zet = a + b
347 xhi = a*b/zet
348 rap = b*rab/zet
349 rbp = -a*rab/zet
350
351 ! [s|s] integral
352 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
353
354 ! Calculate the recurrence relation
355 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
356 IF (PRESENT(dab) .OR. PRESENT(ddab)) THEN
357 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 IF (PRESENT(ddab)) THEN
360 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 DO cob = ofb + 1, ncoset(lb_max)
364 bvec = indco(:, cob)
365 ib = mb + cob - ofb
366 DO coa = ofa + 1, ncoset(la_max)
367 avec = indco(:, coa)
368 ia = ma + coa - ofa
369 sval = [(rr(avec(idir), bvec(idir), idir), idir=1, 3)]
370 ! integrals
371 IF (PRESENT(sab)) THEN
372 sab(ia, ib) = f0*sval(1)*sval(2)*sval(3)
373 END IF
374 ! first derivatives
375 IF (PRESENT(dab)) THEN
376 DO idir = 1, 3
377 jdir = mod(idir, 3) + 1
378 kdir = mod(idir + 1, 3) + 1
379 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 IF (PRESENT(ddab)) THEN
384 DO icomp = 1, 6
385 idir = pairs(1, icomp)
386 jdir = pairs(2, icomp)
387 factor = sval
388 IF (idir == jdir) THEN
389 factor(idir) = -drr_ab(avec(idir), bvec(idir), idir)
390 ELSE
391 factor(idir) = drr_a(avec(idir), bvec(idir), idir)
392 factor(jdir) = drr_a(avec(jdir), bvec(jdir), jdir)
393 END IF
394 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 mb = mb + nb
402 END DO
403 ma = ma + na
404 END DO
405
406 IF (.NOT. PRESENT(rr_work)) DEALLOCATE (rr)
407 NULLIFY (rr)
408
409 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 SUBROUTINE overlap_aab(la1_max, la1_min, npgfa1, rpgfa1, zeta1, &
438 la2_max, la2_min, npgfa2, rpgfa2, zeta2, &
439 lb_max, lb_min, npgfb, rpgfb, zetb, &
440 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
469 tab = sqrt(rab2)
470
471 ! Maximum l for auxiliary integrals
472 cpassert(PRESENT(saab) .OR. PRESENT(daab) .OR. PRESENT(saba) .OR. PRESENT(daba))
473 IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
474 lma = la1_max + la2_max
475 lmb = lb_max
476 END IF
477 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
478 lma = la1_max + la2_max + 1
479 lmb = lb_max
480 END IF
481 ldrr = max(lma, lmb) + 1
482
483 ! Allocate space for auxiliary integrals
484 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
485 IF (PRESENT(daab) .OR. PRESENT(daba)) ALLOCATE (drr(0:lma - 1, 0:lmb, 3))
486
487 ! Number of integrals, check size of arrays
488 ofa1 = ncoset(la1_min - 1)
489 ofa2 = ncoset(la2_min - 1)
490 ofb = ncoset(lb_min - 1)
491 na1 = ncoset(la1_max) - ofa1
492 na2 = ncoset(la2_max) - ofa2
493 nb = ncoset(lb_max) - ofb
494 IF (PRESENT(saab)) THEN
495 cpassert((SIZE(saab, 1) >= na1*npgfa1))
496 cpassert((SIZE(saab, 2) >= na2*npgfa2))
497 cpassert((SIZE(saab, 3) >= nb*npgfb))
498 END IF
499 IF (PRESENT(daab)) THEN
500 cpassert((SIZE(daab, 1) >= na1*npgfa1))
501 cpassert((SIZE(daab, 2) >= na2*npgfa2))
502 cpassert((SIZE(daab, 3) >= nb*npgfb))
503 cpassert((SIZE(daab, 4) >= 3))
504 END IF
505 IF (PRESENT(saba)) THEN
506 cpassert((SIZE(saba, 1) >= na1*npgfa1))
507 cpassert((SIZE(saba, 2) >= nb*npgfb))
508 cpassert((SIZE(saba, 3) >= na2*npgfa2))
509 END IF
510 IF (PRESENT(daba)) THEN
511 cpassert((SIZE(daba, 1) >= na1*npgfa1))
512 cpassert((SIZE(daba, 2) >= nb*npgfb))
513 cpassert((SIZE(daba, 3) >= na2*npgfa2))
514 cpassert((SIZE(daba, 4) >= 3))
515 END IF
516
517 ! Loops over all primitive Gaussian-type functions
518 ma1 = 0
519 DO i1pgf = 1, npgfa1
520 ma2 = 0
521 DO i2pgf = 1, npgfa2
522 rpgfa = min(rpgfa1(i1pgf), rpgfa2(i2pgf))
523 mb = 0
524 DO jpgf = 1, npgfb
525 ! Distance Screening
526 IF (rpgfa + rpgfb(jpgf) < tab) THEN
527 IF (PRESENT(saab)) saab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb) = 0.0_dp
528 IF (PRESENT(daab)) daab(ma1 + 1:ma1 + na1, ma2 + 1:ma2 + na2, mb + 1:mb + nb, 1:3) = 0.0_dp
529 IF (PRESENT(saba)) saba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2) = 0.0_dp
530 IF (PRESENT(daba)) daba(ma1 + 1:ma1 + na1, mb + 1:mb + nb, ma2 + 1:ma2 + na2, 1:3) = 0.0_dp
531 mb = mb + nb
532 cycle
533 END IF
534
535 ! Calculate some prefactors
536 a = zeta1(i1pgf) + zeta2(i2pgf)
537 b = zetb(jpgf)
538 zet = a + b
539 xhi = a*b/zet
540 rap = b*rab/zet
541 rbp = -a*rab/zet
542
543 ! [ss|s] integral
544 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
545
546 ! Calculate the recurrence relation
547 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
548 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
549 CALL apply_cartesian_1d(rr, drr, la1_max + la2_max, lmb, .false., .true., a, &
550 [0.0_dp, 0.0_dp, 0.0_dp])
551 END IF
552
553 DO cob = ofb + 1, ncoset(lb_max)
554 bvec = indco(:, cob)
555 ib = mb + cob - ofb
556 DO coa2 = ofa2 + 1, ncoset(la2_max)
557 a2vec = indco(:, coa2)
558 ia2 = ma2 + coa2 - ofa2
559 DO coa1 = ofa1 + 1, ncoset(la1_max)
560 avec = indco(:, coa1) + a2vec
561 ia1 = ma1 + coa1 - ofa1
562 ! integrals
563 IF (PRESENT(saab) .OR. PRESENT(saba)) THEN
564 value = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
565 IF (PRESENT(saab)) saab(ia1, ia2, ib) = value
566 IF (PRESENT(saba)) saba(ia1, ib, ia2) = value
567 END IF
568 ! first derivatives
569 IF (PRESENT(daab) .OR. PRESENT(daba)) THEN
570 dum(1) = f0*drr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
571 dum(2) = f0*rr(avec(1), bvec(1), 1)*drr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
572 dum(3) = f0*rr(avec(1), bvec(1), 1)*rr(avec(2), bvec(2), 2)*drr(avec(3), bvec(3), 3)
573 IF (PRESENT(daab)) daab(ia1, ia2, ib, 1:3) = dum
574 IF (PRESENT(daba)) daba(ia1, ib, ia2, 1:3) = dum
575 END IF
576 END DO
577 END DO
578 END DO
579
580 mb = mb + nb
581 END DO
582 ma2 = ma2 + na2
583 END DO
584 ma1 = ma1 + na1
585 END DO
586
587 DEALLOCATE (rr)
588
589 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 SUBROUTINE overlap_abb(la_max, la_min, npgfa, rpgfa, zeta, &
616 lb1_max, lb1_min, npgfb1, rpgfb1, zetb1, &
617 lb2_max, lb2_min, npgfb2, rpgfb2, zetb2, &
618 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
642 tab = sqrt(rab2)
643
644 ! Maximum l for auxiliary integrals
645 cpassert(PRESENT(sabb) .OR. PRESENT(dabb))
646 IF (PRESENT(sabb)) THEN
647 lma = la_max
648 lmb = lb1_max + lb2_max
649 END IF
650 IF (PRESENT(dabb)) THEN
651 lma = la_max + 1
652 lmb = lb1_max + lb2_max
653 END IF
654 ldrr = max(lma, lmb) + 1
655
656 ! Allocate space for auxiliary integrals
657 ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
658 IF (PRESENT(dabb)) ALLOCATE (drr(0:lma - 1, 0:lmb, 3))
659
660 ! Number of integrals, check size of arrays
661 ofa = ncoset(la_min - 1)
662 ofb1 = ncoset(lb1_min - 1)
663 ofb2 = ncoset(lb2_min - 1)
664 na = ncoset(la_max) - ofa
665 nb1 = ncoset(lb1_max) - ofb1
666 nb2 = ncoset(lb2_max) - ofb2
667 IF (PRESENT(sabb)) THEN
668 cpassert((SIZE(sabb, 1) >= na*npgfa))
669 cpassert((SIZE(sabb, 2) >= nb1*npgfb1))
670 cpassert((SIZE(sabb, 3) >= nb2*npgfb2))
671 END IF
672 IF (PRESENT(dabb)) THEN
673 cpassert((SIZE(dabb, 1) >= na*npgfa))
674 cpassert((SIZE(dabb, 2) >= nb1*npgfb1))
675 cpassert((SIZE(dabb, 3) >= nb2*npgfb2))
676 cpassert((SIZE(dabb, 4) >= 3))
677 END IF
678
679 ! Loops over all pairs of primitive Gaussian-type functions
680 ma = 0
681 DO ipgf = 1, npgfa
682 mb1 = 0
683 DO j1pgf = 1, npgfb1
684 mb2 = 0
685 DO j2pgf = 1, npgfb2
686 ! Distance Screening
687 rpgfb = min(rpgfb1(j1pgf), rpgfb2(j2pgf))
688 IF (rpgfa(ipgf) + rpgfb < tab) THEN
689 IF (PRESENT(sabb)) sabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2) = 0.0_dp
690 IF (PRESENT(dabb)) dabb(ma + 1:ma + na, mb1 + 1:mb1 + nb1, mb2 + 1:mb2 + nb2, 1:3) = 0.0_dp
691 mb2 = mb2 + nb2
692 cycle
693 END IF
694
695 ! Calculate some prefactors
696 a = zeta(ipgf)
697 b = zetb1(j1pgf) + zetb2(j2pgf)
698 zet = a + b
699 xhi = a*b/zet
700 rap = b*rab/zet
701 rbp = -a*rab/zet
702
703 ! [s|s] integral
704 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
705
706 ! Calculate the recurrence relation
707 CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
708 IF (PRESENT(dabb)) THEN
709 CALL apply_cartesian_1d(rr, drr, la_max, lmb, .false., .true., a, &
710 [0.0_dp, 0.0_dp, 0.0_dp])
711 END IF
712
713 DO cob2 = ofb2 + 1, ncoset(lb2_max)
714 b2vec = indco(:, cob2)
715 ib2 = mb2 + cob2 - ofb2
716 DO cob1 = ofb1 + 1, ncoset(lb1_max)
717 bvec = indco(:, cob1) + b2vec
718 ib1 = mb1 + cob1 - ofb1
719 DO coa = ofa + 1, ncoset(la_max)
720 avec = indco(:, coa)
721 ia = ma + coa - ofa
722 ! integrals
723 IF (PRESENT(sabb)) THEN
724 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 IF (PRESENT(dabb)) THEN
728 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 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 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 mb2 = mb2 + nb2
737 END DO
738 mb1 = mb1 + nb1
739 END DO
740 ma = ma + na
741 END DO
742
743 DEALLOCATE (rr)
744
745 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
775
776 rpgf(1) = huge4
777 za(1) = zeta
778 zb(1) = zetb
779
780 nca = nco(la)
781 ncb = nco(lb)
782 ALLOCATE (cab(nca, ncb))
783 nsa = nso(la)
784 nsb = nso(lb)
785
786 CALL overlap_ab(la, la, 1, rpgf, za, lb, lb, 1, rpgf, zb, rab, cab)
787
788 c2sa => orbtramat(la)%c2s
789 c2sb => orbtramat(lb)%c2s
790 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
791 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
792
793 DEALLOCATE (cab)
794
795 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
823 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cab
824 REAL(kind=dp), DIMENSION(0:3, 0:3) :: fgsum
825 REAL(kind=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
826
827 nca = nco(la)
828 ncb = nco(lb)
829 ALLOCATE (cab(nca, ncb))
830 cab = 0.0_dp
831 nsa = nso(la)
832 nsb = nso(lb)
833
834 zm = min(zeta, zetb)
835 nmax = nint(1.81_dp*alat*sqrt(zm) + 1.0_dp)
836 ALLOCATE (fun(-nmax:nmax, 0:la), gun(-nmax:nmax, 0:lb), &
837 fexp(-nmax:nmax), gexp(-nmax:nmax), gval(-nmax:nmax))
838
839 oa = 1._dp/zeta
840 ob = 1._dp/zetb
841 DO i = -nmax, nmax
842 gval(i) = twopi/alat*real(i, kind=dp)
843 fexp(i) = sqrt(oa*pi)*exp(-0.25_dp*oa*gval(i)**2)
844 gexp(i) = sqrt(ob*pi)*exp(-0.25_dp*ob*gval(i)**2)
845 END DO
846 DO l = 0, la
847 IF (l == 0) THEN
848 fun(:, l) = z_one
849 ELSE IF (l == 1) THEN
850 fun(:, l) = cmplx(0.0_dp, 0.5_dp*oa*gval(:), kind=dp)
851 ELSE IF (l == 2) THEN
852 fun(:, l) = cmplx(-(0.5_dp*oa*gval(:))**2, 0.0_dp, kind=dp)
853 fun(:, l) = fun(:, l) + cmplx(0.5_dp*oa, 0.0_dp, kind=dp)
854 ELSE IF (l == 3) THEN
855 fun(:, l) = cmplx(0.0_dp, -(0.5_dp*oa*gval(:))**3, kind=dp)
856 fun(:, l) = fun(:, l) + cmplx(0.0_dp, 0.75_dp*oa*oa*gval(:), kind=dp)
857 ELSE
858 cpabort("l value too high")
859 END IF
860 END DO
861 DO l = 0, lb
862 IF (l == 0) THEN
863 gun(:, l) = z_one
864 ELSE IF (l == 1) THEN
865 gun(:, l) = cmplx(0.0_dp, 0.5_dp*ob*gval(:), kind=dp)
866 ELSE IF (l == 2) THEN
867 gun(:, l) = cmplx(-(0.5_dp*ob*gval(:))**2, 0.0_dp, kind=dp)
868 gun(:, l) = gun(:, l) + cmplx(0.5_dp*ob, 0.0_dp, kind=dp)
869 ELSE IF (l == 3) THEN
870 gun(:, l) = cmplx(0.0_dp, -(0.5_dp*ob*gval(:))**3, kind=dp)
871 gun(:, l) = gun(:, l) + cmplx(0.0_dp, 0.75_dp*ob*ob*gval(:), kind=dp)
872 ELSE
873 cpabort("l value too high")
874 END IF
875 END DO
876
877 fgsum = 0.0_dp
878 DO l1 = 0, la
879 DO l2 = 0, lb
880 zfg = sum(conjg(fun(:, l1))*fexp(:)*gun(:, l2)*gexp(:))
881 fgsum(l1, l2) = real(zfg, kind=dp)
882 END DO
883 END DO
884
885 na = ncoset(la - 1)
886 nb = ncoset(lb - 1)
887 DO ax = 0, la
888 DO ay = 0, la - ax
889 az = la - ax - ay
890 ia = coset(ax, ay, az) - na
891 DO bx = 0, lb
892 DO by = 0, lb - bx
893 bz = lb - bx - by
894 ib = coset(bx, by, bz) - nb
895 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 c2sa => orbtramat(la)%c2s
902 c2sb => orbtramat(lb)%c2s
903 sab(1:nsa, 1:nsb) = matmul(c2sa(1:nsa, 1:nca), &
904 matmul(cab(1:nca, 1:ncb), transpose(c2sb(1:nsb, 1:ncb))))
905 ovol = 1._dp/(alat**3)
906 sab(1:nsa, 1:nsb) = ovol*sab(1:nsa, 1:nsb)
907
908 DEALLOCATE (cab, fun, gun, fexp, gexp, gval)
909
910 END SUBROUTINE overlap_ab_sp
911
912END MODULE ai_overlap
ai_cartesian_operators::apply_cartesian_1d
Definition ai_cartesian_operators.F:31
ai_cartesian_operators::cartesian_ladder
Definition ai_cartesian_operators.F:27
ai_os_rr::os_rr_ovlp
Definition ai_os_rr.F:23
ai_cartesian_operators
Coordinate and basis-center derivative ladders for Cartesian Gaussian integrals. Derivatives act on a...
Definition ai_cartesian_operators.F:13
ai_cartesian_operators::cartesian_ladder_map
subroutine, public cartesian_ladder_map(lmax, lmin, axis, map, source_lmin)
Construct base, raised, lowered and degree indices for a Cartesian ladder.
Definition ai_cartesian_operators.F:46
ai_os_rr
Definition ai_os_rr.F:8
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:74
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:762
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:252
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:810
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:441
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:619
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:144
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::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
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