Journal Articles

Numeric Solution Of The Non-Linear
Poisson-Boltzmann Equation


Paul Labute and Martin Santavy


Chemical Computing Group Inc.
1010 Sherbrooke Street W, Suite 910, Montreal, Quebec, Canada H3A2R7
 

Abstract. This paper describes a software program that solves the non-linear Poisson-Boltzmann equation in a form that allows for multiple ion classes with different concentrations, van der Waals radii and partial charges. The results of several computational experiments show the method to be efficient and accurate.

 

INTRODUCTION

The Poisson-Boltzmann (PB) equation has been studied extensively. The equation, in terms of the electric potential energy field, is

where d is the relative dielectric, u is the electric potential energy field (with a test charge of 1), f is the molecular partial charge density, and fi is the counter-ion partial charge density given by:

where e is the charge of an electron, N is Avogadro's number, I is the ionic charge concentration (mol/L), k is Boltzmann's constant, and T is the temperature (K). Performing the substitution we find that:

The hyperbolic sine term comes into play because of a few assumptions regarding the ionic charge concentration. Consider a molecule, M, in a salt solution. Let f+ denote the partial charge density of the positive ions, f- the partial charge density of the negative ions, q+ the partial charge of a positive ion, q- the partial charge of a negative ion, n+ the number of positive ions, and n- the number of negative ions. Let W denote the region of space containing the ions.

Suppose that the total number of particles, the volume and the temperature T remain fixed. If u denotes the total electric potential energy, we can then take:

by assuming the Boltzmann distribution of states. Suppose that there exists a point x0 such that u(x0) = 0 (e.g., far enough from the solute molecule), then we can write:

If we assume that q+ = 1, q- = 1 and that there are equal numbers of ions then the ion density will be such that f+(x0) = -f-(x0) and so

from which we can recover the Poisson-Boltzmann equation above by converting I (given in mol/L) to the number of particles per Angstrom cubed.

Procedures for numerically solving the Poisson-Boltzmann equation generally proceed as follows:

  1. Equal concentrations of positive and negative ions are assumed with identical radii and charge.
  2. The hyperbolic sine term is linearized (truncated Taylor expansion at the linear term).
  3. The fields are represented on a grid with finite differencing replacing the continuous operators.
  4. An iterative sparse linear equation solver is used to solve for the electric field.

This article describes the software program, MOE-Electrostatics, that does not rely on the first two assumptions; that is, we solve the full non-linear Poisson-Boltzmann equation allowing for different classes of ions, radii and partial charges. The results of several computational experiments are presented below, proving the method to be efficient and accurate.

METHODS

We shall now present the methods used to solve the Poisson-Boltzmann equation in the following form:

where m denotes the number of ion classes, qi is the partial charge of each ion in class i, Ci is the concentration (mol/L) in ion class i, and Wi is the accessible region of the ions in class i. Our fundamental method is to represent each field in the equation by a uniform grid of function values and use a non-linear partial differential equation solver to solve the discretized Poisson-Boltzmann equation (or equivalent) and thus obtain the potential field u.

Discrete Fields and Operators. All fields are represented by a finite collection of function values sampled on a regular lattice in three dimensions. More precisely, a field v is represented by the collection of function values:

where (x0, y0, z0) is the origin of the grid, > 0 is the grid spacing and nx, ny, and nz are the dimensions of the grid in the x, y and z directions. In this article, we use the symbol v to represent both the continuous function v and the grid representation of v. It will be clear from the context which of the quantities the symbol v represents. We also make use of the notation:

We define the discrete (finite difference) Laplace operator, L, of a discrete field u at a point (i,j,k), denoted by (Lu)ijk, to be:

which has an error of O(h4). The discrete derivative of a field u is defined to be (along one axis):

which has an error of O(h3).

Charge Density. With each grid point (i,j,k) we associate a cube, B, of volume h3 centered at the grid point. Consider a charge density, s. The contribution to B is defined to be:

We shall assume that the density s is separable (i.e., the foregoing integral can be written as a product of three integrals along each axis). With a separable density we need only reason in one dimension since the other dimensions will be similar and the final contribution is the product of the contributions of each dimension.

We model a delta function centered at a point x with an interval of length h centered at x. In this case the contribution to B becomes:

which distributes the point charge to each of the two nearest neighboring grid points. For a collection of point charges, the solute density under this "delta" model becomes the sum of charge scaled discrete delta functions. These contributions can be efficiently calculated with:

