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 : !> \brief Three-center overlap integrals over Cartesian Gaussian functions.
9 : !> \author Matthias Krack (26.06.2001)
10 : !> \par Literature
11 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
12 : ! **************************************************************************************************
13 : MODULE ai_overlap3
14 :
15 : USE ai_cartesian_operators, ONLY: cartesian_ladder,&
16 : cartesian_ladder_map
17 : USE ai_os_rr, ONLY: os_rr_ovlp3
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: pi
20 : USE orbital_pointers, ONLY: indco,&
21 : ncoset
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3'
28 :
29 : PUBLIC :: overlap3
30 :
31 : CONTAINS
32 :
33 : ! **************************************************************************************************
34 : !> \brief Calculation of three-center overlap integrals [a|b|c] over primitive
35 : !> Cartesian Gaussian functions using separable one-dimensional tables.
36 : !> \param la_max_set ...
37 : !> \param npgfa ...
38 : !> \param zeta ...
39 : !> \param rpgfa ...
40 : !> \param la_min_set ...
41 : !> \param lb_max_set ...
42 : !> \param npgfb ...
43 : !> \param zetb ...
44 : !> \param rpgfb ...
45 : !> \param lb_min_set ...
46 : !> \param lc_max_set ...
47 : !> \param npgfc ...
48 : !> \param zetc ...
49 : !> \param rpgfc ...
50 : !> \param lc_min_set ...
51 : !> \param rab displacement B - A
52 : !> \param dab distance between A and B
53 : !> \param rac displacement C - A
54 : !> \param dac distance between A and C
55 : !> \param rbc displacement C - B
56 : !> \param dbc distance between B and C
57 : !> \param sabc integrals [a|b|c]
58 : !> \param sdabc derivative [da/dAi|b|c]
59 : !> \param sabdc derivative [a|b|dc/dCi]
60 : !> \param int_abc_ext largest absolute integral in the requested shell ranges
61 : !> \par History
62 : !> 05.2014 created (Dorothea Golze)
63 : !> \author Dorothea Golze
64 : ! **************************************************************************************************
65 24 : SUBROUTINE overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
66 48 : lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
67 24 : lc_max_set, npgfc, zetc, rpgfc, lc_min_set, &
68 24 : rab, dab, rac, dac, rbc, dbc, sabc, &
69 24 : sdabc, sabdc, int_abc_ext)
70 :
71 : INTEGER, INTENT(IN) :: la_max_set, npgfa
72 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
73 : INTEGER, INTENT(IN) :: la_min_set, lb_max_set, npgfb
74 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
75 : INTEGER, INTENT(IN) :: lb_min_set, lc_max_set, npgfc
76 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
77 : INTEGER, INTENT(IN) :: lc_min_set
78 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
79 : REAL(KIND=dp), INTENT(IN) :: dab
80 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
81 : REAL(KIND=dp), INTENT(IN) :: dac
82 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
83 : REAL(KIND=dp), INTENT(IN) :: dbc
84 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabc
85 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
86 : INTENT(INOUT), OPTIONAL :: sdabc, sabdc
87 : REAL(dp), INTENT(OUT), OPTIONAL :: int_abc_ext
88 :
89 : CHARACTER(len=*), PARAMETER :: routineN = 'overlap3'
90 :
91 : INTEGER :: a_first, a_last, avec(3), axis, b_first, b_last, bvec(3), c_first, c_last, coa, &
92 : cob, coc, cvec(3), handle, ia, ib, ic, ipgf, jaxis, jpgf, kaxis, kpgf, la_max, lc_max, &
93 48 : na, nb, nc, map_c(4, lc_max_set + 1), map_a(4, la_max_set + 1)
94 : REAL(KIND=dp) :: f0, f1, factor, rcp2, value, zet, zetg, &
95 : zetp
96 24 : REAL(KIND=dp), ALLOCATABLE :: derivative_a(:, :, :, :), &
97 24 : derivative_c(:, :, :, :), &
98 : rr(:, :, :, :)
99 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp, rcg, rcp
100 :
101 24 : CALL timeset(routineN, handle)
102 :
103 24 : la_max = la_max_set
104 24 : lc_max = lc_max_set
105 24 : IF (PRESENT(sdabc)) la_max = la_max + 1
106 24 : IF (PRESENT(sabdc)) lc_max = lc_max + 1
107 144 : ALLOCATE (rr(0:la_max, 0:lb_max_set, 0:lc_max, 3))
108 :
109 24 : a_first = ncoset(la_min_set - 1) + 1
110 24 : a_last = ncoset(la_max_set)
111 24 : b_first = ncoset(lb_min_set - 1) + 1
112 24 : b_last = ncoset(lb_max_set)
113 24 : c_first = ncoset(lc_min_set - 1) + 1
114 24 : c_last = ncoset(lc_max_set)
115 :
116 : ! Clear the requested derivative extents, including screened blocks and lower-shell padding.
117 24 : IF (PRESENT(sdabc)) THEN
118 121296 : sdabc(1:npgfa*a_last, 1:npgfb*b_last, 1:npgfc*c_last, 1:3) = 0.0_dp
119 72 : ALLOCATE (derivative_a(0:la_max_set, 0:lb_max_set, 0:lc_max_set, 3))
120 12 : CALL cartesian_ladder_map(la_max_set, 0, 0, map_a)
121 : END IF
122 24 : IF (PRESENT(sabdc)) THEN
123 121296 : sabdc(1:npgfa*a_last, 1:npgfb*b_last, 1:npgfc*c_last, 1:3) = 0.0_dp
124 72 : ALLOCATE (derivative_c(0:la_max_set, 0:lb_max_set, 0:lc_max_set, 3))
125 12 : CALL cartesian_ladder_map(lc_max_set, 0, 0, map_c)
126 : END IF
127 24 : IF (PRESENT(int_abc_ext)) int_abc_ext = 0.0_dp
128 :
129 144 : DO ipgf = 1, npgfa
130 120 : na = (ipgf - 1)*a_last
131 744 : DO jpgf = 1, npgfb
132 600 : nb = (jpgf - 1)*b_last
133 600 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
134 0 : sabc(na + 1:na + a_last, nb + 1:nb + b_last, 1:npgfc*c_last) = 0.0_dp
135 : CYCLE
136 : END IF
137 600 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
138 600 : f1 = zetb(jpgf)*zetp
139 2400 : rcp = f1*rab - rac
140 2400 : rcp2 = DOT_PRODUCT(rcp, rcp)
141 2820 : DO kpgf = 1, npgfc
142 2100 : nc = (kpgf - 1)*c_last
143 2100 : IF ((rpgfb(jpgf) + rpgfc(kpgf) < dbc) .OR. &
144 : (rpgfa(ipgf) + rpgfc(kpgf) < dac)) THEN
145 0 : sabc(na + 1:na + a_last, nb + 1:nb + b_last, nc + 1:nc + c_last) = 0.0_dp
146 : CYCLE
147 : END IF
148 :
149 2100 : zet = zeta(ipgf) + zetb(jpgf) + zetc(kpgf)
150 2100 : zetg = 1.0_dp/zet
151 2100 : f0 = (pi*zetg)**1.5_dp
152 2100 : f0 = f0*EXP(-(zeta(ipgf)*f1*dab*dab + zetc(kpgf)*zetg*rcp2/zetp))
153 :
154 : ! G is the product center of all three Gaussians. The 1D tables are normalized;
155 : ! the three-dimensional [s|s|s] prefactor f0 is applied only on reconstruction.
156 8400 : rap = zetg*(zetb(jpgf)*rab + zetc(kpgf)*rac)
157 8400 : rbp = zetg*(-zeta(ipgf)*rab + zetc(kpgf)*rbc)
158 8400 : rcg = -zetg*(zeta(ipgf)*rac + zetb(jpgf)*rbc)
159 2100 : CALL os_rr_ovlp3(rap, la_max, rbp, lb_max_set, rcg, lc_max, zet, rr)
160 :
161 : ! Differentiate the original A and C basis functions after three-center table construction.
162 2100 : IF (PRESENT(sdabc)) THEN
163 4200 : DO axis = 1, 3
164 8250 : DO ic = 0, lc_max_set
165 : CALL cartesian_ladder(rr(:, :, ic, axis), derivative_a(:, :, ic, axis), &
166 7200 : map_a, .TRUE., zeta(ipgf), 0.0_dp, .FALSE.)
167 : END DO
168 : END DO
169 : END IF
170 2100 : IF (PRESENT(sabdc)) THEN
171 4200 : DO axis = 1, 3
172 10500 : DO ib = 0, lb_max_set
173 : CALL cartesian_ladder(rr(0:la_max_set, ib, :, axis), derivative_c(:, ib, :, axis), &
174 9450 : map_c, .TRUE., zetc(kpgf), 0.0_dp, .TRUE.)
175 : END DO
176 : END DO
177 : END IF
178 :
179 7200 : DO coc = c_first, c_last
180 18000 : cvec = indco(:, coc)
181 4500 : ic = nc + coc
182 24600 : DO cob = b_first, b_last
183 72000 : bvec = indco(:, cob)
184 18000 : ib = nb + cob
185 94500 : DO coa = a_first, a_last
186 288000 : avec = indco(:, coa)
187 72000 : ia = na + coa
188 : value = f0*rr(avec(1), bvec(1), cvec(1), 1)* &
189 72000 : rr(avec(2), bvec(2), cvec(2), 2)*rr(avec(3), bvec(3), cvec(3), 3)
190 72000 : sabc(ia, ib, ic) = value
191 72000 : IF (PRESENT(int_abc_ext)) int_abc_ext = MAX(int_abc_ext, ABS(value))
192 90000 : IF (PRESENT(sdabc) .OR. PRESENT(sabdc)) THEN
193 144000 : DO axis = 1, 3
194 108000 : jaxis = MOD(axis, 3) + 1
195 108000 : kaxis = MOD(axis + 1, 3) + 1
196 : factor = f0*rr(avec(jaxis), bvec(jaxis), cvec(jaxis), jaxis)* &
197 108000 : rr(avec(kaxis), bvec(kaxis), cvec(kaxis), kaxis)
198 108000 : IF (PRESENT(sdabc)) THEN
199 108000 : sdabc(ia, ib, ic, axis) = factor*derivative_a(avec(axis), bvec(axis), cvec(axis), axis)
200 : END IF
201 144000 : IF (PRESENT(sabdc)) THEN
202 108000 : sabdc(ia, ib, ic, axis) = factor*derivative_c(avec(axis), bvec(axis), cvec(axis), axis)
203 : END IF
204 : END DO
205 : END IF
206 : END DO
207 : END DO
208 : END DO
209 : END DO
210 : END DO
211 : END DO
212 :
213 24 : CALL timestop(handle)
214 :
215 48 : END SUBROUTINE overlap3
216 :
217 : END MODULE ai_overlap3
|
