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