Cartesian Grid-Based Poisson-Boltzmann Equation Solver

Continuum Dynamics, Inc. (CDI) has recently developed a Poisson-Boltzmann equation (PBE) solver with full support from the NIH SBIR program (NIH grant numbers 5 R44 GM57764-03 and 5R44GM073391-03) that is based on the adaptive Cartesian grid concept (see Figure 1 below). The Cartesian grid employs an octree data structure to provide a powerful variable mesh spacing capability while simultaneously retaining a simple cell shape and a hierarchical ordering of grid elements that is intrinsically suited to multigrid. The mesh structure is applied to the implicit solvent-based PBE and its generalization to finite non-uniform ion size descriptions for modeling biomolecular electrostatics. The resulting Cartesian grid-based PBE solver, CPB, solves for the electrostatic potential in and around the (bio)molecule and post-processes the solution to evaluate various energy contributions, energy derived properties such as polar solvation and binding free energies, salt sensitivities, pKas, atomic forces, bound ion numbers, and detailed visualization of surface and volume maps of the potential and derivative quantities such as ion concentrations, ionic and dielectric pressures, and induced surface charge. CPB calculates binding energies, with or without conformational change upon binding. Such metrics are useful for drug design applications, understanding the role of nonspecific electrostatics underlying folding, and biomolecular binding (recognition) events under varying environmental conditions (e.g., salt concentration and pH).

