1 Introduction
The PoissonBoltzmann theory has been widely accepted as a meanfield continuum approximation for electrostatic interactions in solvated biomolecular systems Honig1995_science . The PoissonBoltzmann theory treats the solute bimolecules as a singularly charged medium of low dielectric constant () immersed in a high dielectric () solvent with a continuum charge distribution that models the dispersed mobile ions. The large contrast of dielectric constant on the highly complicated biomolecular surfaces poses a significant computational challenge and affords delicate mathematical and numerical treatments. The importance of PoissonBoltzmann theory in biochemistry and biophysics has motivated extensive mathematical and computational investigations, see, for instance PBEQ ; LuB2007e ; RenP2012a ; XieD2013a ; ChenJ2018a ; ZhongY2018a , and references therein for the development of the subjects in the past decades and the latest overviews.
This paper is concerned with the accurate calculation of the gradient of the electrostatic potential near the molecular surface from the numerical solutions of the PoissonBoltzmann equation (PBE) for modeling solvated biomolecules. Traditionally the potential solutions of the PoissonBoltzmann equation are used merely for computing the solvation free energy (and other derived energetic quantities such as pKa value, binding affinity, etc.), and therefore less attention had been paid on the calculation of the potential gradient. The success of PoissonBoltzmann theory in these energetic evaluations promotes the exploration of PBE based electrostatic force calculations for molecular dynamics (MD) simulations MacKerellA1998a ; LuQ2003a ; Prabhu2004a ; GengW2011a and coarsegrained TangY2006a ; MaL2009 or hybrid models ChenX2010a ; ZhouY2010b ; ChengL2013a . There are three types of electrostatic forces defined by using the free energy functional derivatives and shape derivatives GilsonM1993a ; ZhouY2008c ; LiB2011a ; MikuckiM2014a , including (i) the body force exerted at each charged atoms of the solution molecules given by where is the local charge density; (ii) the dielectric boundary force exerted on the dielectric interface that is usually defined by the molecular surface; and (iii) the ionic pressure, exerted also on the molecular surface.
Perpendicular to the dielectric interface, the dielectric boundary force
(1) 
poses significant challenges to the numerical computation because it is defined on the dielectric interface where one usually observes the peaks of the numerical error in the solutions of the PoissonBoltzmann equation Stillinger1961a ; Xiang1995a ; Holst1995a ; Holst2000a ; BoschitschA2002a ; Xie2007a ; LinH2014a ; Geng2017a . For boundary element methods, while can be obtained directly along with from wellconditioned boundary integral formulations, computing involves additional integrations of the numerical solution of the potential with a supersingular kernel over the entire surface LuB2006a ; BajajC2011a ; GengW2013a ; ZhongY2018a , usually causing a larger error. For finite difference interface methods, one usually extrapolates the potential solution across the interface by using the interface conditions to compute potential gradients GengW2011a . For interface finite elements methods, gradient recovery techniques are recently developed so that potential gradient could be computed in a subdomain without using the solution values across the internal interface GuoH2018a ; GuoH2017a ; GuoH2018b .
We will develop an interface gradient recovery method to approximate the screened electrostatic potential from the numerical solution of the PoissonBoltzmann equation. For elliptic problems with smooth coefficients, gradient recovery techniques have been well established on structured or unstructed grids ZienkiewiczO1992a ; ZienkiewiczO1992b ; BankR2003a ; GuoH2015a ; NagaA2005a ; XuJ2003a ; ZhangZ2005a
. By comparison, there is only a small handful of numerical techniques concerning the gradient recovery from the primarily computed solutions of the elliptic problems with piecewise smooth coefficients or singular sources. For 1D interface problems, special interpolation formulas are constructed to recover flux with high order accuracy from numerical solutions of linear and quadratic interface finite element methods
ChouS2012a ; ChouS2015a . For 2D elliptic interface problems on a specially constructed bodyfitted mesh, gradient of the linear finite element solution is shown to be superclose to the gradient of the linear interpolation of the exact solution WeiH2014a . By introducing the jump in the normal derivative as an augmented variable, Li et. al. showed that a secondorder convergence can be obtained for the solution and its gradient LiZ2017a . Recently developed gradient recovery methods for interface problems are based on the secondorder least square reconstruction of the solutions in an individual subdomain GuoH2018a . If the mesh is not conforming to the interface, a practice common in many immersed finite element methods LiZ1998a ; LiuX2000a ; HouS2013a ; JiH2014a , the solutions on the interface are first obtained on approperiately subdivided interface elements then supplied to the least square reconstruction GuoH2017a . When Nitsche’s method HanshoP2005a ; AnnavarapuC2012a is used for solving the elliptic interface problems, this subdivision of the interface element is not necessary, and one can use the solution on the extended fictitious subdomain to carry out the least square reconstruction GuoH2018b . Nevertheless, as we shall demonstrate below, these interface gradient recovery methods are not able to accurately compute the interface potential gradient on a spatial (domain or surface) discretization that are affordable for realistic biomolecular simulations.Our method is based on the least square reconstruction of the numerical solution, with the classical polynomial basis supplemented with where is the distance between a mesh node and the selected charges in the solute biomolecule. These additional functions have been indicated in the Green’s function for the Poisson equation and particularly in the generalized Born (GB) models for the biomolecular electrostatics. In GB models, the electrostatics potential at a charge inside the dielectric boundary is given by
(2) 
where is the singular charge at spatial position , and is the reaction field due to the polarization charges induced at the dielectric boundary ImW2003a ; FeigM2004c . A large variety of GB models have been developed to implicitly account for the effects of solvation through various parameterizations of effective Born radius, solvation accessible surfaces, ionic screening, among others, see OnufrievA2019a for a latest review and references therein. In the solvent domain , although an analytical approximation similar to (2) does not exist, the approximation
(3) 
is widely used to compute the boundary conditions for the PoissonBoltzmann at a boundary point far away from the solute molecules, where is the DebyeHuckel inverse screening length. The approximations (2)(3) motivate us to choose functions of the form to enrich the classical polynomial basis function in the least square reconstruction. Numerical experiments demonstrate that our enriched recovery technique is highly accurate and stable for complicated biomolecules. By construction our recovery technique can be integrated with general interface numerical methods for the PoissonBoltzmann equation, regardless of whether the discretization mesh is interface conforming or not.
The rest of the article is organized as follows. In Section 2.1, the PoissonBoltzmann equation and its numerical treatment are briefly reviewed, followed by the introduction of dielectric boundary force and its calculation. We will present the our gradient recovery technique in Section 3 and discuss its extension to enforcing the interface conditions. Extensive numerical experiments will be conducted in Section 4 to verify the robustness and accuracy of our methods on biomolecules of different complexity. The article concludes with a summary in Section 5.
2 PoissonBoltzmann Equation and Dielectric Boundary Force
This section presents the energetic theory of biomolecular electrostatics. The electrostatic free energy introduced by Sharp, Honig, Gilson et. al. GilsonM1993a ; GilsonM1985a ; SharpK1990a ; SharpK1990b provides a unified theory to connect the electrostatic solvation energy and the PoissonBoltzmann equation. The inconsistency between their definition of the electrostatic force and the sharp dielectric interface model was recently reconciled through the introduction of a shape derivative of the solute biomolecules and the application of the HadamardZolésio structure theorem LiB2011a ; MikuckiM2014a . We will review the regularization techniques for the numerical solution of the PBE, which also sheds light to our choice of the enriching basis functions in the least square reconstruction for the gradient recovery.
2.1 Biomolecular electrostatics theory and the PoissonBoltzmann equation
Consider be a domain that encapsulates the solute biomolecules and the aqueous solution in which species of mobile ions disperse. The occupation domains of solute and solvent are respectively denoted by and , with distinct dielectric constants and . Denote by the dielectric interface separating these two subdomains. Assume dispersion of mobile ions follows the Boltzmann distribution, then the electrostatic free energy of this solvated system is given by
(4) 
where is the bulk concentration of species of ions, is the collection of singular charges of biomolecule atoms, and is the inverse thermal energy at temperature
. The characteristic function
in and vanishes elsewhere. Equation (4) highlights the dependence of the electrostatic free energy on the location of the dielectric interface . Minimization of with respect to leads to the nonlinear PoissonBoltzmann equation(5) 
with interface conditions on :
(6) 
where denotes the jump of the enclosed function and is the normal derivative pointing to .
To derive the dielectric boundary force we introduce a surface transformation corresponding to the slight change of domain due to the motion of interface in its normal direction. The dielectric boundary force is derived by taking the variational derivative of the free energy with respect to a velocity field normal to . Associated with this velocity fields is a mapping from the original surface coordinate (material position) to the deformed coordinate (physical position):
(7) 
Introducing the Jacobian of the surface transformation :
(8) 
one can show that
(9) 
and the shape derivative of
(10) 
where is the dielectric boundary force given in (1) MikuckiM2014a .
Numerical techniques for the PBE are mostly focused on the treatment of the discontinuous dielectric function and the singular charge density . There are three major regularization schemes with which a direct approximation of the singular delta functions can be avoided:

