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
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_coul
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 : CONTAINS
28 :
29 : #:for suffix, integral_type in [('real', 'REAL'), ('complex', 'COMPLEX')]
30 : ! **************************************************************************************************
31 : !> \brief Calculation of the basic Obara-Saika recurrence relation
32 : !> for real or complex product-center displacements. Each Cartesian
33 : !> table is normalized to a unit zeroth-degree integral.
34 : !> \param rap product-center displacement relative to center A
35 : !> \param la_max ...
36 : !> \param rbp product-center displacement relative to center B
37 : !> \param lb_max ...
38 : !> \param zet real sum of the Gaussian exponents
39 : !> \param ldrr ...
40 : !> \param rr ...
41 : !> \date 02.03.2009
42 : !> \author VW
43 : !> \version 1.0
44 : ! **************************************************************************************************
45 155719587 : SUBROUTINE os_rr_ovlp_${suffix}$ (rap, la_max, rbp, lb_max, zet, ldrr, rr)
46 : ${integral_type}$ (dp), DIMENSION(3), INTENT(IN) :: rap
47 : INTEGER, INTENT(IN) :: la_max
48 : ${integral_type}$ (dp), DIMENSION(3), INTENT(IN) :: rbp
49 : INTEGER, INTENT(IN) :: lb_max
50 : REAL(dp), INTENT(IN) :: zet
51 : INTEGER, INTENT(IN) :: ldrr
52 : ${integral_type}$ (dp), DIMENSION(0:ldrr - 1, 0:ldrr - 1, 3) :: rr
53 :
54 : INTEGER :: idir, la, lam1, lap1, lb, lbm1, lbp1
55 : REAL(dp) :: g
56 :
57 155719587 : g = 0.5_dp/zet
58 622878348 : DO idir = 1, 3
59 467158761 : rr(0, 0, idir) = 1.0_dp
60 : ! recursion along la for lb=0
61 467158761 : IF (la_max > 0) THEN
62 312752196 : rr(1, 0, idir) = rap(idir)
63 484382286 : DO la = 1, la_max - 1
64 171630090 : lap1 = la + 1
65 171630090 : lam1 = la - 1
66 484382286 : rr(lap1, 0, idir) = REAL(la, dp)*g*rr(lam1, 0, idir) + rap(idir)*rr(la, 0, idir)
67 : END DO
68 : END IF
69 : ! recursion along lb for all la
70 622878348 : IF (lb_max > 0) THEN
71 254241465 : rr(0, 1, idir) = rbp(idir)
72 580713990 : DO la = 1, la_max
73 326472525 : lam1 = la - 1
74 580713990 : rr(la, 1, idir) = REAL(la, dp)*g*rr(lam1, 0, idir) + rbp(idir)*rr(la, 0, idir)
75 : END DO
76 367332789 : DO lb = 1, lb_max - 1
77 113091324 : lbp1 = lb + 1
78 113091324 : lbm1 = lb - 1
79 113091324 : rr(0, lbp1, idir) = REAL(lb, dp)*g*rr(0, lbm1, idir) + rbp(idir)*rr(0, lb, idir)
80 552695019 : DO la = 1, la_max
81 185362230 : lam1 = la - 1
82 : rr(la, lbp1, idir) = g*(REAL(la, dp)*rr(lam1, lb, idir) + REAL(lb, dp)*rr(la, lbm1, idir)) + &
83 298453554 : rbp(idir)*rr(la, lb, idir)
84 : END DO
85 : END DO
86 : END IF
87 : END DO
88 155719587 : END SUBROUTINE os_rr_ovlp_${suffix}$
89 : #:endfor
90 :
91 : ! **************************************************************************************************
92 : !> \brief Calculation of the Obara-Saika recurrence relation for 1/r_C
93 : !> \param rap ...
94 : !> \param la_max ...
95 : !> \param rbp ...
96 : !> \param lb_max ...
97 : !> \param rcp ...
98 : !> \param zet ...
99 : !> \param ldrr1 ...
100 : !> \param ldrr2 ...
101 : !> \param rr ...
102 : !> \date 02.03.2009
103 : !> \author VW
104 : !> \version 1.0
105 : ! **************************************************************************************************
106 66912 : SUBROUTINE os_rr_coul(rap, la_max, rbp, lb_max, rcp, zet, ldrr1, ldrr2, rr)
107 : REAL(dp), DIMENSION(3), INTENT(IN) :: rap
108 : INTEGER, INTENT(IN) :: la_max
109 : REAL(dp), DIMENSION(3), INTENT(IN) :: rbp
110 : INTEGER, INTENT(IN) :: lb_max
111 : REAL(dp), DIMENSION(3), INTENT(IN) :: rcp
112 : REAL(dp), INTENT(IN) :: zet
113 : INTEGER, INTENT(IN) :: ldrr1, ldrr2
114 : REAL(dp), DIMENSION(0:ldrr1 - 1, ldrr2, *), &
115 : INTENT(INOUT) :: rr
116 :
117 : INTEGER :: ax, ay, az, bx, by, bz, coa, coa1x, &
118 : coa1y, coa1z, coa2x, coa2y, coa2z, &
119 : cob, cob1x, cob1y, cob1z, cob2x, &
120 : cob2y, cob2z, la, lb, m, mmax
121 : REAL(dp) :: g, rcp2, t
122 :
123 66912 : mmax = la_max + lb_max
124 66912 : g = 0.5_dp/zet
125 : !
126 : ! rr(0:mmax) should be initialized before
127 : !
128 66912 : rcp2 = rcp(1)**2 + rcp(2)**2 + rcp(3)**2
129 66912 : t = zet*rcp2
130 66912 : CALL fgamma(mmax, t, rr(0:mmax, 1, 1))
131 : !
132 : ! recursion in la with lb=0
133 : !
134 160924 : DO la = 1, la_max
135 378748 : DO ax = 0, la
136 685972 : DO ay = 0, la - ax
137 374136 : az = la - ax - ay
138 374136 : coa = coset(ax, ay, az)
139 374136 : coa1x = coset(MAX(ax - 1, 0), ay, az)
140 374136 : coa1y = coset(ax, MAX(ay - 1, 0), az)
141 374136 : coa1z = coset(ax, ay, MAX(az - 1, 0))
142 374136 : coa2x = coset(MAX(ax - 2, 0), ay, az)
143 374136 : coa2y = coset(ax, MAX(ay - 2, 0), az)
144 374136 : coa2z = coset(ax, ay, MAX(az - 2, 0))
145 591960 : IF (az > 0) THEN
146 612831 : DO m = 0, mmax - la
147 612831 : rr(m, coa, 1) = rap(3)*rr(m, coa1z, 1) - rcp(3)*rr(m + 1, coa1z, 1)
148 : END DO
149 156312 : IF (az > 1) THEN
150 130919 : DO m = 0, mmax - la
151 130919 : rr(m, coa, 1) = rr(m, coa, 1) + g*REAL(az - 1, dp)*(rr(m, coa2z, 1) - rr(m + 1, coa2z, 1))
152 : END DO
153 : END IF
154 217824 : ELSE IF (ay > 0) THEN
155 481912 : DO m = 0, mmax - la
156 481912 : rr(m, coa, 1) = rap(2)*rr(m, coa1y, 1) - rcp(2)*rr(m + 1, coa1y, 1)
157 : END DO
158 123812 : IF (ay > 1) THEN
159 119115 : DO m = 0, mmax - la
160 119115 : rr(m, coa, 1) = rr(m, coa, 1) + g*REAL(ay - 1, dp)*(rr(m, coa2y, 1) - rr(m + 1, coa2y, 1))
161 : END DO
162 : END IF
163 94012 : ELSE IF (ax > 0) THEN
164 362797 : DO m = 0, mmax - la
165 362797 : rr(m, coa, 1) = rap(1)*rr(m, coa1x, 1) - rcp(1)*rr(m + 1, coa1x, 1)
166 : END DO
167 94012 : IF (ax > 1) THEN
168 107311 : DO m = 0, mmax - la
169 107311 : rr(m, coa, 1) = rr(m, coa, 1) + g*REAL(ax - 1, dp)*(rr(m, coa2x, 1) - rr(m + 1, coa2x, 1))
170 : END DO
171 : END IF
172 : ELSE
173 0 : CPABORT("None of ax, ay, az is positive")
174 : END IF
175 : END DO
176 : END DO
177 : END DO
178 : !
179 : ! recursion in lb with all possible la
180 : !
181 227836 : DO la = 0, la_max
182 512572 : DO ax = 0, la
183 886708 : DO ay = 0, la - ax
184 441048 : az = la - ax - ay
185 441048 : coa = coset(ax, ay, az)
186 441048 : coa1x = coset(MAX(ax - 1, 0), ay, az)
187 441048 : coa1y = coset(ax, MAX(ay - 1, 0), az)
188 441048 : coa1z = coset(ax, ay, MAX(az - 1, 0))
189 441048 : coa2x = coset(MAX(ax - 2, 0), ay, az)
190 441048 : coa2y = coset(ax, MAX(ay - 2, 0), az)
191 441048 : coa2z = coset(ax, ay, MAX(az - 2, 0))
192 1432114 : DO lb = 1, lb_max
193 2864230 : DO bx = 0, lb
194 5493658 : DO by = 0, lb - bx
195 3070476 : bz = lb - bx - by
196 3070476 : cob = coset(bx, by, bz)
197 3070476 : cob1x = coset(MAX(bx - 1, 0), by, bz)
198 3070476 : cob1y = coset(bx, MAX(by - 1, 0), bz)
199 3070476 : cob1z = coset(bx, by, MAX(bz - 1, 0))
200 3070476 : cob2x = coset(MAX(bx - 2, 0), by, bz)
201 3070476 : cob2y = coset(bx, MAX(by - 2, 0), bz)
202 3070476 : cob2z = coset(bx, by, MAX(bz - 2, 0))
203 4787328 : IF (bz > 0) THEN
204 3783551 : DO m = 0, mmax - la - lb
205 3783551 : rr(m, coa, cob) = rbp(3)*rr(m, coa, cob1z) - rcp(3)*rr(m + 1, coa, cob1z)
206 : END DO
207 1353624 : IF (bz > 1) THEN
208 917182 : DO m = 0, mmax - la - lb
209 917182 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(bz - 1, dp)*(rr(m, coa, cob2z) - rr(m + 1, coa, cob2z))
210 : END DO
211 : END IF
212 1353624 : IF (az > 0) THEN
213 1390815 : DO m = 0, mmax - la - lb
214 1390815 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(az, dp)*(rr(m, coa1z, cob1z) - rr(m + 1, coa1z, cob1z))
215 : END DO
216 : END IF
217 1716852 : ELSE IF (by > 0) THEN
218 2866369 : DO m = 0, mmax - la - lb
219 2866369 : rr(m, coa, cob) = rbp(2)*rr(m, coa, cob1y) - rcp(2)*rr(m + 1, coa, cob1y)
220 : END DO
221 1010522 : IF (by > 1) THEN
222 814887 : DO m = 0, mmax - la - lb
223 814887 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(by - 1, dp)*(rr(m, coa, cob2y) - rr(m + 1, coa, cob2y))
224 : END DO
225 : END IF
226 1010522 : IF (ay > 0) THEN
227 1037336 : DO m = 0, mmax - la - lb
228 1037336 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(ay, dp)*(rr(m, coa1y, cob1y) - rr(m + 1, coa1y, cob1y))
229 : END DO
230 : END IF
231 706330 : ELSE IF (bx > 0) THEN
232 2051482 : DO m = 0, mmax - la - lb
233 2051482 : rr(m, coa, cob) = rbp(1)*rr(m, coa, cob1x) - rcp(1)*rr(m + 1, coa, cob1x)
234 : END DO
235 706330 : IF (bx > 1) THEN
236 712592 : DO m = 0, mmax - la - lb
237 712592 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(bx - 1, dp)*(rr(m, coa, cob2x) - rr(m + 1, coa, cob2x))
238 : END DO
239 : END IF
240 706330 : IF (ax > 0) THEN
241 725904 : DO m = 0, mmax - la - lb
242 725904 : rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(ax, dp)*(rr(m, coa1x, cob1x) - rr(m + 1, coa1x, cob1x))
243 : END DO
244 : END IF
245 : ELSE
246 0 : CPABORT("None of bx, by, bz is positive")
247 : END IF
248 : END DO
249 : END DO
250 : END DO
251 : END DO
252 : END DO
253 : END DO
254 : !
255 66912 : END SUBROUTINE os_rr_coul
256 :
257 : END MODULE ai_os_rr
|