As described in the literature (1-3, 7), CPB embodies a number of numerical and computational features that facilitate set up and operation by the user, eliminate singular behavior, promote accurate solutions especially on the surface, ensure reliable convergence behavior and allow rapid evaluation using readily accessible computer resources. Among these features are: a decomposition of the interior potential into singular (Coulombic) and non-singular (reaction field) components and appropriate analytical and numerical treatments of each; an optionally invoked least squares reconstruction of the solution near the surface that respects the local continuity conditions; built in support and generation capability for multiple surface definitions; an extensively tested and optimized multigrid solver capable of handling highly nonlinear problems (e.g., those involving nucleic acid systems)(; and a general outer boundary treatment applicable to nonlinear PBE systems. More details of CPB's unique capabilities are listed in the literature cited below.

  • The adaptive Cartesian grid (or octree) structure to solve the governing equations. As shown in Figure 1, the mesh provides high mesh resolution at the molecular surface since this is where the solution varies most rapidly. Away from the surface (both inside and outside the molecule) the mesh size increases to reduce the number of required mesh points and associated storage cost and computation time.
  • A singularity-free decomposition of the electrostatic solution amenable to grid-based representation. The full potential field contains singularities at the fixed point charges, which introduces a variety of numerical difficulties including increased error, increased mesh resolution requirements to resolve the rapidly varying field and the need for two calculations to compute reaction field energies. CPB decomposes the solution into a Coulombic part that contains the singular contributions and can be calculated analytically, and the reaction field potential which is free of singularities and thus accurately represented on the mesh.
  • Accommodation of multiple ion species and salt mixtures using the linear, nonlinear and size-modified PBE models (3). In the size-modified treatment, the radii and valencies of the ion species can be specified individually (the non-uniform ion size model). A finite size can also be attributed to the solvent. The nonlinear PBE omits ion size effects and the linear PBE further approximates the solvent model by assuming small potentials. The computational costs of the various PBE models are comparable, with the nonlinear and size modified PBE methods generally being 20-30% higher due to the evaluation of exponential functions when computing ion densities.
  • Provision of several surface descriptions used to define the dielectric boundary, including van der Waals, solvent accessible, solvent excluded (Connolly), and a variety of Gaussian surface definitions. Each of these surfaces is generated internally using the analytical definitions. Thus, there is no need to interface the software with external surface generation software. Also, relying on analytical surface definitions ensures accurate problem representation as the mesh is refined.
  • Access to a least squares fitting algorithm (2) for enhanced accuracy and representation of surface potentials and gradients (induced charge). Like other lattice mesh-based PBE codes, the molecular surface cuts through or intersects the grid cells. This makes it difficult to obtain accurate estimates of the potential at the surface. The normal gradients at the surface are even more difficult to obtain since they change discontinuously across the surface. CPB addresses these challenges by formally incorporating the matching conditions between the solution and its gradients up to second order at the surface, and then uses least squares principles to match this representation to surrounding nodes. In practice, this capability increases computation time (for performing the least squares reconstructions), but greatly improves estimates of surface potential, induced charge, pressure and surface-induced atomic forces.
  • Application of outer boundary corrections (7) to account for the electrostatic contributions from outside the computational domain. The computational domain is necessarily of finite extent so that the question of imposing appropriate boundary conditions on the domain boundaries must be addressed. Common options include simply setting the boundary potential to zero or adopting a Debye-Hckel approximation, though neither option is formally applicable to nonlinear problems. Instead, CPB employs charge conservation principles to impose the boundary conditions as described in (7). These principles are valid for both the linear/nonlinear PBE and size-modified PBE models. They are also used to account for the contributions of the region outside the computational domain to various energies and salt sensitivities.
  • Adaptation options to refine the grid near selected sites. Means are provided to enforce higher mesh refinement about selected locations such as candidate binding sites.
  • Compatibility with several structural definition input formats including: APBS (*.pqr files), UHBD (*.qcd files), Delphi's modified PDB Format (*.pdb files) and a simple generic format consisting of atomic positions, partial charges and radii.
  • Multiple output formats ranging from simple text output of energies and atomic forces, surface maps and full 3D solution sets compatible with the TecPlot plotting package. Output in DX format is also available for input to the VMD, Pymol, or Chimera visualization codes.
  • User-friendly set up and operation on readily accessible platforms (LINUX, PCs, OSX devices).

  • Barnase-Barstar complex (PDBid: 1b27).
  • 30S ribosome subunit (PDBid: 1FJG).
  • 50S ribosome subunit (PDBid: 3CC4).
  • Cowpea chlorotic mottle virus (PDBid: 1cwp).


  1. Boschitsch, Alexander H. and Marcia O. Fenley, A Fast and Robust Poisson-Boltzmann Solver Based on Adaptive Cartesian Grids. Journal of Chemical Theory and Computation, 2011. 7(5): p. 1524-1540. (cover PDF)
  2. Boschitsch, Alexander H. and Marcia O. Fenley, The Adaptive Grid-Based Poisson-Boltzmann Solver: Energy and Surface Electrostatic Properties. Computational Electrostatics for Biological Applications, W. Rocchia & M. Spagnuolo (Eds.), Springer, to appear in 2014.
  3. Boschitsch, Alexander H. and Pavel Danilov, Formulation of a New and Simple Non-Uniform Size-Modified Poisson-Boltzmann Description. Journal of Computational Chemistry, 2012. 33(11). (cover PDF)
  4. Fenley, Marcia O., et al., Revisiting the Association of Cationic Groove-Binding Drugs to DNA using a Poisson-Boltzmann Approach. Biophysical Journal, 2010. 99(3): p. 1-8. (cover PDF)
  5. Bredenberg, Johan H., Alexander H. Boschitsch, and Marcia O. Fenley, The Role of Anionic Protein Residues on the Salt Dependence of the Binding of Aminoacyl-tRNA Synthetases to tRNA: A Poisson-Boltzmann Analysis. Communications in Computational Physics, 2008. 3(5): p. 1051-1070. (cover PDF)
  6. Harris, Robert C., et al., Understanding the physical basis of the salt dependence of the electrostatic binding free energy of mutated charged ligand-nucleic acid complexes. Biophysical Chemistry, 2011. 156: p. 79-87. (cover PDF)
  7. Boschitsch, Alexander H. and Marcia O. Fenley, A New Outer Boundary Formulation and Energy Corrections for the Nonlinear Poisson-Boltzmann Equation. Journal of Computational Chemistry, 2007. 28(5): p. 909-921. (cover PDF)
  8. Harris, Robert C., Alexander H. Boschitsch and Marcia O. Fenley, Influence of Grid Spacing in Poisson-Boltzmann Equation Binding Energy Estimation, Journal of Chemical Theory and Computation, 2013. 9(8): p. 3677-3685. (cover PDF)

Barnase-Barstar Complex (PDBid: 1b27)

Reaction field (RF) potential, in kcal/mol/e, and ACG grid structure inside a slice of the barnase-barstar complex (PDBid: 1B27) immersed in a 0.1 M NaCl solution. Note the large concentration of RF contours near the molecular surface and corresponding local concentration of mesh points in the ACG.

Nonlinear PBE potential (in kcal/mol) and (c+)=0.5M iso-surface (yellow).

30S Ribosome Subunit (PDBid: 1FJG, 88.6K atoms)

Electrostatic potential on molecular surface (solvent-excluding) computed using the nonlinear PBE for a 0.1M 1:1 salt environment. The solid arrow indicates the paromomycin binding site. The complex is resolved using a 0.3A grid spacing at the surface resulting in a total of 24 million grid points. With least squares reconstruction of the surface properties, 350 iterations and 3 hours CPU on Linux serial processor workstation are required to generate the mesh, converge the potential solution and post-process the output.

50S Ribosomal Subunit (PDB id: 3cc4, 150K atoms)

Electrostatic potential obtained with nonlinear PBE at 0.1 M NaCl and mapped on the solvent excluded surface of the 50S ribosomal subunit (PDBid: 3CC4). The ACG mesh contains 52 million nodes with a grid spacing about the surface of 0.3. The solution is converged in 170 iterations. The yellow box indicates the approximate region expanded in detail below.

Close up of potential surface map produced by CPB for the 50S ribosomal subunit.

Cowpea chlorotic mottle virus (PDB id: 1cwp, 227K atoms)
Electrostatic potential surface potential (in kcal/mol/e) obtained with the nonlinear PBE at 0.1 M NaCl and mapped on the solvent excluded surface of the cowpea chlorotic mottle virus (PDBid: 1CWP) . The ACG mesh contains 123 million nodes with a grid spacing about the surface of 0.3. The solution is converged in 170 iterations and requires 18.5 hrs on a 24 core (2 AMD Opteron 6234 12-core) workstation with 128 GB memory. The yellow box indicates the approximate region expanded in detail below.

Close up of surface potential (in kcal/mol/e) map produced by CPB.