Subtract the potential
(11) induced by the collection of singular charges from the potential in the entire domain , and thus one only needs to numerically solve the remaining regular potential ChenL2007a . In solvent domain , this potential is much larger than the true potential , so is the regular potential there. Consequently a small relative error in the numerical solution of will present a larger relative error in , constituting an unstable numerical algorithm HolstM2012a .

Subtract the singular potential (11) only in the solute domain . This will create additional jumps in the remaining regular components. Therefore a harmonic potential defined by
(12) is added back to partially compensate the jump and to maintain continuity of the remaining regular potential at the interface ChernY2003a . This strategy greatly improves the stability of the numerical method. On the other hand, one needs to compute on to supply the interface conditions for :
(13) Numerical evaluation of is nontrivial, and might introduce considerable error to the interface condition for . A numerical DirichletNeumann mapping based on boundary integral formulations will make this calculation accurate and efficient ZhouY2007b .

One could choose not to add back the harmonic potential so there are only two components in the solute domain. The absence of this harmonic component leads to the jumps in the potential and its normal derivative on the dielectric interface
(14) On the other hand, this regularization scheme helps reduce the complexity and improve the accuracy as one does not need to supply a numerical computed to the related interface condition for the regularized PoissonBoltzmann equation Geng2017a .
Our enriched gradient recovery method is compatible with all these regularization schemes. We shall present the gradient recovery and force calculation along with the following general form of interface conditions:
(15) 
with known jumps .
2.2 Dielectric boundary force and its calculation
To get the dielectric boundary force , one shall first compute and their respective gradients at selected positions (the centroids of surface triangles, for example) on the interface . There are there different procedures to complete this calculation:

