Harmonic surface mapping algorithm for electrostatic potentials in an atomistic/continuum hybrid model for electrolyte solutions

TL;DR

This paper introduces a harmonic surface mapping algorithm that efficiently computes electrostatic potentials in multi-scale electrolyte models, significantly improving accuracy near dielectric interfaces.

Contribution

The paper presents a novel harmonic surface mapping algorithm combining image charges and charge density discretization for fast electrostatic potential evaluation in hybrid models.

Findings

Improved accuracy by two orders of magnitude over Kirkwood series expansion.

Efficient evaluation of electrostatic reaction potentials in multi-scale electrolyte models.

Demonstrated effectiveness using the fast multipole method for Coulomb summation.

Abstract

Simulating charged many-body systems has been a computational demanding task due to the long-range nature of electrostatic interaction. For the multi-scale model of electrolytes which combines the strengths of atomistic/continuum electrolyte representations, a harmonic surface mapping algorithm is developed for fast and accurate evaluation of the electrostatic reaction potentials. Our method reformulates the reaction potential into a sum of image charges for the near-field, and a charge density on an auxiliary spherical surface for the far-field, which can be further discretized into point charges. Fast multipole method is used to accelerate the pairwise Coulomb summation. The accuracy and efficiency of our algorithm, as well as the choice of relevant numerical parameters are demonstrated in detail. As a concrete example, for charges close to the dielectric interface, our method can improve the accuracy by two orders of magnitudes compared to the Kirkwood series expansion method.

Equations

$$-\nabla \cdot {\epsilon }_1\nabla {\Phi }_1(\mathbf{x})=4\pi \sum _{i=1}^N{q}_i\delta (\mathbf{x}-{\mathbf{x}}_i),\text{in }{\Omega }_1,$$

$$-{\nabla }^2{\Phi }_2(\mathbf{x})+{\kappa }^2{\Phi }_2(\mathbf{x})=0,\text{in }{\Omega }_2,$$

$${\Phi }_1={\Phi }_2,{\epsilon }_1\frac{\partial {\Phi }_1}{\partial r}={\epsilon }_2\frac{\partial {\Phi }_2}{\partial r},\text{for}\mathbf{x}\in \partial {\Omega }_1.$$

$${\Phi }_1(\mathbf{x})={\Phi }_{\mathrm{Coul}}(\mathbf{x})+{\Phi }_{\mathrm{RF}}(\mathbf{x}),$$

$${\Phi }_{\mathrm{Coul}}(\mathbf{x})=\sum _{i=1}^N\frac{q_i}{4\pi {\epsilon }_1|\mathbf{x}-{\mathbf{x}}_i|}.$$

$${\Phi }_{\mathrm{Coul}}(\mathbf{x})=\sum _{i=1}^N\sum _{n,m}^{\infty }\frac{q_i}{(2n+1){\epsilon }_1}\frac{r_{<}^n}{r_{>}^{n+1}}{Y}_n^{m\ast }({\theta }_i,{\varphi }_i)Y_n^m(\theta ,\varphi ),$$

$${\Phi }_{\mathrm{RF}}(\mathbf{x})=\sum _{n,m}^{\infty }{A}_n^m{r}^n{Y}_n^m(\theta ,\varphi ),$$

$${\Phi }_2(\mathbf{x})=\sum _{n,m}^{\infty }{B}_n^m{k}_n(\kappa r)Y_n^m(\theta ,\varphi ),$$

$$k_n(\mu )=\frac{\pi e^{-\mu }}{2\mu }\sum _{l=0}^n\frac{(n+l)!}{l!(n-l)!}\frac{1}{(2\mu {)}^l}.$$

$$A_n^m(\mu )$$

$$B_n^m(\mu )$$

$$S_n(\mu )=-\frac{1}{n+1+\mu }+\mathcal{O}\left(\frac{1}{{\mu }^2}\right),$$

$$S_n(\mu )=-\frac{1}{n+1+\mu }+\mathcal{O}\left(\mu \right),$$

$$M_n(\mu )=\frac{\epsilon (n+1)S_n(\mu )+1}{\epsilon n{S}_n(\mu )-1}.$$

$$M_n(\mu )=\gamma +\frac{\bar{\delta }}{n+\sigma }+{\hat{M}}_n(\mu ),$$

$$\gamma =\frac{\epsilon -1}{\epsilon +1},\sigma =\frac{1+\mu }{1+\epsilon },\bar{\delta }=\gamma (1-\sigma )-\frac{\mu }{1+\epsilon }.$$

