Line data Source code
1 : MODULE m_util
2 : USE m_juDFT
3 : ! error and warning codes for intgrf function
4 : INTEGER, PARAMETER :: NO_ERROR = 0
5 : INTEGER, PARAMETER :: NEGATIVE_EXPONENT_ERROR = 2
6 : ! return type of the pure intgrf function
7 :
8 : INTERFACE derivative
9 : MODULE PROCEDURE derivative_t, derivative_nt
10 : END INTERFACE
11 :
12 : CONTAINS
13 :
14 : ! Calculates Gaunt coefficients, i.e. the integrals of three spherical harmonics
15 : ! integral ( conjg(Y(l1,m1)) * Y(l2,m2) * conjg(Y(l3,m3)) )
16 : ! They are also the coefficients C(l1,l2,l3,m1,m2,m3) in
17 : ! conjg(Y(l1,m1)) * Y(l2,m2) = sum(l3,m3) C(l1,l2,l3,m1,m2,m3) Y(l3,m3)
18 : ! fac contains factorial up to maxfac, i.e. fac(i)= i!
19 : ! sfac contains square root of fac, i.e. sfac(i)= sqrt(i!)
20 :
21 1392704 : FUNCTION gaunt(l1, l2, l3, m1, m2, m3, maxfac, fac, sfac)
22 :
23 : USE m_constants, ONLY: pimach
24 :
25 : IMPLICIT NONE
26 :
27 : REAL :: gaunt
28 : INTEGER, INTENT(IN) :: l1, l2, l3, m1, m2, m3, maxfac
29 : REAL , INTENT(IN) :: fac(0:maxfac)
30 : REAL , INTENT(IN) :: sfac(0:maxfac)
31 :
32 1392704 : gaunt = 0
33 1392704 : IF (m3 /= m2 - m1) RETURN
34 1392704 : IF (abs(m1) > l1) RETURN
35 1392704 : IF (abs(m2) > l2) RETURN
36 1392704 : IF (abs(m3) > l3) RETURN
37 1392704 : IF (l3 < abs(l1 - l2) .OR. l3 > l1 + l2) RETURN
38 : gaunt = (-1)**(m1 + m3)* &
39 : sqrt((2*l1 + 1)*(2*l2 + 1)*(2*l3 + 1)/pimach()/4)* &
40 : wigner3j(l1, l2, l3, -m1, m2, -m3, maxfac, fac, sfac)* &
41 1392704 : wigner3j(l1, l2, l3, 0, 0, 0, maxfac, fac, sfac)
42 1392704 : END FUNCTION gaunt
43 :
44 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45 :
46 : ! Calculates the Wigner 3j symbols using Racah's formula
47 :
48 2785408 : FUNCTION wigner3j(l1, l2, l3, m1, m2, m3, maxfac, fac, sfac)
49 :
50 : IMPLICIT NONE
51 :
52 : REAL :: wigner3j
53 : INTEGER, INTENT(IN) :: l1, l2, l3, m1, m2, m3, maxfac
54 : REAL , INTENT(IN) :: fac(0:maxfac)
55 : REAL , INTENT(IN) :: sfac(0:maxfac)
56 : ! - local -
57 : INTEGER :: tmin, tmax, t, f1, f2, f3, f4, f5
58 :
59 2785408 : wigner3j = 0
60 :
61 : ! The following IF clauses should be in the calling routine and commented here.
62 : ! if(-m3.ne.m1+m2) return
63 : ! if(abs(m1).gt.l1) return
64 : ! if(abs(m2).gt.l2) return
65 : ! if(abs(m3).gt.l3) return
66 : ! if(l3.lt.abs(l1-l2).or.l3.gt.l1+l2) return
67 :
68 2785408 : f1 = l3 - l2 + m1
69 2785408 : f2 = l3 - l1 - m2
70 2785408 : f3 = l1 + l2 - l3
71 2785408 : f4 = l1 - m1
72 2785408 : f5 = l2 + m2
73 2785408 : tmin = max(0, -f1, -f2) ! The arguments to fac (see below)
74 2785408 : tmax = min(f3, f4, f5) ! must not be negative.
75 :
76 : ! The following line is only for testing and should be removed at a later time.
77 2785408 : IF (tmax - tmin /= min(l1 + m1, l1 - m1, l2 + m2, l2 - m2, l3 + m3, l3 - m3, &
78 : l1 + l2 - l3, l1 - l2 + l3, -l1 + l2 + l3)) &
79 0 : call judft_error("wigner3j: Number of terms incorrect.")
80 :
81 2785408 : IF (tmin <= tmax) THEN
82 6464112 : DO t = tmin, tmax
83 : wigner3j = wigner3j + (-1)**t/ &
84 : (fac(t)*fac(f1 + t)*fac(f2 + t) &
85 6464112 : *fac(f3 - t)*fac(f4 - t)*fac(f5 - t))
86 : END DO
87 : wigner3j = wigner3j*(-1)**(l1 - l2 - m3)*sfac(l1 + l2 - l3) &
88 : *sfac(l1 - l2 + l3)*sfac(-l1 + l2 + l3) &
89 : /sfac(l1 + l2 + l3 + 1)* &
90 : sfac(l1 + m1)*sfac(l1 - m1)* &
91 : sfac(l2 + m2)*sfac(l2 - m2)* &
92 2785408 : sfac(l3 + m3)*sfac(l3 - m3)
93 : END IF
94 2785408 : END FUNCTION wigner3j
95 :
96 :
97 :
98 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99 :
100 : ! Calculates the primitive of f, on grid(itypein):
101 :
102 : ! r
103 : ! primf(r) = integral f(r') dr' ( r = grid point )
104 : ! 0
105 :
106 : ! If itypein is negative, the primitive
107 :
108 : ! R
109 : ! primf(r) = integral f(r') dr' ( R = MT sphere radius )
110 : ! r
111 :
112 : ! is calculated instead.
113 :
114 : ! -----------------------------
115 :
116 : ! Fast calculation of primitive.
117 : ! Only Lagrange integration is used
118 :
119 9776 : SUBROUTINE primitivef(primf, fin, rmsh, dx, jri, jmtd, itypein, ntype)
120 :
121 : USE m_constants
122 :
123 : IMPLICIT NONE
124 :
125 : ! - scalars -
126 : INTEGER, INTENT(IN) :: itypein, jmtd, ntype
127 :
128 : ! - arrays -
129 : INTEGER, INTENT(IN) :: jri(ntype)
130 : REAL, INTENT(OUT) :: primf(jri(abs(itypein)))
131 : REAL, INTENT(IN) :: fin(jri(abs(itypein)))
132 : REAL, INTENT(IN) :: rmsh(jmtd, ntype), dx(ntype)
133 :
134 : ! - local scalars -
135 : INTEGER :: itype, n, i, n0
136 : REAL :: h, x, h1
137 : REAL :: intgr, r1, a, dr
138 :
139 : ! - local arrays -
140 : REAL :: fr(7)
141 9776 : REAL :: f(jri(abs(itypein)))
142 9776 : REAL :: r(jri(abs(itypein)))
143 : REAL, PARAMETER :: lagrange(7, 6) = reshape( &
144 : (/19087., 65112., -46461., 37504., -20211., 6312., -863., &
145 : -863., 25128., 46989., -16256., 7299., -2088., 271., &
146 : & 271., -2760., 30819., 37504., -6771., 1608., -191., &
147 : -191., 1608., -6771., 37504., 30819., -2760., 271., &
148 : & 271., -2088., 7299., -16256., 46989., 25128., -863., &
149 : -863., 6312., -20211., 37504., -46461., 65112., 19087./), &
150 : (/7, 6/))
151 :
152 9776 : itype = abs(itypein)
153 :
154 7619760 : primf = 0
155 :
156 9776 : n = jri(itype)
157 9776 : h = dx(itype)
158 :
159 9776 : IF (itypein > 0) THEN
160 4888 : r1 = rmsh(1, itype) ! perform outward integration
161 3809880 : f = fin ! (from 0 to r)
162 : ELSE
163 4888 : r1 = rmsh(jri(itype), itype) ! perform inward integration
164 4888 : h = -h ! (from MT sphere radius to r)
165 3809880 : f = fin(n:1:-1) !
166 : END IF
167 :
168 : ! integral from 0 to r1 approximated by leading term in power series expansion (only if h>0)
169 9776 : IF (h > 0 .AND. f(1)*f(2) > 1e-10) THEN
170 24 : IF (f(2) == f(1)) THEN
171 0 : intgr = r1*f(1)
172 : ELSE
173 24 : x = (f(3) - f(2))/(f(2) - f(1))
174 24 : a = (f(2) - x*f(1))/(1 - x)
175 24 : x = log(x)/h
176 24 : IF (x < 0) THEN
177 0 : IF (x > -1) WRITE (oUnit, '(A,ES9.1)') &
178 : '+intgr: Warning! Negative &exponent x in'// &
179 0 : 'extrapolation a+c*r**x:', x
180 0 : IF (x <= -1) WRITE (oUnit, '(A,ES9.1)') 'intgr: Negative exponent,'// &
181 0 : 'x in extrapolation a+c*r**x:', x
182 0 : IF (x <= -1) call judft_error("intgr:Negative exponent x in extrapolation")
183 : END IF
184 24 : intgr = r1*(f(1) + x*a)/(x + 1)
185 : END IF
186 : ELSE
187 : intgr = 0
188 : END IF
189 :
190 9776 : primf(1) = intgr
191 9776 : dr = exp(h)
192 9776 : r(1) = r1
193 9776 : n0 = 1
194 9776 : h1 = h/60480
195 :
196 : ! Lagrange integration from r(n0) to r(n0+5)
197 8901144 : 1 DO i = 2, 7
198 8901144 : r(i) = r(i - 1)*dr
199 : END DO
200 10172736 : fr = f(n0:n0 + 6)*r(:7)
201 8901144 : DO i = 1, 6
202 61036416 : intgr = intgr + h1*dot_product(lagrange(:, i), fr)
203 8901144 : IF (primf(n0 + i) == 0) primf(n0 + i) = intgr ! avoid double-definition
204 : END DO
205 1271592 : IF (n0 + 12 <= n) THEN
206 1254480 : r(1) = r(7)
207 1254480 : n0 = n0 + 6
208 1254480 : GOTO 1
209 17112 : ELSE IF (n0 + 6 < n) THEN
210 7336 : r(1) = r(n - 5 - n0)
211 7336 : n0 = n - 6
212 7336 : intgr = primf(n - 6)
213 7336 : GOTO 1
214 : END IF
215 :
216 9776 : IF (itypein < 0) THEN !
217 7619760 : primf = -primf(n:1:-1) ! Inward integration
218 : END IF !
219 :
220 9776 : END SUBROUTINE primitivef
221 :
222 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
223 :
224 : ! unction modulo1 maps kpoint into first BZ
225 0 : FUNCTION modulo1(kpoint, nkpt3)
226 :
227 : USE m_constants
228 :
229 : IMPLICIT NONE
230 :
231 : INTEGER, INTENT(IN) :: nkpt3(3)
232 : REAL , INTENT(IN) :: kpoint(3)
233 : REAL :: modulo1(3)
234 : INTEGER :: help(3)
235 :
236 0 : modulo1 = kpoint*nkpt3
237 0 : help = nint(modulo1)
238 0 : IF (any(abs(help - modulo1) > 1e-10)) THEN
239 0 : WRITE (oUnit, '(A,F5.3,2('','',F5.3),A)') 'modulo1: argument (', kpoint, &
240 0 : ') is not an element of the k-point set.'
241 : CALL juDFT_error( &
242 : 'modulo1: argument not an element of k-point set.', &
243 0 : calledby='util:modulo1')
244 : END IF
245 0 : modulo1 = modulo(help, nkpt3)*1.0/nkpt3
246 :
247 : END FUNCTION modulo1
248 :
249 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250 :
251 : ! Returns derivative of f in df.
252 :
253 0 : SUBROUTINE derivative_t(df, f, atoms, itype)
254 : USE m_types_atoms
255 : IMPLICIT NONE
256 : REAL, INTENT(IN) :: f(:)
257 : REAL, INTENT(OUT) :: df(:)
258 : TYPE(t_atoms), INTENT(IN)::atoms
259 : INTEGER, INTENT(IN) :: itype
260 :
261 : call derivative_nt(df, f, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, &
262 0 : atoms%ntype, itype)
263 0 : END SUBROUTINE
264 :
265 0 : SUBROUTINE derivative_nt(df, f, jmtd, jri, dx, rmsh, ntype, itype)
266 :
267 : IMPLICIT NONE
268 :
269 : INTEGER, INTENT(IN) :: ntype, itype, jmtd
270 : INTEGER, INTENT(IN) :: jri(ntype)
271 : REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
272 : REAL, INTENT(IN) :: f(jri(itype))
273 : REAL, INTENT(OUT) :: df(jri(itype))
274 : REAL :: x, h, r, d21, d32, d43, d31, d42, d41, df1, df2
275 : INTEGER :: i, n
276 :
277 0 : n = jri(itype)
278 0 : h = dx(itype)
279 0 : r = rmsh(1, itype)
280 : ! use power series expansion a+c**x for first point
281 0 : IF (f(2) == f(1)) THEN
282 0 : df(1) = 0.0
283 : ELSE
284 0 : x = (f(3) - f(2))/(f(2) - f(1))
285 0 : df(1) = (f(2) - f(1))/(x - 1)*log(x)/h/r
286 : END IF
287 : ! use Lagrange interpolation of 3rd order for all other points (and averaging)
288 0 : d21 = r*(exp(h) - 1); d32 = d21*exp(h); d43 = d32*exp(h)
289 0 : d31 = d21 + d32; d42 = d32 + d43
290 0 : d41 = d31 + d43
291 : df(2) = -d32*d42/(d21*d31*d41)*f(1) &
292 : + (1.0/d21 - 1.0/d32 - 1.0/d42)*f(2) &
293 : + d21*d42/(d31*d32*d43)*f(3) &
294 0 : - d21*d32/(d41*d42*d43)*f(4)
295 : df1 = d32*d43/(d21*d31*d41)*f(1) &
296 : - d31*d43/(d21*d32*d42)*f(2) &
297 : + (1.0/d31 + 1.0/d32 - 1.0/d43)*f(3) &
298 0 : + d31*d32/(d41*d42*d43)*f(4)
299 0 : DO i = 3, n - 2
300 0 : d21 = d32; d32 = d43; d43 = d43*exp(h)
301 0 : d31 = d42; d42 = d42*exp(h)
302 0 : d41 = d41*exp(h)
303 : df2 = -d32*d42/(d21*d31*d41)*f(i - 1) &
304 : + (1.0/d21 - 1.0/d32 - 1.0/d42)*f(i) &
305 : + d21*d42/(d31*d32*d43)*f(i + 1) &
306 0 : - d21*d32/(d41*d42*d43)*f(i + 2)
307 0 : df(i) = (df1 + df2)/2
308 : df1 = d32*d43/(d21*d31*d41)*f(i - 1) &
309 : - d31*d43/(d21*d32*d42)*f(i) &
310 : + (1.0/d31 + 1.0/d32 - 1.0/d43)*f(i + 1) &
311 0 : + d31*d32/(d41*d42*d43)*f(i + 2)
312 : END DO
313 0 : df(n - 1) = df1
314 : df(n) = -d42*d43/(d21*d31*d41)*f(n - 3) &
315 : + d41*d43/(d21*d32*d42)*f(n - 2) &
316 : - d41*d42/(d31*d32*d43)*f(n - 1) &
317 0 : + (1.0/d41 + 1.0/d42 + 1.0/d43)*f(n)
318 :
319 0 : END SUBROUTINE derivative_nt
320 :
321 :
322 : ! Calculates the spherical Bessel functions of orders 0 to l at x
323 : ! by backward recurrence using j_l(x) = (2l+3)/x j_l+1(x) - j_l+2(x) .
324 : ! (Starting points are calculated according to Zhang, Min,
325 : ! "Computation of Special Functions".)
326 0 : pure SUBROUTINE sphbessel(sphbes, x, l)
327 :
328 : IMPLICIT NONE
329 :
330 : INTEGER, INTENT(IN) :: l
331 : REAL , INTENT(IN) :: x
332 : REAL , INTENT(INOUT) :: sphbes(0:l)
333 : REAL :: s0, s1, f, f0, f1, cs
334 : INTEGER :: ll, lsta, lmax, msta2
335 :
336 : ! IF( x.lt.0 ) THEN
337 : ! call judft_error("sphbes: negative argument (bug?).")
338 : ! ELSE
339 0 : IF (x == 0) THEN
340 0 : sphbes(0) = 1.0
341 0 : DO ll = 1, l
342 0 : sphbes(ll) = 0.0
343 : END DO
344 : RETURN
345 : ENDIF
346 0 : sphbes(0) = sin(x)/x
347 0 : IF (l == 0) RETURN
348 0 : sphbes(1) = (sphbes(0) - cos(x))/x
349 : ! IF(l.le.1) RETURN
350 0 : s0 = sphbes(0)
351 0 : s1 = sphbes(1)
352 0 : lsta = lsta1(x, 200) !
353 0 : lmax = l !
354 0 : IF (lsta < l) THEN !
355 0 : lmax = lsta ! determine starting point lsta
356 0 : sphbes(lmax + 1:) = 0.0 ! for backward recurrence
357 : ELSE !
358 0 : lsta = lsta2(x, l, 15) !
359 : END IF !
360 0 : f0 = 0.0 !
361 0 : f1 = 1e-100 !
362 0 : DO ll = lsta, 0, -1 ! backward recurrence
363 0 : f = f1/x*(2*ll + 3) - f0; IF (ll <= lmax) sphbes(ll) = f ! with arbitrary start values
364 0 : f0 = f1 !
365 0 : f1 = f !
366 : END DO !
367 0 : IF (abs(s0) > abs(s1)) THEN; cs = s0/f !
368 0 : ELSE; cs = s1/f0 ! scale to correct values
369 : END IF !
370 0 : sphbes = cs*sphbes !
371 :
372 : CONTAINS
373 :
374 : ! Test starting point
375 : PURE FUNCTION lsta0(x, mp)
376 :
377 : IMPLICIT NONE
378 :
379 : INTEGER :: lsta0
380 : INTEGER, INTENT(IN) :: mp
381 : REAL , INTENT(IN) :: x
382 : REAL :: f, lgx
383 :
384 : lgx = log10(x)
385 : lsta0 = 0
386 : f = lgx
387 : DO WHILE (f > -mp)
388 : lsta0 = lsta0 + 1
389 : f = f + lgx - log10(2.0*lsta0 + 1)
390 : END DO
391 : END FUNCTION lsta0
392 :
393 : ! Returns starting point lsta1 for backward recurrence such that sphbes(lsta1) approx. 10^(-mp).
394 0 : PURE FUNCTION lsta1(x, mp)
395 :
396 : IMPLICIT NONE
397 :
398 : INTEGER :: lsta1
399 : INTEGER, INTENT(IN) :: mp
400 : REAL , INTENT(IN) :: x
401 : REAL :: f0, f1, f
402 : INTEGER :: n0, n1, nn, it
403 :
404 0 : n0 = int(1.1*x) + 1
405 0 : f0 = envj(n0, x) - mp
406 0 : n1 = n0 + 5
407 0 : f1 = envj(n1, x) - mp
408 0 : DO it = 1, 20
409 0 : nn = n1 - (n1 - n0)/(1.0 - f0/f1)
410 0 : f = envj(nn, x) - mp
411 0 : IF (abs(nn - n1) < 1) EXIT
412 0 : n0 = n1
413 0 : f0 = f1
414 0 : n1 = nn
415 0 : f1 = f
416 : END DO
417 0 : lsta1 = nn
418 0 : END FUNCTION lsta1
419 :
420 : ! Returns the starting point lsta2 for backward recurrence such that all sphbes(l) have mp significant digits.
421 0 : PURE FUNCTION lsta2(x, l, mp)
422 :
423 : IMPLICIT NONE
424 :
425 : INTEGER :: lsta2
426 : INTEGER, INTENT(IN) :: l, mp
427 : REAL , INTENT(IN) :: x
428 : REAL :: f0, f1, f, hmp, ejn, obj
429 : INTEGER :: n0, n1, nn, it
430 :
431 0 : hmp = 0.5*mp
432 0 : ejn = envj(l, x)
433 0 : IF (ejn <= hmp) THEN
434 0 : obj = mp
435 0 : n0 = int(1.1*x) + 1
436 : ELSE
437 0 : obj = hmp + ejn
438 0 : n0 = l
439 : END IF
440 0 : f0 = envj(n0, x) - obj
441 0 : n1 = n0 + 5
442 0 : f1 = envj(n1, x) - obj
443 0 : DO it = 1, 20
444 0 : nn = n1 - (n1 - n0)/(1.0 - f0/f1)
445 0 : f = envj(nn, x) - obj
446 0 : IF (abs(nn - n1) < 1) EXIT
447 0 : n0 = n1
448 0 : f0 = f1
449 0 : n1 = nn
450 0 : f1 = f
451 : END DO
452 0 : lsta2 = nn + 10
453 0 : END FUNCTION lsta2
454 :
455 0 : PURE FUNCTION envj(l, x)
456 :
457 : IMPLICIT NONE
458 :
459 : REAL :: envj
460 : REAL , INTENT(IN) :: x
461 : INTEGER, INTENT(IN) :: l
462 :
463 0 : envj = 0.5*log10(6.28*l) - l*log10(1.36*x/l)
464 :
465 0 : END FUNCTION envj
466 :
467 : END SUBROUTINE sphbessel
468 :
469 : ! Returns the spherical harmonics Y_lm(^rvec)
470 : ! for l = 0,...,ll in Y(1,...,(ll+1)**2).
471 :
472 0 : PURE SUBROUTINE harmonicsr(Y, rvec, ll)
473 :
474 : IMPLICIT NONE
475 :
476 : integer , intent(in) :: ll
477 : real , intent(in) :: rvec(3)
478 : complex, intent(out) :: Y((ll + 1)**2)
479 : complex :: c
480 : real :: stheta, ctheta, sphi, cphi, r, rvec1(3)
481 : integer :: l, m, lm
482 : complex, parameter :: img = (0.0, 1.0)
483 :
484 0 : Y(1) = 0.282094791773878
485 0 : IF (ll == 0) RETURN
486 :
487 0 : stheta = 0
488 0 : ctheta = 0
489 0 : sphi = 0
490 0 : cphi = 0
491 0 : r = norm2(rvec)
492 0 : IF (r > 1e-16) THEN
493 0 : rvec1 = rvec/r
494 0 : ctheta = rvec1(3)
495 0 : stheta = sqrt(rvec1(1)**2 + rvec1(2)**2)
496 0 : IF (stheta > 1e-16) THEN
497 0 : cphi = rvec1(1)/stheta
498 0 : sphi = rvec1(2)/stheta
499 : END IF
500 : ELSE
501 0 : Y(2:) = 0.0
502 : RETURN
503 : END IF
504 :
505 : ! define Y,l,-l and Y,l,l
506 0 : r = Y(1)
507 0 : c = 1
508 0 : DO l = 1, ll
509 0 : r = r*stheta*sqrt(1.0 + 1.0/(2*l))
510 0 : c = c*(cphi + img*sphi)
511 0 : Y(l**2 + 1) = r*conjg(c) ! l,-l
512 0 : Y((l + 1)**2) = r*c*(-1)**l ! l,l
513 : END DO
514 :
515 : ! define Y,l,-l+1 and Y,l,l-1
516 0 : Y(3) = 0.48860251190292*ctheta
517 0 : DO l = 2, ll
518 0 : r = sqrt(2.0*l + 1)*ctheta
519 0 : Y(l**2 + 2) = r*Y((l - 1)**2 + 1) ! l,-l+1
520 0 : Y(l*(l + 2)) = r*Y(l**2) ! l,l-1
521 : END DO
522 :
523 : ! define Y,l,m, |m|<l-1
524 0 : DO l = 2, ll
525 0 : lm = l**2 + 2
526 0 : DO m = -l + 2, l - 2
527 0 : lm = lm + 1
528 : Y(lm) = sqrt((2.0*l + 1)/(l + m)/(l - m))*( &
529 : sqrt(2.0*l - 1)*ctheta*Y(lm - 2*l) - &
530 0 : sqrt((l + m - 1.0)*(l - m - 1)/(2*l - 3))*Y(lm - 4*l + 2))
531 : END DO
532 : END DO
533 :
534 : END SUBROUTINE harmonicsr
535 :
536 : ! Returns the complex error function.
537 6 : FUNCTION cerf(z)
538 :
539 : USE m_constants, ONLY: pimach
540 :
541 : IMPLICIT NONE
542 : COMPLEX, INTENT(IN) :: z
543 : COMPLEX :: cerf
544 : COMPLEX :: z1, z2, c, d, delta
545 : REAL :: pi
546 : INTEGER :: i
547 :
548 6 : pi = pimach()
549 6 : z1 = z
550 6 : IF (real(z) < 0) z1 = -z1
551 6 : IF (real(z1) < 2.0) THEN ! McLaurin series
552 0 : z2 = z1**2
553 0 : i = 0
554 0 : c = z1
555 0 : cerf = z1
556 : DO
557 0 : i = i + 1
558 0 : c = -c*z2/i
559 0 : cerf = cerf + c/(2*i + 1)
560 0 : IF (abs(c/(2*i + 1)) < 1e-20) EXIT
561 : END DO
562 0 : cerf = cerf*2/sqrt(pi)
563 : ELSE ! continued fraction using Lentz's method
564 : d = 0.0
565 : c = z1
566 : cerf = z1
567 : i = 0
568 : DO
569 126 : i = i + 1
570 126 : c = 2*z1 + i/c
571 126 : d = (2*z1 + i*d)**(-1)
572 126 : delta = c*d
573 126 : cerf = cerf*delta
574 126 : IF (abs(1 - delta) < 1e-15) EXIT
575 126 : i = i + 1
576 126 : c = z1 + i/c
577 126 : d = (z1 + i*d)**(-1)
578 126 : delta = c*d
579 126 : cerf = cerf*delta
580 126 : IF (abs(1 - delta) < 1e-15) EXIT
581 120 : IF (i == 10000) &
582 6 : call judft_error("cerf: Lentz method not converged after 10000 steps.")
583 : END DO
584 6 : cerf = 1 - exp(-z1**2)/cerf/sqrt(pi)
585 : END IF
586 6 : IF (real(z) < 0) cerf = -cerf
587 6 : END FUNCTION cerf
588 :
589 0 : FUNCTION chr(int)
590 :
591 : IMPLICIT NONE
592 :
593 : CHARACTER(5) :: chr
594 : INTEGER :: int
595 :
596 0 : WRITE (chr, '(I5)') int
597 0 : END FUNCTION chr
598 :
599 : END MODULE m_util
|
