PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Phys Commun. Author manuscript; available in PMC 2010 December 6.
Published in final edited form as:
Comput Phys Commun. 2008 February 1; 178(3): 171–185.
doi:  10.1016/j.cpc.2007.08.015
PMCID: PMC2997574
NIHMSID: NIHMS243075

New Versions of Image Approximations to the Ionic Solvent Induced Reaction Field

Abstract

A recent article by Deng and Cai (Extending the fast multipole method for charges inside a dielectric sphere in an ionic solvent: High-order image approximations for reaction fields, to appear in J. Comput. Phys.) introduced two 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 this Kelvin image point along the radial direction to infinity, with one decaying to zero and the other growing to infinity. In this paper, alternative versions of the fourth-order image approximations are presented, using the same point charge but three different line charges, all decaying to zero along the radial direction. Similar discussions on how to approximate the line charges by discrete image charges and how to apply the resulting multiple discrete image approximations together with the fast multipole method are also included.

Keywords: Method of images, reaction field, ionic solvent, hybrid solvent model

1 Introduction

Electrostatic interactions play a major role in determining the structure, dynamics and function of biological macromolecules; as such, they remain as major objects of theoretical and computational studies of macromolecules. Electrostatic interactions are long-range, and strongly dependent on the solvent and the ions surrounding the biomolecule under study. When modeling biological systems numerically, it has been challenging, however, to account for such solvent environment in a manner that is computationally efficient and physically accurate at the same time.

Explicit representation of solvent molecules [13] offers a detailed and accurate description of a biological macromolecule, yet all-atom simulations are expensive to perform due to the long-range nature of the electric forces. As a result, large explicitly solvated systems typically cannot be simulated for biologically relevant timescales. Alternatively, in implicit solvent models [4,5], an aqueous solvent is modeled as a continuum medium with a large dielectric constant (60 to 85) outside the macromolecule. The macromolecule atoms themselves are explicitly modeled with assigned partial charges embedded in a dissimilar continuum medium of a low dielectric constant (1 to 4) inside the macromolecule volume. The neglect of explicit solvent molecules can significantly reduce the computational cost. Despite the marked success of implicit solvent models, however, they also have fundamental limitations due to the fact that the important atomic details of how the solvent molecules interact with the surface of the macromolecule are ignored. For these reasons, there has been considerable recent interest in developing hybrid explicit/implicit solvent models [68], in which the macromolecule together with a boundary layer or shell of the solvent molecules are considered explicitly within a cavity, and outside the cavity, the solvent is treated implicitly as a dielectric continuum.

Electric charges within the cavity will polarize the surrounding solvent medium, which in turn makes a contribution, called the reaction field, to the electric field throughout the cavity. The electric potential inside the cavity is thus expressed as Φ = ΦS + ΦRF where ΦS is the potential due to direct Coulomb interactions between source charges within the cavity, and ΦRF is the reaction field. Fast and accurate calculation of such a reaction field will have a far-reaching impact on computational simulations for chemical and biological systems involving electrostatic interactions within a solvent. In case of a spherical cavity, a popular approach 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 [9] and references therein. For an ionic solvent, in [10], by assuming that the ionic strength of the solvent is low enough so the product of the inverse Debye screening length of the solvent and the radius of the spherical cavity is less than one, a first- and a second-order image approximation to the ionic solvent induced reaction field were developed. Both approaches employ the same point and line charges as those obtained for the pure water solvent, while the ionic strength effect is included in additional correction terms. Two fourth-order image approximations were further developed in [11], using a point charge at the classical Kelvin image point and two line charges that extend from this Kelvin image point along the radial direction to infinity, with one decaying to zero and the other growing to infinity. In this paper, new versions of the fourth-order image approximations to the ionic solvent induced reaction field shall be developed so that only line charges decaying to zero are utilized.

The paper is organized as follows. In Section 2, after briefly reviewing the series solution to the ionic solvent induced reaction field due to a point charge inside a dielectric sphere, we present the new versions of the fourth-order image approximations to such the reaction field. Then in Section 3, the discretization of the line charges by point image charges is summarized. Numerical results are given next in Section 4 to validate the convergence property and investigate the efficiency of the new versions of the fourth-order image approximations. In addition, how to apply the proposed multiple discrete image approximations together with the fast multipole methods is discussed in Appendix A.

