PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Sci Comput. Author manuscript; available in PMC 2010 December 7.
Published in final edited form as:
J Sci Comput. 2009 December 1; 41(3): 411–435.
doi:  10.1007/s10915-009-9307-z
PMCID: PMC2998241
NIHMSID: NIHMS201927

A Sixth-order Image Approximation to the Ionic Solvent Induced Reaction Field

Abstract

A recent article by Deng and Cai introduced fourth-order image approximations to the reaction field for a charge inside a dielectric sphere immersed in a solvent of low ionic strength. To represent such a reaction field, the image approximations employ a point charge at the classical Kelvin image point and two line charges that extend from the Kelvin image point along the radial direction to infinity. In this paper, a sixth-order image approximation is developed, using the same point charge with three different line charges. Procedures on how to discretize the line charges by point image charges and how to implement the resulting point image approximation in O(N) complexity for potential and force field calculations are included. Numerical results demonstrate the sixth-order convergence rate of the image approximation and the O(N) complexity of the fast implementation of the point image approximation.

Keywords: Method of images, Reaction field, Ionic solvent, Hybrid solvation model

1 Introduction

Following [5, 6], this paper concerns fast and accurate calculation of electrostatic interactions among point charges inside a spherical dielectric cavity embedded in an ionic solvent of dissimilar dielectric constant. Such a problem could be encountered in many applications such as hybrid explicit/implicit solvent biomolecular dynamics simulations [15, 16, 22], in which biomolecules and a part of solvent molecules within a dielectric cavity are explicitly modeled while a surrounding dielectric continuum is used to model bulk effects of the solvent outside the cavity.

The point charges in the dielectric cavity will polarize the surrounding dielectric medium, which in turn generates a reaction field to the electric field throughout the cavity. The electric potential field inside the cavity is thus expressed as Ψin = ΨS + ΨRF, where ΨS is the Coulomb potential given by the Coulomb's Law, and ΨRF is the reaction field which will dominate the computational cost for calculating the electric field inside the cavity. Therefore, fast and accurate calculation of such a reaction field will have a wide impact on computational simulations for chemical and biological systems involving electrostatic interactions within a solvent.

In case of a spherical cavity, a popular approach to calculate the reaction field is the method of images, in which the reaction field is represented in terms of potentials of discrete image charges. For the pure water solvent, namely, with no ions present in the solvent, a variety of approaches exist for calculating the reaction field for charges inside the spherical cavity, for example, the high-order accurate multiple image approximation [2] and references therein. For an ionic solvent, by assuming that the ionic strength of the solvent is low enough so that the product of the inverse Debye screening length of the solvent and the radius of the spherical cavity is less than one, image approximations of various orders (up to fourth-order) to the ionic solvent induced reaction field have been developed. More specifically, in [5], a first- and a second-order image approximations are presented using a point image charge at the classical Kelvin image point and a line image charge that extends from this Kelvin image point along the radial direction to infinity. Later in [6], a fourth-order image approximation and its improved version are obtained using the same point image charge together with two line image charges that extend from the Kelvin image point along the radial direction to infinity. Following the same procedure as used in [5, 6], we shall develop in this paper a sixth-order accurate image approximation to the ionic solvent induced reaction field, which will provide higher accuracy in the evaluation of potential and force fields inside the cavity.

The structure of the paper is as follows. In Sect. 2, the exact series solution to the ionic solvent induced reaction field due to a point charge inside a dielectric sphere is briefly reviewed. Then, in Sect. 3, after a short description of the previous fourth-order image approximations, a sixth-order image approximation to the reaction field by a point image charge at the classical Kelvin image point and three line image charges is developed. How to discretize the line image charges by point image charges is then discussed in Sect. 4. Next in Sects. 5 and 6, details are given on how to calculate in O(N) complexity the potential and the force fields among N charges inside the spherical cavity, respectively. Numerical results are presented in Sect. 7 to validate the convergence property and investigate the efficiency of the sixth-order image approximation. Finally a conclusion is made in Sect. 8.

2 Exact Series Solution to the Ionic Solvent Induced Reaction Field

By the principle of linear superposition, the reaction field due to a single point charge q inside a spherical cavity of radius a centered at the origin only needs to be considered. The sphere has a dielectric constant εi, and the surrounding ionic solvent is represented as a homogeneous medium of a dielectric constant εo. The point charge q is located on the x-axis inside the sphere at a distance rS < a from the spherical center, as shown in Fig. 1.

Fig. 1
A point charge and a dielectric sphere immersed in an ionic solvent

It is well-known that the total electric potential Ψin(r) inside the cavity is given by the solution of the Poisson equation

(iΨin(r))=4πqδ(rrs),
(2.1)

where δ(r) denotes the Dirac delta function. Outside the cavity, on the other hand, by assuming that the mobile ion concentration follows the Debye-Hückel theory, namely, the mobile ion charges follow a Boltzmann distribution in the mean field approximation, for a solvent of low ionic strength, the electric potential Ψout(r) is then given by the solution of the linearized Poisson-Boltzmann equation (LPBE) [14]

2Ψout(r)λ2Ψout(r)=0,
(2.2)

where λ is the inverse Debye screening length determined by the ionic concentration and the exterior dielectric constant εo (for the pure water solvent, λ 0). On the interface Γ of the dielectric cavity and its surrounding dielectric medium, the following two boundary conditions are to be satisfied for the continuities of the potential and the fluxes along the normal direction of the interface

ΨoutΓ=ΨinΓandoΨoutnΓ=iΨinnΓ,
(2.3)

where n is the unit outward vector normal to the surface of the cavity.

Using the classical electrostatic theory, the reaction field of the spherical dielectric can be solved analytically [5]. More precisely, with respect to a spherical coordinate system (r, θ, [var phi]) originating in the center of the sphere (the pole is denoted by the x-axis in this paper), due to the azimuthal symmetry, the reaction field at a point r = (r, θ, [var phi]) inside the sphere takes on the form

