Go to the previous, next section.
The Poisson-Boltzmann equation (PBE) provides a "fast" and approximate method for calculating the effects of solvent on electrostatic interactions in the modeled system. The current solution of the PBE in CONGEN has one unique feature, it can spread the charge of each atom over the van der Waals volume which results in superior convergence properties and makes the calculation less dependent on the position of the system with respect to the grid(3). In addition, CONGEN can display the partition of the Poisson-Boltzmann energy across all the atoms in the system, see the EPBE property described in section Static Properties of Atoms.
The implementation provides two capabilities. First, one can perform individual electrostatic calculations on any set of atoms within the system. Second, one can replace the Coulomb electrostatic energy with the Poisson-Boltzmann energy in the context of a conformational search, see section Poisson-Boltzmann Options. We hope to increase the applicability of this methodology in the future.
One of the most successful implicit models for the treatment of electrostatic effects is the Poisson-Boltzmann equation(4) which is given below where
The equation can be solved over a volume enclosing a molecule of interest using finite difference methods(5) (6) (7) 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,(8) The atomic Debye rule is the most accurate.
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.(9) When the linearized form of the equation is solved, overrelaxation is used. The default spectral radius used for the overrelaxation parameter is given by
In order to understand the function of the commands which use the Poisson-Boltzmann equation, it is important to know the underlying data structures.
The primary data structure is the pbe_info data structure, defined
in the file, `$CGS/pbe.h'. This data
structures contains the grids for the potential, charge density, dielectric
constant, and Debye-Huckel parameter,
There are three grids stored for the dielectric constant because values
for the dielectric constant are stored for the lines connecting grid
points. Since there are three perpendicular lines connecting grid
points to their neighbors, there are three dielectric constant grids,
for the x, y, z directions. The epsx grid stores the dielectric
constant for the midpoint of the line pointing in the x direction
between grid points. The value of epsx[i,j,k] corresponds to the
dielectric constants between grid points [i,j,k] and
[i+1,j,k]. The value of epsy[i,j,k] corresponds to the
dielectric constants between grid points [i,j,k] and
[i,j+1,k], and the value of epsz is
analagous for the the third indices, respectively. One two-dimensional cross
section of each grid is unused.
The pbe_info data structure also contains the spatial origin
of the grid, dimensions, and grid size so that the grid
can be laid out over the space of the molecule of interest. It also contains
a list of atoms included in the grid as well as their positions, radii, and
charges. Finally, it contains a number of parameters which control
the calculation.
The PBE SETUP command is used to setup this data structure. After the first PBE SETUP, the REUSE option controls whether the data structure is initialized or updated, see section PBE SETUP Command.
It is important to remember that this data structure is not automatically updated when the coordinates change (except in the context of a conformational search).
It is possible to generate multiple pbe_info data structures,
and to do operations on them. This is useful for comparing results
generated by different methods or parameters.
All the Poisson-Boltzmann commands begin with the keyword, PBE. The second word on each command line specifies one of the command below. In brief, the commands perform the following operations.
pbe_info data structure under a variable name.
pbe_info data structure into the main data structure.
pbe_info data structure.
PBE SETUP repeat(pbe-setup-options) atom-selection
[SOLVent real]
[interior-options]
[charge-options]
[NWARN int]
[boundary-options]
[IONStr real]
[GRID real]
[SUBDivisions int]
[SPILl]
[TEMPerature real]
pbe-setup-options ::= [WATEr real]
[STERn real]
[XDIM int]
[YDIM int]
[ZDIM int]
[MARGin int]
[REUSE]
[CENTER]
[MOLSURF]
[MINRadius real]
[charge-edit-options]
[EPSAve averaging-option]
[SMOOTH repeat(smooth-option) END]
interior-options ::= repeat([INTErior real exposure-opts atom-selection END])
exposure-opts ::= [ABSQ real] [EXPOsure real] [RESAve]
charge-options ::= [UNIForm ]
[TRILinear]
[ZERO ]
boundary-options ::= BOUNdary [ADEBye ]
[SDEBye ]
[PREVious]
charge-edit-options ::= repeat([CHARge edit-type real atom-selection END])
[SET]
edit-type ::= [SCALe]
[SHIFt]
averaging-option ::= [HARMonic ]
[ARIThmetic]
[TYPE type-option]
[WEIGht weight-option]
smooth-options ::= [OFFGrid offgrid-option]
[RADIus real]
[POINts int]
[ALPHa real]
[NONE ]
type-option ::= [VOLUME]
[GRID ]
weight-option ::= [CONStant]
[GAUSsian]
offgrid-option ::= [SOLVent]
[EDGE ]
The PBE SETUP command creates the pbe_info data structure,
see section PBE Data Structures, which stores the grids upon which the
electrostatic potential is computed. The command has a large number of
options and an atom-selection, see section Atom Selection. The atom
selection allows the user to select any portion of the system for the
calculation. For example, if one is interested in the calculation of the
binding energy of a complex, the atom selection can be used to select
each component and then all the atoms for separate electrostatic
calculations. The options are described in the following table.
The charging options, UNIFORM and TRILINEAR, along with the options SUBDIVISIONS and SPILL, control how the charge grid is set. With the TRILINEAR option, the charge of each atom is divided among the eight nearest points to the atom center. With the UNIFORM option, the charge of each atom is divided evenly among all the points within the van der Waals radii of the atom, except in the case where the number of such points is less than eight. In that case, trilinear interpolation is used. The option, NWARN, controls how many warning of this conversion are displayed.
The uniform charging option can be further improved using the SUBDIVISIONS and SPILL options. One problem with the use of grid based methods to represent molecular structure is the quantization of space. The same problem has arisen in the development of computer graphics using raster devices, and using the techniques of anti-aliasing, CONGEN can smooth the charge distribution of atoms. When the SUBDIVISIONS option is used, the space around each atom is subdivided into virtual cubes, and the charge distribution is calculated over these virtual cubes, and then added onto the real cubes in the grid. If the SPILL option is turned on (recommended!), then charge can spill onto cubes outside the van der Waals radius of each atom. If SPILL is not specified, then the charge will not extend beyond the van der Waals radius.
The boundary-options are used to specify how the potential is set for the boundary of the grid. The ZERO option specifies that the boundary potential should be zero. The ADEBYE option specifies that the atomic Debye method should be used, see section Introduction to Poisson-Boltzmann Electrostatics. This option is very slow. The SDEBYE options specifies that the system Debye method should be used. The keyword, PREVIOUS, is not implemented.
The IONSTR parameter is used to set the ionic strength. The units are in molarity. The default value is 0.0.
The GRID option sets the physical spacing between grid points. This is a critical parameter in the calculation of electrostatics using the Poisson-Boltzmann equation. The smaller the grid, the better the accuracy, but computation time scales as the grid size to the 4th power. The default value is 1.25 Angstroms.
pbe_info data structure. For the
keyword options described here, the default values are taken from the current
values rather than the values described in this documentation.
This option is important for doing binding energy calculations. In order to cancel errors due to grid positioning, it is essential to first calculate the energy of the complex, and then, using the REUSE option, calculate the energies of the components by selecting the atoms of the components, but leaving the atoms and grids in the same spatial positions.
When this option is used, you should not specify any grid dimensions or the MARGIN option. If you do, and they are different than the actual values, then the REUSE option will be turned off.
The SMOOTH command "smooths" the dielectric grids. After the dielectric grids have been assigned internal and external values "smoothing" averages each point with the other points in its neighborhood to produce a dielectric grid with less abrupt changes in value. One hope of the process is to reduce the position dependence of the PBE electrostatic calculations. Its success in this endeavor is still a topic of research.
Several types of smoothing are supported. The smoothing TYPE can either be NONE, VOLUME, or GRID. Of course, no smoothing is done when NONE (the default) is specified. For VOLUME based smoothing the dielectric is averaged over a fixed volume in real space defined by the RADIUS keyword. For GRID based smoothing the dielectric is averaged over a set number of points specified by the POINTS keyword (9 and 15 point averaging are currently supported). In both case all three dielectric grids are averaged togeher. In 9-point averaging each point is averaged with the eight grid points from the other two grids which surround this point in space. In 15-point averaging the six nearest neighbors from the same grid are also included. (Note that volume smoothing with a radius just smaller (larger) than the grid spacing is equivalent to 9-point (15-point) averaging.)
Which ever TYPE is selected, The type of averaging is determined by the EPSAVE keyword above. The WEIGHTING can either be CONSTANT (all points in average counted equally), or GAUSSIAN weighted (points are weighted by specified by the ALPHA keyword).
For points at the edge of the grid the averaging extends to points beyond the grid. Using the OFFGRID keyword, these offgrid points can either be assumed to be SOLVENT or equal to the value at the nearest EDGE of the grid. (Note that the SOLVENT assumption will run faster.)
PBE POTENTIAL repeat(pbe-potential-options)
[NODIelectric]
[TOLerance real]
[RTOLerance real]
pbe-potential-options ::= [MAXITS int]
[HISTsize int]
[RADIUS real]
[problem-option]
problem-option ::= {LINEar }
{NONLinear}
The PBE POTENTIAL command invokes the Gauss-Siedel PBE solver, and the above options affects how the solution is performed. If the PARALLEL LOOPS command is used, see section PARALLEL Command, the potential will be calculated in parallel across available CPU's on a Silicon Graphic multiprocessor workstation. Parallel efficiency is high. N.B., for large memory jobs running on SGI workstations, it is critically important to set the STACK resource limit down from the default value. See section RLIMIT Command -- Set Resource Limits, for more information.
After the potential is calculated, the PBE energy per atom is calculated and stored in an array for this purpose.
The options have the following meaning:
If this option is specified, the potential is calculated assuming the dielectric is equal to the solvent value, keyword SOLVENT in section PBE SETUP Command. The Poisson-Boltzmann equation is not used, but rather the potential contribution from each atom is summed onto each grid point. The boundary is also modified.
If the maximum change in any potential value during a cycle over all the grid points is less than the TOLERANCE option, then the solver terminates. Values smaller than are generally too small to be achieved. The default setting is 0.001.
If the RMS of the changes in potential value during a cycle over all the grid points (maximum residuals) is less than this option, the solver terminates. This option is less restrictive than the TOLERANCE option above, and therefore, terminates sooner. The maximum change is typically an order of magnitude bigger than the RMS change. The default setting is 0.0, in which case, the test is not used.
The specifies the maximum number of iterations allowed for the solver. Because of the odd-even checkboard cycling, two iterations are necessary to cover all the grid points. The default is 100.
When the PBE solver is working toward a small tolerance, it is possible for a potentially infinite loop to be entered. In order to detect this, CONGEN maintains an array of maximum residuals for twice the number of cycles given by HISTSIZE. At the end of each iteration over all the atoms, the minimum value for the maximum residual for the last HISTSIZE cycles is compared against the minimum for the previous HISTSIZE cycles. If it has not improved, then the PBE solver terminates. Effectively, this mechanism checks for any improvement in the solution of HISTSIZE cycles. The default value is 30.
Overrides the default setting, see section Introduction to Poisson-Boltzmann Electrostatics, for the spectral radius. If this value is set to zero, then overrelaxation is not used. If the value is negative, then the absolute value is used as the overrelaxation parameter.
These two options control whether the linearized form of the Poisson-Boltzmann equation is solved. Generally, only the linearized form should be used. The non-linear form has a modest radius of convergence, and large charges will result in divergence of the solver. When the non-linear form is used, execution time is much slower, and overrelaxation is turned off. The default is LINEAR.
PBE ENERGY [ATOM]
This command calculates the electrostatic energy of the system using the formula,
If the option, ATOM, is specified, the PBE energy per atom will be printed. It is usually more convenient to get these energies using the ANALYSIS commands, see section Static Properties of Atoms.
PBE READ UNIT unit
This command reads a Poisson-Boltzmann data structure which has been written to a disk by the PBE WRITE command, see section PBE WRITE Command. The only option, which is required, is the unit number of the file containing the data structure.
PBE WRITE [UNIT unit] [TITLE string del] { ALL }
{ PHI }
{ RHO }
{ EPSX }
{ EPSY }
{ EPSZ }
{ KAPPA }
[X range] [Y range] [Z range] [LINEsz int] [CARD]
[FILE]
range ::= [[int]:[int]]
This command is used to write the entire Poisson-Boltzmann data structure to a binary file, or it can be used to write one of the grids to an output file in either binary or text format. The binary format for the grids is compatible with the PLT2 program, so that contour maps of sections of potentials or other grids can be made.
N.B. the behavior of this command changes depending on whether ALL is specified, so please take note of these confusing changes as described below.
The options have the following meaning:
At least one of RHO, PHI, EPSX, EPSY, EPSZ, KAPPA must be specified.
PBE TEST { CV cv-options } repeat(pbe-test-options)
{ DH dh-options }
cv-options ::= [CHARge real ]
[RADIus real ]
[CAVIty-eps real]
[QPOS real ]
[CHARge real ]
dh-options ::= [RADIus real ]
[CAVIty-eps real ]
[IONStr real ]
[TEMPERATURE real ]
[XDIM int ]
[YDIM int ]
pbe-test-options ::= [ZDIM int ]
[GRID real ]
[SOLVent real]
The PBE TEST command is used to generate test potentials for comparing the finite difference solution against two particular problems for which analytic solutions are available. The two problems are a charge in a spherical cavity of one dielectric with a different surrounding dielectric (the CV option), and a charged centered in a sphere surrounded by a Debye-Huckel fluid (the DH option). These calculations are done using analytic series representations, but they are not optimized for numerical accuracy, so beware of artefactual peaks. CONGEN will report the number of failures of the series convergence. Also, a setting of 3 for the PBE DEBUG variable will display the terms in each series calculation, see section Set Debugging Variables -- DEBUG. A debug setting greater than 10 will cause CONGEN to substitute the number of iterations in place of the potential.
The analytic potentials used in this code were derived by Malcolm Davis.
The command options have following meaning:
PBE SAVE name
PBE RECALL name
PBE DIFF name
The PBE DIFF command calculates the difference in the potential grid and in the PBE energies for each atom between the current pbe_info data structure, and the one stored under name. The results are left in the current data structure.
PBE DESTROY
pbe_info data structure.
The command has no options.
PBE STATISTICS data-name [MIN real] [MAX real]
[LINES int] [COLUMNS int] [IGNORE real]
{ PHI }
{ RHO }
data-name ::= { KAPPA }
{ EPSX }
{ EPSY }
{ EPSZ }
The PBE STATISTICS command prints a histogram of the values in any of the six arrays used in the pbe_info data structure. The options are as follows:
PBE DUMP
The PBE DUMP command displays detailed information about the Poisson-Boltzmann data structures. It is useful primarily for debugging and regression testing.
PBE MASK data-name relop real [ABS] [USING data-name] [SETTO real]
[RECALL name]
{ PHI }
{ RHO }
data-name ::= { KAPPA }
{ EPSX }
{ EPSY }
{ EPSZ }
{ GT }
{ GE }
relop ::= { LT }
{ LE }
{ EQ }
{ NE }
The PBE MASK command allows the user to identify or mask values in any of the grids. This is useful for identifying interesting regions prior to plotting, or for gathering statistics about regions of the system. For example, it is possible to mask all regions of the potential that lie within the van der Waals surface, and then calculate statistics on the outer regions.
The options have the following meanings:
The following commands illustrates setting the potential to zero for all points within and including the Stern layer. Note that this command depends on the ionic strength being positive.
PBE MASK PHI EQ 0.0 USING KAPPA
Two examples of using the Poisson-Boltzmann equation are presented. One example shows the calculation of electrostatic binding energy, and the other shows the basics of using the PBE in a conformational search.
In this example, we illustrate how to calculate the electrostatic binding energy of a complex between two molecules. Assume that the molecules have segment identifiers of A and B.
! Calculate energy of the complex.
pbe setup solvent 78.0 interior 2 all end -
uniform boundary zero ionstr 0.15 -
grid 0.4 temp 300.0 water 0.0 stern 2.0 -
xdim 190 ydim 170 zdim 135 -
all
pbe potential maxiter 2000 rtol 1.0e-4
pbe energy
! Calculate energy of segment A
pbe setup solvent 78.0 interior 2 all end reuse -
uniform boundary zero ionstr 0.15 -
grid 0.4 temp 300.0 water 0.0 stern 2.0 -
xdim 190 ydim 170 zdim 135 -
clear atom A * *
pbe potential maxiter 2000 rtol 1.0e-4
pbe energy
!Calculate energy of segment B
pbe setup solvent 78.0 interior 2 all end reuse -
uniform boundary zero ionstr 0.15 -
grid 0.4 temp 300.0 water 0.0 stern 2.0 -
xdim 190 ydim 170 zdim 135 -
clear atom B * *
pbe potential maxiter 2000 rtol 1.0e-4
pbe energy
The electrostatic binding energy would be the difference between the first energy and the sum of the last two.
In this example, it should be noted that the calculations of each component are done using the same grid as the complex. This results in the clean subtraction of the self-energies.
This example is taken from `$CGT/cgpbe2.inp'. It illustrates using the substitution of the Poisson-Boltzmann electrostatic energy for the Coulomb energy for evaluating all the conformers generated by search, where the sidechains are optimized using the CHARMM potential. This is necessary to keep the execution time reasonable, but this is not a good example for a real application.
Conformational search over NEW
Test of PBE evaluation in parallel
*
OPEN NAME CGDATA:RTOPH8.MOD UNIT 01 READ UNFORM
OPEN NAME CGDATA:PARAM5.MOD UNIT 03 READ UNFORM
OPEN NAME CGTD:newpsf.mod UNIT 12 READ UNFORM
OPEN NAME CGTD:NEWV.MOD UNIT 14 READ UNFORM
READ RTF UNIT 1
READ PARAMETER UNIT 3
READ PSF FILE UNIT 12
READ COOR FILE UNIT 14
parallel loops ncpu 4
!
! Setup the conditions for using Poisson-Boltzmann equation.
! The grid size is too large for a real calculation, though.
!
PBE SETUP IONSTR 0.15 BOUNDARY SDEBYE CENTER SOLVENT 78 GRID 1.5 -
WATER 0.0 STERN 2.0 XDIM 50 YDIM 50 ZDIM 50 -
TRILINEAR INTERIOR 2.0 ALL END
PBE POTENTIAL RTOLER 1.0E-3 MAXITER 1000
PBE ENERGY
!
! Generate conformations for the loop H 31-35 using the CHARMM
! potential energy.
!
OPEN UNIT 60 NAME CGPBE2.CG1 UNFORM WRITE
OPEN UNIT 70 NAME CGDATA:TOPCGEN2.INP FORM READ
OPEN UNIT 51 NAME CGDATA:EMAPGLY30.OMP FORM READ
OPEN UNIT 52 NAME CGDATA:EMAPALA30.OMP FORM READ
OPEN UNIT 53 NAME CGDATA:EMAPPRO30.OMP FORM READ
OPEN UNIT 55 NAME CGDATA:PRO.CNS FORM READ
CGEN -
SEARCH DEPTH END -
NBCG CUTNB 5.0 CTONNB 98.0 CTOFNB 99.0 END -
HBCG CUTHB 4.5 CTONHB 98.0 CTOFHB 99.0 -
CUTHBA 90.0 CTONHA 90.0 CTOFHA 90.0 END -
BACK CISTRANS STARTRES H 31 LASTRES H 32 MAXEVDW 100 CLSA H 35 CA $ -
CHAIN CISTRANS STARTRES H 33 MAXEVDW 100 $ -
SIDE STARTRES H 31 LASTRES H 35 MAXEVDW 100 SIDEOPT ITER SGRID MIN EVAL E $ -
WRITE CUNIT 60 $ -
ERINGPRO 50 GLYMAP 51 ALAMAP 52 PROMAP 53 PROCONS 55 -
GLYEMAX 2 ALAEMAX 2 PROEMAX 2 STUNIT 70
CGPBE2.CG1
CGPBE2 test case
Conformations of H1 loop in NEW.
*
!
! Now evaluate the conformations using the Poisson-Boltzmann energy
! substituting for the Coulomb energy.
!
CLOSE UNIT 60
OPEN UNIT 59 NAME CGPBE2.CG1 UNFORM READ
OPEN UNIT 60 NAME CGPBE2.CG UNFORM WRITE
OPEN UNIT 70 NAME CGDATA:TOPCGEN2.INP FORM READ
OPEN UNIT 51 NAME CGDATA:EMAPGLY30.OMP FORM READ
OPEN UNIT 52 NAME CGDATA:EMAPALA30.OMP FORM READ
OPEN UNIT 53 NAME CGDATA:EMAPPRO30.OMP FORM READ
OPEN UNIT 55 NAME CGDATA:PRO.CNS FORM READ
CGEN -
PBE NEWB END -
SEARCH DEPTH END -
NBCG CUTNB 5.0 CTONNB 98.0 CTOFNB 99.0 END -
HBCG CUTHB 4.5 CTONHB 98.0 CTOFHB 99.0 -
CUTHBA 90.0 CTONHA 90.0 CTOFHA 90.0 END -
RBEST UNIT 59 NBEST 9999 $ -
EVL MINI ENERGY END $ -
WRITE CUNIT 60 $ -
ERINGPRO 50 GLYMAP 51 ALAMAP 52 PROMAP 53 PROCONS 55 -
GLYEMAX 2 ALAEMAX 2 PROEMAX 2 STUNIT 70
CGPBE2.CG
CGPBE2 test case
Conformations of H1 loop in NEW.
*
!
! Extract the best conformation
!
OPEN UNIT 60 NAME CGPBE2.CG READ UNFORM
XCONF 60 BEST 1
Go to the previous, next section.