Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2972553

Formats

Article sections

Authors

Related links

Comput Phys Commun. Author manuscript; available in PMC 2010 November 4.

Published in final edited form as:

Comput Phys Commun. 2007 November 1; 177(9): 689–699.

doi: 10.1016/j.cpc.2007.06.004PMCID: PMC2972553

NIHMSID: NIHMS243074

See other articles in PMC that cite the published article.

The recently developed high-order accurate multiple image approximation to the reaction field for a charge inside a dielectric sphere [J. Comput. Phys., 223 (2007) 846-864] is compared favorably to other commonly employed reaction field schemes. These methods are of particular interest because they are useful in the study of biological macromolecules by the Monte Carlo and Molecular Dynamics methods.

Electrostatics plays a major role in the structure and function of a biomolecule. Electrostatic interactions are long-range, and strongly dependent on the solvent environment surrounding the biomolecule under study. When modeling a biological system numerically, it has been challenging, however, to account for such environment in a manner that is computationally efficient and physically accurate at the same time. As such, theoretical modeling of electrostatic interactions has been and remains an important subject of theoretical and computational studies of biomolecules.

Explicit representation of solvent molecules [1–3] offers a detailed and accurate description of a macromolecule, yet all-atom simulations are expensive to perform due to the long-range nature of the electric forces. 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 the explicit solvent molecules can significantly reduce the computational cost. However, the implicit solvent models 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. In order to benefit the efficiency of the implicit solvent models for replacing the solvent that gives rise to much of the computational cost, while also directly modeling structural effects of the solvent in the proximity of the macromolecule, there has been considerable recent interest in developing hybrid explicit/implicit models [6–8]. In these hybrid models the macromolecule together with a boundary layer or shell of the solvent molecules are considered explicitly within a cavity. 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 potential within the cavity. The electric potential inside the cavity is 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, the reaction field can be calculated using the classical Kirkwood series expansion [9,10]. Although in theory any desired degree of accuracy can be obtained using the series expansion, its convergence rate is slow near the cavity boundary. For this reason, a few image charge approaches have been proposed in which the reaction field is represented in terms of the potential of a single image charge, including the Friedman image approximation [11] and the Abagyan-Totrov modified image approximation [12]. However, by using only one image charge these methods were limited in accuracy. Recently, a high-order accurate approximation using multiple image charges was proposed [13] which was found to perform about 20-30 times faster than the Kirkwood expansion in typical high accuracy calculations. Moreover, combined with the fast multipole methods [14,15], the multiple image approximation has the potential to calculate electrostatic interactions among *N* charges inside the spherical cavity in *O*(*N*) operations.

In this paper, we shall compare this high-order accurate multiple image approximation with other commonly employed reaction field schemes. A brief description of these image methods to the reaction field is first given in Section 2, and a comparison of the numerical results is then given in Section 3.

By linear superposition, the reaction field due to a single source charge *q* inside a spherical cavity only needs to be considered. Without loss of generality, let us consider a sphere of dielectric constant *ε*_{i} embedded in a homogeneous medium of dielectric constant *ε*_{o}. As shown in Fig. 1, the sphere has radius *a* and is centered at the origin of the coordinates. The point charge *q* is located on the *x*-axis at a distance *r*_{s} < *a* from the center of the sphere.

The electrostatic potential Ψ(**r**) at a point **r** is given by the Poisson’s equation

$$\nabla \cdot \left(\u220a\right(\mathbf{r})\nabla \Phi (\mathbf{r}\left)\right)=-4\pi q\delta (\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{s}}\mid ),$$