ΦRF(r)=n=0AnrnPn(cosθ),
(2.4)

where Pn(x) represent the Legendre polynomials and An are the expansion coefficients given by

An=qia1rKni(n+1)kn(u)+oukn(u)inkn(u)oukn(u),n0,
(2.5)

where u = λa, rK = a2/rS with rK = (rK, 0, 0) denoting the classical Kelvin image point, and kn(r) are the modified spherical Hankel functions [1, 9]

kn(r)=π2rerk=0n(n+k)!k!(nk)!1(2r)k,n0.
(2.6)

In theory, any desired degree of accuracy can be obtained using the direct series expansion (2.4). In the case that the point charge is close to the spherical boundary, when calculating the reaction field at an observation point also close to the boundary, the convergence of the series expansion is slow due to the fact that r/rK = rrS/a2 ≈ 1, requiring a great number of terms in the series expansion to achieve satisfactory accuracy in the reaction field.

3 Line Image Approximations to the Reaction Field

Let us now turn ourselves to the problem of finding image charges outside the spherical region giving rise to the reaction potential inside the sphere. For the pure water solvent, such image charges include a point charge at the classical Kelvin image point and a continuous line charge extending along the radial direction from this Kelvin image point to infinity [17, 20, 21]. For an ionic solvent, as mentioned earlier, by assuming that the ionic strength of the solvent is low enough such that the product of the inverse Debye screening length λ (which is proportional to the square root of the ionic concentration) and the radius of the spherical cavity a is less than one (u = λa < 1), several image approximations of various orders (in terms of u = λa) to the ionic solvent induced reaction field have been developed [5, 6]. It should be emphasized again that the assumption of u = λa < 1 is physically justifiable since the condition of low ionic strength is required for the linearization of the Poisson-Boltzmann equation to be meaningful.

3.1 Previous Fourth-Order Image Approximations

The key idea in the development of the foregoing image charge methods is to approximate kn(u)ukn(u) with simple rational functions of n when u is small. In fact, by applying the expansion of the modified spherical Hankel function

kn(r)=π(2n2)!(n1)!(2(2n1)(2r)n+114(2r)n1)+O(1rn3),n2,
(3.1)

which essentially implies

kn(r)rkn(r)=(2+r2)+4n(r22)+(2r2)n+4n2+O(r4),n2,

a fourth-order image approximation to the ionic solvent induced reaction field can be obtained as follows [6]

ΦRF(r)=qKirrK+rKq^L1(x)irxdx+rKq^L2(x)i(1rx1x)dx+Φ^C0+Φ^C1(r)+O(u4),
(3.2)

where the constant position-independent correction potential Φ^C0 and the position-dependent correction potential Φ^C1(r) are defined as

Φ^C0=qia(C0(u)γδ^1σ^1),
(3.3)
Φ^C1(r)=qia(C1(u)γδ^11+σ^1δ^21σ^2)(rrK)cosθ,
(3.4)

respectively. Here we denote

C0(u)=i(1+u)o(1+u)o,
(3.5)
C1(u)=2(1+u)i(2+2u+u2)o(1+u)i+(2+2u+u2)o.
(3.6)

The point image charge qK and the two line image charges q^L1(x) and q^L2(x) are defined by

qK=γarSq,q^L1(x)=δ^1qa(xrK)σ^1,q^L2(x)=δ^2qa(xrK)σ^2,rKx,

respectively, where σ^1, δ^1, σ^2, δ^2 are defined in [6], and γ (εiεo)/(εi + εo).

Moreover, the accuracy of the fourth-order image approximation (3.2) can be further improved by including another position-dependent correction potential, yielding the improved fourth-order image approximation [6]

ΦRF(r)=qKirrK+rKq^L1(x)irxdx+rKq^L2(x)i(1rx1x)dx+Φ^C0+Φ^C1(r)+Φ^C2(r)+O(u4),
(3.7)

with the position-dependent correction potential Φ^C2(r) defined by

Φ^C2(r)=qia(C2(u)γδ^12+σ^1δ^22σ^2)(rrK)2P2(cosθ),

where

C2(u)=3(3+3u+u2)i(9+9u+4u2+u3)o2(3+3u+u2)i+(9+9u+4u2+u3)o.
(3.8)

3.2 Formulation of the Sixth-Order Image Approximation

To construct a sixth-order image approximation for better accuracy, by using the Taylor expansion

er=1r+r22r36+r424r5120+O(r6),

and the truncation

k=0n(n+k)!k!(nk)!1(2r)k=k=05(2nk)!k!(nk)!1(2r)nk+O(1rn6),n5,

one can arrive at the expansion of the modified spherical Hankel function in terms of 1/r as follows

kn(r)=π(2n4)!(n2)![4(2n1)(2n3)(2r)n+12n32(2r)n1+132(2r)n3]+O(1rn5),n5.
(3.9)

However, after expanding k4(r) directly in terms of 1/r, it turns out that (3.9) is also true for n = 4. Correspondingly, the derivative of the modified spherical Hankel function can be expanded as

kn(r)=π(2n4)!(n2)![8(2n3)(2n1)(n+1)(2r)n+2(2n3)(n1)(2r)n+n316(2r)n2]+O(1rn4),n4.
(3.10)

Consequently, for n ≥ 4 we obtain

kn(r)rkn(r)=r4+12r2+248(r2+8)n+32n23(r4+4r28)+(r4+20r240)n8(r2+4)n2+32n3+O(r6).

Inserting this approximation into (2.5) leads to

An=qia1rKna1+a2n+a3n2+a4n3b1+b2n+b3n2+b4n3+O(u6),n4,

where

a1=24(io)+12iu2+12ou2+iu4+3ou4,a2=40(io)+4iu220ou2+iu4ou4,a3=32(io)8iu2+8ou2,a4=32(io),b1=3o(u4+4u28),b2=64o+24(i+o)+12iu2+20ou2+iu4+ou4,b3=8(8i+4o+(i+o)u2),b4=32(i+o).