$$\frac{1}{n+\sigma }=r_{\mathrm{K},i}^{n+\sigma }{\int }_{r_{\mathrm{K},i}}^{\infty }\frac{1}{x^{n+\sigma +1}}dx,$$

$${\Phi }_{\mathrm{RF}}(\mathbf{x})$$

$$+\underset{\text{Higher-order correction term}}{\underbrace{\sum _{i=1}^N\sum _{n,m}^{\infty }\frac{q_i{\hat{M}}_n(\mu )}{(2n+1){\epsilon }_{1}R}{\left(\frac{r}{r_{\mathrm{K},i}}\right)}^n{Y}_n^{m\ast }({\theta }_i,{\varphi }_i)Y_n^m(\theta ,\varphi )}}.$$

$$q_{\mathrm{K},i}=\frac{q_i}{R}\left[\gamma r_{\mathrm{K},i}+\bar{\delta }\frac{(R_{\tau }-r_{\mathrm{K},i})}{2}\right].$$

$${\hat{M}}_n^{{}^{\prime }}(\mu )=M_n(\mu )-\gamma -\bar{\delta }\frac{(R_{\tau }-r_{\mathrm{K},i})}{2r_{\mathrm{K},i}}.$$

$${\Phi }_{\mathrm{RF}}(\mathbf{x})=\sum _{i=1}^N\frac{q_{\mathrm{K},i}}{4\pi {\epsilon }_1|\mathbf{x}-{\mathbf{r}}_{\mathrm{K},i}|}+\sum _{i=1}^N\sum _{n,m}^{\infty }\frac{q_i{\hat{M}}_n^{{}^{\prime }}(\mu )}{(2n+1){\epsilon }_{1}R}{\left(\frac{r}{r_{\mathrm{K},i}}\right)}^n{Y}_n^{m\ast }({\theta }_i,{\varphi }_i)Y_n^m(\theta ,\varphi ).$$

$${\Phi }_{\mathrm{corr}}(\mathbf{x})\approx \sum _{n,m}^p{\hat{A}}_n^m{r}^n{Y}_n^m(\theta ,\varphi ),$$

$${\hat{A}}_n^m=\sum _{i=1}^N\frac{q_i{\hat{M}}_n^{{}^{\prime }}(\mu )}{(2n+1){\epsilon }_{1}R}\frac{Y_n^{m\ast }({\theta }_i,{\varphi }_i)}{r_{\mathrm{K},i}^n}.$$

$${\Phi }_{\mathrm{corr}}(\mathbf{x})=\iint \left[G(\mathbf{x},\mathbf{y})\frac{\partial {\Phi }_{\mathrm{corr}}(\mathbf{y})}{\partial {\mathbf{\nu }}_y}-{\Phi }_{\mathrm{corr}}(\mathbf{y})\frac{\partial G(\mathbf{x},\mathbf{y})}{\partial {\mathbf{\nu }}_y}\right]d{S}_y,$$

$$\frac{\partial G}{\partial {\mathbf{\nu }}_y}\approx \frac{1}{\Delta y}[G(\mathbf{x},{\mathbf{y}}^{+})-G(\mathbf{x},{\mathbf{y}}^{-})],$$

$${\Phi }_{\mathrm{corr}}(\mathbf{x})\approx \iint \frac{{\partial }_{\mathbf{y}}{\Phi }_{\mathrm{corr}}(\mathbf{y})}{4\pi |\mathbf{x}-\mathbf{y}|}d{S}_y+\iint \frac{{\Phi }_{\mathrm{corr}}(\mathbf{y})/△y}{4\pi |\mathbf{x}-{\mathbf{y}}^{-}|}d{S}_y-\iint \frac{{\Phi }_{\mathrm{corr}}(\mathbf{y})/△y}{4\pi |\mathbf{x}-{\mathbf{y}}^{+}|}d{S}_y.$$

$$\iint f(\mathbf{y})d{S}_y\approx \frac{2\pi R_{\tau }^2}{F_2}\sum _{j=0}^{F_2}[1+\cos (\pi z_j)]\times [f({\mathbf{y}}_{2j+1})+f({\mathbf{y}}_{2j+2})],$$

$${\Phi }_{\mathrm{corr}}(\mathbf{x})\approx \sum _{j=1}^{N_{\tau }}\frac{{\bar{q}}_j}{|{\bar{\mathbf{y}}}_j-\mathbf{x}|},$$

