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 : ! **************************************************************************************************
9 : !> \brief Calculation of the G function G_n(t) for 1/R^2 operators
10 : !
11 : ! (1) J. 0. Jensen, A. H. Cameri, C. P. Vlahacos, D. Zeroka, H. F. Hameka, C. N. Merrow,
12 : ! Evaluation of one-electron integrals for arbitrary operators V(r) over cartesian Gaussians:
13 : ! Application to inverse-square distance and Yukawa operators.
14 : ! J. Comput. Chem. 14(8), 986 (1993).
15 : ! doi: 10.1002/jcc.540140814
16 : ! (2) B. Gao, A. J. Thorvaldsen, K. Ruud,
17 : ! GEN1INT: A unified procedure for the evaluation of one-electron integrals over Gaussian
18 : ! basis functions and their geometric derivatives.
19 : ! Int. J. Quantum Chem. 111(4), 858 (2011).
20 : ! doi: 10.1002/qua.22886
21 : ! (3) libgrpp : specfun_gfun.c
22 : ! (4) William Cody, Kathleen Paciorek, Henry Thacher,
23 : ! Chebyshev Approximations for Dawson's Integral,
24 : ! Mathematics of Computation,
25 : ! Volume 24, Number 109, January 1970, pages 171-178.
26 : !
27 : !> \author JHU
28 : ! **************************************************************************************************
29 : MODULE gfun
30 :
31 : USE kinds, ONLY: dp
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : ! *** Global parameters ***
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gfun'
41 :
42 : PUBLIC :: gfun_values
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param nmax ...
49 : !> \param t ...
50 : !> \param g ...
51 : ! **************************************************************************************************
52 540 : SUBROUTINE gfun_values(nmax, t, g)
53 :
54 : INTEGER, INTENT(IN) :: nmax
55 : REAL(KIND=dp), INTENT(IN) :: t
56 : REAL(KIND=dp), DIMENSION(0:nmax), INTENT(OUT) :: g
57 :
58 : INTEGER :: i
59 : REAL(KIND=dp) :: st
60 :
61 1836 : g = 0.0_dp
62 :
63 540 : IF (t <= 12.0_dp) THEN
64 : ! downward recursion
65 373 : g(nmax) = gfun_taylor(nmax, t);
66 891 : DO i = nmax, 1, -1
67 891 : g(i - 1) = (1.0_dp - 2.0_dp*t*g(i))/(2.0_dp*i - 1.0_dp)
68 : END DO
69 : ELSE
70 : ! upward recursion
71 167 : st = SQRT(t)
72 167 : g(0) = daw(st)/st
73 405 : DO i = 0, nmax - 1
74 405 : g(i + 1) = (1.0_dp - (2.0_dp*i + 1.0_dp)*g(i))/(2.0_dp*t)
75 : END DO
76 : END IF
77 :
78 540 : END SUBROUTINE gfun_values
79 :
80 : ! **************************************************************************************************
81 : !> \brief ...
82 : !> \param n ...
83 : !> \param x ...
84 : !> \return ...
85 : ! **************************************************************************************************
86 373 : FUNCTION gfun_taylor(n, x) RESULT(g)
87 : INTEGER, INTENT(IN) :: n
88 : REAL(KIND=dp), INTENT(IN) :: x
89 : REAL(KIND=dp) :: g
90 :
91 : REAL(KIND=dp), PARAMETER :: eps = 1.E-15_dp
92 :
93 : INTEGER :: k
94 : REAL(KIND=dp) :: ex, gk, tkk
95 :
96 373 : ex = EXP(-x)
97 373 : tkk = 1.0_dp
98 373 : g = ex/REAL(2*n + 1, KIND=dp)
99 4641 : DO k = 1, 100
100 4641 : tkk = tkk*x/REAL(k, KIND=dp)
101 4641 : gk = ex*tkk/REAL(2*n + 2*k + 1, KIND=dp)
102 4641 : g = g + gk
103 4641 : IF (gk < eps) EXIT
104 : END DO
105 373 : IF (gk > eps) THEN
106 0 : CPWARN("gfun_taylor did not converge")
107 : END IF
108 :
109 373 : END FUNCTION gfun_taylor
110 :
111 : !*****************************************************************************80
112 : !
113 : ! DAW evaluates Dawson's integral function.
114 : !
115 : ! Discussion:
116 : !
117 : ! This routine evaluates Dawson's integral,
118 : !
119 : ! F(x) = exp ( - x * x ) * Integral ( 0 <= t <= x ) exp ( t * t ) dt
120 : !
121 : ! for a real argument x.
122 : !
123 : ! Licensing:
124 : !
125 : ! This code is distributed under the GNU LGPL license.
126 : !
127 : ! Modified:
128 : !
129 : ! 03 April 2007
130 : !
131 : ! Author:
132 : !
133 : ! Original FORTRAN77 version by William Cody.
134 : ! FORTRAN90 version by John Burkardt.
135 : !
136 : ! Reference:
137 : !
138 : ! William Cody, Kathleen Paciorek, Henry Thacher,
139 : ! Chebyshev Approximations for Dawson's Integral,
140 : ! Mathematics of Computation,
141 : ! Volume 24, Number 109, January 1970, pages 171-178.
142 : !
143 : ! Parameters:
144 : !
145 : ! Input, real ( kind = dp ) XX, the argument of the function.
146 : !
147 : ! Output, real ( kind = dp ) DAW, the value of the function.
148 : !
149 : ! **************************************************************************************************
150 : !> \brief ...
151 : !> \param xx ...
152 : !> \return ...
153 : ! **************************************************************************************************
154 167 : FUNCTION daw(xx)
155 :
156 : REAL(kind=dp) :: xx, daw
157 :
158 : INTEGER :: i
159 : REAL(kind=dp) :: frac, one225, p1(10), p2(10), p3(10), &
160 : p4(10), q1(10), q2(9), q3(9), q4(9), &
161 : six25, sump, sumq, two5, w2, x, &
162 : xlarge, xmax, xsmall, y
163 :
164 : !
165 : ! Mathematical constants.
166 : !
167 : DATA six25/6.25D+00/
168 : DATA one225/12.25d0/
169 : DATA two5/25.0d0/
170 : !
171 : ! Machine-dependent constants
172 : !
173 : DATA xsmall/1.05d-08/
174 : DATA xlarge/9.49d+07/
175 : DATA xmax/2.24d+307/
176 : !
177 : ! Coefficients for R(9,9) approximation for |x| < 2.5
178 : !
179 : DATA p1/-2.69020398788704782410d-12, 4.18572065374337710778d-10, &
180 : -1.34848304455939419963d-08, 9.28264872583444852976d-07, &
181 : -1.23877783329049120592d-05, 4.07205792429155826266d-04, &
182 : -2.84388121441008500446d-03, 4.70139022887204722217d-02, &
183 : -1.38868086253931995101d-01, 1.00000000000000000004d+00/
184 : DATA q1/1.71257170854690554214d-10, 1.19266846372297253797d-08, &
185 : 4.32287827678631772231d-07, 1.03867633767414421898d-05, &
186 : 1.78910965284246249340d-04, 2.26061077235076703171d-03, &
187 : 2.07422774641447644725d-02, 1.32212955897210128811d-01, &
188 : 5.27798580412734677256d-01, 1.00000000000000000000d+00/
189 : !
190 : ! Coefficients for R(9,9) approximation in J-fraction form
191 : ! for x in [2.5, 3.5)
192 : !
193 : DATA p2/-1.70953804700855494930d+00, -3.79258977271042880786d+01, &
194 : 2.61935631268825992835d+01, 1.25808703738951251885d+01, &
195 : -2.27571829525075891337d+01, 4.56604250725163310122d+00, &
196 : -7.33080089896402870750d+00, 4.65842087940015295573d+01, &
197 : -1.73717177843672791149d+01, 5.00260183622027967838d-01/
198 : DATA q2/1.82180093313514478378d+00, 1.10067081034515532891d+03, &
199 : -7.08465686676573000364d+00, 4.53642111102577727153d+02, &
200 : 4.06209742218935689922d+01, 3.02890110610122663923d+02, &
201 : 1.70641269745236227356d+02, 9.51190923960381458747d+02, &
202 : 2.06522691539642105009d-01/
203 : !
204 : ! Coefficients for R(9,9) approximation in J-fraction form
205 : ! for x in [3.5, 5.0]
206 : !
207 : DATA p3/-4.55169503255094815112d+00, -1.86647123338493852582d+01, &
208 : -7.36315669126830526754d+00, -6.68407240337696756838d+01, &
209 : 4.84507265081491452130d+01, 2.69790586735467649969d+01, &
210 : -3.35044149820592449072d+01, 7.50964459838919612289d+00, &
211 : -1.48432341823343965307d+00, 4.99999810924858824981d-01/
212 : DATA q3/4.47820908025971749852d+01, 9.98607198039452081913d+01, &
213 : 1.40238373126149385228d+01, 3.48817758822286353588d+03, &
214 : -9.18871385293215873406d+00, 1.24018500009917163023d+03, &
215 : -6.88024952504512254535d+01, -2.31251575385145143070d+00, &
216 : 2.50041492369922381761d-01/
217 : !
218 : ! Coefficients for R(9,9) approximation in J-fraction form
219 : ! for 5.0 < |x|.
220 : !
221 : DATA p4/-8.11753647558432685797d+00, -3.84043882477454453430d+01, &
222 : -2.23787669028751886675d+01, -2.88301992467056105854d+01, &
223 : -5.99085540418222002197d+00, -1.13867365736066102577d+01, &
224 : -6.52828727526980741590d+00, -4.50002293000355585708d+00, &
225 : -2.50000000088955834952d+00, 5.00000000000000488400d-01/
226 : DATA q4/2.69382300417238816428d+02, 5.04198958742465752861d+01, &
227 : 6.11539671480115846173d+01, 2.08210246935564547889d+02, &
228 : 1.97325365692316183531d+01, -1.22097010558934838708d+01, &
229 : -6.99732735041547247161d+00, -2.49999970104184464568d+00, &
230 : 7.49999999999027092188d-01/
231 :
232 167 : x = xx
233 :
234 167 : IF (xlarge < ABS(x)) THEN
235 :
236 0 : IF (ABS(x) <= xmax) THEN
237 0 : daw = 0.5D+00/x
238 : ELSE
239 : daw = 0.0D+00
240 : END IF
241 :
242 167 : ELSE IF (ABS(x) < xsmall) THEN
243 :
244 : daw = x
245 :
246 : ELSE
247 :
248 167 : y = x*x
249 : !
250 : ! ABS(X) < 2.5.
251 : !
252 167 : IF (y < six25) THEN
253 :
254 0 : sump = p1(1)
255 0 : sumq = q1(1)
256 0 : DO i = 2, 10
257 0 : sump = sump*y + p1(i)
258 0 : sumq = sumq*y + q1(i)
259 : END DO
260 :
261 0 : daw = x*sump/sumq
262 : !
263 : ! 2.5 <= ABS(X) < 3.5.
264 : !
265 167 : ELSE IF (y < one225) THEN
266 :
267 : frac = 0.0D+00
268 0 : DO i = 1, 9
269 0 : frac = q2(i)/(p2(i) + y + frac)
270 : END DO
271 :
272 0 : daw = (p2(10) + frac)/x
273 : !
274 : ! 3.5 <= ABS(X) < 5.0.
275 : !
276 167 : ELSE IF (y < two5) THEN
277 :
278 : frac = 0.0D+00
279 120 : DO i = 1, 9
280 120 : frac = q3(i)/(p3(i) + y + frac)
281 : END DO
282 :
283 12 : daw = (p3(10) + frac)/x
284 :
285 : ELSE
286 : !
287 : ! 5.0 <= ABS(X) <= XLARGE.
288 : !
289 155 : w2 = 1.0D+00/x/x
290 :
291 155 : frac = 0.0D+00
292 1550 : DO i = 1, 9
293 1550 : frac = q4(i)/(p4(i) + y + frac)
294 : END DO
295 155 : frac = p4(10) + frac
296 :
297 155 : daw = (0.5D+00 + 0.5D+00*w2*frac)/x
298 :
299 : END IF
300 :
301 : END IF
302 :
303 167 : END FUNCTION daw
304 :
305 : END MODULE gfun
|