This model generalizes to other densities. The delta function can be replaced with a Gaussian density of a specified variance a2 centered at z. In such a case the contribution (in one dimension) becomes:

By a change of variables we can simplify this to:

This "gauss" model distributes the integrated portions of a Gaussian density onto the grid volumes.

Solute Region. The solute region is used to define both the dielectric field and the ion accessible regions. Analytically, the solute density is a kind of step function that is 1 in the interior of the solute and 0 outside. In this way, the integral of the solute density over all of space is the volume of the solute. We define the solute to be all points that are within a distance D of the van der Waals surface of an atom in the solute.

We approximate a solute volume with a Gaussian form:

where pi is the position of atom i, ri is its radius, w is a probe radius and a is a smoothing factor. A contour at the value 1 of the sum is used to define a molecular surface, the molecular volume consists of all points for which the sum exceeds 1. These Gaussian volumes can be used to approximate solvent accessible volumes and Connolly-type volumes. For example, a water accessible surface uses the parameters w=1.4 and a=2.5 while a Connolly surface can be approximated with w=0 and a=1.2. To construct a solute grid we average several shifted gaussian forms:

where vm is a shift. In the present work we use R=8 shifts representing the corners of a box each of which is 0.25h from the center of the grid. More shifts will result in a more accurate estimate of the portion of each grid cell covered by the solute. Alternatively, random shifts can be used in a Monte-Carlo integration for yet more accuracy.

The dielectric field is constructed from a parameter D, the dielectric offset, i.e., the distance from the solute surface at which the dielectric starts becoming that of the solvent. The dielectric field is defined to be:

where Din is the interior dielectric constant (e.g., 1 or 4) and Dout is the solvent dielectric (e.g., 80 for water).

If R is the radius of an ion class of concentration C with partial charge q then we define the ion concentration field to be

This corresponds to an ion-accessible region multiplied by the concentration data.

Equation Solver. Consider the solution of the following partial differential equation:

Given an initial estimate u, there is some correction v to u that solves the equation; that is, there exists some v such that:

which has the same form as the original equation, namely:

This relation was used as the basis for the nonlinear multi-grid algorithm because we can solve for v on a coarser grid and then apply an interpolated correction to u. The recursion was stopped at some suitably coarse grid where a Newton-Gauss-Seidel method was used to solve for u.

Algorithm Summary. With all of the above quantities defined, the procedure to solve the non-linear Poisson-Boltzmann equation is summarized as follows:

  1. Determine grid spacing h and the bounding box of atoms at positions pi, with radii ri and partial charges qi.
  2. Compute the charge density f.
  3. Compute the dielectric field d.
  4. Compute the ion accessible fields Ki.
  5. Assign boundary values (Coulomb or Debye-Huckel).
  6. Solve the nonlinear Poisson-Boltzmann equation for u with a non-linear multigrid partial differential equation solver.

Materials and Software. The foregoing techniques were implemented in the SVL programming language of the Molecular Operating Environment (MOE) version 1998.10 from Chemical Computing Group Inc. All calculations were performed on a 133 MHz Intel Pentium processor running Windows95. Nonlinear optimization was carried out with the MOE Truncated Newton optimizer preceded by one step of Steepest Descent and terminated when the RMS gradient fell below 0.0001. The MOE implementation of the MMFF94 forcefield was used where structure refinement was performed.

RESULTS AND DISCUSSION

We shall describe a number of computational experiments designed to test the validity of the computed electrostatic potential fields. We began our experiments by restricting our attention to the Poisson equation, i.e., without counter-ions. These tests were designed to test the program's ability to solve the Poisson equation in cases where the analytical solution is known. We then proceeded to more complex cases, which use all of the terms.

1. To test the program's ability to solve the Poisson equation in a system with a homogeneous dielectric, we used the fact that the function r-1 with the dielectric 1 is the potential for the delta function density centered at the origin (up to a constant). We set the center of the charge density at (../2,../2,../2) on a 17x17x17 grid with spacing h = 0.5. Five iterations were required to solve the equation so that the maximum absolute residual was less than 10-6. The calculation required 5 seconds of CPU time. The maximum relative error between the analytic potential and the calculated potential was 2%, and the average relative error was 0.14%. Not surprisingly, the largest errors were nearest to the center of the charge (the potential is infinite at that point).

