8.1 Introduction to Poisson-Boltzmann Electrostatics

One of the most successful implicit models for the treatment of electrostatic effects is the Poisson-Boltzmann equation¹ which is given below

grad (epsilon(x) grad phi(x)) - kappa(x)^2*sinh(phi(x)e/kT) = -4 pi k_e rho(x) where

     phi is the electrostatic potential,
     rho is the charge density,
     epsilon(x) is the dielectric constant in units where the vacuum
     dielectric is exactly 1,
     k_e is the electrostatic force constant,
     e is the charge on the electron,
     k is the Boltzmann constant,
     T is the temperature,
     and
     kappa is a dielectric independent Debye-Huckel parameter,
     kappa = sqrt(epsilon) * kappa

The equation can be solved over a volume enclosing a molecule of interest using finite difference methods^{2 3 4} appropriate for the solution of boundary value problems. In these methods, a grid is laid down over the volume, and each grid point is assigned a value for the ionic strength and charge density. The lines between each grid point are assigned values for dielectric constant. Boundary values for the potential are set according to a variety of models and then, interior values for the electrostatic potential are calculated. Typically, the equation is linearized by substituting x for sinh(x) in the above equation.

The implementation of the finite difference solution of the Poisson-Boltzmann equation in CONGEN is written in C and uses dynamic storage allocation throughout so any size grid can be accommodated. The process of a potential calculation begins with the determination of the grid position. This can be centered on the spatial origin, the center of a rectangular box enclosing the molecule, or the geometric center of the molecule. Next, the dielectric constant stored on all of the grid edges is set to the solvent value, and the Debye-Huckel parameters are likewise set. Next, the program loops over all atoms. Grid points of the protein excluded space are assigned grid charges as described above, their dielectric constants are set to interior values, and the Debye-Huckel parameters are set to zero. The van der Waals, accessible, or molecular surface can be used to delimit the interior. Upon the direction of the user, the program can set the boundary using one of three possible rules. The first possibility is to set it to zero. This is very fast, but is only appropriate when the boundary is very far from the atoms in the system. The second possibility, which is denoted as the system Debye rule, is to set it to the potential caused by a single charge equal to the total charge of the system and whose radius is equal to that of a sphere whose volume equals the solvent excluded volume of the system. This method is also efficient, but it can capture some of the electrostatic screening due to the solvent if there is a thick solvent boundary around the system. The final boundary setting method, which is denoted as the atomic Debye rule, sets the boundary points to the sum of the potentials due to all of the atoms scaled by the Debye-Huckel rule,⁵ The atomic Debye rule is the most accurate.

phi = k_e q_i e^(-kappa r_s) / (\e_s r (1 + \kappa R))

where r is the distance from boundary point to the atom, R is the radius of the atom plus the water radius plus the length of the ion exclusion (Stern) layer and r_s is r-R.

This last option does a more detailed calculation of the boundary, but is more expensive computationally.

The potential is calculated using Gauss-Seidel iteration with odd-even ordering of processing grid points.⁶ When the linearized form of the equation is solved, overrelaxation is used. The default spectral radius used for the overrelaxation parameter is given by

rho = 1/3 (cos(pi/n_x) +cos(pi/n_y) + cos(pi/n_z))

where n_x, n_y, and n_z are the grid dimensions in the x, y, and z directions.