|Home | About | Journals | Submit | Contact Us | Français|
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.
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  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 , 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 , 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.
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.
It is well-known that the total electric potential Ψin(r) inside the cavity is given by the solution of the Poisson equation
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) 
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
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 . 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 a point r = (r, θ, ) inside the sphere takes on the form
where Pn(x) represent the Legendre polynomials and An are the expansion coefficients given by
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.
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.
The key idea in the development of the foregoing image charge methods is to approximate with simple rational functions of n when u is small. In fact, by applying the expansion of the modified spherical Hankel function
which essentially implies
a fourth-order image approximation to the ionic solvent induced reaction field can be obtained as follows 
where the constant position-independent correction potential and the position-dependent correction potential are defined as
respectively. Here we denote
The point image charge qK and the two line image charges and are defined by
respectively, where , , , are defined in , and γ (εi — εo)/(εi + εo).
To construct a sixth-order image approximation for better accuracy, by using the Taylor expansion
and the truncation
one can arrive at the expansion of the modified spherical Hankel function in terms of 1/r as follows
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
Consequently, for n ≥ 4 we obtain
Inserting this approximation into (2.5) leads to
After some algebraic manipulations, the expansion coefficients An can be further expressed as
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
Theorem 1 If 0 ≤ u ≤ 1, then the three roots n1, n2 and n3 of the cubic equation β1 + β2n + β3n2 + n3 = 0 satisfy
Then, we have
where εr = εi/εo > 0. For 0 ≤ u ≤ 1, it can be seen that
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 ≤ ≤ π, we have
and that n1, n2 and n3 are also the three roots of the equation f(n) = 0. Therefore, we must have
Then, using partial fractions gives us
where σ1 = −n1, σ2 = n2, σ3 = n3, and
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
and accordingly, we have
where S1, S2, S3 and S4 represent the following four series
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 
First, in order to obtain an image representation for the series S1, we note that it can be written as
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
Next, to find an image representation for the second series S2, we need the integral identity which can be regarded as a Mellin transformation
and is valid for all n ≥ 0 when σ > 0. Inserting (3.17) with σ = σ1 into S2 yields
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
To find image representations for the series S3 and S4, we use a similar integral identity
which is valid for all n ≥ 1 when σ < 1 and valid for all n ≥ 2 when σ < 2. Inserting this into S3 yields
Then, applying the identity (expansion (3.15) with rS = x, q = 1 and εi = 1)
and noting that
Similarly, we can obtain
where ΦC0 is a position-independent correction potential defined as
and on the other hand, ΦC1(r), ΦC2(r), and ΦC3(r) are position-dependent correction potentials given by
For convenience, in terms of the Kronecker delta δ1n we define
Then, we have
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
where σ > 0. First, by introducing the change of variables rK/x = ((1 − s)/2)τ with τ > 0, we have
where α = τσ − 1 and
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 . Then, the numerical quadrature for the integral in (4.2) is
where for m = 1, 2, ..., M,
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.
Recall that 0 < σ1 < 1. Then using (4.4) with σ = σ1 leads to
where for m = 1, 2, ..., M,
Here and are the quadrature points and weights with α = τσ11 − 1.
Similarly, note that
Recall that 0 < σ2 < 1. Then, applying (4.4) with σ = (1 − σ2) > 0, we get
where for m = 1, 2, ..., M,
Here and are the quadrature points and weights with α = τ(1 − σ2) − 1.
Also, it can be noted that
Recalling that 1 < σ3 < 2, and applying (4.4) with σ = (2 − σ3) > 0, we get
where for m = 1, 2, ..., M,
Here and 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
Likewise, adding the last summation in the right-hand side of (4.10) to ΦC1(r) gives us a modified position-dependent correction potential
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
Then, for m = 1, 2, ..., M, we have
Here the quadrature points and weights are obtained using α = τσc − 1.
Similarly, the line image charges qL2(x) and qL3(x) are discretized into
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.
For the sake of convenience, we simply write (4.17) as
where the summation over m includes the modified Kelvin image charge at the corresponding Kelvin image point rK and all discrete image point charges , , and at , , and , respectively, and Nim is the total number of all point image charges.
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  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 , i = 1, 2, ..., N, be N observation points and , 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 , in the case that the sixth-order image approximation (4.18) is employed, becomes
Here and in the sequel a quantity with a second subscript j designates a quantity associated with the source charge , such as
Obviously, the position-independent correction potential
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
which in component form can be written as
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
which can be written as
and represents the outer product of two vectors.
In component form, (5.4) can be written as
Similarly, the fourth correction potential can be evaluated in O(N) operations as well. By expressing the cosine of the angle θ between and in terms of their rectangular coordinates, from (3.27) one can arrive at
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)
For instance, the adaptive FMM of  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 . 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 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.
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 l, l = 1, 2, ..., L, the calculation of the potential at r = (r, θ, ) inside the dielectric sphere generated by these image charges can be described by a by a local expansion
where are the spherical harmonics, and are the local expansion coefficients given by
Furthermore, for any p ≥ 1,
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 is expressed as
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 )
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
On the other hand, the FMM can be used to calculate
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, θ, ) inside the dielectric sphere can also be described by local expansions. Passing the details, we have
where Re(...) and Im(...) represent the real and the imaginary parts of a complex number, respectively, and
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.
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.
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.
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.
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 . 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.
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
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.
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.
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.