2. An another test of the program's ability to solve the Poisson equation in a system with a homogeneous dielectric, we used the fact that the function erf(../sqrt(2))r-1 with the dielectric 1 is the potential for a Gaussian density centered at the origin (up to a constant). We set the center of the charge density at (../2,../2,../2) on a 17x17x17 grid with spacing h = 0.5. Five iterations were required to solve the equation so that the maximum absolute residual was less than 10-6. The calculation required 5 seconds of CPU time. The maximum relative error between the analytic potential and the calculated potential was 0.4%, and the average relative error was 0.04%.

3. To test the program's ability to solve the Poisson equation in a system with an inhomogeneous dielectric, we noted that the function (1+r)-1 with the dielectric (1+r) is the potential for the density (1+r)-2(r+2)r-1 (up to a constant). We set the center of the charge density at (../2,../2,../2) on a 17x17x17 grid with spacing h = 0.5. Seven iterations were required to solve the equation so that the maximum absolute residual was less than 10-6. The calculation required 6 seconds of CPU time. The maximum relative error between the analytic potential and the calculated potential was 0.65%, and the average relative error was 0.055%.

4. For a more complex check of the homogeneous dielectric solution, we used MOE to build and energy minimize the aspirin molecule:

Using a dielectric constant of d = 4, MMFF94 partial charges modeled with discrete delta functions, a grid of dimensions 33x33x33 at a grid spacing of h = 0.5 Angstroms, six iterations were required to bring the maximum absolute residual below 10-6. The calculation required 46 seconds of CPU time. Figure 1 depicts a +7, -7 kcal/mol isosurface of the computed potential and analytical potential.

Figure 1. Isosurfaces of calculated electrostatic potentials (solid) and analytic potentials (mesh) for the aspirin molecule at 0.5A grid spacing.

We observed that the surfaces for both the calculated and analytical potentials were in remarkable agreement. These results show that the discretization of the partial charge density is not expected to be a large source of errors in such electrostatic calculations. Moreover, reasonably accurate electrostatic potentials can be calculated from a grid spacing as large as 0.5 Angstroms.

5. A spherical charge density of radius a in a homogeneous dielectric d has a total electrostatic energy equal to:

where q is the total partial charge on the surface of the sphere. The difference between two such energies evaluated with different dielectrics results in the Born equation, which states that the free energy difference upon transfer of an ion from one dielectric medium to another can be written as:

We used a 33x33x33 grid with spacing h = 0.25. A spherical charge density was created centered at the origin with radius 2 by performing a Monte Carlo volume integration with radius 2 and subtracting a Monte Carlo volume integration with radius 1.9 (both used 100 probes and resulted in a relative error of less than 0.5% on the volume). We chose q = 1 and d = 1. The calculation required 6 steps to reduce the maximum absolute gradient to less than 10-6 and took 50 seconds of CPU time. The resulting electrostatic energy was -82.1583 as compared to the analytical result of -83.0147 (a relative error of 1%). The same calculation was performed for a dielectric constant of d = 80. The computed energy was -1.0496 as compared to the analytical result of -1.0377 (a relative error of 1.1%). The resulting computed free energy transfer was therefore (-82.1583 + 1.0377) = -81.1087 as compared to the Born equation result of -81.977 (a relative error of 1%).

6. To test the ability of the program to solve for the potential in which counter-ions are present, we note that the potential function u = (1+r)-1 (up to a constant) satisfies the Poisson-Boltzmann equation with an ion concentration of K = 0, dielectric constant d = 1 and partial charge density 2(1+r)-3r-1 - Kexp{-../kT}. We used a 17x17x17 grid with spacing h = 0.1. The calculation required 8 iterations to reduce the residual below 10-6 and took 9 seconds of CPU time. The maximum absolute relative error between the calculated potential and the analytical potential was 0.55%, and the average relative error was 0.04%.

CONCLUSIONS

We have presented a numerical method to solve the full non-linear Poisson-Boltzmann equation. The MOE-Electrostatics program was implemented in the SVL programming language and made use of the SVL multi-dimensional grid technology. Unlike other programs, the method supports different ion classes each of which can have different radii, partial charge and concentration. The results of computational experiments suggest that the method is efficient and accurate. The software was written in the SVL programming language and is available under license from Chemical Computing Group Inc. in the Molecular Operating Environment (MOE).