After some algebraic manipulations, the expansion coefficients An can be further expressed as

An=qia1rKn(γ+α1+α2n+α3n2β1+β2n+β3n2+n3)+O(u6),n4,

where

α1=a1b4γβ1,β1=b1b4,α2=a2b4γβ2,β2=b2b4,α3=a3b4γβ3,β3=b3b4.

Now using the root-finding formula for a cubic equation,1 one can get the three roots of the equation β1 + β2n + β3n2 + n3 = 0 as

n1=13[β3A12(cos(φ3)+3sin(φ3))]=13[β32A12sin(π6+φ3)],n2=13[β3A12(cos(φ3)3sin(φ3))]=13[β32A12sin(π6φ3)],n3=13[β32A12cos(φ3)],

where

A=β323β2,B=9β2β327β12β3354,ρ=(A9)3,φ=arccos(Bρ).

Theorem 1 If 0 ≤ u ≤ 1, then the three roots n1, n2 and n3 of the cubic equation β1 + β2n + β3n2 + n3 = 0 satisfy

1<n1<0<n2<1<n3<2.

Proof Let

f(n)=(b1+b2n+b3n2+b4n3)o.

Then, we have

f(1)=4u440u220ru2120rru4,f(0)=3u4+2412u2,f(1)=2u416+4ru28r+ru4,f(2)=u4+724u28ru2+48r+2ru4,

where εr = εi/εo > 0. For 0 ≤ u ≤ 1, it can be seen that

f(1)<0,f(0)>0,f(1)<0,f(2)>0,

indicating that the equation f(n) = 0 has exactly one root in each of the three open intervals (−1, 0), (0, 1) and (1, 2).

On the other hand, note that for 0 ≤ [var phi]π, we have

n113(β3A12)n213(β3+A12)n3,

and that n1, n2 and n3 are also the three roots of the equation f(n) = 0. Therefore, we must have

1<n1<0<n2<1<n3<2.

Then, using partial fractions gives us

An=qia1rKn(γ+δ1n+σ1+δ2nσ2+δ3nσ3)+O(u6),n4,
(3.11)

where σ1 = −n1, σ2 = n2, σ3 = n3, and

δ1=σ12α3σ1α2+α1(σ1+σ2)(σ1+σ3),δ2=σ22α3+σ2α2+α1(σ1+σ2)(σ2σ3),δ3=σ32α3+σ3α2+α1(σ1+σ3)(σ2+σ3).

On the other hand, for n ≤ 3, directly applying the exact expressions of k0(r), k1(r), k2(r) and k3(r), in fact we can arrive at

k0(r)rk0(r)=11+r,k1(r)rk1(r)=1+r2+2r+r2,k2(r)rk2(r)=3+3r+r29+9r+4r2+r3,k3(r)rk3(r)=90+30r+6r2+r3270+150r+36r2+6r3+r4,

and accordingly, we have

Aj=Cj(u)qia(1rKj),j=0,1,2,3,
(3.12)

where C0(u), C1(u) and C2(u) are given in (3.5), (3.6) and (3.8), respectively, and C3(u) is defined below

C3(u)=4(15+15u+6u2+u3)i(60+60u+27u2+7u3+u4)o3(15+15u+6u2+u3)i+(60+60u+27u2+7u3+u4)o.
(3.13)

Inserting now the approximation of An given by (3.11) into (2.4), the reaction field inside the sphere can be expressed as

ΦRF(r)=n=03AnrnPn(cosθ)+S1+S2+S3+S4+O(u6),
(3.14)

where S1, S2, S3 and S4 represent the following four series

S1=γqian=4(rrK)nPn(cosθ),S2=δ1qian=41n+σ1(rrK)nPn(cosθ),S3=δ2qian=41nσ2(rrK)nPn(cosθ),S4=δ3qian=41nσ3(rrK)nPn(cosθ).

Now, the problem reduces to how to represent each of the above four series by a line image charge. To this end, we recall that the Coulomb potential ΦS(r), the potential at r due to a point charge q at rS , can be expanded in terms of the Legendre polynomials of cos θ as follows [19]

ΦS(r)=qirrS={qirn=0(rSr)nPn(cosθ),rSra,qirSn=0(rrS)nPn(cosθ),0rrS.}
(3.15)

First, in order to obtain an image representation for the series S1, we note that it can be written as

S1=γqirKarSn=0(rrK)nPn(cosθ)γqian=03(rrK)nPn(cosθ),

where the first part is exactly the expansion of the potential at r due to a point charge of magnitude qK outside the sphere at the classical Kelvin image point rK = (rK, 0, 0). Therefore, we have

S1=qKirrKγqian=03(rrK)nPn(cosθ).
(3.16)

Next, to find an image representation for the second series S2, we need the integral identity which can be regarded as a Mellin transformation

1n+σ=rKn+σrK1xn+σ+1dx,
(3.17)

and is valid for all n ≥ 0 when σ > 0. Inserting (3.17) with σ = σ1 into S2 yields

S2=rK[qL1(x)ixn=0(rx)nPn(cosθ)]dxδ1qian=031n+σ1(rrK)nPn(cosθ),

where

qL1(x)=δ1qa(xrK)σ1,rKx.
(3.18)

Note that the integrand in the above integral again becomes the expansion given by (3.15) for a charge of magnitude qL1(x) outside the sphere at the point x = (x, 0, 0). Therefore, qL1(x) can be regarded as a continuous line charge which stretches from the classical Kelvin image point rK along the radial direction to infinity. Thus, the second series S2 becomes

S2=rKqL1(x)irxdxδ1qian=031n+σ1(rrK)nPn(cosθ).
(3.19)

To find image representations for the series S3 and S4, we use a similar integral identity

1nσ=rKnσrK1xnσ+1dx,
(3.20)

which is valid for all n ≥ 1 when σ < 1 and valid for all n ≥ 2 when σ < 2. Inserting this into S3 yields