$${\Phi }_{\mathrm{RF}}(\mathbf{x})\approx \sum _{j=1}^{N+N_{\tau }}\frac{Q_j}{|{\mathbf{Y}}_j-\mathbf{x}|},$$

$$E_{\mathrm{self}}=\frac{1}{2}{\Phi }_{\mathrm{RF}}({\mathbf{x}}_s,{\mathbf{x}}_s).$$

Full text

1 Introduction

Electrostatic effect is ubiquitous in nature, and have caught broad attention in theoretical and numerical investigations, such as the criticality in electrolytes [60, 47, 42], stability of colloid suspensions [38, 28, 17, 43], and charged biomolecular systems [30, 41, 49]. For all these studies, an accurate model of the electrolyte solvent is essential, which have aroused widespread concern up to the present [12, 7, 9, 52, 22, 16]. The explicit solvent model [35, 34], where the solvent is represented explicitly with discrete ions and water molecules, provides an accurate description of the solvent. However, its application becomes limited due to the expensive computational cost. The implicit solvent model [48, 18, 5, 46] replaces atomic details of the solvent with a dielectric continuum, by taking the so-called mean-field approximation of the electrolyte solvent. Such model can dramatically save the computational cost, but the detailed electrostatic interaction between water molecules/ions and the biomolecule is ignored.

An alternative that taking advantage of both models is the multi-scale hybrid model [44, 45]. The hybrid model introduces a spherical cavity within which a microscopic atomistic model is used, while outside the cavity continuum theory is used to describe the (same) electrolyte solvent. In this study, one assumes that the ionic strength lies in the weak coupling regime, where the well-known linearized Poisson–Boltzmann (LPB) equation can be used to approximate the bulk electrolyte solvent accurately [13, 19, 4, 59]. Now the whole simulation system is splitted into two coupled atomistic/continuum regions, one further needs to decide the parameters in the hybrid model, namely the Debye length $\kappa$ and the inside/outside dielectric constants to self-consistently couple the two regions, i.e., minimize the artificial boundary effect near the spherical cavity. For dilute electrolytes considered here, the Debye length $\kappa$ can be accurately determined as a function of the ionic densities, while the choice of inside/outside dielectric constants $\varepsilon_{1}$/$\varepsilon_{2}$ depends on different levels of microscopic descriptions inside the cavity. There are two types of model for the inside region: i) in the hybrid $explicit/implicit$ model [44], both the ions and solvent molecules are treated explicitly, in that case the dielectric constant inside $\varepsilon_{1}$ should be taken as vacuum permittivity while the outside $\varepsilon_{2}$ takes the permittivity of the solvent; ii) by contrast, in the hybrid $primitive/implicit$ model [56] the solvent inside the cavity is modeled implicitly as a dielectric continuum but the ions are treated explicitly, in which case the inside dielectric $\varepsilon_{1}$ should also be taken to be that of the solvent. In recent years, the multi-scale model has been applied in Monte Carlo simulations of $1:1$ electrolytes and compared with the periodic boundary condition (PBC) using Ewald-based methods [37, 36]. The hybrid model shows its advantage in capturing the correct charge density profile with a smaller simulation domain, while the PBC was found to give artifacts [37]. However, after introducing the multi-scale hybrid model, one needs to solve for the reaction potential inside the cavity due to the implicit solvent outside. Thus it becomes very important to improve the performance in solving the LPB equation in the presence of a spherical dielectric interface.

To solve the electrostatic reaction potential for the hybrid model of electrolytes, a variety of approaches have been proposed. For water solvent, Friedman [20] developed the image charge approximation methods, and later Abagyan and Totrov [2] proposed a modified approximation based on Friedman’s approach. These image methods have been extensively used in molecular dynamics or Monte Carlo simulations. The multiple image charge method has also been proposed [10], which can be further accelerated using the fast multipole method (FMM) [25, 26, 11, 57] with $\mathcal{O}(N)$ complexity. And later, the high-order accurate image charge method has been developed [15, 55, 54, 56]. Futhermore, for an ionic solvent, Kirkwood derived the analytical solution of the reaction potential, i.e., the Kirkwood series expansion [33, 50]. However, for large-scale simulations the performance of these methods is still not satisfactory, limiting the applications of the hybrid model.