Independently compute and their gradients using proper gradient recovery techniques.

Compute only one set of or , and use the interface conditions (15) to get the other.

Compute both sets of and by enforcing the interface conditions (15) in the gradient recovery.
Procedure I is inferior to the others because the resulting dielectric boundary force is not consistent with the interface conditions. Procedures II and III will be developed in Section 3 and implemented in Section 4
3 Gradient Recovery Techniques Enriched with Green’s Functions
Our gradient recovery techniques for PBE are based on the least square reconstruction of the numerical solution, which is also the basis of polynomial preserving recovery (PPR) techniques NagaA2005a ; ZhangZ2005a and their recent variations for elliptic interface problems GuoH2017a ; GuoH2018a ; GuoH2018b . In this section we will first summarize the standard PPR method and then present its enrichment with Green’s functions. This enriched gradient recovery technique will be integrated with different regularization schemes described in Section 2 to enforce the respective interface conditions.
3.1 Polynomial preserving recovery and enrichment with Green’s functions
We consider the molecular surface discretized by triangles and assume that the electrostatic potential has been solved from the PoissonBoltzmann equation using a finite difference or finite element interface method. The mesh for the interface method may not be fitting the interface triangulation. The electrostatic force is usually computed at the centroids of triangles which are not necessarily the mesh nodes. Our method can be used to calculate the solution derivatives at vertices of surface triangles but one needs a delicate definition of the normal directions at the corner vertices to finally compute the surface force. For a centroid with coordinate we define be the collection of mesh nodes in solvent domain that are closest to :
(16) 
where are the mesh nodes in or on the molecular surface . The numerical solution at this collection of nodes will be fitted by a polynomial in the least square sense sampled at these nodes:
(17) 
We then define the recovered gradient at surface position as
(18) 
Alternatively, one can choose the collection of nodes in the solute domain to fit a polynomial and to compute the recovered gradient
(19) 
As in most analysis and applications of PPR, we choose to be the space of quadratical polynomials, i.e., .
Remark 1.
The choice of sampling nodes (16) is for general interface methods. For finite element interface methods one can take advantage of mesh structure to choose the mesh nodes in the first layers around the point .
As we will demonstrated in Section 4, the polynomial preserving recovery technique is not able to deliver an accurate surface gradient on grid sizes affordable to practical biomolecular simulations, mostly because polynomial approximation of singular potential components have large truncation errors. Since these singular components are induced by charges at known positions, we will introduce the Green’s function as basis functions in addition to the quadratical polynomials in the least square reconstruction, c.f. Fig. 1. We anticipate that the singular component of the solution will be approximated by these enriching basis functions, while the regular component of the solution will be approximated in the original polynomial space. For the surface point given above, we shall choose charged atoms of the solute molecule that are closest to :
(20) 
where is the set of all charged atoms of the solute molecule. Alternatively, one may choose all charged atoms within a predetermined distance to :
(21) 
Since the singular components decay like , contributions to the potential and force from charged atoms far away from are secondary, and can be well approximated in the original polynomial space. By enriching the polynomial basis of degree with Green’s functions centered at a total of selected charged atoms we define a new space of functions
(22) 
where are the Cartesian coordinates with respect to and is the distance to selected charge. We then fit the numerical solution at mesh nodes in with a function such that
(23) 
This amounts to a least square solution where each row of matrix is corresponding to a sampling mesh node in with coordinate and is given by
(24)  
(25) 
where and is the position of the selected charge. Since the functions in (25) are linearly independent with each other and with the polynomial basis in (24), the least square problem (23) is solvable.
Remark 2.
The potential in does not have a decomposition with an analytical reaction field similar to (2). However, the approximate solution (3) suggests that it is possible to define a component of the form . The space is constructed to approximate this component using and the remaining component using the regular polynomial basis.
Remark 3.
The derivative of the Green’s funtion can also be written in terms of the same Green’s function. Thus the derivative of the basis functions in is also in , allowing us to approximate both the potential and its derivatives in the same space.
Remark 4.
The approximation of the electrostatic potential in for the Born ion RouxB1990a is exact.
3.2 Coupled gradient recovery with interface condition enforcement
We need potentials and potential gradients on both sides of the dielectric interface to compute the dielectric boundary force. These quantities independently recovered on either side using the method described above do not satisfy the interface conditions, although these conditions have been enforced in the interface methods for the numerical solution of the PBE. Procedure II (c.f. Section (2.2)) recovers the solution on one side and generate the surface potential and surface gradient on the other side using the interface conditions. This requires only one least square solution for the gradient recovery at an interface point , and thus is simpler to implement. It is worth noting that if one chooses to compute instead of , then the enrichment using Green’s function is not necessary, because the singular potential component has been removed from in all three regularization schemes here. Coupled gradient recovery on both sides of the interface is also possible by enforcing the interface conditions in least square reconstruction. To this end we shall define in addition to another collection of mesh nodes in the solute domain, , that are closest to :
(26) 
Two functions, one polynomial and the other , respectively sampled on and , will be sought from the following problem:
(27) 
subject to the interface conditions:
(28) 
The solution of (27) leads to a least square problem where the matrix whose first rows correspond to the approximation of the solution respectively at sampling nodes in , i.e.,
(29)  
(30) 
for . The next rows of are for the approximation at sampling nodes in , i.e.,
(31)  
(32)  
(33) 
for . The last two rows of come from the approximation of the interface conditions (28):
(34)  
(35)  
(36)  
(37)  
(38)  
(39)  
(40)  
(41) 
where is the unit outer normal on the interface. Solution of this least square problem reconstructs functions that approximate the potential solutions at sampling mesh nodes on two sides of the interface and the jump conditions at the interface point .
4 Numerical Experiments
In this section we validate the accuracy and stability of our enriched gradient recovery method on a set of biomolecules. Since there do not exist exact solutions of the electrostatic surface potential and gradient for general molecules with complex surface geometries, we make up potential fields for accuracy validation of Procedure II by setting the potential in the solvent domain to be
(42) 
following the approximation boundary condition (3). The exponential screen term is neglected because the ion strength is alway small and makes very slight changes to the surface potential and gradient. Procedure II will be used in enriched recovery with analytical fields. We will compare the classical polynomial preserving (PPR) and enriched gradient recoveries on molecules with the potential field (42). We will further solve the PoissonBoltzmann equation with an interface finite difference method MIB ZhouY2006a ; GengW2007a , and use the solution at a fine mesh as the reference to check the accuracy of both PPR and our enriched gradient recoveries. Procedure III will be used since we have solutions available on both sides of the interface. We will use the following absolute and relative errors for quantifying the accuracy of recovered surface potential , gradient , and force on the interface:
4.0.1 Born ion
The potential field (42) is exact for Born ion so our enriched gradient recovery method will give result with machine error, for that is one of the basis function of . For comparison, we present in Table (1) the errors in surface potential and gradients generated by PPR. It is seen that PPR has proved optimal rate of convergence and can give fairly good results when the mesh is sufficiently refined. In particular, a mesh space smaller than Å is needed to generate a surface potential gradient with a relative error less than , an accuracy barely sufficient for the real biomolecular simulations. The surface potential and magnitude of the surface gradient recovered using PPR at Å are shown in Fig. 2.
0.4  1.33  9.95e2  7.97  1.18 