S3=rKqL2(x)i[n=01x(rx)nPn(cosθ)n=031x(rx)nPn(cosθ)]dx,

where

qL2(x)=δ2qa(xrK)σ2,rKx.
(3.21)

Then, applying the identity (expansion (3.15) with rS = x, q = 1 and εi = 1)

1rx=n=01x(rx)nPn(cosθ),r<x,

and noting that

rKqL2(x)i1x(rx)nPn(cosθ)dx=δ2qia(nσ2)(rrK)nPn(cosθ),1n3,

we obtain

S3=rKqL2(x)i(1rx1x)dxδ2qian=131nσ2(rrK)nPn(cosθ).
(3.22)

Similarly, we can obtain

S4=rKqL3(x)i(1rx1xrx2cosθ)dxδ3qian=231nσ3(rrK)nPn(cosθ),
(3.23)

where

qL3(x)=δ3qa(xrK)σ3,rKx.
(3.24)

Finally, inserting (3.16), (3.19), (3.22) and (3.23) into (3.14) and then combining like terms, we have the sixth-order line image approximation to the ionic solvent induced reaction potential

ΦRF(r)=qKirrK+rKqL1(x)irxdx+rKqL2(x)i(1rx1x)dx+rKqL3(x)i(1rx1xrx2cosθ)dx+ΦC0+ΦC1(r)+ΦC2(r)+ΦC3(r)+O(u6),
(3.25)

where ΦC0 is a position-independent correction potential defined as

ΦC0=qia(C0(u)γδ1σ1),
(3.26)

and on the other hand, ΦC1(r), ΦC2(r), and ΦC3(r) are position-dependent correction potentials given by

ΦC1(r)=qia(C1(u)γδ11+σ1δ21σ2)(rrK)cosθ,ΦC2(r)=qia(C2(u)γδ12+σ1δ22σ2δ32σ3)(rrK)2P2(cosθ),ΦC3(r)=qia(C3(u)γδ13+σ1δ23σ2δ33σ3)(rrK)3P3(cosθ).

For convenience, in terms of the Kronecker delta δ1n we define

cn=1ia(Cn(u)γδ1n+σ1δ2nσ2(1δ1n)δ3nσ3),n=1,2,3.

Then, we have

ΦCn(r)=qcn(rrK)nPn(cosθ),n=1,2,3.
(3.27)

4 Point Image Approximations to the Reaction Field

In this section, we will discuss how to approximate each line image charge involved in the image approximation in the previous section by a set of point image charges. Without losing of generality, we consider

I=rKf(x)(xrK)σdx,
(4.1)

where σ > 0. First, by introducing the change of variables rK/x = ((1 − s)/2)τ with τ > 0, we have

I=2τστ11(1s)αh(s,τ)ds,
(4.2)

where α = τσ − 1 and

h(s,τ)=2τrK(1s)τf(2τrK(1s)τ).
(4.3)

Next, a numerical quadrature will be employed to approximate the integral in (4.2). Although in principle any numerical quadrature can be used, considering that s = −1 corresponds to the Kelvin image point x = rK and that α > −1 because σ > 0 and τ > 0, the Jacobi-Gauss-Radau quadrature is particularly used in this paper. More precisely, let sm, ωm, m = 1, 2, ..., M, be the Jacobi-Gauss-Radau points and weights on the interval [−1, 1] with α = τσ − 1 and β = 0, which can be obtained with the program ORTHPOL [8]. Then, the numerical quadrature for the integral in (4.2) is

I=rKf(x)(xrK)σdx2τστm=1Mωmxmf(xm)
(4.4)

where for m = 1, 2, ..., M,

xm=rK(21sm)τ.
(4.5)

The parameter τ > 0 in the change of variables rK/x = ((1 − s)/2)τ can be used as a parameter to control the accuracy of numerical approximations. When τ = 1/σ we have α = 0, and in this case the quadrature given by (4.4) simply reduces to the usual Gauss-Radau quadrature.

4.1 Discretization by Point Image Charges at Different Locations

Note that

qL1(x)irx=δ1qiarx(xrK)σ1.

Recall that 0 < σ1 < 1. Then using (4.4) with σ = σ1 leads to

rKqL1(x)irxdxm=1MqmL1irxmL1,
(4.6)

where for m = 1, 2, ..., M,

qmL1=2τσ1τδ1ωmL1xmL1aq,xmL1=rK(21smL1)τ.
(4.7)

Here smL1 and ωmL1 are the quadrature points and weights with α = τσ11 − 1.

Similarly, note that

qL2(x)i(1rx1x)=δ2qia(1rx1x)(xrK)(xrK)(1σ2).

Recall that 0 < σ2 < 1. Then, applying (4.4) with σ = (1 − σ2) > 0, we get

rKqL2(x)i(1rx1x)dxm=1MqmL2irxmL2m=1MqmL2ixmL2,
(4.8)

where for m = 1, 2, ..., M,

qmL2=2τ(1σ2)τδ2ωmL2(xmL2rK)xmL2aq,xmL2=rK(21smL2)τ.
(4.9)

Here smL2 and ωmL2 are the quadrature points and weights with α = τ(1 − σ2) − 1.

Also, it can be noted that

qL3(x)i(1rx1xrx2cosθ)=δ3qia(1rx1xrx2cosθ)(xrK)2(xrK)(2σ3).

Recalling that 1 < σ3 < 2, and applying (4.4) with σ = (2 − σ3) > 0, we get

rKqL3(x)i(1rx1xrx2cosθ)dxm=1MqmL3irxmL3m=1MqmL3ixmL3m=1MqmL3rcosθi(xmL3)2,
(4.10)

where for m = 1, 2, ..., M,

qmL3=2τ(2σ3)τδ3ωmL3(xmL3rK)2xmL3aq,xmL3=rK(21smL3)τ.
(4.11)

Here smL3 and ωmL3 are the quadrature points and weights with α = τ(2 − σ3) − 1.