In this paper, we develop a fast algorithm for the multi-scale hybrid model, which combines the image charge method for the near-field contribution, and the harmonic surface mapping algorithm (HSMA) [58] for the far-field. The HSMA is a recently proposed algorithm for fast Coulomb summation as an alternative to periodic boundary condition, but its application is so far restricted to a dielectric homogeneous system. Here we extend the HSMA to the multi-scale hybrid model, which is non-trivial due to the extra dielectric interface condition and the LPB equation instead of Poisson. Particularly, direct-forward applying the HSMA to the hybrid model will suffer from slow convergence problem as the charges approach the interface. To overcome this issue, we further combine it with the method of images to resolve the near-field singularity. Numerical results demonstrate that our algorithm combining HSMA and the method of images can achieve much higher accuracy than the Kirkwood series expansion given the same truncated expansion order $p$. We also examine the challenging case of a source charge very close to the dielectric interface, where our method improves the accuracy by two orders of magnitudes, owing to the analytical resolving of singularity with the image charges. In practice, this improvement may significantly help weaken the artificial boundary effect and reduce the size of the simulation box.

The rest of the paper is organized as follows. We describe the model and derivation of the HSMA in Sec. 2. The resulting algorithm, its computational complexity, and error analysis are summarized in Sec. 3. Then numerical results are given in Sec. 4, where the accuracy and efficiency performance are demonstrated through a few concrete examples. Finally, our conclusion and future works are summarized in Sec. 5.

2 Method

2.1 Model and mathematical formulations

Consider a set of $N$ point sources located at $\bm{x}_{i}=(r_{i},\theta_{i},\phi_{i})$ inside a spherical domain $\Omega_{1}\in\mathbb{R}^{3}$, each carrying charge $q_{i}$, while the outside solvent domain $\Omega_{2}$ is modeled as a dielectric continuum. The dielectric sphere $\Omega_{1}$ centered at the origin with radius $R$, as is illustrated in Fig. 1. The potential inside $\Omega_{1}$ satisfies the Poisson equation, and the linearized Poisson-Boltzmann equation (LPB) can be used to approximate the potential in $\Omega_{2}$ when the ionic strength of the electrolyte solution is week, according to the Debye-Hückel theory [14]. The electrostatic potential at $\bm{x}=(r,\theta,\phi)$ satisfies

[TABLE]

where $\Phi_{\mathrm{1}}(\bm{x})$ and $\Phi_{\mathrm{2}}(\bm{x})$ are electrostatic potentials in $\Omega_{1}$ and $\Omega_{2}$, which satisfy the electrostatic interface conditions at $\partial\Omega_{1}$:

[TABLE]

Note that $\varepsilon_{1}$ and $\varepsilon_{2}$ are dielectric constants in $\Omega_{1}$ and $\Omega_{2}$, respectively, $\varepsilon_{1}=\varepsilon_{2}$ if primitive model is used in $\Omega_{1}$, and $\varepsilon_{1}\approx 1$ if all the molecules are treated explicitly in $\Omega_{1}$. $\kappa$ is the inverse Debye length, $\kappa=(4\pi l_{\mathrm{B}}\sum_{j}\lambda_{j}z_{j}^{2})^{1/2}$, where index $j$ runs over all the ion species, and $\lambda_{j}$ and $z_{j}$ are the bulk concentration and valence of the $j$th ion species. $l_{\mathrm{B}}$ is the solvent Bjerrum length, which equals $7.14$Å for water at room temperature. Finally, the far-field boundary condition is $\Phi_{\mathrm{2}}\to 0$ as $r\to\infty$.

Note that in solving the electrostatics for an arbitrary shaped dielectric interface, boundary integral equation methods have been developed [24, 6, 8]. However, for spherical geometry, more efficient methods can be developed in solving the LPB as was summarized in the introduction, or even for the nonlinear PB equation [51]. One may argue that the spherical geometry as depicted in Fig. 1 is very specialized, has rather limited applications. However, in the multi-scale modeling of electrolytes, the spherical interface is artificially introduced for its simplicity, and serves as an alternative to the popular periodic boundary condition.

2.2 Kirkwood series expansion revisited

We first revisit the Kirkwood series expansion [33, 50] for solving Eqs. (1)–(3). The potential $\Phi_{1}$ inside $\Omega_{1}$ can be written as the sum of two contributions,

[TABLE]

where $\Phi_{\mathrm{Coul}}$ is the Coulomb potential due to the point source charges,

[TABLE]

One can further expands the Coulomb potential using spherical harmonics [31], obtaining the following expression for $\Phi_{\mathrm{Coul}}$:

[TABLE]

