A robust and accurate formulation of molecular and colloidal electrostatics

TL;DR

This paper introduces a boundary integral reformulation for molecular and colloidal electrostatics that eliminates singularities, enabling more accurate and stable numerical computations of potentials, fields, and forces near boundaries and contact points.

Contribution

It presents a boundary regularized integral equation formulation that removes singularities from the conventional BIM, improving accuracy and computational stability in electrostatics modeling.

Findings

Accurate evaluation of surface integrals without special numerical procedures.

Enhanced precision in calculating electric fields and forces at boundaries.

Reliable handling of close boundary configurations without numerical instability.

Abstract

This paper presents a re-formulation of the boundary integral method (BIM) for the Debye-Huckel model of molecular and colloidal electrostatics that removes the mathematical singularities that have been accepted as an intrinsic part of the conventional BIM. The essence of the present boundary regularized integral equation formulation consists of subtracting a known solution from the conventional BIM in such a way as to cancel out the singularities associated with the Green's function. This approach better reflects the non-singular physical behavior of the systems on boundaries with the benefits of (i) the surface integrals can be evaluated accurately using quadrature without any need to devise special numerical integration procedures, (ii) being able to use quadratic or spline function surface elements to represent the surface more accurately and the variation of the functions within each element is represented to a consistent level of precision by appropriate interpolation functions, (iii) being able to calculate electric fields, even at boundaries, accurately and directly from the potential without having to solve hypersingular integral equations and this imparts high precision in calculating the Maxwell stress tensor and consequently, intermolecular or colloidal forces, (iv) a reliable way to handle geometric configurations in which different parts of the boundary can be very close together without being affected by numerical instabilities, therefore potentials, fields and forces between surfaces can be found accurately at surface separations down to near contact, and (v) having the simplicity of a formulation that does not require complex algorithms to handle singularities will result in significant savings in coding effort and in the reduction of opportunities for coding errors. These advantages are illustrated using examples drawn from molecular and colloidal electrostatics.

Equations

$${\nabla }^2\varphi (\mathbf{x})-{\kappa }^2\varphi (\mathbf{x})=0$$

$$c_0\varphi ({\mathbf{x}}_0)+{\int }_S\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\text{ d}S(\mathbf{x})={\int }_{S}G(\mathbf{x},{\mathbf{x}}_0)\frac{\partial \varphi (\mathbf{x})}{\partial n}\text{ d}S(\mathbf{x})$$

$$G(\mathbf{x},{\mathbf{x}}_0)=\frac{\exp (-\kappa |\mathbf{x}-{\mathbf{x}}_0|)}{|\mathbf{x}-{\mathbf{x}}_0|}$$

$${\nabla }^2\varphi (\mathbf{x})=0$$

$$c_0\varphi ({\mathbf{x}}_0)+{\int }_S\varphi (\mathbf{x})\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\text{ d}S(\mathbf{x})={\int }_S{G}_0(\mathbf{x},{\mathbf{x}}_0)\frac{\partial \varphi (\mathbf{x})}{\partial n}\text{ d}S(\mathbf{x})$$

$$G_0(\mathbf{x},{\mathbf{x}}_0)=\frac{1}{|\mathbf{x}-{\mathbf{x}}_0|}.$$

$${\int }_S[\varphi (\mathbf{x})-\varphi ({\mathbf{x}}_0)]\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\text{ d}S(\mathbf{x})={\int }_S{G}_0(\mathbf{x},{\mathbf{x}}_0)\frac{\partial \varphi (\mathbf{x})}{\partial n}\text{ d}S(\mathbf{x}).$$

$${\nabla }^2\psi (\mathbf{x})=0.$$

$$c_0\psi ({\mathbf{x}}_0)+{\int }_S\psi (\mathbf{x})\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\text{ d}S(\mathbf{x})={\int }_S{G}_0(\mathbf{x},{\mathbf{x}}_0)\frac{\partial \psi (\mathbf{x})}{\partial n}\text{ d}S(\mathbf{x})$$

$$\psi ({\mathbf{x}}_0)$$

$$\frac{\partial \psi ({\mathbf{x}}_0)}{\partial n}$$

$${\int }_S\left[\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0)}{\partial n}-\psi (\mathbf{x})\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\right]$$

$${\int }_S$$

$$\psi (\mathbf{x})=\varphi ({\mathbf{x}}_0)g(\mathbf{x})+\frac{\partial \varphi ({\mathbf{x}}_0)}{\partial n}f(\mathbf{x})$$

$${\nabla }^{2}g(\mathbf{x})$$

$${\nabla }^{2}f(\mathbf{x})$$

$${\int }_S\left[\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0)}{\partial n}-\varphi ({\mathbf{x}}_0)g(\mathbf{x})\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}+\varphi ({\mathbf{x}}_0)\frac{\partial g(\mathbf{x})}{\partial n}{G}_0(\mathbf{x},{\mathbf{x}}_0)\right]\text{ d}S(\mathbf{x})=$$

