(git:374b731)
Loading...
Searching...
No Matches
ai_os_rr.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8MODULE 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
23CONTAINS
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 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 g = 0.5_dp/zet
51 rr(0, 0, 1) = 1.0_dp
52 rr(0, 0, 2) = 1.0_dp
53 rr(0, 0, 3) = 1.0_dp
54 !
55 ! recursion along la for lb=0
56 !
57 IF (la_max .GT. 0) THEN
58 rr(1, 0, 1) = rap(1)
59 rr(1, 0, 2) = rap(2)
60 rr(1, 0, 3) = rap(3)
61 !
62 DO la = 1, la_max - 1
63 lap1 = la + 1
64 lam1 = la - 1
65 rr(lap1, 0, 1) = real(la, dp)*g*rr(lam1, 0, 1) + rap(1)*rr(la, 0, 1)
66 rr(lap1, 0, 2) = real(la, dp)*g*rr(lam1, 0, 2) + rap(2)*rr(la, 0, 2)
67 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 IF (lb_max .GT. 0) THEN
74 rr(0, 1, 1) = rbp(1)
75 rr(0, 1, 2) = rbp(2)
76 rr(0, 1, 3) = rbp(3)
77 !
78 DO la = 1, la_max
79 lam1 = la - 1
80 rr(la, 1, 1) = real(la, dp)*g*rr(lam1, 0, 1) + rbp(1)*rr(la, 0, 1)
81 rr(la, 1, 2) = real(la, dp)*g*rr(lam1, 0, 2) + rbp(2)*rr(la, 0, 2)
82 rr(la, 1, 3) = real(la, dp)*g*rr(lam1, 0, 3) + rbp(3)*rr(la, 0, 3)
83 END DO
84 !
85 DO lb = 1, lb_max - 1
86 lbp1 = lb + 1
87 lbm1 = lb - 1
88 rr(0, lbp1, 1) = real(lb, dp)*g*rr(0, lbm1, 1) + rbp(1)*rr(0, lb, 1)
89 rr(0, lbp1, 2) = real(lb, dp)*g*rr(0, lbm1, 2) + rbp(2)*rr(0, lb, 2)
90 rr(0, lbp1, 3) = real(lb, dp)*g*rr(0, lbm1, 3) + rbp(3)*rr(0, lb, 3)
91 DO la = 1, la_max
92 lam1 = la - 1
93 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 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 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 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 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 mmax = la_max + lb_max
135 g = 0.5_dp/zet
136 !
137 ! rr(0:mmax) should be initialized before
138 !
139 rcp2 = rcp(1)**2 + rcp(2)**2 + rcp(3)**2
140 t = zet*rcp2
141 CALL fgamma(mmax, t, rr(0:mmax, 1, 1))
142 !
143 ! recursion in la with lb=0
144 !
145 DO la = 1, la_max
146 DO ax = 0, la
147 DO ay = 0, la - ax
148 az = la - ax - ay
149 coa = coset(ax, ay, az)
150 coa1x = coset(max(ax - 1, 0), ay, az)
151 coa1y = coset(ax, max(ay - 1, 0), az)
152 coa1z = coset(ax, ay, max(az - 1, 0))
153 coa2x = coset(max(ax - 2, 0), ay, az)
154 coa2y = coset(ax, max(ay - 2, 0), az)
155 coa2z = coset(ax, ay, max(az - 2, 0))
156 IF (az .GT. 0) THEN
157 DO m = 0, mmax - la
158 rr(m, coa, 1) = rap(3)*rr(m, coa1z, 1) - rcp(3)*rr(m + 1, coa1z, 1)
159 END DO
160 IF (az .GT. 1) THEN
161 DO m = 0, mmax - la
162 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 ELSEIF (ay .GT. 0) THEN
166 DO m = 0, mmax - la
167 rr(m, coa, 1) = rap(2)*rr(m, coa1y, 1) - rcp(2)*rr(m + 1, coa1y, 1)
168 END DO
169 IF (ay .GT. 1) THEN
170 DO m = 0, mmax - la
171 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 ELSEIF (ax .GT. 0) THEN
175 DO m = 0, mmax - la
176 rr(m, coa, 1) = rap(1)*rr(m, coa1x, 1) - rcp(1)*rr(m + 1, coa1x, 1)
177 END DO
178 IF (ax .GT. 1) THEN
179 DO m = 0, mmax - la
180 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 cpabort("")
185 END IF
186 END DO
187 END DO
188 END DO
189 !
190 ! recursion in lb with all possible la
191 !
192 DO la = 0, la_max
193 DO ax = 0, la
194 DO ay = 0, la - ax
195 az = la - ax - ay
196 coa = coset(ax, ay, az)
197 coa1x = coset(max(ax - 1, 0), ay, az)
198 coa1y = coset(ax, max(ay - 1, 0), az)
199 coa1z = coset(ax, ay, max(az - 1, 0))
200 coa2x = coset(max(ax - 2, 0), ay, az)
201 coa2y = coset(ax, max(ay - 2, 0), az)
202 coa2z = coset(ax, ay, max(az - 2, 0))
203 DO lb = 1, lb_max
204 DO bx = 0, lb
205 DO by = 0, lb - bx
206 bz = lb - bx - by
207 cob = coset(bx, by, bz)
208 cob1x = coset(max(bx - 1, 0), by, bz)
209 cob1y = coset(bx, max(by - 1, 0), bz)
210 cob1z = coset(bx, by, max(bz - 1, 0))
211 cob2x = coset(max(bx - 2, 0), by, bz)
212 cob2y = coset(bx, max(by - 2, 0), bz)
213 cob2z = coset(bx, by, max(bz - 2, 0))
214 IF (bz .GT. 0) THEN
215 DO m = 0, mmax - la - lb
216 rr(m, coa, cob) = rbp(3)*rr(m, coa, cob1z) - rcp(3)*rr(m + 1, coa, cob1z)
217 END DO
218 IF (bz .GT. 1) THEN
219 DO m = 0, mmax - la - lb
220 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 IF (az .GT. 0) THEN
224 DO m = 0, mmax - la - lb
225 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 ELSEIF (by .GT. 0) THEN
229 DO m = 0, mmax - la - lb
230 rr(m, coa, cob) = rbp(2)*rr(m, coa, cob1y) - rcp(2)*rr(m + 1, coa, cob1y)
231 END DO
232 IF (by .GT. 1) THEN
233 DO m = 0, mmax - la - lb
234 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 IF (ay .GT. 0) THEN
238 DO m = 0, mmax - la - lb
239 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 ELSEIF (bx .GT. 0) THEN
243 DO m = 0, mmax - la - lb
244 rr(m, coa, cob) = rbp(1)*rr(m, coa, cob1x) - rcp(1)*rr(m + 1, coa, cob1x)
245 END DO
246 IF (bx .GT. 1) THEN
247 DO m = 0, mmax - la - lb
248 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 IF (ax .GT. 0) THEN
252 DO m = 0, mmax - la - lb
253 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 cpabort("")
258 END IF
259 END DO
260 END DO
261 END DO
262 END DO
263 END DO
264 END DO
265 !
266 END SUBROUTINE os_rr_coul
267
268END MODULE ai_os_rr
ai_os_rr
Definition ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition ai_os_rr.F:39
ai_os_rr::os_rr_coul
subroutine, public os_rr_coul(rap, la_max, rbp, lb_max, rcp, zet, ldrr1, ldrr2, rr)
Calculation of the Obara-Saika recurrence relation for 1/r_C.
Definition ai_os_rr.F:118
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45