2 Image approximations to ionic solvent induced reaction fields

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

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

Inside the sphere, the electrostatic potential Φ is given by the solution to the Poisson's equation

(εiΦ(r))=qδ(|rrs|),
(1)

where δ(r) is the Dirac delta function. Outside the sphere, on the other hand, by assuming that the mobile ion concentration in the ionic solvent is given by the Debye-Hückel theory, for a solvent of low ionic strength, the linearized Poisson-Boltzmann equation [12,13]

2Φ(r)λ2Φ(r)=0
(2)

can be used to approximate the screened Coulomb potential in the solution. Here, λ is the inverse Debye screening length defined by

λ2=8πNAe2I1000εokBT,

where NA is Avogadro's number, e the electron charge, kB the Boltzmann constant, T the absolute temperature, and I the ionic strength of the solvent.

Using the classical electrostatic theory, the reaction field of the spherical dielectric can be solved analytically [10]. More precisely, with respect to a spherical coordinate system (r, θ, ϕ) 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 an observation point r=(r, θ, ϕ) inside the sphere is

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

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

An=qεia1rKnεi(n+1)kn(u)+εoukn(u)εinkn(u)εoukn(u),n0.
(4)

Here u=λa, rK=a2/rS with rK=(rK, 0, 0) denoting the conventional Kelvin image point, and kn(r) are the modified spherical Hankel functions defined as [14,15]

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

In theory, any desired degree of accuracy can be obtained using the direct series expansion (3) to calculate the ionic solvent induced reaction field. 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, however, the convergence by the series expansion is slow due to r/rK=rrS/a2 ≈ 1, requiring a great number of terms in the series expansion to achieve high accuracy in the reaction field. Figure 2 plots the reaction field on a disk inside a spherical cavity of radius 1 in the case that rS=0.8, λ=0.8, εi=2, and εo=80.

Fig. 2
Illustration of the reaction field inside a spherical cavity of radius 1 in the case that rS=0.8, λ=0.8, εi=2, and εo=80.

Let us now consider the problem of finding images outside the spherical region giving rise to the reaction potential inside the sphere. As mentioned earlier, in [11], by assuming that the ionic strength of the solvent is low enough so the product of the inverse Debye screening length and the radius of the spherical cavity is less than one (u=λa < 1), two fourth-order image approximations (in terms of u=λa) to the ionic solvent induced reaction field were developed. Both approximations employ a point charge and two line charges with one decaying to zero and the other growing to infinity. Next, we shall follow similar methodology to develop alternative versions of the fourth-order image approximations. As discussed in detail in [11], 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 [16,17] to be meaningful.

In [11], the construction of the fourth-order image approximations is based on the fact that the modified spherical Hankel functions and their derivatives can be expanded, respectively, in terms of 1/r as

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

and

kn(r)=π(2n2)!(n1)!(4(2n1)(n+1)(2r)n+2n12(2r)n)+O(1rn2),n2.

Correspondingly, we have

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

Inserting this approximation in (4) leads to

An=qεia1rKna1+a2n+a3n2b1+b2n+b3n2+O(u4),n2,

where we denote

a1=(εi+εo)u22(εiεo),a2=(εiεo)(2u2),a3=4(εiεo),b1=εo(2u2),b2=(εi+εo)u22(εiεo),b3=4(εi+εo).

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

An=qεia1rKn(γ+α1+α2nβ1+β2n+n2)+O(u4),n2,

where

γ=εiεoεi+εo,α1=γ(a1a3b1b3),α2=γ(a2a3b2b3),β1=b1b3,β2=b2b3.

Note that β1 < 0 when 0 ≤ u < 1. Then using partial fractions gives us

An=qεia1rKn(γ+δ1n+σ1+δ2nσ2)+O(u4),n2,
(7)

where

σ1=β224β1+β22,δ1=α2σ1α1σ1+σ2,σ2=β224β1β22,δ2=α2σ2+α1σ1+σ2.

It can be shown that 0 < σ1 < 1 and 0 < σ2 < 1. Moreover, 1/2 < σ1 when εi < εo/5, and σ2 < 1/2 when εi < εo, σ2=1/2 when εi=εo, and σ2 > 1/2 when εi > εo. Therefore, for our target applications in hybrid solvent models for biomolecular simulations where the typical dielectric constants are 1 to 4 for εi and 60 to 85 for εo, respectively, we have 1/2 < σ1 < 1 and 0 < σ2 < 1/2.

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

