(git:b195825)
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 
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  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 
268 END 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