where $\sum_{n,m}^{\infty}$,is short for the multipole expansion summation $\sum_{n=0}^{\infty}\sum_{m=-n}^{n}$, $r_{<}(r_{>})$ is the smaller (larger) value between $r_{i}$ and $r$, and $Y_{n}^{m}(\theta,\phi)$ is the spherical harmonic function of degree $n$ and order $m$. Note that the superscript * denotes complex conjugate. The second contribution $\Phi_{\mathrm{RF}}$ is the reaction potential due to the exterior continuum solvent. Since $\Phi_{\mathrm{RF}}$ is a harmonic function, it can also be expanded in terms of spherical harmonics,

[TABLE]

where $A_{n}^{m}$ are the undetermined expansion coefficients.

Analogously, the potential $\Phi_{2}$ in the exterior region $\Omega_{2}$ can also be expanded as,

[TABLE]

where $B_{n}^{m}$ are unknown coefficients and $k_{n}(\cdot)$ is the modified spherical Hankel function of order $n$, defined as [3]:

[TABLE]

Now since both $\Phi_{1}$ and $\Phi_{2}$ are expanded using spherical harmonics, one can further substitute Eqs. (4)–(8) into the interface conditions (3) to solve for the unknown coefficients $A_{n}^{m}$ and $B_{n}^{m}$. By the orthogonality of the spherical harmonics, one obtains the following expressions for $A_{n}^{m}$ and $B_{n}^{m}$,

[TABLE]

where $\bm{r}_{\mathrm{K},i}$ is the so-called Kelvin image point, defined as $\bm{r}_{\mathrm{K},i}=(R/r_{i})^{2}\bm{x}_{i}$, and $\varepsilon=\varepsilon_{1}/\varepsilon_{2}$, $\mu=\kappa R$ and $S_{n}(\mu)=\frac{k_{n}(\mu)}{\mu k_{n}^{\prime}(\mu)}$.

Finally, it is worth noting that $S_{n}(\mu)$ has the following asymptotic approximations, as $\mu\to\infty$ [55, 53],

[TABLE]

and as $\mu\to 0$,

[TABLE]

thus the formula gives correct leading-order asymptotics for both low and high concentrations of the electrolytes. With the asymptotic formulas, the reaction potential can be further simplified into an image charge expression, which will be described in Sec. 2.3.

2.3 Image charge representation

In this section, one derives an image charge representation for the reaction potential. First, consider the expansion coefficients of $A_{n}^{m}(\mu)$ in Eq. (10), denoted here as $M_{n}(\mu)$:

[TABLE]

By substituting the asymptotic formula (i.e. Eq. (12)) into Eq. (14), $M_{n}(\mu)$ can be decomposed into three parts:

[TABLE]

where

[TABLE]

Note that $\widehat{M}_{n}(\mu)$ denotes the $\mathcal{O}(1/\mu^{2})$ higher-order terms in the asymptotic expansion (12). We do not attempt to find the explicit expression of $\widehat{M}_{n}(\mu)$, its contribution will be mapped onto an auxiliary surface in the HSMA method.

Next, substitute Eq. (15) into Eqs. (7) and (10), and further apply the following identity,

[TABLE]

which is valid for all $n\geqslant 0$ with $\sigma>0$ is a constant. The following expression is obtained for the reaction potential $\Phi_{\mathrm{RF}}$ as the sum of Kelvin images, line image densities, and a spherical harmonic expansion which can be regarded as a higher-order correction term:

[TABLE]

If $\varepsilon=1$ (i.e., no dielectric jump at the interface), the Kelvin images will vanish ($\gamma=0$), but the line images and the higher-order correction term will still be non-zero, due to the jump in the inverse Debye length $\kappa$ at the interface.

In [55], numerical results show that keeping the first two leading-order terms can approximate the reaction potential accurately when the source charges are not very close to the dielectric interface. In this work, one aims to keep the third correction term, which will allow us to obtain high-order accuracy even when source charges are close to the interface. And since the direct calculation of the correction term is time consuming, one introduces the harmonic surface mapping technique below to simplify it as a sum of images further.

2.4 Harmonic surface mapping

The harmonic surface mapping algorithm is a recently proposed fast algorithm [58] for solving the Poisson equation with either periodic/non-periodic boundary conditions. But it has not yet been applied to such hybrid model, where one needs to solve the LPB equation, and there exists a spherical dielectric interface.

