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
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(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.
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.(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
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.
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,
kappa_bar.
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.
Syntax
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]
[OLDGRID]
[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.
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. The OLDGRID
option is similar, except it maintains only the grid, but not the values
within
These options
are 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 OLDGRID 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 or OLDGRID option will be turned off.
Syntax
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:
Syntax
PBE ENERGY [ATOM]
This command calculates the electrostatic energy of the system using the formula, E = 1/2 sum over i,j,k rho(i,j,k)*phi(i,j,k)
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.
Syntax
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.
Syntax
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.
Syntax
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:
Syntax
PBE SAVE name
The entire pbe_info data structure is saved under the given name. Any other pbe_info data structure saved under the given name will be destroyed. The PBE RECALL can be used to recover it.
Syntax
PBE RECALL name
The entire pbe_info data structure is replaced with the one stored under name. Memory associated with the previous data structure is freed.
Syntax
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
The PBE DESTROY command destroys the current 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:
Syntax
PBE DUMP
The PBE DUMP command displays detailed information about the Poisson-Boltzmann data structures. It is useful primarily for debugging and regression testing.
Syntax
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 first, previous, next, last section, table of contents.