|Home | About | Journals | Submit | Contact Us | Français|
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.
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 [1–3] 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 [6–8], 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  and references therein. For an ionic solvent, in , 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 , 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.
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.
Inside the sphere, the electrostatic potential Φ is given by the solution to the Poisson's equation
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]
can be used to approximate the screened Coulomb potential in the solution. Here, λ is the inverse Debye screening length defined by
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 . 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
where Pn(x) are 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 (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.
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 , 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 , 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 , 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
Correspondingly, we have
Inserting this approximation in (4) leads to
where we denote
After some algebraic manipulation, the expansion coefficients An can be further expressed as
Note that β1 < 0 when 0 ≤ u < 1. Then using partial fractions gives us
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
with S1, S2 and S3 representing three series, defined respectively by
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 .
Then the first series S1 can be written as
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,
Next, using the integral identity
which is valid for all n ≥ 0 when σ > 0, the second series S2 can be written as
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
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 , 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
to rewrite the third series as S3=S4 − S5 with the two series S4 and S5 being defined as
Now, using partial fractions again leads to
where we denote
Then the series S4 can be written as S4=S6 + S7, where
Now applying the integral identity (11) with σ=1 and σ=1 − σ2, respectively, gives us
Analogously, the series S5 can be written as S5=S8 + S9, where
Applying the integral identity (11) with σ=2 and σ=2 − σ2, respectively, gives us
where the functions qL3,1(x) and qL3,2(x) are
where ΦC1 is a position-independent correction potential given by
and on the other hand, ΦC2(r) defines a position-dependent correction potential
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.
As pointed out in , 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
Now, instead of (8), we express the reaction field inside the sphere as
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.
where the second position-dependent correction potential is
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
where σ > 0.
First, by introducing the change of variables rK/x=((1 − s)/2)τ with τ > 0, we have
where α=τσ − 1 and
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 . Then, we have
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 (28) simply reduces to the usual Gauss or Gauss-Radau quadrature.
Recall that 0 < σ1 < 1, and furthermore, when εi < εo/5, we have σ1 > 1/2. Then using (28) with σ=σ1 leads to
where for m=1, 2, …, M,
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
Then using (28) with σ=1 or σ=(1 − σ2) > 1/2 leads to
where for m=1, 2, …, M,
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
Then using (28) with σ=2 or σ=(2 − σ2) > 3/2 leads to
where for m=1, 2, …, M,
In conclusion, in general we have the following fourth-order multiple discrete image approximations to the reaction potential inside the sphere.
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.
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.
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
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.
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.
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.
The main purpose for discrete image approximations to reaction fields is to apply existing fast algorithms, such as the pre-corrected FFT  or the fast multipole methods (FMMs) [23–27], 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 , 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 improved fourth-order image approximation (39) is employed, becomes
where θij is the angle between and , and 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 the second correction potential
in O(N) operations, however, needs special treatment due to its position-dependence. Expressing the cosine of the angle θij between and in terms of their rectangular coordinates, we get
Rearranging terms in (A.2) leads to
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
where ek, k=1, …, 7, dependent just on the source charges, are gievn by
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
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  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. 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 of qL1(x) at xm,j, to calculate
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
Then each double summation involved such as
can be evaluated in O(N) operations with a FMM run by including in the FMM box all discrete image charges , or at . 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
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 at in the FMM box. After that, the fourth term can be evaluated simply in an additional O(N) operations.
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.