where *δ*(*r*) is the Dirac delta function, and *ε*(**r**) is the dielectric constant. Using the classical electrostatic theory, the reaction field of the spherical dielectric can be solved analytically [9,16]. More precisely, with respect to a spherical coordinate system (*r*, θ, ϕ), due to the azimuthal symmetry, the reaction field Ψ_{RF}(**r**) at an observation point **r**=(*r*, θ, ϕ) inside the sphere can be expressed in terms of the Legendre polynomials of cos θ, namely,

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)=\frac{q({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi {\u220a}_{\mathrm{i}}a}\sum _{n=0}^{\infty}\left(\frac{n+1}{{\u220a}_{\mathrm{i}}n+{\u220a}_{\mathrm{o}}(n+1)}\right){\left(\frac{r{r}_{\mathrm{s}}}{{a}^{2}}\right)}^{n}\phantom{\rule{thinmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right),$$

(1)

where *P _{n}*(

In the case that the point charge is close to the boundary of the sphere, when calculating the reaction field at a point also close to the boundary, the convergence rate by the Kirkwood series expansion is slow due to *rr*_{s}/*a*^{2} ≈ 1, requiring a great number of terms in the series expansion to achieve high accuracy in the reaction field.

Assume that *ε*_{i} < *ε*_{o}. Then expanding the term (*n* +1)/(*ε*_{i}*n* + *ε*_{o}(*n* + 1)) in (1) in terms of *ε*_{i}*n*/ *ε*_{o}(*n* +1) < 1 yields

$$\frac{n+1}{{\u220a}_{\mathrm{i}}n+{\u220a}_{\mathrm{o}}(n+1)}=\frac{1}{{\u220a}_{\mathrm{o}}}\left[1-\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{o}}}\frac{n}{n+1}+{\left(\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{o}}}\right)}^{2}{\left(\frac{n}{n+1}\right)}^{2}-\cdots \right],$$

which enables us to write the reaction field given in (1) as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)={B}^{\left(0\right)}\left(\mathbf{r}\right)+{B}^{\left(1\right)}\left(\mathbf{r}\right)+{B}^{\left(2\right)}\left(\mathbf{r}\right)+\cdots ,$$

(2)

where for *k*=0, 1, 2, , we have

$${B}^{\left(k\right)}\left(\mathbf{r}\right)={(-1)}^{k}\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi {\u220a}_{\mathrm{i}}{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}{\left(\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{o}}}\right)}^{k}\sum _{n=0}^{\infty}{\left(\frac{n}{n+1}\right)}^{k}{\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right).$$

In particular, the first term in (2) is

$${B}^{\left(0\right)}\left(\mathbf{r}\right)=\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi {\u220a}_{\mathrm{i}}{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\sum _{n=0}^{\infty}{\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right),$$

(3)

which is exactly the Legendre polynomial expansion of the Coulomb potential at the point **r** inside the sphere due to a point charge of strength *q*_{K} outside the sphere at the conventional Kelvin image point **r**_{k}=(*r*_{k}, 0, 0) [16], namely,

$${B}^{\left(0\right)}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{K}}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid},$$

where

$${q}_{\mathrm{K}}=\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}q,\phantom{\rule{1em}{0ex}}{r}_{\mathrm{k}}=\frac{{a}^{2}}{{r}_{\mathrm{s}}}.$$

The Kirkwood image approximation to the reaction field is then defined as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{K}}\left(\mathbf{r}\right)={B}^{\left(0\right)}\left(\mathbf{r}\right).$$

(4)

Now let us consider the second term in (2) which can be written as

$$\begin{array}{cc}\hfill {B}^{\left(1\right)}\left(\mathbf{r}\right)=& -\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi {\u220a}_{\mathrm{i}}{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\left(\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{o}}}\right)\sum _{n=0}^{\infty}{\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right),\hfill \\ \hfill & +\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi {\u220a}_{\mathrm{i}}{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\left(\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{o}}}\right)\sum _{n=0}^{\infty}\left(\frac{1}{n+1}\right){\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right).\hfill \end{array}$$

(5)

Similarly, the first series in (5) is exactly the Legendre polynomial expansion of the Coulomb potential at the point **r** inside the sphere due to a point charge of strength ${q}_{\mathrm{K}}^{\prime}$ outside the sphere at the Kelvin image point **r**_{k}, where