0.2  2.10e1  1.45e2  3.02  2.58e1 
0.1  4.15e2  2.83e3  1.05  7.69e2 
0.05  6.69e3  4.55e4  3.16e1  2.19e2 
4.0.2 Diatomic model
Our diatomic model consists of two unit spheres centered respectively at , each carrying a unit positive charge at its center. We first set an analytical potential field (42) in and recover the potential and potential gradient on surface using PPR and enriched recovery. If both charges are included for enrichment then recovered results will be accurate with machine error. Instead we choose the charge nearest to a surface point for the recovery there. It is seen from Table (2) that the enrichment with a single Green’s function significantly improves the accuracy of the gradient recovery. At Å. the relative error in surface gradient from enriched recovery is less than while the traditional PPR method gives a surface gradient with an relative error more than . This significant difference is also highlighted in Fig. 3, where a smoother and symmetric surface gradient is delivered by the enriched recovery.
PPR  Enriched Recovery  

0.4  1.35  4.80e2  9.60  8.17e1  1.91e1  6.44e3  1.46  7.78e2 
0.2  2.77e1  9.51e3  3.98  2.31e1  1.73e2  5.82e4  2.68e1  1.24e2 
0.1  5.57e2  1.90e3  1.39  7.08e2  3.49e3  1.18e4  8.90e2  4.13e3 
0.05  9.12e3  3.11e4  4.23e1  2.06e2  5.48e4  1.85e5  2.69e2  1.25e3 
Th enriched recovery is also tested and compared to PPR on the numerical solutions of the PBE for this diatomic molecule. Since there does not exist analytical values of the surface potential and potential gradient, we would use the results of the enriched recovery at a fine mesh as the reference solution for error measurement. Both atoms are used in the enrichment for optimal recovery. The improvement of the accuracy in the recovered variables due to the enrichment is evidenced by the results collected in Table 3. This table and the plots in Fig. 4 show that both the PPR and enriched recovery can deliver accurate surface potential (top row) from the numerical solutions of the PBE but the accuracies of recovered gradients differ significantly. For all mesh sizes we tested the enriched recovery generates surface gradients with a consistent range. In contrast, the results of PPR vary considerably as the mesh is refined. At a fine mesh with Å , notable variation is still observed in the surface gradient at the center neck of this diatomic molecule. The coupled enriched recovery is also tested on the numerical solutions of the PBE and the results are rather similar to the recovery based on the solution in only.
PPR  Enriched Recovery  

