Numeric Solution Of The NonLinear PoissonBoltzmann 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 nonlinear
PoissonBoltzmann 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 PoissonBoltzmann (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 f_{i} is the counterion 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 x_{0} such that
u(x_{0}) = 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_{+}(x_{0}) = f_{}(x_{0})
and so
from which we can recover the PoissonBoltzmann equation above by
converting I (given in mol/L) to the number of particles per Angstrom
cubed.
Procedures for numerically solving the PoissonBoltzmann
equation generally proceed as follows:
 Equal concentrations of positive and negative ions are assumed with
identical radii and charge.
 The hyperbolic sine term is linearized (truncated Taylor expansion at the
linear term).
 The fields are represented on a grid with finite differencing replacing
the continuous operators.
 An iterative sparse linear equation solver is used to solve for the
electric field.
This article describes the software program, MOEElectrostatics,
that does not rely on the first two
assumptions; that is, we solve the full nonlinear PoissonBoltzmann 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 PoissonBoltzmann
equation in the following form:
where m denotes the number of ion classes, q_{i} is
the partial charge of each ion in class i, C_{i} is the
concentration (mol/L) in ion class i, and W_{i} 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 nonlinear partial differential equation solver to solve the discretized
PoissonBoltzmann 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 (x_{0}, y_{0}, z_{0})
is the origin of the grid, h > 0 is the grid spacing and
n_{x}, n_{y}, and n_{z} 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(h^{4}). The discrete derivative of
a field u is defined to be (along one axis):
which has an error of O(h^{3}).
Charge Density. With each grid point (i,j,k)
we associate a cube, B, of volume h^{3} 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 a^{2}
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 p_{i} is the position of atom i,
r_{i} 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 Connollytype 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 v_{m} 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 MonteCarlo 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 D_{in} is the interior dielectric constant (e.g., 1 or
4) and D_{out} 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 ionaccessible 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 multigrid 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 NewtonGaussSeidel method was used to solve for
u.
Algorithm Summary. With all of the above quantities defined, the
procedure to solve the nonlinear PoissonBoltzmann equation is summarized as
follows:
 Determine grid spacing h and the bounding box of atoms at positions
p_{i}, with radii r_{i} and partial charges
q_{i}.
 Compute the charge density f.
 Compute the dielectric field d.
 Compute the ion accessible fields K_{i}.
 Assign boundary values (Coulomb or DebyeHuckel).
 Solve the nonlinear PoissonBoltzmann equation for u with a
nonlinear 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 counterions. 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 counterions are present, we note that the potential function
u = (1+r)^{1} (up to a constant) satisfies
the PoissonBoltzmann equation with an ion concentration of
K = 0, dielectric constant d = 1 and
partial charge density
2(1+r)^{3}r^{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 nonlinear
PoissonBoltzmann equation. The MOEElectrostatics program was implemented
in the SVL programming language and made use of the SVL multidimensional
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).
