- Email: [email protected]
Calculation of the generalized leaky aquifer integral
Computer Physics Communications 173 (2005) 1–8 www.elsevier.com/locate/cpc
Calculation of the generalized leaky aquifer integral J.A. Alford Quantum Theory Project, Department of Physics, University of Florida, Gainesville, FL 32611, USA Received 13 April 2005; accepted 30 June 2005 Available online 8 August 2005
Abstract We provide a convenient method for the evaluation of the generalized leaky aquifer integral Wβ (x, y), for all real x, y > 0 and integers and half odd integers β 12 . The key ingredient to obtaining values for higher β is a new asymptotic expansion of the integral that is applicable for a wide range of these parameters, particularly large β. An asymptotic value and a value, calculated via previously published methods, for the same x and y but small β, form a well defined two-point boundary value problem that is easily solved. 2005 Elsevier B.V. All rights reserved. PACS: 71.15.-m Keywords: Leaky aquifer; One-dimensional systems; Fourier space methods; Asymptotic expansions; Ewald summation
1. Introduction The practical relevance of the leaky aquifer integral, defined as ∞ W0 (x, y) =
e−xt−y/t dt t
(1)
1
seems to have arisen first in the study of radial fluid flow from artesian aquifers [1], whence the origin of its name. A more general form of the integral ∞ Wβ (x, y) =
e−xt−y/t dt t β+1
1
E-mail address: [email protected] (J.A. Alford). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.06.012
(2)
2
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
for integer and half odd integer β has recently appeared in the Fourier space formulation of electronic structure methods used to calculate properties of extended periodic systems [2]. Before moving on to the technical aspects of how to compute this integral efficiently, we briefly review the context which has motivated this necessity and thus the current work. Common to many electronic structure methods is the need to evaluate lattice sums of Coulomb integrals involving local basis functions, for instance f (r − A)f (r − B − R) V= (3) dr dr , |r − r + R| R
where the basis functions f are centered about the positions A and B + R, where R is a direct space lattice vector. Such integrals are required for example for the calculation of the Hartree term of standard Hamiltonians and for variational Coulomb fitting of the electronic density to basis functions [3,4]. Although various approaches have been introduced to truncate the lattice sum in direct space properly, maintaining precise control over the accuracy has remained a difficult problem [2]. For example, the popular multipolar expansion used for obtaining the long range contributions to the sum converges prohibitively slowly for metallic or near-metallic systems. An alternative to direct space methods is the Ewald summation technique, a mixed space approach in which the short range contributions are evaluated in real space and the long range part is converted into a reciprocal space sum that is also relatively short range and can be computed efficiently [4,5]. To carry out the latter sum, one needs the Fourier transform of the Coulomb integral, the form of which will depend on the specific selection of basis orbitals f . The Hermite Gaussian orbitals provide a rather convenient choice, as the higher angular momenta orbitals of that set can be simply generated via derivatives with respect to an orbital center f (i, a, A, r) = (∂Ax )ix (∂Ay )iy (∂Az )iz e−a|r−A| . 2
(4)
Using these orbitals, it may be shown [4,5] that the Fourier space analog of Eq. (3) is given by ij (G), Vij = V G
ij (G) = (−1) (∂r0x ) (∂r0y ) (∂r0z ) e V |j|
nx
ny
nz iG·r0
e
−iG·r
RN
1 dr
e−γ |r+r0⊥|
2 u2
du,
(5)
(6)
0
where N is the number of dimensions for which the system is periodic, x, y, and z as subscripts refer to the usual Cartesian directions, r0 = A − B, G is a reciprocal lattice vector, and γ = a + b. i and j are multi-indices, the Cartesian components of which are nonnegative integers, n = i + j, |j| = jx + jy + jz , r0 is the component of r0 in the space spanned by RN , and r0⊥ is the component orthogonal to this space. (Note that, irrespective of N , the √ system is three-dimensional.) With the substitution u = 1/ t , the dependence on N becomes more clear ij (G) = (−1)|j| (∂r0x )nx (∂r0y )ny (∂r0z )nz e−iG·r0 V
∞
e
−|G|2 /(4γ )t−γ r02 /t ⊥
2t (3−N )/2
dt.
(7)
1
Note that the derivatives with respect to components of r0 generate higher powers of t in the denominator of the integrand. These integrals have the general form of Eq. (2) if we set x = |G|2 /4γ , y = γ r02⊥ , and β + 1 = (3 − N )/2 + k, where k is an integer resulting from differentiation. They also obey a recursion relation first discussed by Thomas [6] xWβ−1 (x, y) + βWβ (x, y) − yWβ+1 (x, y) = e−(x+y) that is easily derived through integration by parts.
(8)
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
3
For a 3D crystal (N = 3), there is no finite orthogonal vector r0⊥ (y = 0) and Eq. (7) simply reduces to a sum 2 over products of components of G and the factor e−|G| /4γ |G|2 . If N = 2, the indefinite form of the integral can be brought into a known form [7] through a substitution like t = v 2 and evaluated directly. The result for n = 0, corresponding to β + 1 = 12 , is √ √ √ √ √ 00 (G) = √π e2 xy erfc( x + √y ) + e−2 xy erfc( x − √y ) . V 4 x
(9)
To obtain the higher integrals in this particular case, β = n/2, where n is an odd integer, two coupled, twoterm recursion relations may be derived, through derivatives with respect to y, that are generally more stable than Eq. (8). The N = 1 case also requires integrals of the form of Eq. (2) where β is a nonnegative integer. Eq. (8) is the only known applicable recursion for N = 1 and is the primary focus of this paper. The analysis presented here also holds for β equal to half odd integers and thus may be applied when N = 2 as well. The difficulty of evaluating the higher integrals Wβ (x, y) for integer β has in fact prevented the general application of the Ewald summation technique to systems with 1D periodicity, with the notable exception of the studies by Flamant et al. [8] and Langridge et al. [9]. Flamant’s treatment was limited to the use of s orbital basis functions only, in effect restricting β to zero. Langridge and coworkers encountered the leaky aquifer integral in their application of the Ewald technique to the Madelung sums connected with the multipolar expansion. In their expressions the parameters that enter the integral (here x and y) are scaled by the Ewald split parameter µ, and the following expansion is used for the calculation Wβ (x/µ, µy) =
∞ (−µy)m Eβ+m+1 (x/µ). m!
(10)
m=0
The Ewald parameter is adjustable; in practice it is modified in some cases in order to retain sufficient accuracy in this expression [10]. The parameter must play a dual role in this respect and thus the generality of the method is somewhat unclear, although it may represent a significant improvement to more common associated techniques. Harris, however, pointed out that their procedure is not optimal in general [11], a point not disputed by Langridge et al. [10]. The bulk of other previous work related to the integral has focused on deriving efficient expansions of W0 (x, y) for various ranges of x and y [12]. Kryachko [13] has shown that the zeroth-order integral may be written as an infinite sum over Legendre functions of the second kind, and Harris [14] has extended this result to similar forms for W1 (x, y) and W2 (x, y). The obvious method of reaching higher derivatives, namely upward recursion using Eq. (8), is highly unstable except in the case that y is much larger than x and β. It is our intention that the present work will remove the barrier to the general application of the Ewald technique to 1D periodic systems. To our knowledge, no asymptotic form for Wβ has been published, preventing any attempt at or analysis of downward recursion. It is somewhat obvious, however, that the downward recursion will be unstable except in the case that x is much larger than y and β. Rather than working with such uncontrolled approaches to obtain the higher integrals, we view their evaluation as the process of solving a two-point boundary value problem, where the end points are given by Wβ0 with β0 = − 12 , 0, 12 , or 1, and Wβ0 +m , m being a large integer. With the end-point values supplied for a given x and y, the problem of determining the sequence Wβ0 , Wβ0 +1 , Wβ0 +2 , . . . , Wβ0 +m−1 , Wβ0 +m can be solved using standard algebraic methods. In the next section, we derive a general asymptotic expansion for large β and describe its range of applicability. The method of solution for a range of integer or half odd integer β along with some numerical results and discussion is given in Section 3.
4
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
2. Asymptotic formula We follow the procedure used by Gautschi [15] for the exponential integrals En (x) to derive a new expression for Wβ (x, y). Consider the integral t2 I = j
f (t) = −xt − y/t − β ln t,
ef (t) dt,
(11)
t1
where the value of j denotes a specific interval over t , 0, t1 = 1, t2 = ∞, j = 1, t1 = 1, t2 = t ∗ − ε, 2, t1 = t ∗ + ε, t2 = ∞
(12)
and I 0 = lim (I 1 + I 2 ).
(13)
ε→0+
Let t ∗ be defined by f (t ∗ ) = 0, 2 β y β ∗ + . t =− + 2x 2x x
(14)
We deal with the case for which 1 < t ∗ < ∞. Except for t ∗ = 1, the following is obviously valid for t ∗ not lying in this interval, for example, when x > y. With the definitions g0 (t) = 1/f (t), integration of
Ij
gk+1 (t) = gk (t)g0 (t),
(15)
by parts yields the expansion
j
j
j
j
j
j
I j = v0 − v1 + v2 − v3 + · · · + (−1)k−1 vk−1 + Rk ,
(16)
where vk = gk (t2 )ef (t2 ) − gk (t1 )ef (t1 )
(17)
and j Rk
t2 = (−1)
k
gk−1 (t)ef (t) dt.
(18)
t1
All gk (t) in addition to g0 (t) are singular at t = t ∗ . Adding I 1 and I 2 , we see that all of the singular terms cancel term by term,
∗ ∗ lim vk1 + vk2 = lim gk (t ∗ + ε)ef (t +ε) − gk (t ∗ − ε)ef (t −ε) + gk (∞)ef (∞) − gk (1)ef (1) ε→0+
ε→0+
= gk (∞)ef (∞) − gk (1)ef (1) ≡ vk0 .
(19)
Thus I 0 may be expressed as I0 =
n−1 (−1)k vk0 + Rn0 . k=0
(20)
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
5
Before analyzing the remainder term, we construct a recursion equation for generating the gk . Calculation of the first few gk using Eq. (15) shows that each takes the form gk (u, α) =
(−1)k hk (u, α) g0 (u, α), (1 + u − αu )2k β k
(21)
where xy xt , α= 2 (22) β β and the hk (u, α) are polynomials in u and α. The recurring pattern that determines the form of the polynomials is seen to be
hk+1 (u, α) = (1 + u − α/u) · ∂u uhk (u, α) − (2k + 1)(1 + α/u2 )uhk (u, α) (23) u=
which gives the recursion
hk+1 (u, α) = u(u + 1) − α
hk (u, α) +
1+k (1 − 2ku) − 2α hk (u, α), u
(24)
where the prime now denotes differentiation with respect to u rather than t . Setting α = 0, equivalent to putting y = 0, Eqs. (21) and (24) manifestly reduce to the forms originally derived by Gautschi [15] for the exponential integral as they should. Evaluation of the gk (u, α) at the proper limits to obtain the vk0 in Eq. (19) and the further substitutions u = x/β,
v = y/β
(25)
finally yield Wβ (x, y) =
n−1 hk (u, v) e−(x+y) + Rn0 . x −y +β (u − v + 1)2k β k
(26)
k=0
The asymptotic prefactor to the expansion has been left in terms of x and y to show that it is in fact the asymptotic value expected from the recursion in Eq. (8), at least for x > y. The first few polynomials hk (u, v) are listed in Table 1. The parameter space over which this expansion converges is at least all u, v > 0 and β > 1 such that x −y +β >1
(27)
as can be seen by rewriting the sum as n−1 k=0
hk · β k . (x − y + β)2k
Table 1 Polynomials hk (u, v) for the asymptotic expansion of Wβ (x, y) generated by Eq. (24) k
hk (u, v)
0 1 2 3 4 5
1 1 − 2v 1 − 4v − 2u + 6uv + 6v 2 1 − 6v − 8u + 40uv + 14v 2 + 6u2 − 24u2 v − 72uv 2 − 24v 3 1 − 8v + 154uv + 28v 2 − 348u2 v − 432uv 2 − 36v 3 + 120u3 v + 720u2 v 2 + 720uv 3 − 22u + 58u2 − 24u3 + 120v 4 1 − 10v + 468uv + 44v 2 − 2624u2 v − 1820uv 2 − 152v 3 − 52u + 328u2 + 3108u3 v + 8436u2 v 2 + 3564uv 3 − 108v 4 − 720u4 v − 7200u3 v 2 − 14400u2 v 3 − 7200uv 4 − 444u3 + 120u4 − 720v 5
6
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
Table 2 Minimum bound L1k , and maximum bound L2k , of hk (u, v)/(u − v + 1)2k for the parameter range |u − v + 1| > 1, 0 < u 1, and v 0 used for bounding the remainder term Rn0 of the asymptotic expansion in Eq. (26) k
L1k
L2k
1 2 3 4 5
−5.0000 −0.0625 −1141.0 −0.3700 −1.070 × 106
1.0000 59.000 1.0000 30799 1.0000
It is straightforward to bound the remainder term by calculating bounds for the ratios hk (u, v) . (u − v + 1)2k
(28)
A simple bound on the range of x is found by noting Wβ (x, y) Wβ (x, 0) = Eβ+1 (x)
(29)
for all y > 0. Using an asymptotic expression [16] for En (x), we are safe in estimating that Wβ < 10−15 for all x > 40. Choosing a reasonable minimum value for β that would be used in Eq. (26), β = 40, bounds on Eq. (28) may be determined for the parameter range |u − v + 1| > 1, where 0 < u 1 and v 0. In Table 2, we have listed the minimum, L1k , and maximum, L2k , values for the first few k. Substituting Eq. (21) into Eq. (18), the remainder term is bounded as Wβ β −k · L1k < Rk0 < Wβ β −k · L2k .
(30)
3. Method of solution Given a starting value Wβ0 (β0 = 0, 12 , 1, etc. . . . ), an ending value Wβasymp , and the recursion of Eq. (8), the problem of determining the set Wβ0 . . . Wβ . . . Wβasymp is equivalent to finding the solution vector of Wβ 1 0 0 ··· ··· ··· ··· 0 0 Wβ0 −y 0 ··· ··· ··· 0 Wβ0 +1 e−(x+y) x β0 + 1 −(x+y) −y ··· ··· ··· 0 x β0 + 2 0 Wβ0 +2 e . ··· 0 0 x β0 + 3 −y · · · 0 .. .. = (31) · . . . . .. .. .. . . . . . . . . · · · · · · · · · . .. .. . .. e−(x+y) ··· ··· ··· · · · x βasymp − 1 −y Wβ −1 asymp Wβasymp 0 ··· ··· ··· ··· ··· 0 1 Wβ asymp
This simple linear tri-diagonal system may be solved using standard methods (e.g., Gaussian elimination or LU decomposition). Results for a variety of parameter values are provided in Table 3. We preferred starting with the value β0 = 1, since W0 is singular at x = 0 whereas W1 is not. It was found that including six terms in the asymptotic expansion of Eq. (26) is sufficient to obtain results to 15 significant figures. The solution vectors were generated by C language versions of the popular LAPACK routines [17] using “long double” precision arithmetic; real*8 or “double” precision tends to limit the number of significant digits roughly to six or seven. The accuracy of the results in Table 3 extends to all digits shown and has been checked via direct numerical integration using the Maple software package [18].
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
7
Table 3 Values of Wβ (x, y) x = 1.000 y = 20.00
x = 10.00 y = 9.000
x = 8.000 y = 0.100
x = 0.00001 y = 250.0
β
βasymp = 60
βasymp = 100
βasymp = 80
βasymp = 300
1 2 4 6 7 8 9 10 11 12 20 25 30 40 50 51 52 60 75 80 90 100
0.000025446720827 0.000006668340198 0.000000624248919 0.000000086616653 0.000000037035507 0.000000017255347 0.000000008716002 0.000000004747055 0.000000002771415 0.000000001723718 0.000000000183145 0.000000000094353 0.000000000060898 0.000000000034638 0.000000000023970 0.000000000023248 0.000000000022568 0.000000000018276
0.000000001150518 0.000000000998764 0.000000000779794 0.000000000632428 0.000000000576130 0.000000000528266 0.000000000487181 0.000000000451610 0.000000000420568 0.000000000393283 0.000000000256506 0.000000000209870 0.000000000177348 0.000000000135156 0.000000000109072 0.000000000107003 0.000000000105011 0.000000000091387 0.000000000073481 0.000000000068972 0.000000000061429 0.000000000055370
0.000031158613748 0.000028442098138 0.000024166924896 0.000020970651323 0.000019660563227 0.000018500150928 0.000017465751596 0.000016538337063 0.000015702453169 0.000014945432766 0.000010765987568 0.000009156607197 0.000007963400162 0.000006314346311 0.000005229684208 0.000005141318681 0.000005055882924 0.000004462459473 0.000003657158329 0.000003449583424
0.003941537912348 0.000015960419059 0.000000001534731 0.000000000000491 0.000000000000012 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000
It might be expected that in some extreme parameter regimes and particularly when Wβ0 and Wβasymp differ by several orders of magnitude, that maintaining precision would be a significant problem. An example might be when y is very large and substantially larger than x. In this case βasymp must also be chosen to be very large, in order that the asymptotic expansion be guaranteed to converge. If β only up to a reasonable value is needed, say β = 10, upward recursion maybe preferred. We have found that as long as β is chosen to be larger than x and y, all values are computed correctly to within 15 decimal places and typically better. This is illustrated by the fourth column of Table 3. Of course, eventually βasymp may become so large that the tri-diagonal matrix defined by the recursion cannot be inverted. Upward recursion is then stable in this situation up to β ∼ y as may be verified by performing a standard recursion analysis [19] of Eq. (8) with x = 0. We have also found that even when x is adjusted so that the starting value is on the order of Wβ0 = 10−14 and the asymptotic value is much smaller, all values in between are computed correctly to within the stated precision. Surprisingly, the use of this boundary value approach and upward recursion is sufficient to cover the entire range of integer and half odd integer β and x, y > 0 for which Wβ > 10−15 and no other expansions for the higher integrals are needed.
Acknowledgements We would like to acknowledge Samuel B. Trickey and Frank E. Harris for useful discussions. In addition, we would like to thank Trickey for helpful editing suggestions and Harris for assistance with the calculation of W0 (x, y) and W1 (x, y). This work was supported by the National Science Foundation under grant ITR DMR0218957.
8
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
References [1] M.S. Hantush, C.E. Jacob, Trans. Amer. Geophys. Union 36 (1955) 95; M.S. Hantush, Trans. Amer. Geophys. Union 37 (1956) 702. [2] I. Flamant, J.G. Fripiat, J. Delhalle, F.E. Harris, Theo. Chem. Acc. 104 (2000) 350. [3] B.I. Dunlap, J.W.D. Connolly, J.R. Sabin, J. Chem. Phys. 71 (1979) 3396; J.W. Mintmire, J.R. Sabin, S.B. Trickey, Phys. Rev. B 26 (1982) 1743 and references therein. [4] S.B. Trickey, J.A. Alford, J.C. Boettger, in: J. Lesczynski (Ed.), Computation Materials Science, in: Theoretical Computational Chemistry, vol. 15, Elsevier, 2004 (Chapter 6). [5] U. Birkenheuer, Quantenmechanische Beschreibung von Oberflachen und geordneten Adsorbatsystemen, Methoden und Andwenduncgen, Ph.D. Thesis, T.U. Munich, 1994. [6] G.W. Thomas, R.G. Keys, A.C. Reynolds Jr., J. Hydrol. 36 (1978) 173. [7] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1964 (Eq. (7.4.33)). [8] I. Flamant, PhD Thesis, Presses Universitaires de Namur, University of Namur; I. Flamant, J.G. Fripiat, J. Delhalle, Int. J. Quantum Chem. 70 (1998) 1045. [9] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 134 (2001) 78. [10] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 146 (2002) 274. [11] F.E. Harris, Comput. Phys. Comm. 146 (2002) 271. [12] F.E. Harris, Int. J. Quantum Chem. 70 (1998) 623; F.E. Harris, Int. J. Quantum Chem. 63 (1997) 913. [13] E.S. Kryachko, Int. J. Quantum Chem. 78 (2001) 303. [14] F.E. Harris, Int. J. Quantum Chem. 81 (2001) 332. [15] W. Gautschi, J. Res. Nat. Bureau of Standards 62 (1959) 123; J.G. van der Corput, J. Franklin, Nederl. Akad. Wetensch. Proc. Ser. A 213 (1951). [16] Eq. (5.1.51) of Ref. [7]. [17] http://www.netlib.org/clapack. [18] Maple 8 is a product of the Symbolic Computation Group at the University of Waterloo, Ontario, Canada. [19] See, for example: W. Gautschi, J. ACM 8 (1961) 21.
Calculation of the generalized leaky aquifer integral J.A. Alford Quantum Theory Project, Department of Physics, University of Florida, Gainesville, FL 32611, USA Received 13 April 2005; accepted 30 June 2005 Available online 8 August 2005
Abstract We provide a convenient method for the evaluation of the generalized leaky aquifer integral Wβ (x, y), for all real x, y > 0 and integers and half odd integers β 12 . The key ingredient to obtaining values for higher β is a new asymptotic expansion of the integral that is applicable for a wide range of these parameters, particularly large β. An asymptotic value and a value, calculated via previously published methods, for the same x and y but small β, form a well defined two-point boundary value problem that is easily solved. 2005 Elsevier B.V. All rights reserved. PACS: 71.15.-m Keywords: Leaky aquifer; One-dimensional systems; Fourier space methods; Asymptotic expansions; Ewald summation
1. Introduction The practical relevance of the leaky aquifer integral, defined as ∞ W0 (x, y) =
e−xt−y/t dt t
(1)
1
seems to have arisen first in the study of radial fluid flow from artesian aquifers [1], whence the origin of its name. A more general form of the integral ∞ Wβ (x, y) =
e−xt−y/t dt t β+1
1
E-mail address: [email protected] (J.A. Alford). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.06.012
(2)
2
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
for integer and half odd integer β has recently appeared in the Fourier space formulation of electronic structure methods used to calculate properties of extended periodic systems [2]. Before moving on to the technical aspects of how to compute this integral efficiently, we briefly review the context which has motivated this necessity and thus the current work. Common to many electronic structure methods is the need to evaluate lattice sums of Coulomb integrals involving local basis functions, for instance f (r − A)f (r − B − R) V= (3) dr dr , |r − r + R| R
where the basis functions f are centered about the positions A and B + R, where R is a direct space lattice vector. Such integrals are required for example for the calculation of the Hartree term of standard Hamiltonians and for variational Coulomb fitting of the electronic density to basis functions [3,4]. Although various approaches have been introduced to truncate the lattice sum in direct space properly, maintaining precise control over the accuracy has remained a difficult problem [2]. For example, the popular multipolar expansion used for obtaining the long range contributions to the sum converges prohibitively slowly for metallic or near-metallic systems. An alternative to direct space methods is the Ewald summation technique, a mixed space approach in which the short range contributions are evaluated in real space and the long range part is converted into a reciprocal space sum that is also relatively short range and can be computed efficiently [4,5]. To carry out the latter sum, one needs the Fourier transform of the Coulomb integral, the form of which will depend on the specific selection of basis orbitals f . The Hermite Gaussian orbitals provide a rather convenient choice, as the higher angular momenta orbitals of that set can be simply generated via derivatives with respect to an orbital center f (i, a, A, r) = (∂Ax )ix (∂Ay )iy (∂Az )iz e−a|r−A| . 2
(4)
Using these orbitals, it may be shown [4,5] that the Fourier space analog of Eq. (3) is given by ij (G), Vij = V G
ij (G) = (−1) (∂r0x ) (∂r0y ) (∂r0z ) e V |j|
nx
ny
nz iG·r0
e
−iG·r
RN
1 dr
e−γ |r+r0⊥|
2 u2
du,
(5)
(6)
0
where N is the number of dimensions for which the system is periodic, x, y, and z as subscripts refer to the usual Cartesian directions, r0 = A − B, G is a reciprocal lattice vector, and γ = a + b. i and j are multi-indices, the Cartesian components of which are nonnegative integers, n = i + j, |j| = jx + jy + jz , r0 is the component of r0 in the space spanned by RN , and r0⊥ is the component orthogonal to this space. (Note that, irrespective of N , the √ system is three-dimensional.) With the substitution u = 1/ t , the dependence on N becomes more clear ij (G) = (−1)|j| (∂r0x )nx (∂r0y )ny (∂r0z )nz e−iG·r0 V
∞
e
−|G|2 /(4γ )t−γ r02 /t ⊥
2t (3−N )/2
dt.
(7)
1
Note that the derivatives with respect to components of r0 generate higher powers of t in the denominator of the integrand. These integrals have the general form of Eq. (2) if we set x = |G|2 /4γ , y = γ r02⊥ , and β + 1 = (3 − N )/2 + k, where k is an integer resulting from differentiation. They also obey a recursion relation first discussed by Thomas [6] xWβ−1 (x, y) + βWβ (x, y) − yWβ+1 (x, y) = e−(x+y) that is easily derived through integration by parts.
(8)
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
3
For a 3D crystal (N = 3), there is no finite orthogonal vector r0⊥ (y = 0) and Eq. (7) simply reduces to a sum 2 over products of components of G and the factor e−|G| /4γ |G|2 . If N = 2, the indefinite form of the integral can be brought into a known form [7] through a substitution like t = v 2 and evaluated directly. The result for n = 0, corresponding to β + 1 = 12 , is √ √ √ √ √ 00 (G) = √π e2 xy erfc( x + √y ) + e−2 xy erfc( x − √y ) . V 4 x
(9)
To obtain the higher integrals in this particular case, β = n/2, where n is an odd integer, two coupled, twoterm recursion relations may be derived, through derivatives with respect to y, that are generally more stable than Eq. (8). The N = 1 case also requires integrals of the form of Eq. (2) where β is a nonnegative integer. Eq. (8) is the only known applicable recursion for N = 1 and is the primary focus of this paper. The analysis presented here also holds for β equal to half odd integers and thus may be applied when N = 2 as well. The difficulty of evaluating the higher integrals Wβ (x, y) for integer β has in fact prevented the general application of the Ewald summation technique to systems with 1D periodicity, with the notable exception of the studies by Flamant et al. [8] and Langridge et al. [9]. Flamant’s treatment was limited to the use of s orbital basis functions only, in effect restricting β to zero. Langridge and coworkers encountered the leaky aquifer integral in their application of the Ewald technique to the Madelung sums connected with the multipolar expansion. In their expressions the parameters that enter the integral (here x and y) are scaled by the Ewald split parameter µ, and the following expansion is used for the calculation Wβ (x/µ, µy) =
∞ (−µy)m Eβ+m+1 (x/µ). m!
(10)
m=0
The Ewald parameter is adjustable; in practice it is modified in some cases in order to retain sufficient accuracy in this expression [10]. The parameter must play a dual role in this respect and thus the generality of the method is somewhat unclear, although it may represent a significant improvement to more common associated techniques. Harris, however, pointed out that their procedure is not optimal in general [11], a point not disputed by Langridge et al. [10]. The bulk of other previous work related to the integral has focused on deriving efficient expansions of W0 (x, y) for various ranges of x and y [12]. Kryachko [13] has shown that the zeroth-order integral may be written as an infinite sum over Legendre functions of the second kind, and Harris [14] has extended this result to similar forms for W1 (x, y) and W2 (x, y). The obvious method of reaching higher derivatives, namely upward recursion using Eq. (8), is highly unstable except in the case that y is much larger than x and β. It is our intention that the present work will remove the barrier to the general application of the Ewald technique to 1D periodic systems. To our knowledge, no asymptotic form for Wβ has been published, preventing any attempt at or analysis of downward recursion. It is somewhat obvious, however, that the downward recursion will be unstable except in the case that x is much larger than y and β. Rather than working with such uncontrolled approaches to obtain the higher integrals, we view their evaluation as the process of solving a two-point boundary value problem, where the end points are given by Wβ0 with β0 = − 12 , 0, 12 , or 1, and Wβ0 +m , m being a large integer. With the end-point values supplied for a given x and y, the problem of determining the sequence Wβ0 , Wβ0 +1 , Wβ0 +2 , . . . , Wβ0 +m−1 , Wβ0 +m can be solved using standard algebraic methods. In the next section, we derive a general asymptotic expansion for large β and describe its range of applicability. The method of solution for a range of integer or half odd integer β along with some numerical results and discussion is given in Section 3.
4
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
2. Asymptotic formula We follow the procedure used by Gautschi [15] for the exponential integrals En (x) to derive a new expression for Wβ (x, y). Consider the integral t2 I = j
f (t) = −xt − y/t − β ln t,
ef (t) dt,
(11)
t1
where the value of j denotes a specific interval over t , 0, t1 = 1, t2 = ∞, j = 1, t1 = 1, t2 = t ∗ − ε, 2, t1 = t ∗ + ε, t2 = ∞
(12)
and I 0 = lim (I 1 + I 2 ).
(13)
ε→0+
Let t ∗ be defined by f (t ∗ ) = 0, 2 β y β ∗ + . t =− + 2x 2x x
(14)
We deal with the case for which 1 < t ∗ < ∞. Except for t ∗ = 1, the following is obviously valid for t ∗ not lying in this interval, for example, when x > y. With the definitions g0 (t) = 1/f (t), integration of
Ij
gk+1 (t) = gk (t)g0 (t),
(15)
by parts yields the expansion
j
j
j
j
j
j
I j = v0 − v1 + v2 − v3 + · · · + (−1)k−1 vk−1 + Rk ,
(16)
where vk = gk (t2 )ef (t2 ) − gk (t1 )ef (t1 )
(17)
and j Rk
t2 = (−1)
k
gk−1 (t)ef (t) dt.
(18)
t1
All gk (t) in addition to g0 (t) are singular at t = t ∗ . Adding I 1 and I 2 , we see that all of the singular terms cancel term by term,
∗ ∗ lim vk1 + vk2 = lim gk (t ∗ + ε)ef (t +ε) − gk (t ∗ − ε)ef (t −ε) + gk (∞)ef (∞) − gk (1)ef (1) ε→0+
ε→0+
= gk (∞)ef (∞) − gk (1)ef (1) ≡ vk0 .
(19)
Thus I 0 may be expressed as I0 =
n−1 (−1)k vk0 + Rn0 . k=0
(20)
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
5
Before analyzing the remainder term, we construct a recursion equation for generating the gk . Calculation of the first few gk using Eq. (15) shows that each takes the form gk (u, α) =
(−1)k hk (u, α) g0 (u, α), (1 + u − αu )2k β k
(21)
where xy xt , α= 2 (22) β β and the hk (u, α) are polynomials in u and α. The recurring pattern that determines the form of the polynomials is seen to be
hk+1 (u, α) = (1 + u − α/u) · ∂u uhk (u, α) − (2k + 1)(1 + α/u2 )uhk (u, α) (23) u=
which gives the recursion
hk+1 (u, α) = u(u + 1) − α
hk (u, α) +
1+k (1 − 2ku) − 2α hk (u, α), u
(24)
where the prime now denotes differentiation with respect to u rather than t . Setting α = 0, equivalent to putting y = 0, Eqs. (21) and (24) manifestly reduce to the forms originally derived by Gautschi [15] for the exponential integral as they should. Evaluation of the gk (u, α) at the proper limits to obtain the vk0 in Eq. (19) and the further substitutions u = x/β,
v = y/β
(25)
finally yield Wβ (x, y) =
n−1 hk (u, v) e−(x+y) + Rn0 . x −y +β (u − v + 1)2k β k
(26)
k=0
The asymptotic prefactor to the expansion has been left in terms of x and y to show that it is in fact the asymptotic value expected from the recursion in Eq. (8), at least for x > y. The first few polynomials hk (u, v) are listed in Table 1. The parameter space over which this expansion converges is at least all u, v > 0 and β > 1 such that x −y +β >1
(27)
as can be seen by rewriting the sum as n−1 k=0
hk · β k . (x − y + β)2k
Table 1 Polynomials hk (u, v) for the asymptotic expansion of Wβ (x, y) generated by Eq. (24) k
hk (u, v)
0 1 2 3 4 5
1 1 − 2v 1 − 4v − 2u + 6uv + 6v 2 1 − 6v − 8u + 40uv + 14v 2 + 6u2 − 24u2 v − 72uv 2 − 24v 3 1 − 8v + 154uv + 28v 2 − 348u2 v − 432uv 2 − 36v 3 + 120u3 v + 720u2 v 2 + 720uv 3 − 22u + 58u2 − 24u3 + 120v 4 1 − 10v + 468uv + 44v 2 − 2624u2 v − 1820uv 2 − 152v 3 − 52u + 328u2 + 3108u3 v + 8436u2 v 2 + 3564uv 3 − 108v 4 − 720u4 v − 7200u3 v 2 − 14400u2 v 3 − 7200uv 4 − 444u3 + 120u4 − 720v 5
6
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
Table 2 Minimum bound L1k , and maximum bound L2k , of hk (u, v)/(u − v + 1)2k for the parameter range |u − v + 1| > 1, 0 < u 1, and v 0 used for bounding the remainder term Rn0 of the asymptotic expansion in Eq. (26) k
L1k
L2k
1 2 3 4 5
−5.0000 −0.0625 −1141.0 −0.3700 −1.070 × 106
1.0000 59.000 1.0000 30799 1.0000
It is straightforward to bound the remainder term by calculating bounds for the ratios hk (u, v) . (u − v + 1)2k
(28)
A simple bound on the range of x is found by noting Wβ (x, y) Wβ (x, 0) = Eβ+1 (x)
(29)
for all y > 0. Using an asymptotic expression [16] for En (x), we are safe in estimating that Wβ < 10−15 for all x > 40. Choosing a reasonable minimum value for β that would be used in Eq. (26), β = 40, bounds on Eq. (28) may be determined for the parameter range |u − v + 1| > 1, where 0 < u 1 and v 0. In Table 2, we have listed the minimum, L1k , and maximum, L2k , values for the first few k. Substituting Eq. (21) into Eq. (18), the remainder term is bounded as Wβ β −k · L1k < Rk0 < Wβ β −k · L2k .
(30)
3. Method of solution Given a starting value Wβ0 (β0 = 0, 12 , 1, etc. . . . ), an ending value Wβasymp , and the recursion of Eq. (8), the problem of determining the set Wβ0 . . . Wβ . . . Wβasymp is equivalent to finding the solution vector of Wβ 1 0 0 ··· ··· ··· ··· 0 0 Wβ0 −y 0 ··· ··· ··· 0 Wβ0 +1 e−(x+y) x β0 + 1 −(x+y) −y ··· ··· ··· 0 x β0 + 2 0 Wβ0 +2 e . ··· 0 0 x β0 + 3 −y · · · 0 .. .. = (31) · . . . . .. .. .. . . . . . . . . · · · · · · · · · . .. .. . .. e−(x+y) ··· ··· ··· · · · x βasymp − 1 −y Wβ −1 asymp Wβasymp 0 ··· ··· ··· ··· ··· 0 1 Wβ asymp
This simple linear tri-diagonal system may be solved using standard methods (e.g., Gaussian elimination or LU decomposition). Results for a variety of parameter values are provided in Table 3. We preferred starting with the value β0 = 1, since W0 is singular at x = 0 whereas W1 is not. It was found that including six terms in the asymptotic expansion of Eq. (26) is sufficient to obtain results to 15 significant figures. The solution vectors were generated by C language versions of the popular LAPACK routines [17] using “long double” precision arithmetic; real*8 or “double” precision tends to limit the number of significant digits roughly to six or seven. The accuracy of the results in Table 3 extends to all digits shown and has been checked via direct numerical integration using the Maple software package [18].
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
7
Table 3 Values of Wβ (x, y) x = 1.000 y = 20.00
x = 10.00 y = 9.000
x = 8.000 y = 0.100
x = 0.00001 y = 250.0
β
βasymp = 60
βasymp = 100
βasymp = 80
βasymp = 300
1 2 4 6 7 8 9 10 11 12 20 25 30 40 50 51 52 60 75 80 90 100
0.000025446720827 0.000006668340198 0.000000624248919 0.000000086616653 0.000000037035507 0.000000017255347 0.000000008716002 0.000000004747055 0.000000002771415 0.000000001723718 0.000000000183145 0.000000000094353 0.000000000060898 0.000000000034638 0.000000000023970 0.000000000023248 0.000000000022568 0.000000000018276
0.000000001150518 0.000000000998764 0.000000000779794 0.000000000632428 0.000000000576130 0.000000000528266 0.000000000487181 0.000000000451610 0.000000000420568 0.000000000393283 0.000000000256506 0.000000000209870 0.000000000177348 0.000000000135156 0.000000000109072 0.000000000107003 0.000000000105011 0.000000000091387 0.000000000073481 0.000000000068972 0.000000000061429 0.000000000055370
0.000031158613748 0.000028442098138 0.000024166924896 0.000020970651323 0.000019660563227 0.000018500150928 0.000017465751596 0.000016538337063 0.000015702453169 0.000014945432766 0.000010765987568 0.000009156607197 0.000007963400162 0.000006314346311 0.000005229684208 0.000005141318681 0.000005055882924 0.000004462459473 0.000003657158329 0.000003449583424
0.003941537912348 0.000015960419059 0.000000001534731 0.000000000000491 0.000000000000012 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000
It might be expected that in some extreme parameter regimes and particularly when Wβ0 and Wβasymp differ by several orders of magnitude, that maintaining precision would be a significant problem. An example might be when y is very large and substantially larger than x. In this case βasymp must also be chosen to be very large, in order that the asymptotic expansion be guaranteed to converge. If β only up to a reasonable value is needed, say β = 10, upward recursion maybe preferred. We have found that as long as β is chosen to be larger than x and y, all values are computed correctly to within 15 decimal places and typically better. This is illustrated by the fourth column of Table 3. Of course, eventually βasymp may become so large that the tri-diagonal matrix defined by the recursion cannot be inverted. Upward recursion is then stable in this situation up to β ∼ y as may be verified by performing a standard recursion analysis [19] of Eq. (8) with x = 0. We have also found that even when x is adjusted so that the starting value is on the order of Wβ0 = 10−14 and the asymptotic value is much smaller, all values in between are computed correctly to within the stated precision. Surprisingly, the use of this boundary value approach and upward recursion is sufficient to cover the entire range of integer and half odd integer β and x, y > 0 for which Wβ > 10−15 and no other expansions for the higher integrals are needed.
Acknowledgements We would like to acknowledge Samuel B. Trickey and Frank E. Harris for useful discussions. In addition, we would like to thank Trickey for helpful editing suggestions and Harris for assistance with the calculation of W0 (x, y) and W1 (x, y). This work was supported by the National Science Foundation under grant ITR DMR0218957.
8
J.A. Alford / Computer Physics Communications 173 (2005) 1–8
References [1] M.S. Hantush, C.E. Jacob, Trans. Amer. Geophys. Union 36 (1955) 95; M.S. Hantush, Trans. Amer. Geophys. Union 37 (1956) 702. [2] I. Flamant, J.G. Fripiat, J. Delhalle, F.E. Harris, Theo. Chem. Acc. 104 (2000) 350. [3] B.I. Dunlap, J.W.D. Connolly, J.R. Sabin, J. Chem. Phys. 71 (1979) 3396; J.W. Mintmire, J.R. Sabin, S.B. Trickey, Phys. Rev. B 26 (1982) 1743 and references therein. [4] S.B. Trickey, J.A. Alford, J.C. Boettger, in: J. Lesczynski (Ed.), Computation Materials Science, in: Theoretical Computational Chemistry, vol. 15, Elsevier, 2004 (Chapter 6). [5] U. Birkenheuer, Quantenmechanische Beschreibung von Oberflachen und geordneten Adsorbatsystemen, Methoden und Andwenduncgen, Ph.D. Thesis, T.U. Munich, 1994. [6] G.W. Thomas, R.G. Keys, A.C. Reynolds Jr., J. Hydrol. 36 (1978) 173. [7] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1964 (Eq. (7.4.33)). [8] I. Flamant, PhD Thesis, Presses Universitaires de Namur, University of Namur; I. Flamant, J.G. Fripiat, J. Delhalle, Int. J. Quantum Chem. 70 (1998) 1045. [9] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 134 (2001) 78. [10] D.J. Langridge, J.F. Hart, S. Crampin, Comput. Phys. Comm. 146 (2002) 274. [11] F.E. Harris, Comput. Phys. Comm. 146 (2002) 271. [12] F.E. Harris, Int. J. Quantum Chem. 70 (1998) 623; F.E. Harris, Int. J. Quantum Chem. 63 (1997) 913. [13] E.S. Kryachko, Int. J. Quantum Chem. 78 (2001) 303. [14] F.E. Harris, Int. J. Quantum Chem. 81 (2001) 332. [15] W. Gautschi, J. Res. Nat. Bureau of Standards 62 (1959) 123; J.G. van der Corput, J. Franklin, Nederl. Akad. Wetensch. Proc. Ser. A 213 (1951). [16] Eq. (5.1.51) of Ref. [7]. [17] http://www.netlib.org/clapack. [18] Maple 8 is a product of the Symbolic Computation Group at the University of Waterloo, Ontario, Canada. [19] See, for example: W. Gautschi, J. ACM 8 (1961) 21.