Note that the second summation in the right-hand side of (4.8) and that in the right-hand side of (4.10) are position-independent. Adding them to ΦC0 leaves us with a modified position-independent correction potential

ΦC0=ΦC0m=1MqmL2ixmL2m=1MqmL3ixmL3.
(4.12)

Likewise, adding the last summation in the right-hand side of (4.10) to ΦC1(r) gives us a modified position-dependent correction potential

ΦC1(r)=ΦC1(r)m=1MqmL3rcosθi(xmL3)2=qc1(rrK)cosθ,
(4.13)

where

c1=c11iam=1M2τ(2σ3)τδ3ωmL3(21smL3)τ.

4.2 Discretization by Point Image Charges at the Same Locations

For computational efficiency, we can choose to discretize the three line image charges qL1(x), qL2(x) and qL3(x) by point charges at the same locations. More specifically, one can choose a common parameter σc > 0 to rewrite the line charge qL1(x) as

qL1(x)=δ1qa(xrK)σcσ1(xrK)σc.

Then, for m = 1, 2, ..., M, we have

qmL1=2τσcτδ1ωm(xmrK)σcσ1xmaq.
(4.14)

Here the quadrature points and weights are obtained using α = τσc − 1.

Similarly, the line image charges qL2(x) and qL3(x) are discretized into

qmL2=2τσcτδ2ωm(xmrK)σc+σ2xmaq,
(4.15)
qmL3=2τσcτδ3ωm(xmrK)σc+σ3xmaq.
(4.16)

The parameter σc > 0 is tunable for optimal computational efficiency. For example, depending on the value of u = λa, any of the three natural choices σc = σ1, σc = 1 − σ2 or σc = 2 − σ3 could perform well.

Furthermore, since s1 = −1 and consequently x1 = rK , the classical Kelvin image charge and the first discrete image charge of each image line charge can be combined. In conclusion, in general we have the following multiple discrete image approximation to the reaction potential inside the sphere in terms of the potentials of 3M − 2 point charges (or M charges if a common parameter σc is utilized) and some correction potentials.

ΦRF(r)1i[qKrrK+m=2M(qmL1rxmL1+qmL2rxmL2+qmL3rxmL3)]+ΦC0+ΦC1(r)+ΦC2(r)+ΦC3(r),
(4.17)

where

qK=qK+q1L1+q1L2+q1L3.

For the sake of convenience, we simply write (4.17) as

ΦRF(r)m=1NimqmImirxm+ΦC0+ΦC1(r)+ΦC2(r)+ΦC3(r),
(4.18)

where the summation over m includes the modified Kelvin image charge qK at the corresponding Kelvin image point rK and all discrete image point charges qmL1, qmL2, and qmL3 at xmL1, xmL2, and xmL3, respectively, and Nim is the total number of all point image charges.

5 O(N) Implementation of the Point Image Approximation

The main purpose for the discrete image approximation to the reaction field is to enable us to apply existing fast algorithms, such as the pre-corrected FFT [23] or the fast multipole method (FMM) [3, 4, 10, 12, 13], directly in calculating the electrostatic interactions among N source point charges inside the spherical cavity in O(N log N) or even O(N) operations. In particular, below we give a straightforward O(N) implementation of the discrete image approximation with using the FMM.

For convenience, let riF=(xiF,yiF,ziF), i = 1, 2, ..., N, be N observation points and rjS=(xjS,yjS,zjS), j = 1, 2, ..., N, be the locations of N source charges with charge strengths q1, q2, ..., qN. By linear superposition, the reaction field at an observation point riF, in the case that the sixth-order image approximation (4.18) is employed, becomes

ΦRF(riF)j=1Nm=1Nimqm,jImiriFxm,j+j=1N(ΦC0,j+ΦC1,j(riF)+ΦC2,j(riF)+ΦC3,j(riF)).
(5.1)

Here and in the sequel a quantity with a second subscript j designates a quantity associated with the source charge rjS, such as

rK,j=a2rjS,xm,j=rK,j(21sm)τ.

5.1 O(N) Calculation of the Correction Potentials

Obviously, the position-independent correction potential

j=1NΦC0,j=j=1N[qjia(C0(u)γδ1σ1)m=1Mqm,jL2ixm,jL2m=1Mqm,jL3ixm,jL3]

can be evaluated in O(N) operations. The evaluation of other correction potentials in O(N) operations, however, needs some manipulations due to their position-dependence. First of all, from (4.13) we have

j=1NΦC1,j(riF)=(c1a2j=1NqjrjS)riF,
(5.2)

which in component form can be written as

j=1NΦC1,j(riF)=d1xiF+d2yiF+d3ziF,
(5.3)

where

d1=c1a2j=1NqjxjS,d2=c1a2j=1NqjyjS,d3=c1a2j=1NqjzjS.

Now it becomes clear that the second correction potential can be evaluated in O(N) operations as well. Analogously, the third correction potential can also be evaluated in O(N) operations. In fact, from (3.27) we have

j=1NΦC2,j(riF)=3c22a4j=1Nqj(rjSriF)2(c22a4j=1Nqj(rjS)2)(riF)2,

which can be written as

j=1NΦC2,j(riF)=e0(riF)2+(MriF)riF,
(5.4)

where

e0=c22a4j=1Nqj(rjS)2,M=3c22a4j=1NqjrjSrjS,

and [multiply sign in circle] represents the outer product of two vectors.

In component form, (5.4) can be written as

j=1NΦC2,j(riF)=e0(riF)2+e1(xiF)2+e2(yiF)2+e3(ziF)2+e4xiFyiF+e5xiFziF+e6yiFziF,
(5.5)

where

e1=3c22a4j=1Nqj(xjS)2,e2=3c22a4j=1Nqj(yjS)2,e3=3c22a4j=1Nqj(zjS)2,e4=3c2a4j=1NqjxjSyjS,e5=3c2a4j=1NqjxjSzjS,e6=3c2a4j=1NqjyjSzjS.