0.4  5.46e1  1.86e2  9.51  7.72e1  8.23e1  3.43e2  1.08  7.22e2 
0.2  1.86e1  6.32e3  4.37  2.53e1  2.76e1  9.54e3  4.72e1  2.05e2 
0.1  5.71e2  1.94e3  1.54  7.74e2  2.97e2  9.89e4  2.09e1  9.18e3 
0.05  9.76e3  3.31e4  4.58e1  2.20e2 
4.0.3 Amino Acids
The numerical experiments above show that the accuracy of recovered gradients can be increased by enriching Green’s functions at all charged atoms. For real biomolecular simulations we might not be able to afford this global enrichment because of the increasing computational cost. Global enrichment might not be necessary because the Green’s functions decay in , which indicates the potential is prone to be dominated by the charges nearby. It is therefore possible only to choose a small set of charged atoms close to the surface point of recovery. Here we will test the enriched recovery on six amino acids: TRY, ASN, ASP, ARG, LYS and GLU dimer. For each amino acid, we set an analytical potential field using (42) assuming a unit positive charge at each atom.
As seen in Tables (4,5), the classical PPR () has to use a very fine mesh Åto generate a surface gradient with a relative error about . Recovery enriched with Green’s functions at nearest atoms is able to deliver results of the similar accuracy at Å. This number is also identified on the results for other residues, which are neglected as they show the very similar behavior as in these two tables. Enrichment with more Green’s functions will further improve the accuracy of the recovery but for the cause of efficiency we will choose in the numerical experiments below.
0  1.16e1  8.00e1  2.53e2  2.46e1  7.12e4  2.10e2 

1  3.66e2  1.62e1  7.06e3  5.68e2  1.86e4  5.49e3 
3  1.05e2  4.35e2  1.41e3  1.17e2  3.87e5  1.17e3 
5  5.20e3  2.22e2  6.87e4  5.44e3  2.00e5  5.71e4 
8  2.28e3  9.67e3  1.81e4  1.47e3  2.02e5  5.19e4 
11  8.80e4  3.62e3  5.90e5  4.99e4  5.70e6  1.39e4 
15  3.13e4  1.35e3  1.49e5  1.23e4  1.08e6  2.84e5 
20  1.70e5  9.04e5  5.80e7  5.17e6  5.25e7  1.32e5 
0  1.05e1  7.89e1  2.34e2  2.45e1  6.56e4  2.09e2 

1  3.94e2  1.79e1  8.25e3  6.94e2  2.07e4  6.44e3 
3  1.12e2  4.86e2  1.88e3  1.50e2  3.63e5  1.16e3 
5  5.20e3  2.22e2  6.85e4  5.57e3  1.41e5  4.28e4 
7  2.56e3  1.11e2  4.27e4  3.19e3  2.02e5  5.97e4 
9  1.09e3  4.72e3  1.34e4  9.83e4  1.03e5  2.89e4 
12  1.80e4  8.52e4  1.73e5  1.38e4  2.26e6  6.23e5 
We shall also compare PPR with enriched recovery on numerical solutions of the PBE for these amino acids. The surface potential and gradients recovered with the global enrichment at Å are taken as the reference solutions to examine the accuracy of recovery scheme enriched with different number of Green’s function and different . The plots for ARG in Fig. 5 and for GLU dimer in Fig. 6 show that surface potential recovered by the classical PPR has an error more than with respect to enriched recovery on a coarse mesh with Å; the surface gradient delivered by PPR is indeed errant on the central neck of the GLU dimer. Despite the salient errors in the surface quantities recovered by PPR on coarse meshes, these results converge to the same limit as the enriched recovery when the mesh of the numerical solutions for the PBE is sufficiently refined. In contrast, the enriched recovery produces very consistent surface potentials and surface gradients on all meshes tested, c.f., the ranges of these surface quantities in these plots.
4.0.4 Proteins
Our last numerical experiments are on a lowdensity lipoprotein receptor 1AJJ which has 37 residues and 519 atoms. The accuracy of the recovered surface potential and gradient from the numerical solutions of the PBE is documented in Tab. 6, where the results of the enriched recovery on the uniform mesh with Å is used as the reference. This table provides evidence of the higher accuracy of the enriched recovery over the classical PPR. For Å, a relative error over is found in the surface gradient recovered by PPR. This error is reduced to about on a very fine mesh with Å, which is still larger than the enriched recovery on a coarse mesh with Å. The accuracy improvement due to the enrichment is also indicated in Fig. 7. Consistent with the observations on model molecules and amino acids above, enriched recovery produces more stable ranges of surface potential and gradient than PPR.
PPR  Enriched Recovery  

