Poisson-Bolztmann Equation Module
The PBEQ module allows the setting up and the numerical solution of
the Poisson-Boltzmann equation on a discretized grid for a solute molecule.
WARNING: The module has not been used in numerous situations and caution
should be the rule. All problems should be reported to
. Benoit Roux at rouxb@PLGCN.Umontreal.CA, phone (514) 343-7105
. Dmitrii Beglov at beglovd@PLGCN.Umontreal.CA
. Wonpil Im at imwo@PLGCN.Umontreal.CA
* Menu:
* Syntax:: Syntax of the PBEQ commands
* Function:: Purpose of each of the commands
* Examples:: Usage examples of the PBEQ module
Syntax
[SYNTAX PBEQ functions]
Syntax:
PBEQ enter the PBEQ module
END exit the PBEQ module
Subcommands:
SOLVe convergence-specifications grid-specifications
physical variable-spec. boundary potential-spec.
membrane-specifications spherical droplet-spec.
cylinder-specifications dielectric boundary-spec.
force-specifications charge interpolation-spec.
atoms-selection
ITERate convergence-specifications
ENPB
CAPAcitance
WRITE property [CARD] [write-range] [UNIT integer]
READ PHI [UNIT integer]
COOR coordinate-manipulation-command
SCALar scalar-manipulation-command
HELP
RESET
convergence-specifications::=[MAXI integer] [DEPS real] [DOMEga real]
[OSOR]
MAXI [2000] : number of iterations
DEPS [0.000002] : parameter of convergence
DOME [1.0] : initial mixing factor
OSOR [.FALSE.] : optimization of the over-relaxation parameter is used?
grid-specifications::= [NCEL integer] [DCEL real]
[NCLX integer] [NCLY integer] [NCLZ integer]
[XCEN real] [YCEN real] [ZCEN real]
NCEL [65] : number of grid point in 1D for a cubic
DCEL [0.1] : size of grid unit cell
NCLX [NCEL] : number of grid point in X for general parallelepiped
NCLY [NCEL] : number of grid point in Y for general parallelepiped
NCLZ [NCEL] : number of grid point in Z for general parallelepiped
XCEN [0.0] : the center of a box in X
YCEN [0.0] : the center of a box in Y
ZCEN [0.0] : the center of a box in Z
physical variable-specifications::= [EPSW real] [EPSP real] [WATR real]
[CONC real] [TEMP real]
EPSW [80.0] : bulk solvent dielectric constant
EPSP [1.0] : protein interior dielectric constant
WATR [0.0] : solvent probe radius
CONC [0.0] : salt concentration [moles/liter]
TEMP [300.0] : Temperature [K]
boundary potential-specifications::= [ZERO] [INTBP] [FOCUS]
INTBP [.FALSE.] : INTerpolation of Boundary Potentials is used?
ZERO [.FALSE.] : boundary potentials are set to ZERO ?
(metalic conductor boundary conditions)
FOCUS [.FALSE.] : previous potentials are used to set up
boundary potentials ?
membrane-specifications:: [TMEMb real] [HTMEmb real] [ZMEMb real] [EPSM real]
[EPSH real] [VMEMB real]
TMEMB [0.0] : thickness of membrane (along Z)
HTMEMB [0.0] : thickness of headgroup region
ZMEMB [0.0] : membrane position (along Z)
EPSM [1.0] : membrane dielectric constant
EPSH [EPSM] : membrane headgroup dielectric constant (optional)
VMEMB [0.0] : potential difference across membrane (entered in [volts])
QIMGB [.FALSE.] : use the nearest image atoms
for boundary potential calculations in the XY plane
spherical droplet-spec.::= [DROPlet real] [EPSD real]
[XDROplet real] [YDROplet real] [ZDROplet real]
DROPLET [0.0] : radius of spherical droplet
EPSD [1.0] : dielectric constant of spherical droplet
XDROPLET [0.0] : position of spherical droplet in X
YDROPLET [0.0] : position of spherical droplet in Y
ZDROPLET [0.0] : position of spherical droplet in Z
cylinder-specifications::= [RCYLN real] [HCYLN real] [EPSC real]
[XCYLN real] [YCYLN real] [ZCYLN real]
[CTOM] [CKAP]
RCYLN [0.0] : radius of cylinder
HCYLN [0.0] : height of cylinder
EPSC [1.0] : dielectric constant inside cylinder
XCYLN [0.0] : position of cylinder in X
YCYLN [0.0] : position of cylinder in Y
ZCYLN [0.0] : position of cylinder in Z
CTOM [.FALSE.] : the dielectric constant of the overlapped region
with membrane is set to EPSM ?
CKAP [.FALSE.] : the Debye-Huckel factor inside the cylinder
is set to KAPPA ?
dielectric boundary-specifications::= [SMOOTH] [SWIN real]
SWIN [0.5] : solute-solvent boundary Smoothing WINdow
force-specifications::= [FORCE] [STEN real] [NPBEQ integer]
STEN [0.0] : surface tension coefficient (in kcal/mol/A^2)
NPBEQ [1] : the frequency for calculating solvation forces
during minimizations and MD simulations
charge interpolation-spec.::= [BSPLine]
BSPLine [.FALSE.] : the Cardinal B-spline method is used?
atoms-selection::= a selection of a group of atoms
range::= XFIRST real XLAST real YFIRST real YLAST real ZFIRST real ZLAST real
property::= [PHI] { KCAL/MOL }
{ VOLTS }
[FKAPPA2]
[CHRG]
[EPSX]
[EPSY]
[EPSZ]
[TITLE]
PHI : Electrostatic Potential [ KCAL/MOL ]
[ VOLTS ]
FKAPPA2 : Debye screening factor
CHRG : charges on the lattice
EPSX : X sets of dielectric function
EPSY : Y sets of dielectric function
EPSZ : Z sets of dielectric function
TITLE : formatted title line
General discussion regarding the PBEQ module
1. SOLVE
Prepare grids and solve PB equation for the selected atoms and return the
electrostatic free energy in ?enpb = (1/2)*Sum Q_i PHI_i over the lattice.
The factor of 1/2 is there for the linear response free energy of charging.
The atomic contributions are returned in WMAIN (destroying the radii).
The maximum number of iterations can be specified. The convergence parameters
DEPS and DOMEga should not be modified. The numerical algorithm is pretty
much the same as DelPhi as described in: Klapper et al. Proteins 1, 47 (1986).
NOTE: The number of grid points in X, Y, and Z (NCEL,NCLX,NCLY,NCLZ) must
**** be odd. Otherwise, the number of grid points will be increased by ONE
without any WARNING message.
NOTE: At the first stage of PBEQ or after "RESET", WMAIN should be set to
the atomic radii for the calculation. After a call to SOLVE the atomic
radii are saved in a special array. The atomic contibution to the
electrostatic free energy are returned in WMAIN (destroying the radii).
To modify the value of the radii, the keyword RESET must be issued.
1) OSOR (Optimized Successive OverRelaxation)
With OSOR keyword, the relaxation parameter will be optimized.
(Reference: A. Nicholls et al; J. Comput. Chem, 12(4),435-445 (1991))
Note that one can compare the iteration time using TIMER 2.
2) Boundary Potential
By default, boundary conditions are calculated using the Debye-Huckel
approximation. However, the computational time increases prohibitively
as the number of grid points and of atoms in the system increases.
ZERO keyword sets boundary conditions at the edge of the grid to zero.
Significant computational time for boundary potentials of a large
system. However, the result will depend upon the grid parameters and
molecular properties.
INTBP keyword will use the bilinear interpolation to construct
boundary potentials in a box with DCEL and (NCLx,NCLy,NCLz) from those
in a box with 2*DCEL and (NCLx/2+1,NCLy/2+1,NCLz/2+1).
FOCUS keyword will use previously calculated potentials to set up
boundary potentials.
(Reference: M.K. Gilson et al; J. Comput. Chem. 9(4),327-335 (1987))
(see also an example below)
3) SMOOTH
SMOOTH changes the attribute of the solute-solvent boundary.
By default (NO SMOOTH), the boundary is defined by the van der Waals
surface or the molecular surface (with WATR). SMOOTH keyword change
the boundary as a region having +/- SWIN (Smoothing WINdow) from the
surface of the solute. Within the solute-solvent boundary,
the dielectric constant and the Debye screening factor will be changed
continuously from EPSP and zero to EPSW and the screening factor
at bulk solvent.
4) FORCE
This keyword invokes the calculation of the solvation free energy and
forces and must be followed by SMOOTH keyword. The solvation energy is
taken as a sum of electrostatic and nonpolar solvatin energy.
The former is calculated from the PB equation and the latter by using
the surface tension coefficient (STEN) that relates free energy with
surface area. Note that the calculated surface is approximately the
van der Waals surface. If membrane is considered, the surface of the
membrane is also approximately included. The corresponding forces are
also calculated and will be used in minimizations and MD simulations
where NPBEQ can be used to specify the frequency for calculating the
solvation forces. Note that SWIN must be equal or greater to DCEL to
get correct solvation free energy and forces.
(Reference: W. Im, D. Beglov and B. Roux
Continuum Solvation Model: computation of electrostatic
forces from numberical solutions to the PB equation,
Comput. Phys. Commun. 109,1-17 (1998))
NOTE:To print out the force of each atom, PRNLEV should be greater
than 6.
4) Charge Distribution Method
The default is the trilinear method to distribute a charge over
nearest 8 gird points. BSPLINE keyword will invoke the 3rd-order
B-splines interpolation over nearest 27 grid points.
B-splines method removes discontinuities in the reaction field forces.
2. ITERATE
Pursue the iteration on the grid. SOLVE must have been called first.
3. ENPB
Compute the electrostatic PB energy Sum Q_i PHI_i over the lattice.
Notice that the electrostatic energy is twice as much as the electrostatic
free energy (see above). The value of the electrostatic energy is passed
through the substitution parameter enpb.
4. CAPACITANCE
Compute the capacitance based on the net induced charge in the double
layer. The induced charge beyond the limits of the box are estimated based on
the analytical solution to a planar membrane.
5. WRITE
The WRITE command is used to write out the grid properties. By default,
a binary file of the property will be written for the whole grid. The keyword
CARD implies that a formatted output will be produced. In that case, the
spatial range can be specified for the output. By default, the electrostatic
potential PHI is given in [UNIT CHARGE]/[ANGS]. If specified, the PHI can be
given in [VOLTS] or in [KCAL/MOL].
6. READ
The READ command is used to read the electrostatic potential PHI in
[UNIT CHARGE]/[ANGS] written in a binary file with the whole grid.
7. RESET
Resets all assignements of the PBEQ module and free the HEAP array.
Destroys all lists and grids. By default, the grids and arrays remain assigned
when exiting and re-entering the PBEQ module. This is to allow multiple call
to PBEQ without having to free the HEAP and other arrays if they are going
to be used again. The RESET keyword must be used to re-assign new values for
the atomic radii.
8. Miscellaneous command manipulations
*note misc: (miscom.doc) are supported within the PBEQ module,
allowing opening and closing of files, streaming of files, label assignments
(e.g., LABEL), and repeated loops (e.g., GOTO), parameter substitutions
(e.g., @1,@2, etc...) control (e.g., IF 1 eq 10.0 GOTO LOOP) and CALC
(e.g., CALC energy = ?enpb).
NOTE: TIMER 2 gives the times of various components in PBEQ module;
the grid parameter preparation (subroutine MAYER),
iterative solution (subroutine PBEQ1), and,
force calculation (subroutine RFORCE and BFORCE).
9. COORMAN and SCALAR commands
*note misc: (corman.doc) and (scalar.doc) are supported within
the PBEQ module, allowing the easy manipulation of charges, radii, rotation
and translations of molecules, etc...
10. A set of "ATOMIC BORN RADII"
Atomic radii derived from solvant electrostatic charge distribution may be
used. (test/data/radius.str) These radii were tested with free energy
perturbation with explicit solvent.
(Reference: M. Nina, D. Beglov and B. Roux.
Atomic Radii for Continuum Electrostatics Calculations based on
Molecular Dynamics Free Energy Simulations.
J. Phys. Chem. 101(26),5239-5248,1997).
NOTE: A typo for residue HSD was present in the original set of radii.
Check with M. Nina for new updated file.
To get the set of appropriate radii when using SWIN,
the commands are as follows;
STREAM RADIUS.STR
SCALAR WMAIN ADD {SWIN}
SCALAR WMAIN MULT {FACTOR}
SCALAR WMAIN SET 0.0 SELE TYPE H* END
The factor has a linear relationship with SWIN.
-----------------------------------------------------------------------------
SWIN 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
FACTOR 0.979 0.965 0.952 0.939 0.927 0.914 0.901 0.888 0.875 0.861
-----------------------------------------------------------------------------
** FACTOR = -0.1296 x SWIN + 0.9914 (a least-square fit)
Examples
This examples are meant to be a partial guide in setting up
an input file for PBEQ. There are two test files, pbeqtest1.inp,
pbeqtest2.inp, and pbeqtest3.inp.
Example (1)
-----------
This example shows how to perform two PB calculations, one for a surrounding
dielectric of 80 (water) and one for a surrounding of 1.0 (vaccum). The
difference between the two energies then corresponds to the electrostatic
contribution to the solvation free energy. The salt concentration was zero
in this calculation.
PBEQ
scalar wmain = radius
SOLVE epsw 80.0 conc 0.0 ncel 30 dcel 0.4
set ener80 = ?ENPB
SOLVE epsw 1.0
set ener1 = ?ENPB
CALC total = @ener80 - @ener1
RESET
END
Example(2)
----------
This example shows how to use a set of atomic Born radii with a smoothing
window.
set sw 0.4
set factor 0.939
PBEQ
stream radius.str
scalar wmain add @sw
scalar wmain mult @factor
scalar wmain set 0.0 sele type H* end
scalar wmain show
SOLVE epsw 80.0 ncel 100 dcel 0.3 -
smooth swin @sw force sten 0.03 npbeq 1
RESET !! If you consider a minimization or dynamics with PB forces,
!! don't use RESET here.
END
Example(3)
----------
This example shows how to set up a membrane potential and how to get
the electrostatic contribution to the solvation free energy in the membrane
environment. Note that a non-zero concentration is required for a sensible
system with a membrane potential.
PBEQ
scalar wmain = radius
SOLVE epsw 80.0 ncel 150 dcel 0.5 conc 0.150 -
Tmemb 25.0 Zmemb 0.0 epsm 2.0 vmemb 0.100
set ener80 = ?ENPB
SOLVE epsw 1.0 conc 0.000 -
Tmemb 25.0 Zmemb 0.0 epsm 1.0 vmemb 0.000
set ener1 = ?ENPB
CALC total = @ener80 - @ener1
RESET
END
Example(4)
----------
This example shows how to set up boundary potentials using FOCUS keyword,
how to read the saved potential, and how to calculate the electrostatic
contribution to the solvation free energy using FOCUS.
PBEQ
scalar wmain = radius
SOLVE epsw 1.0 ncel 60 dcel 0.4
open write file unit 40 name phi.dat
write phi unit 40
SOLVE epsw 1.0 dcel 0.2 focus ! boundary potentials from DCEL 0.4 potentials
! NOTE: YOU CAN CHANGE NCEL IN THE FOCUSSED SYSTEM AS FOLLOWS;
! SOLVE epsw 1.0 ncel 80 dcel 0.2 focus
SOLVE epsw 1.0 dcel 0.1 focus ! boundary potentials from DCEL 0.2 potentials
open read file unit 41 name phi.dat
read phi unit 41
SOLVE epsw 1.0 dcel 0.1 focus ! boundary potentials from DCEL 0.4 potentials
RESET
END
PBEQ
scalar wmain = radius
SOLVE epsw 80.0 ncel 60 dcel 0.4
set ener81 = ?ENPB
SOLVE epsw 80.0 dcel 0.2 focus
set ener82 = ?ENPB
SOLVE epsw 80.0 dcel 0.1 focus
set ener83 = ?ENPB
SOLVE epsw 80.0 dcel 0.05 focus
set ener84 = ?ENPB
SOLVE epsw 1.0 dcel 0.4
set ener11 = ?ENPB
SOLVE epsw 1.0 dcel 0.2 focus
set ener12 = ?ENPB
SOLVE epsw 1.0 dcel 0.1 focus
set ener13 = ?ENPB
SOLVE epsw 1.0 dcel 0.05 focus
set ener14 = ?ENPB
calc total = @ener81 - @ener11
calc total = @ener82 - @ener12
calc total = @ener83 - @ener13
calc total = @ener84 - @ener14
SOLVE epsw 80.0 ncel 120 dcel 0.2
set ener80 = ?ENPB
SOLVE epsw 1.0
set ener1 = ?ENPB
calc total = @ener80 - @ener1
RESET
END
CHARMM Documentation / Rick_Venable@nih.gov