Similarly, the fourth correction potential can be evaluated in O(N) operations as well. By expressing the cosine of the angle θ between riF and rjS in terms of their rectangular coordinates, from (3.27) one can arrive at

j=1NΦC3,j(riF)=f1(xiF)3+f2(yiF)3+f3(ziF)3+f4(xiF)2yiF+f5(xiF)2ziF+f6xiF(yiF)2+f7(yiF)2ziF+f8xiF(ziF)2+f9yiF(ziF)2+f10xiFyiFziF+f11xiF(riF)2+f12yiF(riF)2+f13ziF(riF)2,
(5.6)

where

f1=5c32a6j=1Nqj(xjS)3,f2=5c32a6j=1Nqj(yjS)3,f3=5c32a6j=1Nqj(zjS)3,f4=15c32a6j=1Nqj(xjS)2yjS,f5=15c32a6j=1Nqj(xjS)2zjS,f6=15c32a6j=1NqjxjS(yjS)2,f7=15c32a6j=1Nqj(yjS)2zjS,f8=15c32a6j=1NqjxjS(zjS)2,f9=15c32a6j=1NqjyjS(zjS)2,f10=15c3a6j=1NqjxjSyjSzjS,f11=3c32a6j=1NqjxjS(rjS)2,f12=3c32a6j=1NqjyjS(rjS)2,f13=3c32a6j=1NqjzjS(rjS)2.

5.2 O(N) Calculation of the Potentials of the Image Charges

The FMM is known to be extremely efficient in the evaluation of pairwise interactions in large ensembles of particles, such as that included in (5.1)

j=1Nm=1Nimqm,jImiriFxm,j,i=1,2,,N.

For instance, the adaptive FMM of [4] requires O(N) work and breaks even with the direct calculation at about N = 750 for three-digit precision, N = 1500 for six-digit precision, and N = 2500 for nine-digit precision, respectively [11]. Using such an adaptive FMM with O(N) computational complexity, the calculation of the potentials of the discrete image charges for all observation points can be evaluated in O(N) operations in a straightforward way.

In the simplest implementation, such evaluation can be carried out with a single FMM run by including all point image charges qm,jIm into the FMM cube. In the case that the total potential is to be calculated, all original source charges are also included in the FMM box. All charges are taken as acting in a homogeneous medium of the dielectric constant εi.

5.3 Local Expansions for Potential Calculations

The introduced discrete image charges outside the sphere will result in a large computational domain and the image charges are highly nonuniformly distributed, particularly because the image charges of those source charges close to the center of the sphere are far away from the spherical boundary. Therefore, direct application of the FMM by including all image charges in this large computational box is not efficient. Instead, a simple but more efficient way would be to calculate the local expansion due to the far field image charges directly inside the sphere. This way, one can achieve not only a smaller FMM box but also a smaller number of total charges in the FMM box.

More specifically, in practice we could introduce a bigger cut-off sphere of radius κa centered at the origin with κ > 1. The calculation of the potentials inside the original dielectric sphere due to those image charges inside this cut-off sphere is still carried out by the chosen fast method. For all image charges outside this cut-off sphere at (ρl, αl, βl), l = 1, 2, ..., L, with charge strengths ql, l = 1, 2, ..., L, the calculation of the potential at r = (r, θ, [var phi]) inside the dielectric sphere generated by these image charges can be described by a by a local expansion

Φ(r)j=0pk=jjLjkYjk(θ,ϕ)rj,
(5.7)

where Yjk(θ,ϕ) are the spherical harmonics, and Ljk are the local expansion coefficients given by

Ljk=l=1Lq^lYjk(αl,βl)ρlj+1.
(5.8)

Furthermore, for any p ≥ 1,

Φ(r)j=0pk=jjLjkYjk(θ,ϕ)rj(l=1Lq^lκar)(rκa)p+1.
(5.9)

6 O(N) Calculation of the Force Field

In practice, most time in molecular dynamics simulations the electric force field needs to be calculated. Although force equations are not difficult to derive, we would like to include them here for completeness. In general, electric forces are computed by taking gradients of electric potentials. Therefore, the electric force exerted on Particle i at the position riF is expressed as

f(riF)=qirΦ(riF).

In the case that the reaction potential field is approximated by the sixth-order image approximation (4.18), the image approximation for the electric force field becomes (note that rΦC0,j=0)

f(riF)=qij=1Nm=1Nimqm,jIm(riFxm,j)iriFxm,j3qirj=1N(ΦC1,j(riF)+ΦC2,j(riF)+ΦC3,j(riF)).

First of all, O(N) calculation of the gradients of the correction potentials can be derived directly from the O(N) calculation of the corresponding correction potentials. For examples, from (5.2) and (5.4) we can obtain

rj=1NΦC1,j(riF)=c1a2j=1NqjrjS,

and

rj=1NΦC2,j(riF)=2(e0I+M)riF=(2(e0+e1)e4e5e42(e0+e2)e6e5e62(e0+e3))riF.

On the other hand, the FMM can be used to calculate

qij=1Nm=1Nimqm,jIm(riFxm,j)iriFxm,j3,i=1,2,,N

in O(N) complexity.

Moreover, the force field inside the sphere due to the far field image charges can also be described in local expansions. Local expansion for force calculations in the FMM can be found in [7, 18]. When using the local expansion (5.7) to calculate the far-field electric potential due to all image charges outside the cut-off sphere, the corresponding far-field force f(r) exerted on a particle q at r = (r, θ, [var phi]) inside the dielectric sphere can also be described by local expansions. Passing the details, we have

fx(r)=qxΦ(r)=qRe(H2H3),
(6.1)
fy(r)=qyΦ(r)=qIm(H2+H3),
(6.2)
fz(r)=qzΦ(r)=q(H0+2Re(H1)),
(6.3)

where Re(...) and Im(...) represent the real and the imaginary parts of a complex number, respectively, and