0.4  1.61  1.87e2  11.0  4.65e1  1.52e1  1.76e3  9.49e1  3.03e2 
0.2  3.59e1  4.15e3  4.50  1.60e1  1.73e2  2.00e4  2.18e1  6.96e3 
0.1  6.31e2  7.30e4  1.48  4.90e2  2.70e3  3.13e5  6.28e2  2.01e3 
5 Concluding Remarks
We present in this paper a novel numerical method for computing the potential and its gradient on the molecular surfaces from the numerical solutions of the PoissonBoltzmann equation, which is widely used in modeling the electrostatic interactions for biomolecules in solvent. Our method reconstructs a solution in the least square sense locally in the solvent region. In contrast to the classical polynomial preserving recovery (PPR) method that is based only on the polynomial basis functions, our reconstruction enriches the basis with Green’s functions modeling the electrostatic potential field induced by selected charges. This enrichment is motivated by the Generalized Born method which decomposes the electrostatic potential inside the solute molecules into a Coulomb potential and a reaction field, and also by the approximate boundary condition for the PoissonBoltzmann equation. Extensive numerical methods on Born atom, model diatomic molecule, amino acids, and proteins demonstrate that our enriched recovery is more accurate than PPR on the same mesh, more stable along the mesh refinement, and indeed provides a stable reference to which the PPR solutions will converge.
Our enriched recovery is easy to implement, and can be readily integrated with other interface finite difference or finite element methods for the PoissonBoltzmann equation. As the reconstruction is independently computed at each surface point of interest, parallelization of this enriched recovery is straightforward. We are currently working to supply the recovered surface gradient to compute the dielectric boundary force and use the force to drive the large scale continuum molecular deformation as induced due to varying solvation states such as binding of molecules or ligands, changes in ion concentration, or protonation.
Acknowledgements
This work has been partially supported by National Institutes of Health through the grant R01GM117593 as part of the joint DMS/NIGMS initiative to support research at the interface of the biological and mathematical sciences.
References
 (1) B. Honig, A. Nicholls, Classical electrostatics in biology and chemistry, Science 268 (5214) (1995) 1144–1149. doi:10.1126/science.7761829.
 (2) W. Im, D. Beglov, B. Roux, Continuum solvation model: Computation of electrostatic forces from numerical solutions to the PoissonBoltzmann equation, Comput. Phys. Commun. 111 (1998) 59–75. doi:10.1016/S00104655(98)000162.
 (3) B. Lu, Y. C. Zhou, M. Holst, J. A. McCammon, Recent progress in numerical solution of the PoissonBoltzmann equation for biophysical applications, Commun. Comput. Phys. 3 (2008) 973–1009.
 (4) P. Ren, J. Chun, D. G. Thomas, M. J. Schnieders, M. Marucho, J. Zhang, N. A. Baker, Biomolecular electrostatics and solvation: a computational perspective 45 (2012) 427–491. doi:10.1017/S003358351200011X.
 (5) D. Xie, Y. Jiang, L. R. Scott, Efficient algorithms for a nonlocal dielectric model for protein in ionic solvent, SIAM J. Sci. Comput. 35 (2013) 1267–1284. doi:10.1137/120899078.
 (6) J. Chen, W. Geng, On preconditioning the treecodeaccelerated boundary integral (tabi) PoissonBoltzmann solver, J. Comput. Phys. 373 (2018) 750–762. doi:10.1016/j.jcp.2018.07.011.
 (7) Y. Zhong, K. Ren, R. Tsai, An implicit boundary integral method for computing electric potential of macromolecules in solvent, J. Comput. Phys. 359 (2018) 199 – 215. doi:10.1016/j.jcp.2018.01.021.
 (8) A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. JosephMcCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. WiórkiewiczKuczera, D. Yin, M. Karplus, Allatom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B 102 (18) (1998) 3586–3616. arXiv:https://doi.org/10.1021/jp973084f, doi:10.1021/jp973084f.
 (9) Q. Lu, R. Luo, A PoissonBoltzmann dynamics method with nonperiodic boundary condition, J. Chem. Phys. 119 (21) (2003) 11035–11047.
 (10) N. V. Prabhu, P. J. Zhu, K. A. Sharp, Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smoothpermittivity finite difference PoissonBoltzmannmethod, J. Comput. Chem. 25 (16) (2004) 2049–2064. doi:10.1002/jcc.20138.
 (11) W. Geng, G. W. Wei, Multiscale molecular dynamics using the matched interface and boundary method, J. Comput. Phys. 230 (2011) 435–457. doi:10.1016/j.jcp.2010.09.031.
 (12) Y. Tang, G. Cao, X. Chen, J. Yoo, A. Yethiraj, Q. Cui, A finite element framework for studying the mechanical response of macromolecules: application to the gating of machanosensitive channel MscL, Biophys. J. 91 (2006) 1248–1263. doi:10.1529/biophysj.106.085985.
 (13) L. Ma, A. Yethiraj, X. Chen, Q. Cui, A computational framework for mechanical response of macromolecules: applications to the salt concentration dependence of DNA bendability, Biophys. J. 96 (2009) 3542–3554.
 (14) X. Chen, Q. Cui, Computational molecular biomechanics: A hierarchical multiscale framework with applications to gating of mechanosensitive channels of large conductance, in: T. Dumitrica (Ed.), Trends in Computational Nanomechanics, Vol. 9 of Challenges and Advances in Computational Chemistry and Physics, Springer Netherlands, 2010, pp. 535–556.
 (15) Y. C. Zhou, B. Lu, A. A. Gorfe, Continuum electromechanical modeling of proteinmembrane interactions, Phys. Rev. E 82 (4) (2010) 041923. doi:10.1103/PhysRevE.82.041923.
 (16) L.T. Cheng, B. Li, M. White, S. Zhou, Motion of a cylindrical dielectric boundary, SIAM J. Appl. Math. 73 (2013) 594–616. doi:10.1137/120867986.
 (17) M. K. Gilson, M. E. Davis, B. A. Luty, J. A. McCammon, Computation of electrostatic forces on solvated molecules using the PoissonBoltzmann equation, J. Phys. Chem. 97 (14) (1993) 3591–3600. doi:10.1021/j100116a025.
 (18) Y. Zhou, M. Holst, J. A. McCammon, Nonlinear elastic modeling of macromolecular conformational change induced by electrostatic forces, J. Math. Anal. Appl. 340 (2008) 135–164. doi:10.1016/j.jmaa.2007.07.084.
 (19) B. Li, X. Cheng, Z. Zhou, Dielectric boundary force in molecular solvation with the PoissonBoltzmann free energy: A shape derivative approach, SIAM J. Appl. Math. 71 (2011) 2093–2111. doi:10.1137/110826436.
 (20) M. Mikucki, Y. C. Zhou, Electrostatic forces on charged surfaces of bilayer lipid membranes, SIAM J. Appl. Math. 74 (2014) 1–21. doi:10.1137/130904600.
 (21) J. Frank H. Stillinger, Interfacial solutions of the poissonboltzmann equation, J. Chem. Phys. 35 (5) (1961) 1584–1589. doi:10.1063/1.1732113.
 (22) Z. X. Xiang, Y. Y. Shi, Y. W. Xu, Solving the finitedifference, nonlinear, PoissonBoltzmann equation under a linearapproach, J. Comput. Chem. 16 (2) (1995) 200–206. doi:10.1002/jcc.540160207.
 (23) M. Holst, F. Saied, Numerical solution of the nonlinear PoissonBoltzmann equation: Developing more robust and efficient methods., J. Comput. Chem. 16 (1995) 337–364. doi:10.1002/jcc.540160308.
 (24) M. Holst, N. A. Baker, F. Wang, Adaptive multilevel finite element solution of the PoissonBoltzmann equation I: algorithms and examples, J. Comput. Chem. 21 (2000) 1319–1342.
 (25) A. H. Boschitsch, M. O. Fenley, H.X. Zhou, Fast boundary element method for the linear PoissonBoltzmann equation, J. Phys. Chem. B 106 (10) (2002) 2741–2754. doi:10.1021/jp013607q.
 (26) D. Xie, S. Zhou, A new minimization protocol for solving nonlinear PoissonBoltzmann mortar finite element equation, BIT 47 (2007) 853–871, in press. doi:10.1007/s1054300701459.
 (27) H. Lin, H. Tang, W. Cai, Accuracy and efficiency in computing electrostatic potential for an ion channel model in layered dielectric/electrolyte media, J. Comput. Phys. 259 (2014) 488 – 512. doi:10.1016/j.jcp.2013.12.017.
 (28) W. Geng, S. Zhao, A twocomponent matched interface and boundary (MIB) regularization for charge singularity in implicit solvation, J. Comput. Phys. 351 (2017) 25 – 39. doi:10.1016/j.jcp.2017.09.026.
 (29) B. Z. Lu, X. L. Cheng, J. F. Huang, J. A. McCammon, Order algorithm for computation of electrostatic interactions in biomolecular systems, Proc. Natl. Acad. Sci. U. S. A. 103 (51) (2006) 19314–19319. doi:10.1073/pnas.0605166103.
 (30) C. Bajaj, S. Chen, A. Rand, An efficient higherorder fast multipole boundary element solution for PoissonBoltzmannbased molecular electrostatics, SIAM J. Sci. Comput. 33 (2) (2011) 826–848. doi:10.1137/090764645.
 (31) W. Geng, R. Krasny, A treecodeaccelerated boundary integral PoissonBoltzmann solver for electrostatics of solvated biomolecules, J. Comput. Phys. 247 (2013) 62–78. doi:10.1016/j.jcp.2013.03.056.
 (32) H. Guo, X. Yang, Gradient recovery for elliptic interface problem: I. bodyfitted mesh, Commun. Comput. Phys. 23 (2018) 1488–1511. doi:10.4208/cicp.OA20170026.
 (33) H. Guo, X. Yang, Gradient recovery for elliptic interface problem: Ii. immersed finite element methods, J. Comput. Phys. 338 (2017) 606 – 619. doi:10.1016/j.jcp.2017.03.003.
 (34) H. Guo, X. Yang, Gradient recovery for elliptic interface problem: III. Nitsche’s method, J. Comput. Phys. 356 (2018) 46 – 63. doi:10.1016/j.jcp.2017.11.031.