Let us introduce the auxiliary spherical surface $\partial\Omega_{\tau}$. As was shown in Fig. 1, it is concentric with $\partial\Omega_{1}$ and encloses the whole interior domain $\Omega_{1}$. The radius of the auxiliary surface is $R_{\tau}=(1+\tau)R$, note that $R$ is the radius of $\Omega_{1}$ and we have an adjustable parameter $\tau>0$ (practically setting the value of $\tau$ will be discussed in Sec. 4). Then the line image integrals in Eq. (2.3) can be divided as $\int_{r_{\mathrm{K},i}}^{\infty}=\int_{r_{\mathrm{K},i}}^{R_{\tau}}+\int_{R_{\tau}}^{\infty}$. One uses the trapezoidal rule to approximate the first line integral on $[r_{\mathrm{K},i},R_{\tau}]$, then the quadrature weight assigned at the Kelvin point $\bm{r}_{\mathrm{K},i}$ is $\bar{\delta}\frac{(R_{\tau}-r_{\mathrm{K},i})}{2}\frac{q_{i}}{R}$. Note that its location overlaps with the original Kelvin image, so one just modifies the original Kelvin image and obtains the new Kelvin image magnitude

[TABLE]

The other trapezoidal point at $\bm{R}_{\tau}$ and the numerical discretization error can be both absorbed into the correction term, thus the modified harmonic coefficient $\widehat{M}_{n}^{{}^{\prime}}(\mu)$ is defined,

[TABLE]

Then the reaction potential $\Phi_{\mathrm{RF}}$ can be rewritten as,

[TABLE]

Note that the trapezoidal rule is used to obtain the modified Kelvin image magnitude $q_{\mathrm{K},i}$, but the above expression for $\Phi_{\mathrm{RF}}$ is still exact. However, directly calculating the second correction term in Eq. (21) will be again time-consuming, so one should discuss below how to handle it computationally through the HSMA approach.

We first define the correction term in Eq. (21) (truncated at order $p$) as

[TABLE]

where

[TABLE]

Through the Green’s second identity and the fact that $\Phi_{\mathrm{corr}}$ is a harmonic function, one can convert $\Phi_{\mathrm{corr}}$ into a surface integral on the auxiliary surface $\partial\Omega_{\tau}$,

[TABLE]

where $G(\bm{x},\bm{y})=\frac{1}{4\pi|\bm{x}-\bm{y}|}$ is the Green’s function for Poisson equation in free space and $\bm{\nu}_{y}$ is the unit outward normal vector at $\bm{y}$. One further uses the central difference scheme to approximate $\frac{\partial G}{\partial\bm{\nu}_{y}}$ in Eq. (24), i.e.,

[TABLE]

where $\bm{y}^{\pm}=(y\pm\Delta y/2,\theta,\phi)$ and $\Delta y$ is the central difference step size. Then $\Phi_{\mathrm{corr}}(\bm{x})$ can be expressed as the sum of three surface integrals,

[TABLE]

The first term represents a surface charge density, while the second and third terms are essential a central difference approximation for a surface dipole density. It is worth noting that since $\bm{x}\in\Omega_{1}$ while $\bm{y}\in\partial\Omega_{\tau}$, all three integrands in Eq. (26) are non-singular. Thus the Fibonacci numerical integration scheme [29] can be applied, which achieves $\mathcal{O}(N_{\tau}^{-6})$ convergence for $N_{\tau}$ discretized grid points.

Suppose $f(\bm{y})$ is a non-singular integrand, the Fibonacci integration method discretizes the surface integral over $f(\bm{y})$ on a sphere $\partial\Omega_{\tau}$ as

[TABLE]

where $z_{j}=(-1+2j/F_{2})$, $\bm{y}_{2j+1}=(R_{\tau},\arccos(z_{j}+\sin(\pi z_{j})/\pi),\pi jF_{1}/F_{2})$, $\bm{y}_{2j+2}=(R_{\tau},\arccos(z_{j}+\sin(\pi z_{j})/\pi),\pi+\pi jF_{1}/F_{2})$, and $F_{1}<F_{2}$ are two successive Fibonacci numbers. Therefore, after discretization using the Fibonacci integration scheme, Eq. (26) can be approximated by a sum of discrete images located on $\Omega_{\tau}$ and $\Omega_{\tau}^{\pm}$, i.e.,

[TABLE]