H0=j=1pjLj0Pj1(cosθ)rj1,H1=j=1pk=1j1(j+k)CjkLjkeikϕPj1k(cosθ)rj1,H2=j=1pLj0eiϕPj11(cosθ)rj1+j=1pk=1jCjkLjkei(k+1)ϕPj1k+1(cosθ)rj1,H3=j=1pCj1Lj1Pj1(cosθ)rj1+j=1pk=2jBjkCjkLjkei(k1)ϕPj1k1(cosθ)rj1,

and

Bjk=(j+k)(j+k1),Cjk=(jk)!(j+k)!.

7 Numerical Examples

For illustration purpose, a unit dielectric sphere is used. The dielectric constants of the sphere and its surrounding medium are assumed to be εi = 2 (normally 1, 2 or 4) and εo = 80 (the dielectric constant of water), respectively. Unless otherwise specified, the results obtained by the direct series expansion with 400 terms are treated as the exact reaction fields to calculate the errors of the sixth-order image approximation.

7.1 Accuracy vs the Ionic Strength

We begin by considering a single point charge located on the x-axis inside the sphere at a distance rS = 0.5 or rS = 0.95 from the center of the sphere. Different σ values are used in discretizing the underlying image line charges, but for simplicity, we always choose τ = 1/σ and M = 40 so that the same Gauss-Radau quadrature points and weights sm and ωm are involved. For each selected value of u = λa, we calculate the relative error of the image approximation in the reaction field, respectively, at 10,000 observation points uniformly distributed (under the polar coordinates) within the unit disk in the plane y = 0. The maximal relative error || E || at these observation points for various u values and the corresponding order of convergence are shown in Tables 1 and and2.2. For sake of comparison, the results obtained using the improved fourth order image approximation are also included. As can be observed, the results clearly demonstrate the O(u6) convergence rate of the sixth-order image approximation.

Table 1
Convergence rate of the proposed image approximation (rS = 0.5)
Table 2
Convergence rate of the proposed image approximation (rS = 0.95)

7.2 Accuracy vs the Number of Discrete Image Charges

One natural concern with the proposed discrete image approximation is the final number of point image charges required to achieve certain order of degree of accuracy. For a desired accuracy, this number depends on the locations of the source charge and the observation point. It should be small if compared to the number of terms needed to achieve the same degree of accuracy in the direct series expansion to make the image approximation useful in practice.

In this test, the source location is fixed at rS = 0.95. For each selected value of u = λa ranging from 0.05 to 0.9, we approximate the reaction fields at the same 10,000 points within the sphere by the sixth-order image approximation with several different numbers of point image charges. Figure 2 shows the maximum relative errors of the point image approximation using M = 2, M = 6 and M = 10, without a common σ being used. For the sake of comparison, the corresponding error analysis results for the improved fourth-order image approximation are also plotted.

Fig. 2
Maximum relative errors in the ionic solvent induced reaction field due to a source charge inside the unit sphere at distance rS = 0.95 from the spherical center. No common σ is used

As can be seen, when M = 2, the sixth-order multiple discrete image approximation can achieve 10−3 accuracy. However, as M is small, the accuracy of the sixth-order image approximation is no better than that of the improved fourth-order image approximation. This is because when M is small, the numerical error stemming from the approximation of the line image charges by point image charges dominates. Nevertheless, when M is relatively large, the advantage of the sixth-order image approximation becomes evident in the sense of accuracy. For example, as shown in Fig. 2, the sixth-order image approximation is clearly shown to be more accurate than the fourth-order one when M = 10 is used, particularly for cases with large u values. It should be pointed out that, since different line images are discretized by point images at different locations, the number of total point images in the sixth-order method is 3M − 2, while that in the fourth-order method is only 2M − 1 in this particular test.

Similar error analysis results with a common σ being used are displayed in Fig. 3. Note that in this case the sixth-order and the fourth-order methods both use M point images for a given M value. As shown in Fig. 3, for low accuracy (such as 10−2), for all u values in the range of [0.05, 0.9], both the sixth-order and the fourth-order methods need 2 point images and thus the same computational cost to achieve the required accuracy. For high accuracy (such as 10−5), on the one hand, when the u value is small, the two methods still have comparable performance, namely, they need to use the same number of point images to achieve the same accuracy. On the other hand, when the u value is large such as 0.9, the error of the sixth-order method using 10 point images can be less than 4 × 10−5, but that of the fourth-order method using the same number of point images is around 2 × 10−4, implying that if it is possible for the fourth-order method to realize the same 4 × 10−5 accuracy, it must use more than 10 image charges and consequently more computational cost. Actually, as shown by Fig. 4, the best accuracy the fourth-order method can achieve is around 1.2 × 10−3 for u = 0.9 and 1.5 × 10−4 for u = 0.5, respectively. Therefore, when the u value is large, high-accuracy can only be accomplished by using the sixth-order method.

Fig. 3
Maximum relative errors in the ionic solvent induced reaction field due to a source charge inside the unit sphere at distance rS = 0.95 from the spherical center. A common σ is used
Fig. 4
Maximum relative errors in the ionic solvent induced reaction field due to a source charge inside the unit sphere at distance rS = 0.95 from the spherical center. A common σ is used

7.3 Computational Complexity of the Point Image Approximation

To demonstrate the O(N) computational complexity of the O(N) implementation of the point image approximation as described in Sects. 5 and 6, the sixth-order image approximation has been implemented with using the free software KIFMM, developed by L. Ying using a so-called kernel-independent adaptive fast multipole method [24]. The experiments are carried out on a Linux-based Dell OptiPlex 745 workstation with a CPU clock rate of 3 GHz and a memory of 4 GB, using GNU Fortran 3.4.6 and C++ 4.1.2 compilers. In Tables 3 and and4,4, timing results of the potential and the force evaluations are reported and compared with those obtained without using the FMM. The Gauss-Radau quadrature with M = 2 and a common parameter σc = σ2 are used to construct discrete image charges, so for each source charge, there are two discrete image charges. As can be seen, the timing scales as O(N2) without using the FMM, and linearly with using the FMM. In the tables, N denotes the number of total source charges, either uniformly or non-uniformly distributed inside the unit sphere, and NC the number of total point charges included in the FMM box, respectively.