(35)
O. C. Zienkiewicz, J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique, Internat. J. Numer. Methods Eng. 33 (1992) 1331–1364.
doi:10.1002/nme.1620330702.  (36) O. C. Zienkiewicz, J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity, Internat. J. Numer. Methods Eng. 33 (1992) 1365–1382. doi:10.1002/nme.1620330703.
 (37) R. Bank, J. Xu, Asymptotically exact a posteriori error estimators, part ii: General unstructured grids, SIAM J. Numer. Anal. 41 (6) (2003) 2313–2332. doi:10.1137/S0036142901398751.
 (38) H. Guo, Z. Zhang, Gradient recovery for the Crouzeix–Raviart element, J. Sci. Comput. 64 (2) (2015) 456–476. doi:10.1007/s1091501499395.
 (39) A. Naga, Z. Zhang, The polynomialpreserving recovery for higher order finite element methods in 2D and 3D 5 (2005) 769–798. doi:10.3934/dcdsb.2005.5.769.
 (40) J. Xu, Z. Zhang, Analysis of recovery type a posteriori error estimators for mildly structured grids, Math. Comp. 247 (2003) 1139–1152. doi:10.1090/S0025571803016004.
 (41) Z. Zhang, A. Naga, A new finite element gradient recovery method: Superconvergence property, SIAM J. Sci. Comput. 26 (4) (2005) 1192–1213. doi:10.1137/S1064827503402837.
 (42) S.H. Chou, An immersed linear finite element method with interface flux capturing recovery 17 (2012) 2343. doi:10.3934/dcdsb.2012.17.2343.
 (43) S.H. Chou, C. Attanayake, Flux recovery and superconvergence of quadratic immersed interface finite elements, Int. J. Numer. Anal. Mod. 14 (2015) 88–102.
 (44) H. Wei, L. Chen, Y. Huang, B. Zheng, Adaptive mesh refinement and superconvergence for twodimensional interface problems, SIAM J. Sci. Comput. 36 (4) (2014) A1478–A1499. doi:10.1137/120866622.
 (45) Z. Li, H. Ji, X. Chen, Accurate solution and gradient computation for elliptic interface problems with variable coefficients, SIAM J. Numer. Anal. 55 (2) (2017) 570–597. doi:10.1137/15M1040244.
 (46) Z. Li, The immersed interface method using a finite element formulation, Appl. Numer. Math. 27 (3) (1998) 253–267. doi:10.1016/S01689274(98)000154.
 (47) X.D. Liu, R. P. Fedkiw, M. Kang, A boundary condition capturing method for poisson’s equation on irregular domains, J. Comput. Phys. 160 (1) (2000) 151–178. doi:10.1006/jcph.2000.6444.
 (48) S. Hou, P. Song, L. Wang, H. Zhao, A weak formulation for solving elliptic interface problems without body fitted grid, J. Comput. Phys. 249 (2013) 80 – 95. doi:10.1016/j.jcp.2013.04.025.
 (49) H. Ji, J. Chen, Z. Li, A symmetric and consistent immersed finite element method for interface problems, J. Sci. Comput. 61 (3) (2014) 533–557. doi:10.1007/s109150149837x.
 (50) P. Hansbo, Nitsche’s method for interface problems in computational mechanics, GAMM  Mitteilungen 28 (2005) 183–206. doi:10.1002/gamm.201490018.
 (51) C. Annavarapu, M. Hautefeuille, J. E. Dolbow, A robust Nitsche’s formulation for interface problems, Comput. Methods Appl. Mech. Eng. 225228 (2012) 44–54. doi:https://doi.org/10.1016/j.cma.2012.03.008.
 (52) W. Im, M. S. Lee, C. L. B. III, Generalized Born model with a simple smoothing function, J. Comput. Chem. 24 (2003) 1691–1702. doi:10.1002/jcc.10321.
 (53) M. Feig, W. Im, C. L. B. III, Implicit solvation based on generalized Born theory in different dielectric environments, J. Comput. Chem. 120 (2004) 903–911. doi:10.1063/1.1631258.
 (54) A. V. Onufriev, D. A. Case, Generalized Born implicit solvent models for biomolecules, Annu. Rev. Biophys. 48 (2019) 275–296. doi:10.1146/annurevbiophys052118115325.
 (55) M. K. Gilson, A. Rashin, R. Fine, B. Honig, On the calculation of electrostatic interactions in proteins, J. Mol. Biol. 184 (3) (1985) 503–516. doi:10.1016/00222836(85)902979.
 (56) K. A. Sharp, B. Honig, Calculating total electrostatic energies with the nonlinear PoissonBoltzmann equation, J. Phys. Chem. 94 (19) (1990) 7684–7692. doi:10.1021/j100382a068.
 (57) K. A. Sharp, B. Honig, Electrostatic interactions in macromolecules  theory and applications, Annu. Rev. Biophys. Biophys. Chem. 19 (1990) 301–332. doi:10.1146/annurev.bb.19.060190.001505.
 (58) L. Chen, M. Holst, J. Xu, The finite element approximation of the nonlinear PoissonBoltzmann equation, SIAM J. Numer. Anal. 45 (2007) 2298–2320. doi:10.1137/060675514.
 (59) M. Holst, J. A. McCammon, Z. Yu, Y. C. Zhou, Y. Zhu, Adaptive finite element modeling techniques for the PoissonBoltzmann equation, Commun. Comput. Phys. 11 (2012) 179–214. doi:10.4208/cicp.081009.130611a.
 (60) I.L. Chern, J.G. Liu, W.C. Wang, Accurate evaluation of electrostatics for macromolecules in solution, Methods Appl. Anal. 10 (2003) 309–328.
 (61) Y. C. Zhou, B. Z. Lu, G. A. Huber, M. J. Holst, J. A. McCammon, Continuum simulations of acetylcholine consumption by acetylcholinesterase  a PoissonNernstPlanck approach, J. Phys. Chem. B 112 (2) (2008) 270–275.
 (62) B. Roux, H.A. Yu, M. Karplus, Molecular basis for the Born model of ion solvation, J. Phys. Chem. 94 (1990) 4683–4688. doi:10.1021/j100374a057.
 (63) Y. C. Zhou, S. Zhao, M. Feig, G. W. Wei, High order matched interface and boundary (MIB) schemes for elliptic equations with discontinuous coefficients and singular sources, J. Comput. Phys. 213 (2006) 1–30.
 (64) W. Geng, S. Yu, G. Wei, Treatment of charge singularities in implicit solvent models, J. Chem. Phys. 127 (2007) 114106.
Comments
There are no comments yet.