where $N_{\tau}=6(F_{2}+1)$ is the total number of images, and $\bar{q}_{j}$ and $\bm{\bar{y}}_{j}$ are the charge and location from the Fibonacci numerical integration. It is a noteworthy fact that $N_{\tau}$ is independent of the total number of sources $N$. Finally, one substitutes Eq. (28) to the reaction potential Eq. (21) and combines it with the Kelvin images, then obtains a simple expression for $\Phi_{\mathrm{RF}}$ as a Coulomb sum of $N+N_{\tau}$ image charges, i.e.,

[TABLE]

where $Q_{j}=q_{\mathrm{K},j},\bm{Y}_{j}=\bm{r}_{\mathrm{K},j}$ for $j=1,\dots,N$, and $Q_{j}=\bar{q}_{j},\bm{Y}_{j}=\bm{\bar{y}}_{j}$ for $j=N+1,\dots,N+N_{\tau}$. For the case of $\varepsilon=1$, although the original Kelvin images vanish, the total number of the image charges of the reaction potential is still $N+N_{\tau}$ due to the trapezoidal rule used to discretize the line images.

3 Algorithm, complexity and error analysis

We now describe the algorithm steps and its computational complexity, (as summarized in Algorithm 1).

The numerical error of HSMA comes from three parts: (a) The $p$-th order truncation error of the spherical harmonic expansion in Eq. (21). An error estimation for this part was given in [27], i.e., if truncated at order $p$, the truncation error $\mathcal{E}_{\mathrm{trunc}}\sim\mathcal{O}\left(\frac{1}{1+\tau}\right)^{p}$; (b) The discretization error from the central difference in Eq. (25) for calculating $\partial G/\partial\nu_{y}$, $\mathcal{E}_{\mathrm{diff}}\sim\mathcal{O}(\Delta y^{2})$; (c) The Fibonacci numerical integration error, $\mathcal{E}_{\mathrm{Fibo}}\sim\mathcal{O}(N_{\tau}^{-6})$. In practice, one finds that the truncation error from part (a) is the dominant part, as long as a reasonable $\Delta y$ and $N_{\tau}$ is chosen, the errors from parts (b) and (c) are minor. However, it should be noted that given the same truncation order $p$, the HSMA can achieve better accuracy than merely using the Kirkwood series, due to the fact that the numerical singularity is mainly caused by the Kelvin images, which has been handled here analytically. Numerical evidence will be shown in Sec. 4.

4 Numerical results

In this section, one tests the performance of the HSMA in terms of both accuracy and efficiency. In all the calculations, one fixes $R=1,\mu=5$, $\varepsilon_{1}=2$ and $\varepsilon_{2}=80$, and varies parameters $\tau$ and $p$ to check the accuracy. Since the error from the central difference is minor, one takes the step size to be $\Delta y=10^{-5}R_{\tau}$. And one takes the results from the Kirkwood series truncated at sufficiently large $p$ (for the worst case here one takes $p=201$) as the reference solution.

4.1 Accuracy tests

One first tests the accuracy of the HSMA by considering a unit source charge inside the spherical dielectric interface. Suppose a unit point source located at $\bm{x}_{s}=(r_{s},\theta_{s},\phi_{s})$ inside $\Omega_{1}$, i.e. $r_{s}\in(0,1)$, one considers the error in its self energy, the self energy $E_{\mathrm{self}}$ is defined as

[TABLE]

One first tests the error dependence on the adjustable parameter $\tau$. In Fig. 2 (a), one shows the numerical error in $E_{\mathrm{self}}$ as a function of the source charge location $r_{s}$ with $\tau=0.05$, $0.1$ and $0.15$, while fixing $p=30$ and $F_{2}=1597$.

One observes that for $r_{s}$ ranging from [math] to $0.95$, the error is not very sensitive about the value of $\tau$. As expected, one sees that the error increases a few orders of magnitudes as the charge approaches the interface. But even when $r_{s}=0.95$, the HSMA can still obtain absolute error $\sim 10^{-4}$. In Fig. 2 (b), one also tests the error in $E_{\mathrm{self}}$ as a function of the truncated order $p$ for the same set of values for $\tau$. Consistently, one finds that for $p$ ranging from $10$ to $65$, the error does not change much for different $\tau$ values. And as expected, the error decays exponentially as a function of $p$, e.g., the method exhibits spectral convergence in $p$. As a result, in the following numerical tests, one will fix $\tau=0.1$.

One now moves to the error dependence on the truncation order $p$ and the Fibonacci number $F_{2}$ used in the numerical integration. Here one will focus on a challenging case by taking $r_{s}=0.95$, i.e., the charge is very close to the interface. The results are shown in Fig. 3.