Table 3
Timing results for potential calculation
Table 4
Timing results for force calculation

7.4 A Simple Application of the Method

Finally, the proposed sixth-order method is tested through calculation of solvation effects. A spherical cavity of radius 15 Å contains 333 TIP3P water molecules and is immersed in an ionic solvent. To calculate the total solvation energy of these water molecules, εi = 1 and εo = 80 are used in this test. Here the solvation energy of a collection of N charges qi located at ri, i = 1, 2, ..., N, is defined as

Usol=12i=1NqiΦRF(ri).

Figure 5 shows the relative error in the solvation energy obtained with the fourth-order and the sixth-order image approximations by using M = 6 and M = 10 point image charges, respectively, where the exact solvation energy is calculated by the series solution with 400 terms. As can be seen, the sixth-order accuracy of the present method has a noticeable effect on the solvation energy.

Fig. 5
Relative error in the solvation energy calculation of water molecules in ionic solvents by two image approximation methods

8 Conclusions

In this paper, we have presented a sixth-order image approximation to the reaction field due to a point charge inside a dielectric sphere immersed in an ionic solvent for small values of u = λa (λ—the inverse Debye screening length of the ionic solvent, a—the radius of the dielectric sphere). O(N) implementations of the image approximation in the electric potential and force field evaluations have been described. Numerical results have demonstrated the sixth-order convergence rate of the image approximation as well as the O(N) complexity of the O(N) implementations of the image approximation. Numerical results have also demonstrated that, compared to the previous fourth-order image approximation, the proposed sixth-order image approximation is more accurate for cases with large u values when a relatively large M value is used.

Acknowledgements

The authors would like to thank the support by the National Institutes of Health (grant number: 1R01GM083600-02) for the work reported in this paper. S. Deng was supported partially by funds provided by the University of North Carolina at Charlotte as well.

References

1. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover; New York: 1972.
2. Cai W, Deng S, Jacobs D. Extending the fast multipole method to charges inside or outside a dielectric sphere. J. Comput. Phys. 2007;223:846–864.
3. Carrier J, Greengard LF, Rokhlin V. A fast adaptive multipole algorithm for particle simulations. SIAM J. Sci. Stat. Comput. 1988;9:669–686.
4. Cheng H, Greengard LF, Roklin V. A fast adaptive multipole algorithm in three dimensions. J. Comput. Phys. 1999;155:468–498.
5. Deng S, Cai W. Discrete image approximations of ionic solvent induced reaction field to charges. Commun. Comput. Phys. 2007;2:1007–1026.
6. Deng S, Cai W. Extending the fast multipole method for charges inside a dielectric sphere in an ionic solvent: High-order image approximations for reaction fields. J. Comput. Phys. 2007;227:1246–1266. [PMC free article] [PubMed]
7. Esselink L. A comparison of algorithms for long-range interactions. Comput. Phys. Commun. 1995;87:375–395.
8. Gautschi W. Algorithm 726: ORTHPOL—a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans. Math. Softw. 1994;20:21–62.
9. Gradshteyn IS, Ayzhik IM. Table of Integrals, Series, and Products. Academic Press; Boston: 1994.
10. Greengard LF. The Rapid Evaluation of Potential Fields in Particle Systems. MIT; Cambridge: 1987.
11. Greengard LF, Huang J. A new version of the fast multipole method for screened Coulomb interactions in three dimensions. J. Comput. Phys. 2002;180:642–658.
12. Greengard LF, Rokhlin V. A fast algorithm for particle simulations. J. Comput. Phys. 1987;73:325–348.
13. Greengard LF, Rokhlin V. A new version of the fast multipole method for the Laplace equation in three dimensions. Acta Numer. 1997;6:229–269.
14. Juffer A, Botta EFF, van Keulen BAM, van der Ploeg A, Berendsen HJC. The electric potential of a macromolecule in a solvent: A fundamental approach. J. Comput. Phys. 1991;97:144–171.
15. Lee MS, Olson MA. Evaluation of Poisson solvation models using a hybrid explicit/implicit solvent method. J. Phys. Chem. B. 2005;109:5223–5236. [PubMed]
16. Lee MS, Salsbury FR, Jr., Olson MA. An efficient hybrid explicit/implicit solvent method for bio-molecular simulations. J. Comput. Chem. 2004;25:1967–1978. [PubMed]
17. Lindell IV. Electrostatic image theory for the dielectric sphere. Radio Sci. 1992;27:1–8.
18. Lu B, Cheng X, McCammon JA. New-version-fast-multipole-method accelerated electrostatic calculations in biomolecular systems. J. Comput. Phys. 2007;226:1348–1366. [PMC free article] [PubMed]
19. Morse PM, Feshbach H. Methods of Theoretical Physics. McGraw-Hill; New York: 1953.
20. Neumann C. Hydrodynamische Untersuchen nebst einem Anhang uber die Probleme der Elecktrostatik und der magnetischen Induktion. Teubner; Leipzig: 1883. pp. 279–282.
21. Norris WT. Charge images in a dielectric sphere. IEE Proc. Sci. Meas. Technol. 1995;142:142–150.
22. Okur A, Simmerling C. Spellmeyer D, editor. Hybrid explicit/implicit solvation methods. Annu. Rep. Comput. Chem. 2006;2 chap. 6.
23. Phillips JR, White JK. A precorrected-FFT method for electrostatic analysis of complicated 3-D structures. IEEE Trans. Comput. Aided Des. Integr. Circ. Syst. 1997;10:1059–1072.
24. Ying L, Biros G, Zorin D. A kernel-independent adaptive fast multipole algorithm in two and three dimensions. J. Comput. Phys. 2004;196:591–626.