$${q}_{\mathrm{K}}^{\prime}=-\frac{{\u220a}_{\mathrm{i}}({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{{\u220a}_{\mathrm{o}}^{2}}\frac{a}{{r}_{\mathrm{s}}}q.$$

On the other hand, using the integral identity

$$\frac{1}{n+1}={r}_{\mathrm{k}}^{n+1}{\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{1}{{x}^{n+2}}\mathrm{d}x,$$

(6)

which is valid for all *n* ≥ 0, the second series in (5) can be written as

$${\int}_{{r}_{\mathrm{k}}}^{\infty}\left[\frac{{\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)}{4\pi {\u220a}_{\mathrm{i}}x}\sum _{n=0}^{\infty}{\left(\frac{r}{x}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right)\right]\mathrm{d}x,$$

(7)

where

$${\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)=\frac{{\u220a}_{\mathrm{i}}({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{{\u220a}_{\mathrm{o}}^{2}}\frac{q}{a}\left(\frac{{r}_{\mathrm{k}}}{x}\right),\phantom{\rule{1em}{0ex}}{r}_{\mathrm{k}}\le x.$$

Note that ${\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)$ can be regarded as the density function of a continuous line charge extending along the radial direction from the Kelvin image point **r**_{k} to infinity. Also, the integrand in (7) is exactly the Legendre polynomial expansion of the Coulomb potential at the point **r** inside the sphere due to a charge of strength ${\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)$ outside the sphere at the point **x**=(*x*, 0, 0). Hence we get

$${B}^{\left(1\right)}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{K}}^{\prime}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid}+{\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{{\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-\mathbf{x}\mid}\mathrm{d}x.$$

The Kirkwood image approximation (4) can then be improved by including *B*^{(1)}(**r**) as a correction potential, and the resulting image approximation is referred to by us as the improved Kirkwood image approximation, namely,

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{IK}}\left(\mathbf{r}\right)={B}^{\left(0\right)}\left(\mathbf{r}\right)+{B}^{\left(1\right)}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{IK}}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid}+{\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{{\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-\mathbf{x}\mid}\mathrm{d}x,$$

(8)

where

$${q}_{\mathrm{IK}}={q}_{\mathrm{K}}+{q}_{\mathrm{K}}^{\prime}=-\frac{{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}^{2}}{{\u220a}_{\mathrm{o}}^{2}}\frac{a}{{r}_{\mathrm{s}}}q.$$

Furthermore, by evaluating the integral in (8) explicitly, a compact analytical form of the improved Kirkwood image approximation can be obtained as

$${\widehat{\Phi}}_{\mathrm{IK}}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{IK}}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid}+\{\begin{array}{cc}\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{4\pi {\u220a}_{\mathrm{o}}^{2}}\frac{q}{a\xi}\mathrm{ln}\left(\frac{\sqrt{1-2\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \xi +{\xi}^{2}}+\xi -\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}{1-\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}\right),\hfill & \text{if}\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \ne 1,\hfill \\ \frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{4\pi {\u220a}_{\mathrm{o}}^{2}}\frac{q}{a\xi}\mathrm{ln}(1-\xi ),\hfill & \text{if}\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta =1,\hfill \end{array}\phantom{\}}$$

where ξ=*rr*_{s}/*a*^{2}.

Alternatively, expanding the term (*n* +1)/(*ε*_{i}*n* + *ε*_{o}(*n* + 1)) in (1) in terms of *ε*_{i}/((*ε*_{i} + *ε*_{o})(*n* +1)) < 1 results in

$$\frac{n+1}{{\u220a}_{\mathrm{i}}n+{\u220a}_{\mathrm{o}}(n+1)}=\frac{1}{{\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}}}\left[1+\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}}}\frac{1}{n+1}+{\left(\frac{{\u220a}_{\mathrm{i}}}{{\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}}}\right)}^{2}{\left(\frac{1}{n+1}\right)}^{2}+\cdots \right],$$

which enables us to write the reaction field given in (1) as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)={R}^{\left(0\right)}\left(\mathbf{r}\right)+{R}^{\left(1\right)}\left(\mathbf{r}\right)+{R}^{\left(2\right)}\left(\mathbf{r}\right)+\cdots ,$$

(9)

where for *k*=0, 1, 2, , we have

$${R}^{\left(k\right)}\left(\mathbf{r}\right)=\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}){\u220a}_{\mathrm{i}}^{k}}{4\pi {\u220a}_{\mathrm{i}}{({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{k+1}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\sum _{n=0}^{\infty}{\left(\frac{1}{n+1}\right)}^{k}{\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right).$$

In particular, the first term in (9) is

$${R}^{\left(0\right)}\left(\mathbf{r}\right)=\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{4\pi {\u220a}_{\mathrm{i}}({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\sum _{n=0}^{\infty}{\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right),$$

which is exactly the Legendre polynomial expansion of the Coulomb potential at the point **r** inside the sphere due to a point charge of strength *q*_{F} outside the sphere at the Kelvin image point **r**_{k}, namely,

$${R}^{\left(0\right)}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{F}}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid},$$

where

$${q}_{\mathrm{F}}=\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{{\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{\mathrm{s}}}q.$$

The Friedman image approximation to the reaction field is thus defined as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{F}}\left(\mathbf{r}\right)={R}^{\left(0\right)}\left(\mathbf{r}\right).$$

(10)

Next, let us consider the second term in (9) which is

$${R}^{\left(1\right)}\left(\mathbf{r}\right)=\frac{({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}){\u220a}_{\mathrm{i}}}{4\pi {\u220a}_{\mathrm{i}}{({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{2}}\frac{a}{{r}_{\mathrm{s}}}\frac{q}{({a}^{2}\u2215{r}_{\mathrm{s}})}\sum _{n=0}^{\infty}\left(\frac{1}{n+1}\right){\left(\frac{r}{{a}^{2}\u2215{r}_{\mathrm{s}}}\right)}^{n}\phantom{\rule{thickmathspace}{0ex}}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \right).$$

Using the integral identity (6) again, *R*^{(1)}(**r**) can be regarded as the potential of a continuous line charge extending along the radial direction from the Kelvin image point **r**_{k} to infinity with the charge density function given by

$${\stackrel{\u2012}{q}}_{\mathrm{F}}\left(x\right)=\frac{{\u220a}_{\mathrm{i}}({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{{({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{2}}\frac{q}{a}\left(\frac{{r}_{\mathrm{k}}}{x}\right),\phantom{\rule{1em}{0ex}}{r}_{\mathrm{k}}\le x,$$

namely,

$${R}^{\left(1\right)}\left(\mathbf{r}\right)={\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{{\stackrel{\u2012}{q}}_{\mathrm{F}}\left(x\right)}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-\mathbf{x}\mid}\mathrm{d}x.$$

(11)

Note that the line charge ${\stackrel{\u2012}{q}}_{\mathrm{F}}\left(x\right)$ is a constant multiple of the line charge ${\stackrel{\u2012}{q}}_{\mathrm{K}}\left(x\right)$.

Also, integrating (11) leads to a compact analytical form for *R*^{(1)}(**r**) as

$${R}^{\left(1\right)}\left(\mathbf{r}\right)=\{\begin{array}{cc}\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{4\pi {({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{2}}\frac{q}{a\xi}\mathrm{ln}\left(\frac{\sqrt{1-2\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \xi +{\xi}^{2}}+\xi -\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}{1-\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}\right),\hfill & \text{if}\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \ne 1,\hfill \\ -\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{4\pi {({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{2}}\frac{q}{a\xi}\mathrm{ln}(1-\xi ),\hfill & \text{if}\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta =1.\hfill \end{array}\phantom{\}}$$

Consequently, the Friedman image approximation can be improved by including *R*^{(1)}(**r**) as a correction potential, and the resulting image approximation is referred to by us as the improved Friedman image approximation, namely,

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{IF}}\left(\mathbf{r}\right)={R}^{\left(0\right)}\left(\mathbf{r}\right)+{R}^{\left(1\right)}\left(\mathbf{r}\right).$$

(12)

Since the Friedman image approximation was proposed in 1975, the paper [11] has been cited more than 140 times (*Source: Web of Science, 2006*), and the method has been applied in many areas including molecular dynamics or Monte Carlo simulations [17–19].

The Friedman image approximation (10) or (12) provides insufficient accuracy. In particular, when *r*_{s} tends to zero so that the point charge is located in the center of the sphere, the reaction field energy based on the Friedman image approximation (10) does not reproduce the Born formula [12].

Based on the Friedman image approximation (10), Abagyan and Totrov proposed a modified image approximation which is more accurate and less computationally intensive than the improved Friedman image approximation (12). Instead of an exact expression for *R*^{(1)}(**r**), a position-independent correction potential *R*^{corr} is added to the Friedman image approximation (10) so that for the particular case of a charge in the center of the sphere one gets the exact solution. The position-independent correction potential *R*^{corr} is defined as

$${R}^{\mathrm{corr}}=\frac{q({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{4\pi a{\u220a}_{\mathrm{o}}({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}.$$

The Abagyan-Totrov modified image approximation to the reaction field is then defined as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{AT}}\left(\mathbf{r}\right)={R}^{\left(0\right)}\left(\mathbf{r}\right)+{R}^{\mathrm{corr}}.$$

(13)

The reaction field energy based on the Kirkwood image approximation (4), however, reproduces the Born formula when *r*_{s} tends to zero. This implies that, when the source is located around the center of the sphere, the Kirkwood image approximation (4) should perform better than the Friedman image approximation (10), which actually has been verified numerically in Section 3. On the other hand, one can show that ${\widehat{\Phi}}_{\mathrm{K}}\left(\mathbf{r}\right)={\widehat{\Phi}}_{\mathrm{AT}}\left(\mathbf{r}\right)$ as *r*_{s} tends to zero, indicating that for relatively small values of *r*_{s}, the Kirkwood and the Abagyan-Totrov image approximations are comparable in terms of their accuracy.

Since its publication in 1994, the paper [12] has been cited more than 269 times (*Source: Web of Science, 2006*). Applications of the Abagyan-Totrov modified image approximation can be found in [12,20,21].

In essence, the image approximations to the reaction field discussed in the previous subsections all employ a single image charge to represent the reaction field with limited accuracy. Recently, a high-order accurate multiple image (MI) approximation to the reaction field was proposed [13] based on a little-known result dating back more than 120 years ago! First published in 1883 by C. Neumann [22] and later rediscovered in the 1990s independently by Lindell and Norris [23,24], a point charge at the Kelvin image point together with a continuous line charge extending from this Kelvin image point along the radial direction to infinity can be used to represent the reaction field *exactly*, namely,

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{MI}}}{4\pi {\u220a}_{i}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid}+{\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{{\stackrel{\u2012}{q}}_{\mathrm{MI}}}{4\pi {\u220a}_{i}\mid \mathbf{r}-\mathbf{x}\mid}\mathrm{d}x,$$

(14)

where

$${q}_{\mathrm{MI}}={q}_{\mathrm{F}}=\frac{{\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}}}{{\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}}}\frac{a}{{r}_{s}}q,{\stackrel{\u2012}{q}}_{\mathrm{MI}}\left(x\right)=\frac{{\u220a}_{\mathrm{i}}({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{{({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}^{2}}\frac{q}{a}{\left(\frac{{r}_{\mathrm{k}}}{x}\right)}^{{\u220a}_{\mathrm{o}}\u2215({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})},\phantom{\rule{1em}{0ex}}{r}_{\mathrm{k}}\le x.$$

The MI approximation to the reaction field is then obtained by representing the continuous line charge ${\stackrel{\u2012}{q}}_{\mathrm{MI}}\left(x\right)$ with discrete charges constructed through an appropriate numerical quadrature [13]. More precisely, without losing any generality, let *s _{m}*,

$${\int}_{{r}_{\mathrm{k}}}^{\infty}\frac{{\stackrel{\u2012}{q}}_{\mathrm{MI}}\left(x\right)}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-\mathbf{x}\mid}\mathrm{d}x\approx \sum _{m=1}^{M}\frac{{q}_{m}}{4\pi {\u220a}_{\mathrm{i}}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid},$$

where **x**_{m}=(*x _{m}*, 0, 0), and for

$${q}_{m}=\frac{{\u220a}_{i}({\u220a}_{\mathrm{i}}-{\u220a}_{\mathrm{o}})}{2{\u220a}_{\mathrm{o}}({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})}\frac{{w}_{m}{x}_{mj}}{a}q,\phantom{\rule{1em}{0ex}}{x}_{m}={r}_{\mathrm{k}}{\left(\frac{2}{1-{s}_{m}}\right)}^{({\u220a}_{\mathrm{i}}+{\u220a}_{\mathrm{o}})\u2215{\u220a}_{\mathrm{o}}}.$$

Accordingly, the high-order accurate MI approximation to the reaction field is defined as

$${\Phi}_{\mathrm{RF}}\left(\mathbf{r}\right)\approx {\widehat{\Phi}}_{\mathrm{MI}}\left(\mathbf{r}\right)=\frac{{q}_{\mathrm{MI}}}{4\pi {\u220a}_{i}\mid \mathbf{r}-{\mathbf{r}}_{\mathrm{k}}\mid}+\sum _{m=1}^{M}\frac{{q}_{m}}{4\pi {\u220a}_{i}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid}.$$

(15)

In this section, the image approximations to the reaction field described in the previous section are numerically tested. For demonstration purpose, a unit dielectric sphere is adopted. Unless otherwise specified, the dielectric constants of the unit dielectric sphere and the surrounding medium are assumed to be *ε*_{i}=2 and *ε*_{o}=80, respectively. Also, the single point charge *q* is located on the *x*-axis inside the sphere at a distance *r*_{s} < *a*=1 from the center of the sphere. In addition, the reaction field obtained by using the high-order accurate MI approximation with 201 discrete charges (*M*=200) is chosen as the exact reaction field to calculate the errors of the image approximations.

For two different source locations with *r*_{s}=0.4 and 0.999, respectively, the image approximations are used to calculate the reaction field at 21 observation points equally spaced either on the *x*-axis or on the *z*-axis, from -1 to 1. For each selected observation point, the relative errors of the image approximations in the reaction field are calculated and displayed in Fig. 2, from which several phenomena can be observed.

- All results clearly demonstrate that, for the particular choice of
*ε*_{i}=2 and*ε*_{o}=80, the reaction field obtained by using the high-order accurate MI approximation has the highest accuracy, even with only one additional image charge being included (*M*=1), regardless of where the source and the observation point are located. - As mentioned in Section 2.4, in cases that the source charge is close to the center of the sphere, the Kirkwood and the Abagyan-Totrov image approximations have comparable accuracy (Fig. 2 (a)(b)). When the source charge is close to the boundary of the sphere, the two approaches in general seem to still have comparable accuracy. However, the Abagyan-Totrov image approach yields better approximation at those observation points around the source charge (Fig. 2 (c)).
- For the particular choice of
*ε*_{i}=2 and*ε*_{o}=80, the Friedman image approximation can achieve only an accuracy of around 1% error except when the source charge is close to the boundary of the sphere and the observation point is close to the source at the same time (Fig. 2 (c)). The improved Friedman image approximation has better accuracy compared to the original one, but overall it is still not as efficient as the Abagyan-Totrov modified image approximation. - Surprisingly, the improved Kirkwood image approximation improves the original scheme only when the source is close to the boundary and the observation point is close to the source (Fig. 2 (c)). In addition, the improved Kirkwood and the improved Friedman image approximations are shown to have comparable accuracy.

For each selected source location with *r*_{s} varying from 0.01 to 0.999, the relative errors of the image approximations in the reaction field are calculated at each of *N*_{F}=8, 000 observation points uniformly distributed within the sphere. Then the maximal relative error and the normalized *L*_{2} error for each approximation scheme are calculated as

$${\Vert E\Vert}_{\infty}=\underset{1\le i\le {N}_{\mathrm{F}}}{\mathrm{max}}\mid \frac{{\Phi}_{\mathrm{RF}}\left({\mathbf{r}}_{i}\right)-\widehat{\Phi}\left({\mathbf{r}}_{i}\right)}{{\Phi}_{\mathrm{RF}}\left({\mathbf{r}}_{i}\right)}\mid ,$$

(16)

and

$${\Vert E\Vert}_{2}=\frac{{\left(\sum _{i=1}^{{N}_{\mathrm{F}}}{\mid {\Phi}_{\mathrm{RF}}\left({\mathbf{r}}_{i}\right)-\widehat{\Phi}\left({\mathbf{r}}_{i}\right)\mid}^{2}\right)}^{1\u22152}}{{\left(\sum _{i=1}^{{N}_{\mathrm{F}}}{\mid {\Phi}_{\mathrm{RF}}\left({\mathbf{r}}_{i}\right)\mid}^{2}\right)}^{1\u22152}},$$

(17)

respectively, where Ψ_{RF}(**r**_{i})and $\widehat{\Phi}\left({\mathbf{r}}_{i}\right)$ denote the actual reaction potential and its image approximation at an observation point **r**_{i}, respectively. The approximation errors using different image methods are shown in Fig. 3. As can be seen, the high-order accurate MI approximation is much more accurate than all other image methods, regardless of the source charge location. For most cases, the Abagyan-Totrov image approximation is much more accurate than both the Friedman image approaches and the improved Kirkwood image approximation, but when the source is near the boundary of the sphere, the improved Kirkwood or Friedman image approximation becomes comparable to the Abagyan-Totrov image approximation. Also, the MI method, the Kirkwood and the Abagyan-Totrov image methods are sensitive to the source location and they are all much more accurate when the source is near the center of the sphere than when the source is near the boundary of the sphere.

In this test, the dielectric constant of the surrounding medium is fixed as *ε*_{o}=80, but that of the dielectric sphere varies between *ε*_{i}=1 and *ε*_{i}=7. For each selected value of *ε*_{i}, the relative errors of the various image approximations in the reaction field are calculated at each of the same *N*_{F}=8, 000 observation points inside the sphere. Then the corresponding maximal relative error for each approximation scheme is calculated like in (16). The test is repeated for four different source locations with *r*_{s}=0.4, 0.6, 0.8 and 0.999, respectively, and the results are plotted in Fig. 4, from which the following facts can be observed.

- For large values of
*ε*_{i}(say*ε*_{i}> 5), , the accuracy that can be achieved by the Friedman, the improved Friedman and the improved Kirkwood image approximations all is around 10%, and the Abagyan-Totrov image approximation is typically accurate to only a few percent error as well, indicating that these approaches can only provide very limited accuracy in practice when the dielectric constant of the spherical cavity is relatively large (say*ε*_{i}≥ 3) while, at the same time, the dielectric constant of the surrounding solvent is relatively small (say*ε*_{o}≤ 70). On the other hand, the high-order accurate MI approximation can still achieve an accuracy of less than 1% error when only one additional image charge is adopted. - For small values of
*ε*_{i}(say*ε*_{i}< 2), the improved Kirkwood or Friedman image approximation is more accurate than other single charge image approximations, particularly when the source charge is close to the boundary of the sphere. - For all cases with
*ε*_{i}≥ 1.25, the results obtained using the high-order accurate MI approximation has the highest accuracy, even with only one additional image charge being included. - When
*ε*_{i}=1, the improved Kirkwood or Friedman image approximation is more accurate than the high-order MI image approximation with*M*=1, but still less accurate than the high-order accurate MI approximation with*M*=2 or*M*=3, depending on where the source charge is located.

As demonstrated already, typically the high-order accurate MI approximation is more accurate than all other commonly employed image approximation approaches, even with only one additional image charge being included. Moreover, like the Kirkwood series expansion, in theory any desired degree of accuracy can be achieved using the MI approximation by including a sufficient number of discrete image charges. One natural concern with this MI approach, however, is the number of discrete image charges required to achieve a certain order of degree of accuracy in the reaction field within the whole dielectric sphere. This number should be small if compared to the number of terms needed to achieve the same order of degree of accuracy in the Kirkwood series expansion to make the MI approximation useful in the practice.

To address this issue, for each of several selected source locations, the relative errors of the MI approximation with different numbers of discrete images in the reaction field are calculated at each of the same *N*_{F}=8, 000 observation points uniformly distributed inside the sphere. And then the maximal relative error for each selected number of images is calculated and displayed in Fig. 5. As can be seen, only a small number of discrete image charges, five or less than five (including the point image charge at the conventional Kelvin image point) for most cases, are needed for a 10^{−4} accuracy.

The dependence of the accuracy of the MI approximation on *M* +1, the number of the total discrete image charges, including the point image charge at the conventional Kelvin image point.

As the final remark, it should be pointed out that the high-order accurate approximation using multiple image charges was found to perform about 20-30 times faster than the Kirkwood series expansion in typical high accuracy calculations. For more details the readers may consult Ref. [13].

In this paper, the recently developed high-order accurate MI approximation and other commonly employed single image approximations were compared to calculate the reaction field due to a point charge inside a dielectric sphere embedded in a homogeneous dielectric medium. Overall, the MI approximation including one additional image charge performed well compared to the single image approximations, regardless of the source location. Only when the dielectric constant of the sphere is very small compared to that of the surrounding medium, two to three additional charges may be needed for the MI approximation to perform better than the improved Kirkwood or Friedman image approximation.

Wei Cai and Shaozhong Deng thank the support of the National Science Foundation (grant numbers: DMS-0408309, CCF-0513179) and the Department of Energy (grant number: DEFG0205ER25678), and Donald Jacobs thanks the support from the National Institutes of Health (grant number: NIGMS R01 GM073082-01A1), for the work reported in this paper.

**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.

[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, Chapter 6. In: Spellmeyer David., editor. Annu. Rep. Comput. Chem. Vol. 2 2006.

[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, An efficient hybrid explicit/implicit solvent method for biomolecular simulations. J. Phys. Chem. B. 2005;109:5223–5236. [PubMed]

[9] Kirkwood JG. Theory of solutions of molecules containing widely separated charges with special applications to awitterions. J. Chem. Phys. 1934;2:351–461.

[10] Kirkwood JG. Statistical mechanics of liquid solutions. Chem. Rev. 1936;19:275–307.

[11] Friedman HL. Image approximation to the reaction field. Mol. Phys. 1975;29:1533–1543.

[12] Abagyan R, Totrov M. Biased probability Monte Carlo conformational searches and electrostatic calculations for peptides and proteins. J. Mol. Biol. 1994;235:983–1002. [PubMed]

[13] 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.

[14] Greengard L, Rokhlin V. A fast algorithm for particle simulations. J. Comput. Phys. 1987;73:325–348.

[15] Greengard L. The Rapid Evaluation of Potential Fields in Particle Systems. MIT; Cambridge: 1987.

[16] Morse PM, Feshbach H. Methods of Theoretical Physics. McGraw-Hill; New York: 1953.

[17] Wallqvist A. On the implementation of Friedman boundary conditions in liquid water simulations. Mol. Simul. 1993;10:13–17.

[18] Wang L, Hermans J. Reaction field molecular dynamics simulation with Friedman’s image method. J. Phys. Chem. 1995;99:12001–12007.

[19] Rullmann JAC, Duijnen P.Th.V. Analysis of discrete and continuum dielectric models: application to the calculation of protonation energies in solution. Mol. Phys. 1987;61:293–311.

[20] Petraglio G. Nonperiodic boundary conditions for solvated systems. J. Chem. Phys. 2005;123:044103. [PubMed]

[21] Havranek JJ, Harbury PB. Tanford-Kirkwood electrostatics for protein modeling. Proc. Natl. Acad. Sci. USA. 1999;96:11145–11150. [PubMed]

[22] Neumann C. Hydrodynamische Untersuchen nebst einem Anhang uber die Probleme der Elecktrostatik und der magnetischen Induktion. Teubner; Leipzig: 1883. pp. 279–282.

[23] Lindell IV. Electrostatic image theory for the dielectric sphere. Radio Sci. 1992;27:1–8.

[24] Norris WT. Charge images in a dielectric sphere. IEE Proc.-Sci. Meas. Technol. 1995;142:142–150.

[25] 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |