Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! Copyright (c) 2016 Anton Shterenlikht, The University of Bristol, UK
9 : !
10 : ! Redistribution and use in source and binary forms, with or without
11 : ! modification, are permitted provided that the following conditions are
12 : ! met:
13 : !
14 : ! 1. Redistributions of source code must retain the above copyright
15 : ! notice, this list of conditions and the following disclaimer.
16 : !
17 : ! 2. Redistributions in binary form must reproduce the above copyright
18 : ! notice, this list of conditions and the following disclaimer in the
19 : ! documentation and/or other materials provided with the distribution.
20 : !
21 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
22 : ! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
23 : ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24 : ! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 : ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 : ! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28 : ! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29 : ! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30 : ! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 : ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 : !
33 : ! Module with error function and related functions.
34 : ! Two algorithms are implemented: Poppe and Wijers (faddeyeva_fast, erfz_fast)
35 : ! and Zaghloul and Ali (faddeyeva_accurate, erfz_accurate). The second algorithm is
36 : ! supposed to be more accurate for some values. However, the first
37 : ! algorithm can be over an order of magnitude faster.
38 : !
39 : ! This code was written before the author became aware of
40 : ! M. R. Zaghloul, Remark on "Algorithm 916: Computing the
41 : ! Faddeyeva and Voight functions": efficiency improvements and
42 : ! Fortran translation, ACM Trans. Math. Software 42, Article 26, 2016.
43 :
44 : ! **************************************************************************************************
45 : !> \brief Module to compute the error function of a complex argument
46 : !> \par History
47 : !> 08.2025 Adapted to use CP2K intrinsics and constants
48 : !> \author Stefano Battaglia
49 : ! **************************************************************************************************
50 : MODULE erf_complex
51 :
52 : USE kinds, ONLY: dp
53 : USE mathconstants, ONLY: half,&
54 : one,&
55 : oorootpi,&
56 : pi,&
57 : two,&
58 : zero
59 :
60 : IMPLICIT NONE
61 :
62 : REAL(kind=dp), PARAMETER :: rmin = TINY(one), eps0 = EPSILON(one), &
63 : sqrt_log_rmin = SQRT(-LOG(rmin)), &
64 : pi2 = pi*pi
65 : COMPLEX(kind=dp), PARAMETER :: &
66 : cmplxj = CMPLX(zero, one, kind=dp), &
67 : cmplx0 = CMPLX(zero, zero, kind=dp)
68 :
69 : PRIVATE
70 : PUBLIC :: faddeyeva_fast, faddeyeva_accurate, erfz_fast, erfz_accurate
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Computes the Faddeyeva function w(z) = exp(-z**2) * erfc(-i*z)
76 : !> \param z complex argument
77 : !> \param err desired accuracy (positive)
78 : !> \return Faddeyeva function w(z)
79 : ! **************************************************************************************************
80 0 : elemental COMPLEX(kind=dp) FUNCTION faddeyeva_accurate(z, err)
81 :
82 : ! This is the Faddeyeva or the plasma dispersion function,
83 : ! w(z) = exp(-z**2) * erfc(-i*z). erfc(z) is the complex complementary
84 : ! error function of z. z is a complex number.
85 : !
86 : ! Adapted from the Matlab code implementing TOMS
87 : ! Algorithm 916: http://www.netlib.org/toms/
88 : !
89 : ! file: 916.zip
90 : ! ref: TOMS 38,2 (Dec 2011) Article: 15
91 : ! for: Computing the Faddeyeva and Voigt Functions
92 : ! by: Mofreh R. Zaghloul and Ahmed N. Ali
93 : !
94 : ! Most of the code is calculation of equations (13)-(19) of the
95 : ! above paper.
96 : !
97 : ! Inputs:
98 : ! z - the argument of the function
99 : ! err - The desired accuracy (positive). err must be (err .le. 1.0e-4)
100 : ! and (err .ge. 0.06447*epsilon). For efficiency,
101 : ! no checks are made!
102 : ! The lowest accuracy and the fastest calculation are obtained
103 : ! with err .eq. 1.0e-4. For higher accuracy use smaller err.
104 :
105 : COMPLEX(kind=dp), INTENT(in) :: z
106 : REAL(kind=dp), INTENT(in) :: err
107 :
108 : INTEGER :: n, n3, n3_3
109 : REAL(kind=dp) :: a, a_pi, a_sqr, aux13, cos_2yx, del2_tmp, del3, del3_3_tmp, del3_tmp, del5, &
110 : delta3, delta5, den1, erfcsy, exp1, exp2, exp3, exp3_3_den, exp3_den, exp_del1, &
111 : exp_x_sqr, four_a_sqr, half_a, l_old, myerr, sigma1, sigma2, sigma3, sigma4, sigma4_5, &
112 : sigma5, sin_2yx, two_a, two_a_pi, two_a_sqr, two_a_x, two_exp_x_sqr_ysqr, two_yx, v_old, &
113 : x, x_sqr, xsign, y, y_sqr, ysign
114 :
115 0 : x = REAL(z)
116 0 : y = AIMAG(z)
117 :
118 : ! For purely imaginary z, use intrinsic scaled complement of
119 : ! the error function, erfc_scaled (F2008 and beyond).
120 : ! Return immediately.
121 0 : IF (ABS(x) == zero) THEN
122 0 : faddeyeva_accurate = erfc_scaled(y)
123 0 : RETURN
124 : END IF
125 :
126 0 : myerr = MAX(err, eps0)
127 0 : a = SQRT(-pi2/LOG(err/2.0_dp))
128 0 : half_a = half*a
129 0 : a_sqr = a**2
130 0 : two_a = 2*a
131 0 : two_a_sqr = 2*a_sqr
132 0 : four_a_sqr = 4*a_sqr
133 0 : a_pi = a/pi
134 0 : two_a_pi = 2*a_pi
135 0 : erfcsy = erfc_scaled(ABS(y))
136 0 : xsign = SIGN(one, x)
137 0 : ysign = SIGN(one, y)
138 0 : x = ABS(x)
139 0 : y = MAX(rmin, ABS(y))
140 0 : x_sqr = x**2
141 0 : y_sqr = y**2
142 0 : two_yx = 2*y*x
143 0 : two_a_x = two_a*x
144 0 : exp_x_sqr = EXP(-x_sqr)
145 0 : cos_2yx = COS(two_yx)
146 0 : sin_2yx = SIN(two_yx)
147 : v_old = exp_x_sqr* &
148 0 : (erfcsy*cos_2yx + two_a_pi*SIN(two_yx/2)**2/y)
149 0 : l_old = -erfcsy + a_pi/y
150 0 : sigma3 = rmin
151 0 : sigma5 = rmin
152 0 : sigma4_5 = zero
153 0 : delta3 = one
154 0 : delta5 = one
155 0 : n = 0
156 0 : n3 = CEILING(x/a)
157 0 : n3_3 = n3 - 1
158 :
159 0 : outer: IF ((sqrt_log_rmin - x) > 0) THEN
160 0 : sigma1 = rmin
161 0 : sigma2 = rmin
162 0 : sigma4 = rmin
163 0 : exp1 = EXP(-two_a_x)
164 0 : exp2 = EXP(four_a_sqr*n3 - 2*two_a_x - two_a_sqr)
165 0 : exp3 = EXP(-(two_a_sqr*n3 - two_a_x - two_a_sqr))
166 0 : del2_tmp = one
167 : del3_tmp = EXP(-(a_sqr*n3**2 - two_a_x*n3 &
168 0 : - two_a_sqr*n3 + x_sqr + two_a_x + a_sqr))
169 0 : del3_3_tmp = EXP(a_sqr - (two_a_sqr*n3 - two_a_x))
170 :
171 0 : loop1: DO
172 0 : IF (delta3 < myerr .AND. delta5 < myerr .AND. &
173 : n > 50) EXIT loop1
174 0 : n = n + 1
175 0 : den1 = a_sqr*n**2 + y_sqr
176 0 : exp_del1 = EXP(-(a_sqr*n**2))/den1
177 0 : del2_tmp = del2_tmp*exp1
178 :
179 0 : minor: IF (n3_3 >= 1) THEN
180 0 : del3_tmp = del3_tmp*exp3
181 : exp3_den = del3_tmp*exp_del1* &
182 0 : (den1/(a_sqr*n3**2 + y_sqr))
183 0 : del3_3_tmp = del3_3_tmp*exp2
184 : exp3_3_den = exp3_den*del3_3_tmp* &
185 : ((a_sqr*n3**2 + y_sqr)/ &
186 0 : (a_sqr*n3_3**2 + y_sqr))
187 0 : del5 = n3_3*exp3_3_den + n3*exp3_den
188 0 : del3 = exp3_3_den + exp3_den
189 : ELSE
190 0 : del3_tmp = del3_tmp*exp3
191 : del3 = del3_tmp*exp_del1* &
192 0 : (den1/(a_sqr*n3**2 + y_sqr))
193 0 : del5 = n3*del3
194 : END IF minor
195 :
196 0 : delta3 = del3/sigma3
197 0 : delta5 = del5/sigma5
198 0 : sigma1 = sigma1 + exp_del1
199 0 : sigma2 = sigma2 + del2_tmp*exp_x_sqr*exp_del1
200 0 : sigma3 = sigma3 + del3
201 0 : sigma4 = sigma4 + n*del2_tmp*exp_x_sqr*exp_del1
202 0 : sigma5 = sigma5 + del5
203 :
204 0 : IF (x >= 5.0e-4_dp) THEN
205 0 : sigma4_5 = -sigma4 + sigma5
206 : ELSE
207 : sigma4_5 = sigma4_5 + 2*n**2*two_a_x*exp_x_sqr &
208 : *exp_del1*(one + 1.666666666666667e-1_dp &
209 : *(two_a_x*n)**2 + 8.333333333333333e-3_dp &
210 0 : *(two_a_x*n)**4)
211 : END IF
212 :
213 0 : n3 = n3 + 1
214 0 : n3_3 = n3_3 - 1
215 :
216 : END DO loop1
217 :
218 : ! Second line of Eqn (13)
219 : aux13 = y*two_a_pi* &
220 0 : (-cos_2yx*exp_x_sqr*sigma1 + half*(sigma2 + sigma3))
221 0 : mumu: IF (y <= 5.0_dp .AND. two_yx > rmin) THEN
222 : faddeyeva_accurate = v_old + aux13 + cmplxj*xsign &
223 : *(sin_2yx*exp_x_sqr*(l_old + two_a_pi*y*sigma1) &
224 0 : + two_a_pi*half_a*sigma4_5)
225 0 : ELSE IF (y <= 5.0_dp .AND. two_yx <= rmin) THEN
226 : faddeyeva_accurate = v_old + aux13 + cmplxj*xsign &
227 : *(two*y*exp_x_sqr*(x*l_old + x*two_a_pi*y &
228 0 : *sigma1) + two_a_pi*half_a*sigma4_5)
229 : ELSE
230 : faddeyeva_accurate = v_old + aux13 + cmplxj*xsign &
231 : *(sin_2yx*exp_x_sqr*MIN(zero, ABS(l_old &
232 : + (two_a_pi*y*sigma1))) + two_a_pi*half_a &
233 0 : *sigma4_5)
234 : END IF mumu
235 :
236 0 : ELSE IF (x >= sqrt_log_rmin .AND. x < 1.0e15_dp) THEN
237 :
238 0 : exp2 = EXP(four_a_sqr*n3 - 2*two_a_x - two_a_sqr)
239 0 : del3_3_tmp = EXP(a_sqr + two_a_x - two_a_sqr*n3)
240 :
241 0 : loop2: DO
242 0 : IF (delta3 < myerr .AND. delta5 < myerr .AND. &
243 : n > 50) EXIT loop2
244 0 : n = n + 1
245 0 : IF (n3_3 >= 1) THEN
246 : exp3_den = EXP(-(a*n3 - x)*(a*n3 - x)) &
247 0 : /(a_sqr*n3**2 + y_sqr)
248 0 : del3_3_tmp = del3_3_tmp*exp2
249 : exp3_3_den = exp3_den*del3_3_tmp*((a_sqr*n3**2 + y_sqr) &
250 0 : /(a_sqr*n3_3**2 + y_sqr))
251 0 : del5 = n3_3*exp3_3_den + n3*exp3_den
252 0 : del3 = exp3_3_den + exp3_den
253 : ELSE
254 0 : del3 = EXP(-(a*n3 - x)**2)/(a_sqr*n3**2 + y_sqr)
255 0 : del5 = n3*del3
256 : END IF
257 :
258 0 : delta3 = del3/sigma3
259 0 : delta5 = del5/sigma5
260 0 : sigma3 = sigma3 + del3
261 0 : sigma5 = sigma5 + del5
262 0 : n3 = n3 + 1
263 0 : n3_3 = n3_3 - 1
264 :
265 : END DO loop2
266 :
267 : faddeyeva_accurate = v_old + y*a_pi*sigma3 + cmplxj*xsign*(sin_2yx &
268 0 : *exp_x_sqr*l_old + two_a_pi*half_a*sigma5)
269 :
270 : ELSE
271 0 : faddeyeva_accurate = oorootpi*((y + cmplxj*xsign*x)/(x_sqr + y_sqr))
272 : END IF outer
273 :
274 0 : IF (ysign < zero) THEN
275 0 : two_exp_x_sqr_ysqr = two*EXP(-x_sqr + y_sqr)
276 : faddeyeva_accurate = two_exp_x_sqr_ysqr*cos_2yx - REAL(faddeyeva_accurate) - cmplxj &
277 0 : *(-xsign*two_exp_x_sqr_ysqr*sin_2yx - AIMAG(faddeyeva_accurate))
278 : END IF
279 :
280 : END FUNCTION faddeyeva_accurate
281 :
282 : ! **************************************************************************************************
283 : !> \brief Computes the error function of a complex argument using the Zaghloul and Ali algorithm
284 : !> \param z complex argument
285 : !> \param err desired accuracy (positive)
286 : !> \return error function of z
287 : ! **************************************************************************************************
288 0 : elemental COMPLEX(kind=dp) FUNCTION erfz_accurate(z, err)
289 :
290 : ! This is an error function of a complex argument, which uses faddeyeva_accurate(z).
291 :
292 : COMPLEX(kind=dp), INTENT(in) :: z
293 : REAL(kind=dp), INTENT(in) :: err
294 :
295 0 : erfz_accurate = one - faddeyeva_accurate(cmplxj*z, err)*EXP(-z**2)
296 :
297 0 : END FUNCTION erfz_accurate
298 :
299 : ! **************************************************************************************************
300 : !> \brief Computes the Faddeyeva function w(z) = exp(-z**2) * erfc(-i*z)
301 : !> \param z complex argument
302 : !> \return Faddeyeva function w(z)
303 : ! **************************************************************************************************
304 842260 : elemental COMPLEX(kind=dp) FUNCTION faddeyeva_fast(z)
305 :
306 : ! A modified version of algorithm 680, rewritten in Fortran 2008.
307 : ! G.P.M. Poppe, C.M.J. Wijers, More efficient computation of
308 : ! the complex error-function, ACM Trans. Math. Software 16:38-46, 1990.
309 : ! and
310 : ! G.P.M. Poppe, C.M.J. Wijers, Algorithm 680, Evaluation of the
311 : ! complex error function, ACM Trans. Math. Software 16:47, 1990.
312 : !
313 : ! Given a complex number z, this function computes
314 : ! the value of the Faddeeva-function w(z) = exp(-z**2)*erfc(-i*z),
315 : ! where erfc is the complex complementary error-function and i
316 : ! means sqrt(-1). The accuracy of the algorithm for z in the 1st
317 : ! and 2nd quadrant is 14 significant digits; in the 3rd and 4th
318 : ! it is 13 significant digits outside a circular region with radius
319 : ! 0.126 around a zero of the function.
320 :
321 : COMPLEX(kind=dp), INTENT(in) :: z
322 :
323 : REAL(kind=dp), PARAMETER :: factor = 1.12837916709551257388_dp
324 :
325 : INTEGER :: i, j, kapn, n, np1, nu
326 : LOGICAL :: a, b
327 : REAL(kind=dp) :: c, daux, h, h2, qlambda, qrho, rx, ry, &
328 : sx, sy, tx, ty, u, u1, u2, v, v1, v2, &
329 : w1, x, xabs, xabsq, xaux, xi, xquad, &
330 : xsum, y, yabs, yi, yquad, ysum
331 :
332 : ! factor is 2/sqrt(pi)
333 :
334 : ! To avoid the complier uninitialised varning
335 842260 : h2 = zero
336 :
337 842260 : xi = REAL(z)
338 842260 : yi = AIMAG(z)
339 842260 : xabs = ABS(xi)
340 842260 : yabs = ABS(yi)
341 842260 : x = xabs/6.3_dp
342 842260 : y = yabs/4.4_dp
343 842260 : qrho = x**2 + y**2
344 842260 : xabsq = xabs**2
345 842260 : xquad = xabsq - yabs**2
346 842260 : yquad = 2*xabs*yabs
347 :
348 842260 : a = qrho < 0.085264_dp
349 :
350 842260 : branch1: IF (a) THEN
351 :
352 : ! If ( qrho .lt. 0.085264 ) then the Faddeeva-function is evaluated
353 : ! using a power-series (abramowitz/stegun, equation (7.1.5), p.297)
354 : ! n is the minimum number of terms needed to obtain the required
355 : ! accuracy
356 :
357 0 : qrho = (one - 0.85_dp*y)*SQRT(qrho)
358 0 : n = NINT(6.0_dp + 72.0_dp*qrho)
359 0 : j = 2*n + 1
360 0 : xsum = one/REAL(j, kind=dp)
361 0 : ysum = zero
362 :
363 0 : DO i = n, 1, -1
364 0 : j = j - 2
365 0 : xaux = (xsum*xquad - ysum*yquad)/REAL(i, kind=dp)
366 0 : ysum = (xsum*yquad + ysum*xquad)/REAL(i, kind=dp)
367 0 : xsum = xaux + one/REAL(j, kind=dp)
368 : END DO
369 :
370 0 : u1 = -factor*(xsum*yabs + ysum*xabs) + one
371 0 : v1 = factor*(xsum*xabs - ysum*yabs)
372 0 : daux = EXP(-xquad)
373 0 : u2 = daux*COS(yquad)
374 0 : v2 = -daux*SIN(yquad)
375 0 : u = u1*u2 - v1*v2
376 0 : v = u1*v2 + v1*u2
377 : ELSE
378 :
379 842260 : bran2: IF (qrho > one) THEN
380 :
381 : ! If ( qrho .gt. 1) then w(z) is evaluated using the laplace
382 : ! continued fraction. nu is the minimum number of terms needed
383 : ! to obtain the required accuracy.
384 :
385 841340 : h = zero
386 841340 : kapn = 0
387 841340 : qrho = SQRT(qrho)
388 : nu = INT(3.0_dp + (1442.0_dp/(26.0_dp*qrho &
389 841340 : + 77.0_dp)))
390 :
391 : ELSE
392 :
393 : ! If ( qrho .ge. 0.085264 .and. qrho .le. one ) then
394 : ! w(z) is evaluated by a truncated Taylor expansion,
395 : ! where the Laplace continued fraction is used to calculate
396 : ! the derivatives of w(z). KAPN is the minimum number of terms
397 : ! in the Taylor expansion needed to obtain the required accuracy.
398 : ! NU is the minimum number of terms of the continued fraction
399 : ! needed to calculate the derivatives with the required accuracy.
400 :
401 920 : qrho = (one - y)*SQRT(one - qrho)
402 920 : h = 1.88_dp*qrho
403 920 : h2 = two*h
404 920 : kapn = NINT(7.0_dp + 34.0_dp*qrho)
405 920 : nu = NINT(16.0_dp + 26.0_dp*qrho)
406 :
407 : END IF bran2
408 :
409 842260 : b = h > zero
410 :
411 : ! To avoid uninitialise compiler warning. qlambda is used
412 : ! only if (b), so can define to any value otherwise.
413 841340 : qlambda = zero
414 920 : IF (b) qlambda = h2**kapn
415 :
416 842260 : rx = zero
417 842260 : ry = zero
418 842260 : sx = zero
419 842260 : sy = zero
420 :
421 13006912 : DO n = nu, 0, -1
422 12164652 : np1 = n + 1
423 12164652 : tx = yabs + h + np1*rx
424 12164652 : ty = xabs - np1*ry
425 12164652 : c = half/(tx**2 + ty**2)
426 12164652 : rx = c*tx
427 12164652 : ry = c*ty
428 13006912 : IF (b .AND. n <= kapn) THEN
429 7360 : tx = qlambda + sx
430 7360 : sx = rx*tx - ry*sy
431 7360 : sy = ry*tx + rx*sy
432 7360 : qlambda = qlambda/h2
433 : END IF
434 : END DO
435 :
436 842260 : IF (h == zero) THEN
437 841340 : u = factor*rx
438 841340 : v = factor*ry
439 : ELSE
440 920 : u = factor*sx
441 920 : v = factor*sy
442 : END IF
443 :
444 842260 : IF (yabs == zero) u = EXP(-xabs**2)
445 :
446 : END IF branch1
447 :
448 : ! Evaluation of w(z) in the other quadrants
449 :
450 842260 : IF (yi < zero) THEN
451 :
452 0 : IF (a) THEN
453 0 : u2 = two*u2
454 0 : v2 = two*v2
455 : ELSE
456 0 : xquad = -xquad
457 0 : w1 = two*EXP(xquad)
458 0 : u2 = w1*COS(yquad)
459 0 : v2 = -w1*SIN(yquad)
460 : END IF
461 :
462 0 : u = u2 - u
463 0 : v = v2 - v
464 0 : IF (xi > zero) v = -v
465 : ELSE
466 842260 : IF (xi < zero) v = -v
467 : END IF
468 :
469 842260 : faddeyeva_fast = CMPLX(u, v, kind=dp)
470 :
471 842260 : END FUNCTION faddeyeva_fast
472 :
473 : ! **************************************************************************************************
474 : !> \brief Computes the error function of a complex argument using the Poppe and Wijers algorithm
475 : !> \param z complex argument
476 : !> \return error function of z
477 : ! **************************************************************************************************
478 842260 : elemental COMPLEX(kind=dp) FUNCTION erfz_fast(z)
479 :
480 : ! This is an error function of a complex argument, which uses faddeyeva_fast(z).
481 :
482 : COMPLEX(kind=dp), INTENT(in) :: z
483 :
484 842260 : erfz_fast = one - faddeyeva_fast(cmplxj*z)*EXP(-z**2)
485 :
486 842260 : END FUNCTION erfz_fast
487 :
488 : !*********************************************************************
489 : END MODULE erf_complex
|