A0=C0(u)qεia,A1=C1(u)qεia1rK,

where

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

Inserting now the approximation of An given by (7) into (3), the reaction field inside the sphere is taken as

ΦRF(r)=A0+A1rcosθ+S1+S2+S3+O(u4),
(8)

with S1, S2 and S3 representing three series, defined respectively by

S1=γqεian=2(rrK)nPn(cosθ),S2=δ1qεian=21n+σ1(rrK)nPn(cosθ),S3=δ2qεian=21nσ2(rrK)nPn(cosθ).

To find image representations of the three series, we first recall that the Coulomb potential ΦS(r), the potential at r due to a point charge at rS, can be expanded in terms of the Legendre polynomials of cosθ as follows [18].

ΦS(r)=qεi|rrS|={qεirn=0(rSr)nPn(cosθ),rSra,qεirSn=0(rrS)nPn(cosθ),0rrS.
(9)

Then the first series S1 can be written as

S1=γqεirKarSn=0(rrK)nPn(cosθ)γqεia(1+rrKcosθ),

where the first part of the right-hand side is exactly the expansion given by (9) for a point charge of magnitude qK=γaq/rS outside the sphere at the classical Kelvin image point rK=(rK, 0, 0), namely,

S1=qKεi|rrK|γqεia(1+rrKcosθ).
(10)

Next, using the integral identity

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

which is valid for all n ≥ 0 when σ > 0, the second series S2 can be written as

S2=rK[qL1(x)εixn=0(rx)nPn(cosθ)]dxδ1qεia(1σ1+11+σ1rrKcosθ),

where

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

Note that the integrand of the above integral again is the expansion given by (9) 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 that extends from the classical Kelvin image point rK along the radial direction to infinity. Thus, the second series S2 becomes

S2=rKqL1(x)εi|rx|dxδ1qεia(1σ1+11+σ1rrKcosθ).
(13)

However, an image representation similar to (13), in which the line charge decays to zero along the radial direction, cannot be obtained for the third series S3 by following the same methodology. Instead, in [11], this series is represented using a line charge that grows to infinity. In this paper, we shall develop a new approach in which only line charges that decay to zero are involved. To this end, we first use the recursion formula

Pn(x)=2n1nxPn1(x)n1nPn2(x),n2

to rewrite the third series as S3=S4S5 with the two series S4 and S5 being defined as

S4=qεian=2(2n1)δ2n(nσ2)(rrK)ncosθPn1(cosθ),S5=qεian=2(n1)δ2n(nσ2)(rrK)nPn2(cosθ).

Now, using partial fractions again leads to

(2n1)δ2n(nσ2)=δ2,1n+δ2,2nσ2

and

(n1)δ2n(nσ2)=δ2,1n+δ2,3nσ2,

where we denote

δ2,1=δ2σ2,δ2,2=(2σ21)δ2σ2,δ2,3=(σ21)δ2σ2.

Then the series S4 can be written as S4=S6 + S7, where

S6=qεian=2δ2,1n(rrK)ncosθPn1(cosθ)=qεia(rcosθrK)m=1δ2,1m+1(rrK)mPm(cosθ),

and

S7=qεian=2δ2,2nσ2(rrK)ncosθPn1(cosθ)=qεia(rcosθrK)m=1δ2,2m+(1σ2)(rrK)mPm(cosθ).

Now applying the integral identity (11) with σ=1 and σ=1 − σ2, respectively, gives us

S6=(rKqL2,1(x)εi|rx|dxδ2,1qεia)(rcosθrK),
(14)

and

S7=(rKqL2,2(x)εi|rx|dxδ2,2qεia(1σ2))(rcosθrK).
(15)

Here,

qL2,1(x)=δ2,1qa(xrK)1,rKx,

and

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

Combining (14) and (15) yields

S4=(rKqL2(x)εi|rx|dxqεiaδ21σ2)(rcosθrK),
(16)

where

qL2(x)=qL2,1(x)+qL2,2(x)=δ2,1qa[(xrK)1+(2σ21)(xrK)(1σ2)].

Analogously, the series S5 can be written as S5=S8 + S9, where

S8=qεian=2δ2,1n(rrK)nPn2(cosθ)=qεia(rrK)2m=0δ2,1m+2(rrK)mPm(cosθ),

and

S9=qεian=2δ2,3nσ2(rrK)nPn2(cosθ)=qεia(rrK)2m=0δ2,3m+(2σ2)(rrK)mPm(cosθ).

Applying the integral identity (11) with σ=2 and σ=2 − σ2, respectively, gives us

S8=rKqL3,1(x)εi|rx|dx(rrK)2,
(17)

and

S9=rKqL3,2(x)εi|rx|dx(rrK)2,
(18)

where the functions qL3,1(x) and qL3,2(x) are

qL3,1(x)=δ2,1qa(xrK)2,rKx,

and

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

Combining (17) and (18) yields

S5=rKqL3(x)εi|rx|dx(rrK)2,
(19)

where

qL3(x)=qL3,1(x)+qL3,2(x)=δ2,1qa[(xrK)2+(σ21)(xrK)(2σ2)].

Combining (10), (13), (16), and (19), we finally have a fourth-order image approximation to the ionic solvent induced reaction field as follows.

ΦRF(r)=qKεi|rrK|+rKqL1(x)εi|rx|dx+rKqL2(x)εi|rx|dx(rcosθrK)rKqL3(x)εi|rx|dx(rrK)2+ΦC1+ΦC2(r)+O(u4),
(20)

where ΦC1 is a position-independent correction potential given by

ΦC1=qεia(C0(u)γδ1σ1),
(21)

and on the other hand, ΦC2(r) defines a position-dependent correction potential

ΦC2(r)=qεia(C1(u)γδ11+σ1δ21σ2)rrKcosθ.
(22)

Likewise, qL2(x) and qL3(x) can be regarded as continuous line charges extending from the classical Kelvin image point rK along the radial direction to infinity. Asymptotically, all three line charges qL1(x), qL2(x) and qL3(x) decay to zero faster than x−1/2 at infinity. More precisely, asymptotically qL1(x) decays to zero as rapidly as xσ1, qL2(x) as x−(1−σ2), and qL3(x) as x−(2−σ2), respectively. Figure 3 shows the density distributions of the image line charges qL1(x), qL2(x) and qL3(x) for a source point charge q=−1 at rs=0.8 on the x-axis in the case of εi=2, εo=80, a=1, and λ=0.8.

Fig. 3
Density distributions of the image line charges qL1(x), qL2(x) and qL3(x) for a source point charge q=−1 at rs=0.8 in the case of εi=2, εo=80, a=1, and λ=0.8.

As pointed out in [11], the accuracy of the fourth-order image approximation (20) can be improved by including more correction terms. In particular, for n=2, from the exact expression of k2(r), in fact we can arrive at

A2=C2(u)qεia1rK2,

where

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

Now, instead of (8), we express the reaction field inside the sphere as

ΦRF(r)=A0+A1rcosθ+A2r2P2(cosθ)+n=3qεia1rKn(γ+δ1n+σ1+δ2nσ2)rnPn(cosθ)+O(u4).
(23)

Note that the only difference between (8) and (23) is in the calculation of the r2P2(cos θ)-terms. Therefore, the fourth-order image approximation (20) can be improved simply by including another position-dependent correction potential, resulting in an improved fourth-order image approximation as follows.

ΦRF(r)=qKεi|rrK|+rKqL1(x)εi|rx|dx+rKqL2(x)εi|rx|dx(rrK)cosθrKqL3(x)εi|rx|dx(rrK)2+ΦC1+ΦC2(r)+ΦC3(r)+O(u4),
(24)

where the second position-dependent correction potential is

ΦC3(r)=qεia(C2(u)γδ12+σ1δ22σ2)(rrK)2P2(cosθ).
(25)

3 Fourth-order multiple discrete image approximations

Now let us consider how to represent approximately each image line charge involved in the image approximations (20) and (24) by a set of discrete image charges. This can be achieved by approximating the involved integrals by numerical quadratures. Without losing any generality, let us consider

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

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,
(27)

where α=τσ − 1 and

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

Note that α > −1 because σ > 0 and τ > 0. Also s=−1 corresponds to the Kelvin image point x=rK. Therefore, we can choose either Jacobi-Gauss or Jacobi-Gauss-Radau quadrature to approximate the integral in (27). More precisely, let sm, ωm, m=1, 2, …, M, be the Jacobi-Gauss or Jacobi-Gauss-Radau points and weights on the interval [−1, 1] with α=τσ − 1 and β=0 (s1=−1 if Jacobi-Gauss-Radau quadrature is used.), which can be obtained with the program ORTHPOL [21]. Then, we have

I2τστm=1Mωmxmf(xm),
(28)

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

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

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 (28) simply reduces to the usual Gauss or Gauss-Radau quadrature.

3.1 Discretization of qL1(x)

Note that

rKqL1(x)εi|rx|dx=rKδ1qεia|rx|(xrK)σ1dx.

Recall that 0 < σ1 < 1, and furthermore, when εi < εo/5, we have σ1 > 1/2. Then using (28) with σ=σ1 leads to

rKqL1(x)εi|rx|dxm=1MqmL1εi|rxm|,
(30)

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

qmL1=2τσ1τδ1ωmxmaq.
(31)

3.2 Discretization of qL2(x)

Discretization of qL2 (x) can be carried out by discretizing qL2,1(x) and qL2,2(x) separately in the same way as how qL1(x) is discretized, which in general will introduce a total of 2M discrete charges even when the same τ value is used to discretize qL2,1(x) and qL2,2(x). In order to obtain an approach that uses only M discrete charges, we reformulate qL2(x) as

qL2(x)=δ2,1qa[1+(2σ21)(xrK)σ2](xrK)1

or

qL2(x)=δ2,1qa[(xrK)σ2+(2σ21)](xrK)(1σ2).

Then using (28) with σ=1 or σ=(1 − σ2) > 1/2 leads to

rKqL2(x)εi|rx|dxm=1MqmL2εi|rxm|,
(32)

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

qmL2=2ττδ2,1ωm[1+(2σ21)(xmrK)σ2]xmaq,
(33)

or

qmL2=2τ(1σ2)τδ2,1ωm[(xmrK)σ2+(2σ21)]xmaq.
(34)

3.3 Discretization of qL3(x)

Discretization of qL3(x) can be carried out in the same way as how qL2(x) is discretized. For instance, to obtain an approach that uses only M discrete charges, we reformulate qL3(x) as

qL3(x)=δ2,1qa[1+(σ21)(xrK)σ2](xrK)2,

or

qL3(x)=δ2,1qa[(xrK)σ2+(σ21)](xrK)(2σ2).

Then using (28) with σ=2 or σ=(2 − σ2) > 3/2 leads to

rKqL3(x)εi|rx|dxm=1MqmL3εi|rxm|,
(35)

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

qmL3=22ττδ2,1ωm[1+(σ21)(xmrK)σ2]xmaq,
(36)

or

qmL3=2τ(2σ2)τδ2,1ωm[(xmrK)σ2+(σ21)]xmaq.
(37)

In conclusion, in general we have the following fourth-order multiple discrete image approximations to the reaction potential inside the sphere.

  1. A fourth-order multiple discrete image approximation:
    ΦRF(r)qKεi|rrK|+m=1MqmL1εi|rxm|+(m=1MqmL2εi|rxm|)(rrK)cosθ(m=1MqmL3εi|rxm|)(rrK)2ΦC1+ΦC2(r).
    (38)
  2. An improved fourth-order multiple discrete image approximation:
    ΦRF(r)qKεi|rrK|+m=1MqmL1εi|rxm|+(m=1MqmL2εi|rxm|)(rrK)cosθ(m=1MqmL3εi|rxm|)(rrK)2+ΦC1+ΦC2(r)+ΦC3(r).
    (39)

4 Numerical results

For demonstration purpose, let us consider a unit dielectric sphere. The dielectric constants of the sphere and its surrounding medium are assumed to be εi=2 and εo=80, respectively. The results obtained by the direct series expansion with 400 terms are treated as the exact reaction fields to calculate the errors of various image approximations. Unless otherwise specified, (31) with τ=1/σ1, (34) with τ=1/(1 − σ2), and (37) with τ=1/(2 − σ2) are used to approximate the image line charges qL1(x), qL2(x) and qL3(x), respectively, and for simplicity, we always choose the same M value so that the same Gauss quadrature points and weights sm and ωm are involved.

We begin by verifying the convergence property of the proposed line image approximations. To this end, let us consider a single point charge located on the x-axis inside the sphere at a distance rS=0.95 from the center of the sphere. A large M value, M=50, is used to discretize the line charges so that the error arising in the discretization of the image line charges by discrete charges is negligible compared to the error arising in the approximation of the reaction field by the point and image line charges. For each selected value of u=λa, we calculate the relative error of the image approximations in the reaction field, respectively, at 10, 000 observation points uniformly distributed (under the polar coordinates) within the unit disk that contains the x-axis. The maximal relative error ‖E‖ at the 10, 000 observation points for various u values and the corresponding order of convergence are shown in Table 1. As can be observed, the results clearly demonstrate the O(u4) convergence property of the new versions of the fourth-order image approximations.

Table 1
Convergence property of the proposed image approximations.

One natural concern with the proposed multiple discrete image approximations is the final number of discrete image charges required to achieve certain order of degree of accuracy. For a desired accuracy, this number depends on both the locations of the source charge and the ionic strength. 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 approximations useful in the practice.

In this test, we shall confine ourselves to the improved fourth-order image approximation, and we consider four different source locations rS=0.8, 0.9, 0.95, and 0.98 as well as different u=λa values ranging from 0.05 to 0.95. For each selected source position and u value, we approximate the reaction fields at the same 10,000 points within the sphere first by the direct series expansion with various numbers of terms, and then by the improved fourth-order image approximation with various numbers of discrete charges. The results are compared to the exact ones obtained by the direct series expansion with 400 terms, and the maximal relative errors are plotted in Fig. 4.

Fig. 4
Accuracy of the improved fourth-order image approximation to the ionic solvent induced reaction field due to a source charge inside the unit sphere at distance rS from the center. The total number of discrete image charges is 3M + 1.

First of all, as can be seen, the approximation error increases as the source charge moves to the spherical boundary. Second, as pointed out earlier, in the case that the source charge is close to the spherical boundary, the Kelvin image point rK is close to the boundary as well. Therefore, when calculating the reaction field at an observation point also close to the spherical boundary, the convergence by the direct series expansion

ΦRF(r)=qεian=0εi(n+1)kn(u)+εoukn(u)εinkn(u)εoukn(u)(rrK)nPn(cosθ)

shall be slow due to r/rK=rrS/a2 ≈ 1, requiring a great number of terms to achieve high accuracy in the reaction field.

As indicated, to achieve an accuracy with the error being less than 10−2, around 20, 40, 90 and 225 terms have to be included in the direct series expansion for the cases of rS=0.8, 0.9, 0.95 and 0.98, respectively. To achieve the same accuracy, however, for all source locations and all u ≤ 0.95 only 4 discrete charges including the point charge at the Kelvin image point (M=1) are needed in the multiple discrete image approximation. Similarly, as can be seen, to achieve an accuracy with the error being less than 10−3, around 30, 65, 130 and 320 terms have to be included in the direct series expansion for the cases of rS=0.8, 0.9, 0.95 and 0.98, respectively. On the other hand, to achieve the same accuracy, for relatively small and relatively large u values (u < 0.5 and u > 0.5), 7 and 10 discrete charges (M=2 and M=3) are needed in the multiple discrete image approximation, respectively. Therefore, to calculate reaction fields at points close to the spherical boundary due to source charges also close to the boundary, the improved fourth-order image approximation is clearly much more efficient than the method of direct series expansion.

5 Conclusions

In this paper, we have presented new versions of the fourth-order image approximations to the reaction field due to a point charge inside a dielectric sphere immersed in an ionic solvent. Compared to the original ones which employ two image line charges, with one decaying to zero and the other growing to infinity, the new versions utilize three image line charges, all decaying to zero along the radial direction. Numerical results have demonstrated the efficiency of the new versions of the fourth-order image approximations compared to the direct series expansion.

Acknowledgments

C. Xue thanks the support of the National Natural Science Foundation of China (grant number: 10632010) and S. Deng thanks the support of the National Institutes of Health (grant number: 1R01GM08300-01) for the work reported in this paper. The authors also thank Dr. Wei Cai for many interesting discussions during this work.

A Fast calculation of the image approximations

The main purpose for discrete image approximations to reaction fields is to apply existing fast algorithms, such as the pre-corrected FFT [22] or the fast multipole methods (FMMs) [2327], directly in calculating the electrostatic interactions among N charges inside the spherical cavity in O(N log N) or even O(N) operations. In particular, below we give a straightforward but far from optimal FMM-based O(N) implementation of the discrete image approximations.

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 improved fourth-order image approximation (39) is employed, becomes

ΦRF(riF)j=1N[qK,jεi|riFrK,j|+m=1Mqm,jL1εi|riFxm,j|+(riFrK,jcosθij)m=1Mqm,jL2εi|riFxm,j|(riFrK,j)2m=1Mqm,jL3εi|riFxm,j|]+j=1NΦC1,j+j=1NΦC2,j(riF)+j=1NΦC3,j(riF),
(A.1)

where θij is the angle between riF and rjS, and a quantity with a second subscript j designates a quantity associated with the source charge rjS, such as

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

A.1 Calculation of the correction potentials in O(N) operations

Obviously, the position-independent correction potential

j=1NΦC1,j=j=1Nqjεia(C0(u)γδ1σ1)

can be evaluated in O(N) operations. The evaluation of the second correction potential

j=1NΦC2,j(riF)=j=1Nqjεia(C1(u)γδ11+σ1δ21σ2)(riFrK,j)cosθij

in O(N) operations, however, needs special treatment due to its position-dependence. Expressing the cosine of the angle θij between riF and rjS in terms of their rectangular coordinates, we get

j=1NΦC2,j(riF)=j=1Nc1qjriFrK,j(xiFxjS+yiFyjS+ziFzjSriFrjS),
(A.2)

where

c1=1εia(C1(u)γδ11+σ1δ21σ2).

Rearranging terms in (A.2) leads to

j=1NΦC2,j(riF)=d1xiF+d2yiF+d3ziF,i=1,2,,N,
(A.3)

where

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

Note that d1, d2 and d3 depend just on the source charges, and each can be calculated in O(N) operations. Then it is clear that, after obtaining d1, d2 and d3, the calculation of (A.3) for all observation points can be evaluated in O(N) operations. Analogously, the third correction term can also be evaluated in O(N) operations. In fact, using P2(cos θ) = (3 cos2 θ − 1)/2, we get

j=1NΦC3,j(riF)=e1(xiF)2+e2(yiF)2+e3(ziF)2+e4xiFyiF+e5xiFziF+e6yiFziF+e7(riF)2,i=1,,N,

where ek, k=1, …, 7, dependent just on the source charges, are gievn by

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

and

c2=1εia(C2(u)γδ12+σ1δ22σ2).

A.2 Calculation of the potentials of the image charges in O(N) operations

The FMMs are known to be extremely efficient in the evaluation of pairwise interactions in large ensembles of particles, such as expressions of the form

Φ(ri)=j=1,jiNqj|rirj|,i=1,2,,N

for the electrostatic potential, where r1,r2, …, rN are points in R3 and q1, q2, …, qN are the corresponding charge strengths. For instance, the adaptive FMM of [27] 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 [28]. 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. Such evaluation can be carried out with five FMM runs.

The first run includes all point charges qK,j at the corresponding Kelvin image points rK,j and all discrete image point charges qm,jL1 of qL1(x) at xm,j, to calculate

j=1N(qK,jεi|riFrK,j|+m=1Mqm,jL1εi|riFxm,j|).

In the case that the total potential is to be calculated, all original source charges are also included in the FMM computational box. And all charges are taken as acting in a homogeneous medium of the dielectric constant εi.

To calculate the third term in (A.1), we write it as

j=1N[(riFrK,jcosθij)m=1Mqm,jL2εi|riFxm,j|]=j=1N[(riFrK,j)(xiFxjS+yiFyjS+ziFzjSriFrjS)m=1Mqm,jL2εi|riFxm,j|]=xiFj=1Nm=1Mqm,jL2xεi|riFxm,j|+yiFj=1Nm=1Mqm,jL2yεi|riFxm,j|+ziFj=1Nm=1Mqm,jL2xεi|riFxm,j|,

where

qm,jL2x=xjSa2qm,jL2,qm,jL2y=yjSa2qm,jL2,qm,jL2z=zjSa2qm,jL2.

Then each double summation involved such as

j=1Nm=1Mqm,jL2xεi|riFxm,j|,i=1,2,,N

can be evaluated in O(N) operations with a FMM run by including in the FMM box all discrete image charges qm,jL2x, qm,jL2y or qm,jL2z at xm,j. That being done, the third term can then be evaluated in another O(N) operations.

Similarly, the fourth term in (A.1) can be written as

j=1N[(riFrK,j)2m=1Mqm,jL3εi|riFxm,j|]=(riF)2j=1Nm=1Mqm,jL3rεi|riFxm,j|,

where

qm,jL3r=(rjSa2)2qm,jL3.

The double summation for i = 1, 2, …, N can be calculated in O(N) operations with the last FMM run by including all discrete image charges qm,jL3r at xm,j in the FMM box. After that, the fourth term can be evaluated simply in an additional O(N) operations.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Koehl P. Electrostatics calculations: latest methodological advances. Curr Opin Struct Biol. 2006;16:142–151. [PubMed]
2. Levy RM, Gallicchio E. Computer simulations with explicit solvent: recent progress in the thermodynamic decomposition of free energies and in modeling electrostatic effects. Annu Rev Phys Chem. 1998;49:531–567. [PubMed]
3. Sagui C, Darden TA. Molecular dynamics simulation of biomolecules: long-range electrostatic effects. Annu Rev Biophys Biomol Struct. 1999;28:155–179. [PubMed]
4. Feig M, Brooks CL., III Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr Opin Struct Biol. 2004;14:217–224. [PubMed]
5. Baker NA. Improving implicit solvent simulations: a Poisson-centric view. Curr Opin Struct Biol. 2005;15:137–143. [PubMed]
6. Okur A, Simmerling C. Hybrid explicit/implicit solvation methods. In: Spellmeyer D, editor. Annu Rep Comput Chem. Vol. 2. 2006. Chapter 6.
7. Lee MS, Salsbury FR, Jr, Olson MA. An efficient hybrid explicit/implicit solvent method for biomolecular simulations. J Comput Chem. 2004;25:1967–1978. [PubMed]
8. 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]
9. 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.
10. Deng S, Cai W. Discrete image approximations of ionic solvent induced reaction field to charges. Commun Comput Phys. 2007;2:1007–1026.
11. 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. in press. [PMC free article] [PubMed]
12. Bordner AJ, Huber GA. Boundary element solution of the linear Poisson-Boltzmann equation and a multipole method for the rapid calculation of forces on macromolecules in solution. J Comput Chem. 2003;24:353–367. [PubMed]
13. 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.
14. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications; New York: 1972.
15. Gradshteyn IS, Ayzhik IM. Table of Integrals, Series, and Products. Academic Press; Boston: 1994.
16. Baker NA. Poisson-Boltzmann methods for biomolecular electrostatics. Methods Enzymol. 2004;383:94–118. [PubMed]
17. Fogolari F, Brigo A, Molinari H. The Poisson-Boltzmann equation for biomolecular electrostatics. J Mol Biol. 2002;15:377–392. [PubMed]
18. Morse PM, Feshbach H. Methods of Theoretical Physics. McGraw-Hill; New York: 1953.
19. Lindell IV. Electrostatic image theory for the dielectric sphere. Radio Sci. 1992;27:1–8.
20. Norris WT. Charge images in a dielectric sphere. IEE Proc-Sci Meas Technol. 1995;142:142–150.
21. 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.
22. Phillips JR, White JK. A precorrected-FFT method for electrostatic analysis of complicated 3-D structures, IEEE Trans. Comput-Aided Des Integr Circuits Syst. 1997;10:1059–1072.
23. Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comput Phys. 1987;73:325–348.
24. Greengard L. The Rapid Evaluation of Potential Fields in Particle Systems. MIT; Cambridge: 1987.
25. Greengard L, Roklin V. A new version of the fast multipole method for the Laplace equation in three dimensions. Acta Numer. 1997;6:229–269.
26. Carrier J, Greengard L, Rokhlin V. A fast adaptive multipole algorithm for particle simulations. SIAM J Sci Stat Comput. 1988;9:669–686.
27. Cheng H, Greengard L, Roklin V. A fast adaptive multipole algorithm in three dimensions. J Comput Phys. 1999;155:468–498.
28. 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.