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 : MODULE ai_os_rr
9 :
10 : USE gamma, ONLY: fgamma => fgamma_0
11 : USE kinds, ONLY: dp
12 : USE orbital_pointers, ONLY: coset, indco, ncoset
13 : #include "../base/base_uses.f90"
14 :
15 : IMPLICIT NONE
16 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_os_rr'
17 : PRIVATE
18 :
19 : ! *** Public subroutines ***
20 :
21 : PUBLIC :: os_rr_ovlp, os_rr_ovlp3, os_rr_coul, os_rr_radial, os_rr_coulomb3
22 :
23 : INTERFACE os_rr_ovlp
24 : MODULE PROCEDURE os_rr_ovlp_real, os_rr_ovlp_complex
25 : END INTERFACE os_rr_ovlp
26 :
27 : INTERFACE os_rr_ovlp3
28 : MODULE PROCEDURE os_rr_ovlp3_real, os_rr_ovlp3_complex
29 : END INTERFACE os_rr_ovlp3
30 :
31 : INTERFACE os_rr_radial
32 : MODULE PROCEDURE os_rr_radial_abm
33 : END INTERFACE os_rr_radial
34 :
35 : CONTAINS
36 :
37 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
38 : ! **************************************************************************************************
39 : !> \brief Calculation of the basic Obara-Saika recurrence relation
40 : !> for real or complex product-center displacements. Each Cartesian
41 : !> table is normalized to a unit zeroth-degree integral.
42 : !> \param rap product-center displacement relative to center A
43 : !> \param la_max ...
44 : !> \param rbp product-center displacement relative to center B
45 : !> \param lb_max ...
46 : !> \param zet real sum of the Gaussian exponents
47 : !> \param ldrr ...
48 : !> \param rr ...
49 : !> \date 02.03.2009
50 : !> \author VW
51 : !> \version 1.0
52 : ! **************************************************************************************************
53 155824387 : SUBROUTINE os_rr_ovlp_${suffix}$ (rap, la_max, rbp, lb_max, zet, ldrr, rr)
54 : ${integral_type}$ (dp), DIMENSION(3), INTENT(IN) :: rap
55 : INTEGER, INTENT(IN) :: la_max
56 : ${integral_type}$ (dp), DIMENSION(3), INTENT(IN) :: rbp
57 : INTEGER, INTENT(IN) :: lb_max
58 : REAL(dp), INTENT(IN) :: zet
59 : INTEGER, INTENT(IN) :: ldrr
60 : ${integral_type}$ (dp), DIMENSION(0:ldrr - 1, 0:ldrr - 1, 3) :: rr
61 :
62 : INTEGER :: idir
63 :
64 623297548 : DO idir = 1, 3
65 623297548 : CALL os_rr_ovlp_1d_${suffix}$ (rap(idir), la_max, rbp(idir), lb_max, zet, rr(:, :, idir))
66 : END DO
67 155824387 : END SUBROUTINE os_rr_ovlp_${suffix}$
68 :
69 : ! **************************************************************************************************
70 : !> \brief Normalized one-dimensional two-center Gaussian moments.
71 : !> \param rap product-center displacement relative to A
72 : !> \param la_max highest A degree
73 : !> \param rbp product-center displacement relative to B
74 : !> \param lb_max highest B degree
75 : !> \param zet sum of Gaussian exponents
76 : !> \param rr normalized moments with zeroth-degree value one
77 : ! **************************************************************************************************
78 467479461 : PURE SUBROUTINE os_rr_ovlp_1d_${suffix}$ (rap, la_max, rbp, lb_max, zet, rr)
79 : ${integral_type}$ (dp), INTENT(IN) :: rap, rbp
80 : INTEGER, INTENT(IN) :: la_max, lb_max
81 : REAL(dp), INTENT(IN) :: zet
82 : ${integral_type}$ (dp), INTENT(INOUT) :: rr(0:, 0:)
83 :
84 : INTEGER :: la, lam1, lap1, lb, lbm1, lbp1
85 : REAL(dp) :: g
86 :
87 467479461 : g = 0.5_dp/zet
88 467479461 : rr(0, 0) = 1.0_dp
89 : ! recursion along la for lb=0
90 467479461 : IF (la_max > 0) THEN
91 313054812 : rr(1, 0) = rap
92 484858527 : DO la = 1, la_max - 1
93 171803715 : lap1 = la + 1
94 171803715 : lam1 = la - 1
95 484858527 : rr(lap1, 0) = REAL(la, dp)*g*rr(lam1, 0) + rap*rr(la, 0)
96 : END DO
97 : END IF
98 : ! recursion along lb for all la
99 467479461 : IF (lb_max > 0) THEN
100 254571438 : rr(0, 1) = rbp
101 581535981 : DO la = 1, la_max
102 326964543 : lam1 = la - 1
103 581535981 : rr(la, 1) = REAL(la, dp)*g*rr(lam1, 0) + rbp*rr(la, 0)
104 : END DO
105 367886766 : DO lb = 1, lb_max - 1
106 113315328 : lbp1 = lb + 1
107 113315328 : lbm1 = lb - 1
108 113315328 : rr(0, lbp1) = REAL(lb, dp)*g*rr(0, lbm1) + rbp*rr(0, lb)
109 553670691 : DO la = 1, la_max
110 185783925 : lam1 = la - 1
111 : rr(la, lbp1) = g*(REAL(la, dp)*rr(lam1, lb) + REAL(lb, dp)*rr(la, lbm1)) + &
112 299099253 : rbp*rr(la, lb)
113 : END DO
114 : END DO
115 : END IF
116 467479461 : END SUBROUTINE os_rr_ovlp_1d_${suffix}$
117 :
118 : ! **************************************************************************************************
119 : !> \brief Normalized separable three-center Gaussian overlap tables.
120 : !> The C recurrence uses the directly supplied G-C displacement.
121 : !> \param rap product-center displacement G-A
122 : !> \param la_max highest A degree
123 : !> \param rbp product-center displacement G-B
124 : !> \param lb_max highest B degree
125 : !> \param rcp product-center displacement G-C
126 : !> \param lc_max highest C degree
127 : !> \param zet sum of all three Gaussian exponents
128 : !> \param rr normalized tables indexed by A/B/C degrees and Cartesian axis
129 : ! **************************************************************************************************
130 2100 : SUBROUTINE os_rr_ovlp3_${suffix}$ (rap, la_max, rbp, lb_max, rcp, lc_max, zet, rr)
131 : ${integral_type}$ (dp), INTENT(IN) :: rap(3), rbp(3), rcp(3)
132 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
133 : REAL(dp), INTENT(IN) :: zet
134 : ${integral_type}$ (dp), INTENT(INOUT) :: rr(0:, 0:, 0:, :)
135 :
136 : INTEGER :: a, b, c, idir
137 : REAL(dp) :: g
138 :
139 2100 : g = 0.5_dp/zet
140 8400 : DO idir = 1, 3
141 6300 : CALL os_rr_ovlp_1d_${suffix}$ (rap(idir), la_max, rbp(idir), lb_max, zet, rr(:, :, 0, idir))
142 13350 : DO c = 0, lc_max - 1
143 21150 : DO b = 0, lb_max
144 42750 : DO a = 0, la_max
145 : rr(a, b, c + 1, idir) = rcp(idir)*rr(a, b, c, idir) + &
146 : g*(REAL(a, dp)*rr(MAX(0, a - 1), b, c, idir) + &
147 : REAL(b, dp)*rr(a, MAX(0, b - 1), c, idir) + &
148 37800 : REAL(c, dp)*rr(a, b, MAX(0, c - 1), idir))
149 : END DO
150 : END DO
151 : END DO
152 : END DO
153 2100 : END SUBROUTINE os_rr_ovlp3_${suffix}$
154 : #:endfor
155 :
156 : #:for layout in ['abm', 'mab']
157 : #:set idx = (lambda a, b, m: a + ', ' + b + ', ' + m) if layout == 'abm' else (lambda a, b, m: m + ', ' + a + ', ' + b)
158 : #:set bounds = ':, :, 0:' if layout == 'abm' else '0:, :, :'
159 : #:set outer, inner = ('m', 'ia') if layout == 'abm' else ('ia', 'm')
160 : #:set starts = {'m': '0', 'ia': 'ncoset(la - 1) + 1'}
161 : #:set stops = {'m': 'mmax - la', 'ia': 'ncoset(la)'}
162 : ! **************************************************************************************************
163 : !> \brief Propagate caller-supplied radial seeds through the A/B angular recurrence.
164 : !> Both supported layouts use the same algebra and require no transpose.
165 : !> \param rap displacement P-A
166 : !> \param la_max highest A shell
167 : !> \param rbp displacement P-B
168 : !> \param lb_max highest B shell
169 : !> \param rpw displacement multiplying the next radial order
170 : !> \param zetp reciprocal exponent sum on A/B
171 : !> \param eta reduction factor of the next radial order (one for point Coulomb)
172 : !> \param mmax highest supplied seed order, including subsequent C recurrence
173 : !> \param rr caller-initialized s/s seeds and resulting angular table
174 : ! **************************************************************************************************
175 12813366 : SUBROUTINE os_rr_radial_${layout}$ (rap, la_max, rbp, lb_max, rpw, zetp, eta, mmax, rr)
176 : REAL(dp), INTENT(IN) :: rap(3), rbp(3), rpw(3), zetp, eta
177 : INTEGER, INTENT(IN) :: la_max, lb_max, mmax
178 : REAL(dp), INTENT(INOUT) :: rr(${bounds}$)
179 :
180 : INTEGER :: axis, degree, ia, ia1, ia2, ib, ib1, ib2, ico, iinner, iouter, la, lb, m, power(3)
181 25626732 : INTEGER :: lowered(3, ncoset(MAX(la_max, lb_max)))
182 : REAL(dp) :: g
183 :
184 12813366 : g = 0.5_dp*zetp
185 94009440 : DO ico = 1, ncoset(MAX(la_max, lb_max))
186 324784296 : power = indco(:, ico)
187 : lowered(:, ico) = [coset(MAX(0, power(1) - 1), power(2), power(3)), &
188 : coset(power(1), MAX(0, power(2) - 1), power(3)), &
189 337597662 : coset(power(1), power(2), MAX(0, power(3) - 1))]
190 : END DO
191 27382296 : DO la = 1, la_max
192 57825250 : DO iouter = ${starts[outer]}$, ${stops[outer]}$
193 164905166 : DO iinner = ${starts[inner]}$, ${stops[inner]}$
194 119893282 : ${outer}$ = iouter
195 119893282 : ${inner}$ = iinner
196 119893282 : axis = 3
197 119893282 : IF (indco(3, ia) == 0) axis = 2
198 119893282 : IF (indco(axis, ia) == 0) axis = 1
199 119893282 : ia1 = lowered(axis, ia)
200 119893282 : ia2 = lowered(axis, ia1)
201 119893282 : degree = indco(axis, ia) - 1
202 : rr(${idx('ia', '1', 'm')}$) = rap(axis)*rr(${idx('ia1', '1', 'm')}$) + &
203 119893282 : rpw(axis)*rr(${idx('ia1', '1', 'm + 1')}$)
204 150336236 : IF (degree > 0) THEN
205 : rr(${idx('ia', '1', 'm')}$) = rr(${idx('ia', '1', 'm')}$) + g*REAL(degree, dp)* &
206 28056045 : (rr(${idx('ia2', '1', 'm')}$) - eta*rr(${idx('ia2', '1', 'm + 1')}$))
207 : END IF
208 : END DO
209 : END DO
210 : END DO
211 21851316 : DO lb = 1, lb_max
212 52241070 : DO ib = ncoset(lb - 1) + 1, ncoset(lb)
213 30389754 : axis = 3
214 30389754 : IF (indco(3, ib) == 0) axis = 2
215 30389754 : IF (indco(axis, ib) == 0) axis = 1
216 30389754 : ib1 = lowered(axis, ib)
217 30389754 : ib2 = lowered(axis, ib1)
218 30389754 : degree = indco(axis, ib) - 1
219 105830440 : DO la = 0, la_max
220 : #:set stops = {'m': 'mmax - la - lb', 'ia': 'ncoset(la)'}
221 220340264 : DO iouter = ${starts[outer]}$, ${stops[outer]}$
222 467049616 : DO iinner = ${starts[inner]}$, ${stops[inner]}$
223 277099106 : ${outer}$ = iouter
224 277099106 : ${inner}$ = iinner
225 277099106 : ia1 = lowered(axis, ia)
226 : rr(${idx('ia', 'ib', 'm')}$) = rbp(axis)*rr(${idx('ia', 'ib1', 'm')}$) + &
227 277099106 : rpw(axis)*rr(${idx('ia', 'ib1', 'm + 1')}$)
228 277099106 : IF (degree > 0) THEN
229 : rr(${idx('ia', 'ib', 'm')}$) = rr(${idx('ia', 'ib', 'm')}$) + g*REAL(degree, dp)* &
230 29220211 : (rr(${idx('ia', 'ib2', 'm')}$) - eta*rr(${idx('ia', 'ib2', 'm + 1')}$))
231 : END IF
232 400646880 : IF (indco(axis, ia) > 0) THEN
233 : rr(${idx('ia', 'ib', 'm')}$) = rr(${idx('ia', 'ib', 'm')}$) + g*REAL(indco(axis, ia), dp)* &
234 85808025 : (rr(${idx('ia1', 'ib1', 'm')}$) - eta*rr(${idx('ia1', 'ib1', 'm + 1')}$))
235 : END IF
236 : END DO
237 : END DO
238 : END DO
239 : END DO
240 : END DO
241 12813366 : END SUBROUTINE os_rr_radial_${layout}$
242 : #:endfor
243 :
244 : ! **************************************************************************************************
245 : !> \brief Angular three-center Coulomb recurrence from caller-supplied s/s/s seeds.
246 : !> \param rap displacement P-A
247 : !> \param la_max highest A shell
248 : !> \param rbp displacement P-B
249 : !> \param lb_max highest B shell
250 : !> \param rcp displacement P-C
251 : !> \param lc_max highest C shell
252 : !> \param zetp reciprocal exponent sum on A/B
253 : !> \param zetq reciprocal exponent on C
254 : !> \param zetw reciprocal exponent sum on A/B/C
255 : !> \param rr caller-initialized seeds and angular table, with radial order last
256 : ! **************************************************************************************************
257 4500 : SUBROUTINE os_rr_coulomb3(rap, la_max, rbp, lb_max, rcp, lc_max, zetp, zetq, zetw, rr)
258 : REAL(dp), INTENT(IN) :: rap(3), rbp(3), rcp(3), zetp, zetq, zetw
259 : INTEGER, INTENT(IN) :: la_max, lb_max, lc_max
260 : REAL(dp), INTENT(INOUT) :: rr(:, :, :, 0:)
261 :
262 : INTEGER :: axis, ia, ia1, ib, ib1, ic, ic1, ic2, ico, la, lb, lc, m, mmax, power(3)
263 4500 : INTEGER :: lowered(3, ncoset(MAX(la_max, lb_max, lc_max)))
264 : REAL(dp) :: eta, eta_c, g, g_c, rcw(3), rpw(3)
265 :
266 4500 : mmax = la_max + lb_max + lc_max
267 4500 : eta = zetw/zetq
268 4500 : eta_c = zetw/zetp
269 18000 : rpw = -eta*rcp
270 : ! Direct evaluation avoids subtracting P-C and its nearly equal Coulomb shift.
271 18000 : rcw = eta_c*rcp
272 4500 : CALL os_rr_radial(rap, la_max, rbp, lb_max, rpw, zetp, eta, mmax, rr(:, :, 1, :))
273 4500 : IF (lc_max == 0) RETURN
274 2700 : g = 0.5_dp*zetw
275 2700 : g_c = 0.5_dp*zetq
276 27126 : DO ico = 1, ncoset(MAX(la_max, lb_max, lc_max))
277 97704 : power = indco(:, ico)
278 : lowered(:, ico) = [coset(MAX(0, power(1) - 1), power(2), power(3)), &
279 : coset(power(1), MAX(0, power(2) - 1), power(3)), &
280 100404 : coset(power(1), power(2), MAX(0, power(3) - 1))]
281 : END DO
282 7200 : DO lc = 1, lc_max
283 27900 : DO ic = ncoset(lc - 1) + 1, ncoset(lc)
284 20700 : axis = 3
285 20700 : IF (indco(3, ic) == 0) axis = 2
286 20700 : IF (indco(axis, ic) == 0) axis = 1
287 20700 : ic1 = lowered(axis, ic)
288 20700 : ic2 = lowered(axis, ic1)
289 56940 : DO lb = 0, lb_max
290 102258 : DO la = 0, la_max
291 191718 : DO m = 0, mmax - la - lb - lc
292 354246 : DO ib = ncoset(lb - 1) + 1, ncoset(lb)
293 194268 : ib1 = lowered(axis, ib)
294 649212 : DO ia = ncoset(la - 1) + 1, ncoset(la)
295 344784 : ia1 = lowered(axis, ia)
296 344784 : rr(ia, ib, ic, m) = rcw(axis)*rr(ia, ib, ic1, m + 1)
297 344784 : IF (indco(axis, ic) > 1) THEN
298 : rr(ia, ib, ic, m) = rr(ia, ib, ic, m) + g_c*REAL(indco(axis, ic) - 1, dp)* &
299 101376 : (rr(ia, ib, ic2, m) - eta_c*rr(ia, ib, ic2, m + 1))
300 : END IF
301 344784 : IF (indco(axis, ia) > 0) THEN
302 : rr(ia, ib, ic, m) = rr(ia, ib, ic, m) + &
303 79212 : g*REAL(indco(axis, ia), dp)*rr(ia1, ib, ic1, m + 1)
304 : END IF
305 539052 : IF (indco(axis, ib) > 0) THEN
306 : rr(ia, ib, ic, m) = rr(ia, ib, ic, m) + &
307 79212 : g*REAL(indco(axis, ib), dp)*rr(ia, ib1, ic1, m + 1)
308 : END IF
309 : END DO
310 : END DO
311 : END DO
312 : END DO
313 : END DO
314 : END DO
315 : END DO
316 : END SUBROUTINE os_rr_coulomb3
317 :
318 : ! **************************************************************************************************
319 : !> \brief Calculation of the Obara-Saika recurrence relation for 1/r_C
320 : !> \param rap ...
321 : !> \param la_max ...
322 : !> \param rbp ...
323 : !> \param lb_max ...
324 : !> \param rcp ...
325 : !> \param zet ...
326 : !> \param ldrr1 ...
327 : !> \param ldrr2 ...
328 : !> \param rr ...
329 : !> \date 02.03.2009
330 : !> \author VW
331 : !> \version 1.0
332 : ! **************************************************************************************************
333 66912 : SUBROUTINE os_rr_coul(rap, la_max, rbp, lb_max, rcp, zet, ldrr1, ldrr2, rr)
334 : REAL(dp), DIMENSION(3), INTENT(IN) :: rap
335 : INTEGER, INTENT(IN) :: la_max
336 : REAL(dp), DIMENSION(3), INTENT(IN) :: rbp
337 : INTEGER, INTENT(IN) :: lb_max
338 : REAL(dp), DIMENSION(3), INTENT(IN) :: rcp
339 : REAL(dp), INTENT(IN) :: zet
340 : INTEGER, INTENT(IN) :: ldrr1, ldrr2
341 : REAL(dp), DIMENSION(0:ldrr1 - 1, ldrr2, *), &
342 : INTENT(INOUT) :: rr
343 :
344 : INTEGER :: mmax
345 : REAL(dp) :: t
346 :
347 66912 : mmax = la_max + lb_max
348 267648 : t = zet*DOT_PRODUCT(rcp, rcp)
349 66912 : CALL fgamma(mmax, t, rr(0:mmax, 1, 1))
350 : CALL os_rr_radial_mab(rap, la_max, rbp, lb_max, -rcp, 1.0_dp/zet, 1.0_dp, mmax, &
351 267648 : rr(0:mmax, 1:ncoset(la_max), 1:ncoset(lb_max)))
352 :
353 66912 : END SUBROUTINE os_rr_coul
354 :
355 : END MODULE ai_os_rr
|