$${\int }_S\left[\frac{\partial \varphi (\mathbf{x})}{\partial n}G(\mathbf{x},{\mathbf{x}}_0)-\frac{\partial \varphi ({\mathbf{x}}_0)}{\partial n}\frac{\partial f(\mathbf{x})}{\partial n}{G}_0(\mathbf{x},{\mathbf{x}}_0)+\frac{\partial \varphi ({\mathbf{x}}_0)}{\partial n}f(\mathbf{x})\frac{\partial G_0(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\right]\text{ d}S(\mathbf{x}).$$

$$g(\mathbf{x})=1,f(\mathbf{x})=\mathbf{n}({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0).$$

$${\nabla }^2\varphi (\mathbf{x})$$

$$=$$

$$\varphi (\mathbf{x})=\frac{q}{4\pi {ϵ}_0{ϵ}_{in}}\frac{1}{|\mathbf{x}-{\mathbf{x}}_s|}+{\varphi }_{react}(\mathbf{x}),0<r<a.$$

$$E_{solv}=\frac{1}{2}q{\varphi }_{react}({\mathbf{x}}_s).$$

$${\left[(Z+c{)}^2+R^2\right]}^{-\frac{3}{2}}+{\left[(Z-c{)}^2+R^2\right]}^{-\frac{3}{2}}=2{\left(c^2+d^2\right)}^{-\frac{3}{2}}.$$

$${\sigma }_{ij}={ϵ}_0{ϵ}_{out}\left[E_i{E}_j-\frac{1}{2}\left[E_k{E}_k+{\left(\kappa \varphi \right)}^2\right]{\delta }_{ij}\right]$$

$$F_i={\int }_S{\sigma }_{ij}{n}_j\text{ d}S$$

$$T_i={\int }_S{\epsilon }_{ijk}{r}_j\left({\sigma }_{km}{n}_m\right)\text{ d}S,$$

$$\varphi (\mathbf{x})$$

$$=$$

$$\frac{1}{{ϵ}_{in}}{\left(\frac{r_s}{a}\right)}^n+\frac{1}{{ϵ}_{in}}{A}_n$$

$$-(n+1){\left(\frac{r_s}{a}\right)}^n+n{A}_n$$

$$\frac{1}{{ϵ}_{st}}\left[B_n{\left(\frac{b}{a}\right)}^n+C_n{\left(\frac{a}{b}\right)}^{n+1}\right]$$

$$n{B}_n{\left(\frac{b}{a}\right)}^n-(n+1)C_n{\left(\frac{a}{b}\right)}^{n+1}$$

Full text

molecular and colloidal electrostatics, boundary integral method, electrostatic forces, solvation energy

††preprint: J. Chem. Phys

I Introduction

Quantitative estimates of the electrostatic interaction between constituents in molecular and colloidal systems are central to understanding the role of structures and functions in areas ranging from biology, engineering to material science. Specific applications include understanding the translocation of DNA molecules in external fields Kim and Doyle (2006) or in the vicinity of membranes Buyukdagli and Ala-Nissila (2016), estimating the energetics of protein complexation Kim, Song, and Song (2010), modeling biomembrane mechanics Miloshevsky et al. (2010), quantifying the effects of induced charges at dielectric boundaries in soft condensed matter Tyagi et al. (2010), modeling charge transfer energies in electric cells Xiao and Song (2014), electrolyte theory Lange and Herbert (2011) and simulation studies of colloidal systems Alvarez and Tellez (2010). In the majority of cases, the focus is on the forces and energetics of systems comprised of charged components in an aqueous electrolyte. For general qualitative understanding and in many instances even sufficient for quantitative precision, the linear Debye-H$\ddot{\text{u}}$ckel or the linearised Poisson-Boltzmann model has served as an informative and tractable starting point that can capture most of the important physical ingredients. This has been recognized as such since the classic contributions of Kirkwood Kirkwood (1934) on the solvation energies of ionic species and of the Derjaguin-Landau-Verwey-Overbeek (DLVO) model Derjaguin and Landau (1941); Verwey and Overbeek (1948) of the interaction between colloidal particles developed in the early part of the last century. Even at present, the Kirkwood model of charged molecules treated as point charges embedded in a continuum dielectric body that is immersed in an implicit continuum solvent model electrolyte is still very much in use for understanding the interaction between complex biological molecules. Furthermore, the available analytic solutions for simple geometries are valuable benchmarks for complex computational methods in molecular electrostatics. Similarly the paradigm introduced by the DLVO theory still underpins current understanding in colloidal interactions.

Understandably, the Kirkwood and DLVO models are based on simple geometries such as spheres, cylinders or planes for which analytical solutions are available for the Debye-H$\ddot{\text{u}}$ckel equation for the electrostatic potential $\phi(\bm{x})$ in an electrolyte:

[TABLE]

where the ionic concentration of the electrolyte is characterized by the Debye length, $1/\kappa$. For general geometries, the partial differential equation, Eq. (1), has to be solved numerically using finite difference methods, finite element methods or boundary element methods.

The finite difference methods and finite element methods are based on discretization of the 3D domain. This approach requires appropriate choice of variable grid resolution that can faithfully represent complex surface shapes and surface spacings as well as handling infinite domains. For complex geometries, this can be challenging. However, these methods generate large but sparse matrix equations for which well-tested computational algorithms are available to handle the numerical task.

In contrast, the boundary element method uses Green’s identity to express the solution of Eq. (1) in terms of $\phi(\bm{x})$ and its normal derivative $\partial\phi(\bm{x})/\partial n\equiv\bm{n}(\bm{x})\cdot\nabla\phi(\bm{x})$ on the 2D boundary, $S$, of the 3D domain where $\bm{n}(\bm{x})$ is the outward unit normal of $S$ at $\bm{x}$. The values of these functions on the boundary are obtained by solving the surface integral equation Becker (1992)

[TABLE]

where the points $\bm{x}$ and $\bm{x}_{0}$ both lie on the surface $S$ that may be the surface of a molecule or a colloidal particle. Here

[TABLE]

is the Green’s function:$\nabla^{2}G-\kappa^{2}G=-4\pi\delta(\bm{x}-\bm{x}_{0})$ of Eq. (1). The constant, $c_{0}$, in Eq. (2) is the solid angle at $\bm{x}_{0}$. It is equal to $2\pi$ if the tangent plane of $S$ at $\bm{x}_{0}$ is defined, otherwise it has to be calculated from the local geometry Mantic̆ (1993); Sun et al. (2015). As we shall see later, our formulation of the boundary integral equation is independent of the value of $c_{0}$. Therefore, once the values of $\phi(\bm{x})$ and $\partial\phi(\bm{x})/\partial n$ on the surface, $S$, have been found, the value of $\phi(\bm{x})$ anywhere in the 3D domain can be determined by an integral over the surface Becker (1992).

For an electrostatic problem, either $\phi(\bm{x})$ or $\partial\phi(\bm{x})/\partial n$ or a relation between them Carnie and Chan (1993) on the surface is known from the boundary conditions appropriate to the problem. Thus the unknown quantity in Eq. (2) can be found by solving a problem in 2D. In contrast, finite difference or finite element methods require solving a problem in 3D. This reduction in dimension and hence the number of unknowns, together with the accurate account of conditions at infinity is a desirable tradeoff for a slightly more complex formulation that favors the boundary element approach.

Traditionally, the boundary element method is regarded to have two major disadvantages. Firstly, discretization of Eq. (2) results in a dense matrix equation that has $O(N^{3})$ complexity, rendering it unsuited for large problems. However, with the recent advent of fast solvers based on Krylov subspace iteration methods together with Fast Fourier Transform and Singular Value Decomposition algorithms Altman et al. (2009) the computational cost can be reduced to a practically competitive level of $O(N\log N)$. Secondly, the Green’s function formulation of Eq. (2) contains singularities as $\bm{x}\rightarrow\bm{x}_{0}$ in the two surface integrals that are inherent in the Green’s function, $G$, and its normal derivative, $\partial G/\partial n$. Such singular behavior is a consequence of the mathematical formulation and does not relate to any physical divergences on the boundary surfaces. Although these integrals are bounded in spite of the singularities in $G$ and $\partial G/\partial n$, that is, these singularities are integrable, nonetheless, the mathematical complexity needed to handle such singular behavior means that usually only the simplest area elements are used to discretize $S$. By far the most common approach is to represent $S$ by a mesh of planar triangular area elements and assume that $\phi$ and $\partial\phi/\partial n$ are constants within each element as the unknowns to be determined.

A more serious consequence of the singular nature of $G$ and $\partial G/\partial n$ is that when problems involve having two parts of the boundary being very close together, for example, when the surfaces of two ions or colloidal particles are nearly in contact, such singular behavior will limit the precision for which the surface integrals can be evaluated numerically because the near singular behavior at one surface can adversely affect the precision of the evaluation of integrals on a nearby surface.

A further problematic consequence of such singular behavior is that in calculating the force between charged entities, it is necessary to determine the electric field, $\bm{E}=-\nabla\phi$, that enters into the Maxwell stress tensor. It has been suggested that a boundary integral equation for $\bm{E}$ can be found by taking the gradient of Eq. (2) Lu, Zhang, and McCammon (2005), resulting in hypersingular integral equations in which the more strongly divergent integrals need to be interpreted as principal value integrals. Numerical evaluation of such integrals requires special care that impacts adversely on the achievable precision. This hypersingular behavior arises from the interchange of integration and differentiation – an ill-advised procedure for integrals that do not converge absolutely.

Given the achieved advances in the development of $O(N\log N)$ fast solvers for dense linear systems that arise from the boundary integral solution of the Debye-H$\ddot{\text{u}}$ckel equation, it is timely to re-examine the foundation of the integral equation, Eq. (2), and seek to eliminate the mathematical singularities that originate from the Green’s function $G(\bm{x},\bm{x}_{0})$. Such singularities, while long accepted in boundary integral equations as unavoidable, do not have any actual physical basis. Thus success in eliminating them will obviously be very beneficial to many areas in chemical physics in which molecular and colloidal electrostatics feature as key components of a larger framework.

In Secs. II and III, we present a new formulation of the boundary integral equation solution of the Debye-H$\ddot{\text{u}}$ckel equation that does not contain the traditional singularities. The consequences of this non-singular formulation are:

the surface integrals can be evaluated accurately using quadrature without any need to devise special numerical integration procedures so that, for example, standard Gauss quadrature can be used, 2. 2.

we can use quadratic surface elements to represent the surface $S$ more accurately and the variation of the functions within each element is represented to a consistent level of precision by quadratic interpolation functions; and opens the possibility to use spline functions or higher order functions to represent surfaces, 3. 3.

electric fields, even at boundaries, can be evaluated accurately and directly from the potential without having to solve hypersingular integral equations and this imparts high precision in calculating the Maxwell stress tensor and consequently, intermolecular or colloidal forces, 4. 4.

geometric configurations in which different parts of the boundary are very close together will not cause numerical instabilities, thus potentials, fields and forces between surfaces can be found accurately at surface separations down to near contact, 5. 5.

the simplicity of the formulation in not requiring complex algorithms to handle singularities means significant savings in coding effort and reduction of opportunities for coding errors, and 6. 6.

multiple domains connected by boundary conditions can be implemented with relative ease.

Since the implementation of fast $O(N\log N)$ algorithms is well-developed and documented, we will only focus on the non-singular formulation of relevant physical problems that will be called the boundary regularized integral equation formulation - BRIEF, and we present numerical examples that will highlight the precision that our approach can furnish.

II Motivation

To provide background to our boundary regularized integral equation formulation (BRIEF) for the Debye-H$\ddot{\text{u}}$ckel equation, we first consider the boundary integral formulation of the solution of the Laplace equation for $\phi$:

[TABLE]

for which the standard boundary integral equation on the surface $S$ that encloses the solution domain is Becker (1992)

[TABLE]

where the Green’s function for the Laplace equation (4) is

[TABLE]

The standard way to ameliorate the singularity in ${\partial G_{0}}/{\partial n}$ on the LHS of Eq. (5) is to note that the function $[\phi(\bm{x})-\phi(\bm{x}_{0})]$ also satisfies the Laplace equation, Eq. (4), and therefore the corresponding standard boundary integral equation for $[\phi(\bm{x})-\phi(\bm{x}_{0})]$ is

[TABLE]

Thus provided the function $\phi(\bm{x})$ is continuous at $\bm{x}_{0}$, the use of Eq. (7) instead of Eq. (5) means that the remaining integrable weak singularities can be accommodated by a local change of variables when $\bm{x}\rightarrow\bm{x}_{0}$.

We now extend this approach to the case when $\phi(\bm{x})$ satisfies the Debye-H$\ddot{\text{u}}$ckel equation.

III Formulation

For ease of future reference, we designate the boundary integral equation in Eq. (2) derived from the conventional boundary integral method as CBIM. Our objective is to remove analytically the singularities in $G$ and $\partial G/\partial n$ in Eq. (2), that occur as $\bm{x}\rightarrow\bm{x}_{0}$. Such singularities are the consequence of the mathematical derivation of Eq. (2) and have no intrinsic physical basis. It is also worth noticing that the singular behavior of $G$ is in fact identical to $G_{0}\equiv 1/|\bm{x}-\bm{x}_{0}|$ since $G\equiv G_{0}+\Delta G$ where $\Delta G\equiv[\exp(-\kappa|\bm{x}-\bm{x}_{0}|)-1]/|\bm{x}-\bm{x}_{0}|$ is finite as $\bm{x}\rightarrow\bm{x}_{0}$. The same analysis can be applied to $\partial G/\partial n$ as well. The approach that we use to deal with these singularities, a simple version of which is given in Sec. II, is adapted from de-singularisation of integral equations that occur in fluid mechanics, elasticity and acoustics Klaseboer, Sun, and Chan (2012); Sun et al. (2014, 2015).

We begin by first considering the function $\psi(\bm{x})$ that satisfies the Laplace equation for the same domain as where Eq. (1) is valid:

[TABLE]

The corresponding conventional boundary integral equation for $\psi(\bm{x})$ for the same surface $S$ is

[TABLE]

where $G_{0}$ is given by Eq. (6). In order to use $\psi(\bm{x})$ to remove the singularities in Eq. (2), we require it to assume the following special values at $\bm{x}=\bm{x}_{0}$:

[TABLE]

Thus subtracting Eq. (9) from Eq. (2) we can eliminate the term $c_{0}\phi(\bm{x}_{0})$ to give

[TABLE]

To satisfy Eqs. (10) and (11), we can choose $\psi(\bm{x})$ to have the form

[TABLE]

where $g(\bm{x})$ and $f(\bm{x})$ satisfy the Laplace equation with the following conditions at $\bm{x}=\bm{x}_{0}$

[TABLE]

and this will ensure both integrands in Eq. (III) will not be singular at $\bm{x}=\bm{x}_{0}$. Note that the conditions on $g(\bm{x})$ and $f(\bm{x})$ in Eqs. (14) and (15) are constraints on at a single position $\bm{x}=\bm{x}_{0}$ and are not general boundary data on a surface.

Finally, we have the key result for our boundary regularized integral equation formulation (BRIEF):

[TABLE]

It can be shown that this integral equation contains no singularities provided $g(\bm{x})$ and $f(\bm{x})$ have simple mathematical smoothness properties,Klaseboer, Sun, and Chan (2012) and so Eq. (III) can be solved by straightforward numerical methods. It replaces Eq. (2), the equation from the conventional boundary integral method (CBIM) that contains singular integrands Altman et al. (2009); Lu, Zhang, and McCammon (2005); Yoon and Lenhoff (1990); Cooper, Bardhan, and Barba (2014); Juffer et al. (1991). This equation connecting $\phi$ and $\partial\phi/\partial n$ on the surface $S$ is the starting point of our boundary regularized integral equation formulation (BRIEF) of molecular and colloidal electrostatics. As we shall see, the absence of mathematical singularities in the BRIEF simplifies numerical implementation and coding effort resulting in very significant improvements in numerical precision even for problems that pose nearly insurmountable difficulties for the CBIM.

There are many possible choices for $g(\bm{x})$ and $f(\bm{x})$ satisfying Eqs. (14) and (15) that will ensure Eq. (III) is non-singular. For instance, we can take

[TABLE]

Although for this simple choice, $\psi(\bm{x})$ does not vanish at infinity for external problems, its integral over the closed surface at infinity can be found analytically and its magnitude is equal to $4\pi\phi(\bm{x}_{0})$ Klaseboer, Sun, and Chan (2012); Sun et al. (2015) while its sign depends on the direction of the normal vector. However, there are other choices that may better suit the problem at hand Sun et al. (2014).

Although we have presented the derivation of the BRIEF by considering an integral equation on a single surface, $S$, the generalisation to more complex surface topologies and multiple domains is straightforward Altman et al. (2009), see also the examples below.

IV Implementation and Examples

To illustrate the utility and precision that can be achieved by our boundary regularized integral equation formulation (BRIEF) we benchmark against results from recent approaches using the conventional boundary integral method (CBIM) to calculate the solvation energy of the venerable Kirkwood model of an ion in electrolyte. For applications in colloidal electrostatics, we compare against analytic results for the interaction between two spheres available as infinite series expansions in terms of orthogonal functions. We also present illustrative examples of systems with surfaces that are almost in contact for which the geometry can be quite challenging for the CBIM because of the presence of singularities in the integral equation but in contrast, pose no numerical problems using BRIEF because of the absence of singularities in the integrands.

IV.1 Molecular electrostatics

IV.1.1 Kirkwood ion

To benchmark against known analytic results we consider calculations of the solvation energy of a Kirkwood ion Kirkwood (1934) that was the first molecular electrostatics model on which the conventional boundary integral method was testedYoon and Lenhoff (1990). The model comprises a point charge, $q$ embedded at a location $\bm{x}_{s}=(0,0,r_{s})$ from the center of a sphere of radius, $a~{}(>r_{s})$ and dielectric constant, $\epsilon_{in}$, immersed in a solvent of dielectric constant, $\epsilon_{out}$. The solvent can also be an electrolyte characterized by a Debye length, $1/\kappa$. In addition, the spherical ion can have a concentric outer shell or Stern layer, of radius, $b$, that excludes ionic species in the thin layer, $a<r<b$. The Stern layer may also have a different dielectric constant, $\epsilon_{st}$. A schematic representation of the Kirkwood ion is given in Fig. 1b.

The potential $\phi$ is determined by the following equations in different spatial domains

[TABLE]

The solution inside the ion $0<r<a$ can be written as the sum of the coulomb potential due to the point ion and a reaction potential

[TABLE]

At the boundaries $r=a$ and $r=b$, we have the continuity conditions of $\phi$ and $\epsilon~{}\partial\phi/\partial r$, whereas $\phi$ vanishes as $r\rightarrow\infty$. The solution expressed in terms of infinite series is given in Appendix A.

The solvation energy, $E_{solv}$ is given in terms of the reaction potential evaluated at the location of the point ion at $\bm{x}=\bm{x}_{s}$ Kirkwood (1934)

[TABLE]

Details of using BRIEF to calculate $\phi_{react}(\bm{x}_{s})$ are given in Appendix B. The surface of the Kirkwood sphere is represented by triangular shaped quadratic elements. The unknowns are potential values at the nodes of the elements. Variation of the potential within each quadratic element is obtained by quadratic interpolation of the nodal values. The surface integral over each element is calculated by quadrature and as the resulting linear system is not large, the linear equations are solved by Gauss elimination.

In Fig. 1a we compare the relative error in $E_{solv}$ obtained by BRIEF for a Kirkwood ion to corresponding results obtained by CBIM in a dielectric solvent Cooper, Bardhan, and Barba (2014) and an electrolyte solvent Altman et al. (2009) for different mesh size or number of nodes used in the evaluation of the surface integrals. It is evident that the relative error obtained using the BRIEF that has no singular integrals can be 2 orders of magnitude smaller than that of the CBIM, or conversely for the CBIM to achieve the same precision as the BRIEF, over a 100-fold increase in the number of nodes will be required. These results clearly demonstrate the superior efficiency of the BRIEF. In Fig. 1b, we quantify how the position of the point ion inside the solute sphere can affect the relative error in the solvation energy. We see that even when the ion is located at 1 Å from the surface of a 10 Å radius solvent sphere, accurate results can still be obtained by the BRIEF with a modest number of nodes. When the distance of the ion from the surface of the ion is comparable to the thickness of the Stern layer, the rate of convergence with respect to the number of nodes is slower.

IV.1.2 Dumbbell zwitterion

To illustrate how our boundary regularized integral equation formulation (BRIEF) can be used to calculate accurately the potential and electric field $\bm{E}=-\nabla\phi$ at the surface of charged molecules without using the hyper-singular boundary integral formulations (see Appendix C for details), we consider an axisymmetric dumbbell shaped zwitterion described in the body axis frame $(X,Y,Z)$ where $Z$ is the axis of rotation and $R^{2}=X^{2}+Y^{2}$ byChwang and Wu (1974), (see Fig. 2)

[TABLE]

A point charge, $q$ or $-q$ is placed at each of the foci at $Z=\pm c$ inside the dielectric dumbbell that has dielectric constant $\epsilon_{in}$. The zwitterion is immersed in a continuum electrolyte characterized by dielectric constant $\epsilon_{out}$ and Debye length $1/\kappa$. The dumbbell has length, $2a$, with a narrow neck of width $2d$ at $Z=0$ in between two lobes of width $2b$.

In Fig. 2 we show a pair of zwitterions at separation 9.514 Å between their points of closest approach with the surface potential $\phi$ on each zwitterion indicated by a color scale. The corresponding external electric vector field on the surface is calculated according to Eq. (59). Also, corresponding to Fig. 2a, we show a video of the variation of the surface potential as the relative orientation of the zwitterions changes (multimedia view).

IV.2 Colloidal electrostatics

IV.2.1 Potential and field

We now illustrate the utility of the boundary regularized integral equation formulation (BRIEF) in colloidal electrostatics. In the Debye-H$\ddot{\text{u}}$ckel model, the colloidal particles of dielectric constant, $\epsilon_{in}$, are assumed to carry a specified uniform surface charge density, $\sigma_{ch}$. They interact across an electrolyte characterized by $\epsilon_{out}$ and $\kappa$. Inside the dielectric particle, the electrostatic potential, $\phi$, obeys the Laplace equation and in the outer electrotrolyte, $\phi$ is governed by Eq. (1). The usual electrostatic boundary conditions at the surface are the continuity of $\phi$ and the condition $[\epsilon~{}\partial\phi/\partial n]_{in}-[\epsilon~{}\partial\phi/\partial n]_{out}=\sigma_{ch}/\epsilon_{0}$. For simple geometries such as for the interaction between two colloidal spheres, it is possible to obtain explicit forms for $\phi$ in terms of infinite series expansions of orthogonal functions Glendinning and Russel (1983). The boundary integral method has been used as the starting point of a perturbation calculation McCartney and Levine (1969).

We have verified that our non-singular boundary regularized integral equation formulation (BRIEF) can reproduce the infinite series solution for the force between two spheres of the same size Glendinning and Russel (1983). Here we highlight the advantages of the BRIEF in being able to calculate the potential very accurately in the region between two very nearly touching dielectric spheres - a problem that is very challenging using the series expansion method or the conventional formulation of the boundary integral method. The dielectric spheres with $\epsilon_{in}=2$, radii $a$ and $3a$ are positioned at a minimum separation $h=a/1000$. The spheres carry equal and opposite uniform surface charge densities, $\pm\sigma_{ch}$ and are immersed in an electrolyte characterized by $\epsilon_{out}=80$ and $\kappa a=1$. The centers of the spheres are located along the $z$-axis and the origin of the Cartesian axes system is midway between the surfaces of the nearly touching spheres with $z=0$ being the median plane (see Fig. 3a).

In Fig. 3b, we show the variation of the potential in the median plane obtained using the conventional boundary integral method (CBIM in Eq. (61)) and by the present non-singular boundary regularized integral equation formulation (BRIEF in Eq. (C)). At this small separation ($a/h=1000$, or $\kappa h=10^{-3}$), the numerical precision of calculating the integral on one surface is adversely affected by the intrinsic near-singular behavior of the CBIM from the nearby surface is evident in the large errors or variations in the region $|x|/a<0.5$. In contrast, the non-singular nature of the BRIEF means that in the same region, the calculation of the potential is not sensitive to the influence of proximal surfaces and that the potential variation is smooth and well-behaved as expected.

IV.2.2 Forces and torques

The force acting on a colloid particle is calculated by integrating the Maxwell stress tensor over the particle surface. In the Debye-H$\ddot{\text{u}}$ckel model, the $ij$-component of the Maxwell stress tensor, $\sigma_{ij}$, is Bell, Levine, and McCartney (1970)

[TABLE]

where $E_{i}$ is the $i$th component of electric field $\bm{E}$, and $\delta_{ij}$ is the Kronecker Delta function. Thus, the $i$th component of the force, $\bm{F}$, acting on the particle is

[TABLE]

in which $n_{j}$ is the $j$th component of the outward surface normal. The torque about the Cartesian axis $i$ is calculated using

[TABLE]

where $\varepsilon_{ijk}$ is the three-dimensional Levi-Civita symbol, and $r_{j}$ is the $j$th component of the vector between the surface node and the center of the particle. For a pair of particles we confirm that the numerical results for the force and torque are consistent whether we integrate over the surface of one particle or the other.

In Fig. 4 we show results for the force and torque between two identically charged dumbbells, described by Eq. (23), that carry a uniform surface charge density, $\sigma_{ch}$. The dumbbells have dielectric constant $\epsilon_{in}=2$ and are immersed in a solvent of dielectric constant $\epsilon_{out}=80$ and Debye length $1/\kappa=7.955$ Å. The relevant dimensionless parameters are: $\kappa a=1$, $b/a=0.4757$, $c/a=0.5303$, $d/a=0.25$. We give results for two separations: (i) at $\kappa l=0.4524$ (see Fig. 4b for the definition) so that the distance of closest approach between surfaces of the dumbbells at their widest part is very small: $\kappa h=10^{-3}$; and (ii) at $\kappa l=0.3$, where the waists of the 2 dumbbells dovetail around each other so that the relative orientation between the axes of the dumbbells is confined within a limited range around $\pi/2$. These results show that present non-singular boundary regularized integral equation formulation (BRIEF) is very robust and being able to furnish stable numerical results when the system parameters such as separation and size ratios are at quite extreme limits.

V Conclusions

The boundary integral method is a very powerful approach for studying molecular and colloidal electrostatics. Its principal advantage from the physical point of view is that the solution is cast in terms of solving for unknowns on the surfaces that define the charged system. This is particularly beneficial if the surfaces have complex and important facets or when the primary interest is in the potentials and electric fields near surfaces or forces at small separations. Unlike volume based methods such as finite difference or finite element methods, one does not have to be concerned with multi-scale meshing issues of the 3D domain or with the behavior of the potential and field at infinity.

The technical challenge of boundary integral methods in giving rise to dense matrix systems that required impractical execution times has been overcome with the development of fast algorithms that are of very acceptable $O(N\log N)$ complexity Bardhan et al. (2007). The remaining hurdle facing the wider adoption of the boundary integral method is the appearance of mathematical singularities in the conventional boundary integral method. While these singularities have no physical origin, they make it very difficult to use higher order methods to represent surfaces more accurately or to develop algorithms to evaluate the surface integrals with higher precision.

In this paper we have focused on a general formulation of the boundary integral method for the Debye-H$\ddot{\text{u}}$ckel model that removes the singular behavior that has to date been accepted as an intrinsic part of the conventional boundary integral equation method. By re-casting the problem in a way that better reflects the non-singular physical behavior of the system, we have shown with examples drawn from molecular and colloidal electrostatics that this is a robust, efficient and accurate approach. The removal of physically irrelevant singularities affords considerable savings in coding effort and results in orders of magnitude improvement in numerical precision for the same problem size. The enhanced accuracy also allows the electric field at boundaries to be calculated easily and accurately without having to solve hypersingular integral equations. As a consequence, physically important quantities such as forces and torques can be calculated easily. These advances should therefore provide an impetus to use the present non-singular boundary regularized integral equation formulation (BRIEF) to tackle complex and important problems. The present framework can easily be accommodated in existing boundary integral codes with the addition of a few lines of new code whereas all existing code that handles the singularities can be discarded.

Acknowledgements.

This work is supported by the Australian Research Council (ARC) through a Discovery Early Career Research Award to QS and a Discovery Project Grant to DYCC.

Appendix A Kirkwood ion potential distribution

By positioning the point ion of charge, $q$, on the $z$-axis at $\bm{x}_{s}=(0,0,r_{s})$, the electrostatic potential, $\phi(\bm{x})$, inside and outside a Kirkwood ion with a Stern layer is governed by Eqs. (18) - (20). The analytic solution of which can be expressed as infinite series in Legendre polynomials, $P_{n}(\cos\theta)$, of order $n$ that depends on the polar angle $\theta$ measured relative to the $z$-axis Kirkwood (1934)

[TABLE]

The first sum in Eq. (27) is the spherical harmonic expansion of the coulomb potential due to the point charge at $\bm{x}_{s}$ with $r_{>}\equiv$ max$(r,r_{s})$ and $r_{<}\equiv$ min$(r,r_{s})$. The second sum in Eq. (27) is the reaction potential, $\phi_{react}(\bm{x})$ – see Eq. (21). In Eq. (29), $k_{n}(z)$ is a modified spherical function of the second kind of order $n$ Abramowitz and Stegun (1970). If the solvent in the region $r>b$ is a dielectric, $k_{n}(\kappa r)$ is to be replaced by $(b/r)^{n+1}$. The 4 coefficients $A_{n},...,D_{n}$ are determined by the continuity of $\phi$ and $\epsilon(\partial\phi/\partial r)$ at $r=a$ and $r=b$.

The continuity of $\phi$ and $\epsilon(\partial\phi/\partial r)$ at $r=a$ gives

[TABLE]

and at at $r=b$ gives

[TABLE]

These 4 equations can readily be solved for $A_{n}$. For a dielectric solvent with $\kappa=0$, we make the replacements: $k_{n}(\kappa b)D_{n}\rightarrow D_{n}$ and $(\kappa b)k^{\prime}_{n}(\kappa b)D_{n}\rightarrow-(n+1)D_{n}$.

From Eqs. (21) and (22), the solvation energy, $E_{solv}$ can be expressed in terms of the reaction potential, $\phi_{react}(\bm{x})$ evaluated at the point ion: $(r,\theta)=(r_{s},0)$ Kirkwood (1934)

[TABLE]

Appendix B Numerical implementation

As an example of how our boundary regularized integral equation formulation (BRIEF) given by Eq. (III) can be discretized to give a linear system to be solved for the unknown potential, we consider the Kirkwood ion defined by Eqs. (18) - (20). This simple model has all the physical features of general problems in molecular electrostatics in which a charged molecule is modeled as having point ions embedded in a dielectric region of dielectric constant, $\epsilon_{in}$, that is immersed in an electrolyte of dielectric constant, $\epsilon_{out}$, and ionic concentration characterized by a Debye length, $1/\kappa$. In addition, a thin dielectric Stern layer of dielectric constant, $\epsilon_{st}$, that excludes ionic species separates the electrolyte from the charged molecule.

Such a system comprises two boundaries that are the inner, $S_{in}$, and outer, $S_{out}$, surfaces of the Stern layer (see Fig. 5). The potential is therefore defined by implementing the BRIEF, Eq. (III), at these two surface that result in a pair of coupled surface integral equations. The usual electrostatic boundary conditions for the continuity of $\phi$ and $(\epsilon\;\partial{\phi}/\partial{n})$ are applied at these two surfaces. On the inner surface, $S_{in}$, we have

[TABLE]

where $\phi_{in}$ and $\partial{\phi_{in}}/\partial{n}$ are the potential and its normal derivative within the particle that are on $S_{in}$ and $\phi_{st,in}$ and $\partial{\phi_{st,in}}/\partial{n}$ are the potential and its normal derivative within the Stern layer that are on $S_{in}$. Similarly, on the outer surface, $S_{out}$, we have

[TABLE]

where $\phi_{out}$ and $\partial{\phi_{out}}/\partial{n}$ are the potential and its normal derivative within the particle that are on $S_{out}$ and $\phi_{st,out}$ and $\partial{\phi_{st,out}}/\partial{n}$ are the potential and its normal derivative within the Stern layer that are on $S_{out}$.

Using Eq. (III) for the potential distribution in the region prescribed by Eq. (18), we can write a relation between the potential $\phi_{in}$ and its normal derivative $\partial{\phi_{in}}/\partial{n}$ on $S_{in}$ that is the surface of the molecule or the inner surface of the Stern layer as

[TABLE]

For the potential in the Stern layer prescribed by Eq. (19), the surface integrals of the BRIEF will be taken on both the inner surface, $S_{st,in}$, and the outer surface, $S_{st,out}$, of the Stern layer as

[TABLE]

where Eqs. (35) to (38) have been used.

Finally, the integrals from the BRIEF for the potential in the external domain are

[TABLE]

The integrals over the surface at infinity, $S_{\infty}$ in the last two terms in Eq. (B) can be easily simplified to $4\pi\phi_{out}(\bm{x}_{0})$ when $g(\bm{x})$ and $f(\bm{x})$ are chosen according to Eq. (17) thus giving

[TABLE]

In Eqs. (B) to (B), $G\equiv G(\bm{x},\bm{x}_{0})$ is given in Eq. (3) and $G_{0}\equiv G_{0}(\bm{x},\bm{x}_{0})$ is given in Eq. (6).

To solve Eqs. (B) to (B) to obtain $\phi_{in}$ $\phi_{out}$, $\partial{\phi_{in}}/\partial{n}$ and $\partial{\phi_{out}}/\partial{n}$ numerically, the surfaces, $S_{in}$ and $S_{out}$, are discretized by using quadratic triangular area elements on which each element is bounded by 3 nodes on the vertices and 3 nodes on the edges, see Fig. 6 for a total of $N$ nodes on the surface. The coordinates of a point in each element and the function values there are obtained by quadratic interpolation from the values at the nodes using the standard quadratic interpolation function ($\nu\equiv 1-\xi-\eta$)

[TABLE]

in terms of the local coordinates $(\xi,\eta)$ (see Fig. 6).

The solution of the potentials and their normal derivatives on the surfaces are expressed in terms of the values at the $N$ surface nodes. As illustrated in Fig. 5, when two observation points, $\bm{x}_{{0}_{in}}$ and $\bm{x}_{{0}_{out}}$, are chosen to be located on the inner surface, $S_{in}$, and outer surface, $S_{out}$, of the Stern layer, respectively, the surface integral equations, Eqs. (B) to (B), can be expressed as a system of linear equations

[TABLE]

where $\bm{I}$ is the identity matrix, the elements of the matrices ${\cal{H}}_{in}^{0}(\bm{x}_{0_{in}})$, ${\cal{H}}_{out}^{0}(\bm{x}_{0_{in}})$, ${\cal{G}}_{in}^{0}(\bm{x}_{0_{in}})$ and ${\cal{G}}_{out}^{0}(\bm{x}_{0_{in}})$ are the results of integrals (influence matrices) over the surface elements involving the unknown $4N$-vector ($\phi_{in}$ $\phi_{out}$, $\partial{\phi_{in}}/\partial{n}$, $\partial{\phi_{out}}/\partial{n}$) corresponding to Eqs. (B) and (B) as $\bm{x}_{{0}_{in}}$ is on $S_{in}$, and ${\cal{H}}_{in}^{0}(\bm{x}_{0_{out}})$, ${\cal{H}}_{out}^{0}(\bm{x}_{0_{out}})$, ${\cal{H}}_{out}(\bm{x}_{0_{out}})$, ${\cal{G}}_{in}^{0}(\bm{x}_{0_{out}})$, ${\cal{G}}_{out}^{0}(\bm{x}_{0_{out}})$ and ${\cal{G}}_{out}(\bm{x}_{0_{out}})$ are the results of integrals over the surface elements corresponding to Eqs. (B) and (B) as $\bm{x}_{{0}_{out}}$ is on $S_{out}$. Since the surface integral equations (B) to (B) do not have any singular behavior, these matrix elements can be calculated accurately using standard Gauss quadrature. The above set of equations is a $4N\times 4N$ linear system for the unknown complex $4N$-vectors: $\phi_{in}$ $\phi_{out}$, $\partial{\phi_{in}}/\partial{n}$, $\partial{\phi_{out}}/\partial{n}$ on the surface in the final form

[TABLE]

Appendix C Calculating the potential and field

The absence of singular integrals in the BRIEF means that the potential on the surface, $S$, can be evaluated accurately without numerical instabilities. Consequently, the electric field on the surface can be calculated without the need to solve hypersingular integral equations Lu, Zhang, and McCammon (2005). Suppose we seek the electric field, $\bm{E}=-\nabla\phi$ at node $k$ on the mesh that represents the surface $S$, shown in Fig. 7. Consider one of the surface elements, $m$, assumed for simplicity to be a planar triangle with vertices at $\bm{x}_{A}$ (node $k$), $\bm{x}_{B}$ and $\bm{x}_{C}$. The normal component of the electric field of node $k$ at $\bm{x}_{A}$, estimated from element $m$ is $\bm{E}_{k,m}$ is given by $\partial\phi/\partial n$ from the boundary integral equation, and tangential components of $\bm{E}_{k,m}$ can be estimated via finite differencing using the potential values at the nodes of this area element $m$

[TABLE]

Now all components of $\bm{E}_{k,m}$ can be found by solving Eq. (58). However, all area elements that share node $k$ also contribute to the estimate of the electric field of node $k$ at $\bm{x}_{A}$, therefore the field $\bm{E}_{k}$ at node $k$ will be the weighted contribution from $N_{k}$ such elements according to

[TABLE]

The weight, $w_{m}$, corresponding to element $m$ is taken to be inversely proportional to its area, $S_{m}$:

[TABLE]

Also, BRIEF provides a robust way to calculate the potentials at field positions with the same level of accuracy within the entire solution domain. In CBIM, the loss of accuracy due to the near singularity when the field position is close to the boundaries is usually more difficult to deal with compared to the singular behavior on the boundaries.

To calculate the potential accurately at position $\bm{x}_{p}$ in the 3D domain, we first use the CBIM to get $\phi(\bm{x}_{p})$, with $c_{0}=4\pi$,

[TABLE]

In the same manner, from Eq. (9), we have

[TABLE]

in which the $\bm{x}_{0}$ used to construct $\psi(\bm{x})$ in Eq. (13) is chosen as the node position that is closest to $\bm{x}_{p}$. Subtracting Eq. (62) from Eq. (61) yields

[TABLE]

The near singular behavior when $\bm{x}_{p}$ is close to the boundary can be eliminated by subtracting the BRIEF in Eq. (III) from the above boundary integral equation (C). This then provides a numerically robust expression for $\phi(\bm{x}_{p})$ whose accuracy is not affected by the distance between $\bm{x}_{p}$ and any boundary:

[TABLE]

This expression for the potential $\phi(\bm{x}_{p})$ at any point $\bm{x}_{p}$ in the solution domain given by Eq. (C) contains no singular or near singular behavior and will give equally good precision irrespective of the location of the field point $\bm{x}_{p}$.