First, one observes that for different values of $p$, the error converges very soon as $F_{2}$ increases due to the high-order convergence of the Fibonacci integration scheme, e.g., for the case $p=65$, the error saturates if one takes $F_{2}\geq 6765$. Second, as long as one uses sufficiently large value for $F_{2}$, the magnitude of the saturated error is decided by the truncation order $p$ one chooses. For the case $r_{s}=0.95$, if one wants to achieve 5-digits accuracy, then one needs to choose $p\sim 55$. Finally, one also compares the HSMA with the original Kirkwood series solution with the same truncation order $p$. As is shown in Fig. 4, one observes that HSMA can improve the accuracy by two orders of magnitudes, compare to the Kirkwood series expansion with the same truncation order $p$.

For example, HSMA with $p=20$ achieves even better accuracy than the Kirkwood series with $p=40$ over the whole range of $r_{s}\in[0.9,0.98]$. Moreover, by choosing $p=20$ for HSMA, it is guaranteed that the absolute error remains less than $10^{-2}$, while the Kirkwood series can not achieve the same goal even if one takes $p=60$. Thus the HSMA has a clear advantage in terms of accuracy compared with the Kirkwood series approach.

4.2 CPU time tests

Here one tests the CPU time performance of our method for a large number of source charges inside the dielectric sphere. One uses FMM (the software package FMM3DLIB [1]) to accelerate the pairwise Coulomb summations, with the FMM precision fixed to be $10^{-6}$. The timing results are obtained on a 64-core workstation(4 AMD operation Processors Model 6272, 2.1 GHz with 16 cores each), and one uses 32 cores for each run. The parameters of the HSMA are chosen to be $\tau=0.1,p=20,F_{1}=987$ and $F_{2}=1597$. As was tested in Sec. 4.1, the parameters chosen here are sufficient to obtain numerical error less than $10^{-5}$ if $r_{s}<0.93R$. In the numerical tests here, one randomly generates $N$ source charges inside the dielectric sphere, with $N$ ranging from $10^{3}$ to $10^{6}$, and one calculates the reaction potential of the system. As is shown in Fig. 5, one finds that HSMA accelerated by FMM can achieve linear $\mathcal{O}(N)$ scaling. And compare with the direct sum, the break-even is around $N=1000$. For large-scale simulations, say if the system contains $10^{6}$ source charges, the CPU time of the HSMA accelerated by FMM is $24$s, while with direct sum the cost becomes $8.9\times 10^{4}$s. Thus the FMM-accelerated HSMA can be very attractive for large-scale simulations of charged particles immersed in electrolytes.

5 Conclusion

We have developed a harmonic surface mapping algorithm for calculating the electrostatic reaction potentials in multi-scale model of electrolytes. Based on the Kirkwood series solution, the asymptotic expansion is first used to rewrite the reaction potential into the sum of Kelvin images, line images, and an extra correction term. Then an auxiliary surface is introduced, by using the Green’s second identity, allowing us to transform the correction term into an integral on the auxiliary surface. Further combined with the Fibonacci integration scheme, we manage to rewrite the reaction potential into a simple pairwise Coulomb summation, which can be further accelerated by FMM in large-scale simulations to achieve linear scaling. Numerical tests demonstrate that HSMA can achieve much better accuracy comparing with the Kirkwood series solution. Particularly, even when a source charge is very close to the dielectric interface, the HSMA can still obtain pretty good accuracy, which will significantly help weaken the artificial boundary effect and reduce the size of the simulation domain. Thus the HSMA can be a useful tool for large-scale simulations of charged systems using the multi-scale hybrid model.

In the future, we plan to combine the hybrid solvent model with solute molecules inside the spherical cavity. One can couple our method with boundary integral equation methods for the simulation of biomolecules [32, 39, 40, 24, 21]; or moment/image method for colloidal suspensions [22, 23]. One advantage of the hybrid model lies in the explicit treatment of the electrolyte solvent inside the cavity, thus the ion specific/electrostatic correlation effect can be investigated, and it also avoids the artifacts of PBC [37]. Another goal is to implement the HSMA on GPUs to speed up its performance, and we shall also try to apply this multi-scale strategy to particle-based simulations.

Acknowledgements

J. F. acknowledge the financial support from the Natural Science Foundation of China (Grant Nos:11571236 and 21773165). Z. G. was supported by NSF grants DMS-1418966 and DMS-1819094. We want to thank Prof. Zhenli Xu for helpful discussions.