Generalized Born Solvation Energy and Forces Module
The GBORN modeule permits the calculation of the Generalized Born
solvation energy and forces of this energy following a formulation
similar to that of Still and co-workers (ref) and as described in the
manuscript from B. Dominy and C.L. Brooks, III (J. Comp. Chem.,
manuscript in preparation). This module implements the following
equation for the polarization energy, Gpol:
q q
N N i j
G = -C (1-1/eps){1/2 sum sum ------------------------------------ }
pol el i=1 j=1 [r^2 + alpha *alpha exp(-D )]^(0.5)
ij i j ij
The gradient of the function is also computed so forces due to solvent
polarization can be utilized in energy minimization and dynamics. In
its current implementation, the calculation of the alpha(i) variables
and the sums over particles indicated in the sums above are done
without cutoffs, therefore for large systems these can be costly
calculations (though still less so than for explicit solvent).
Questions and comments regarding implementation of these equations or
there parameterization for the CHARMM forcefields (param19/toph19,
param22 for proteins and amino acids) should be directed to
Charles L. Brooks, III at brooks@scripps.edu.
* Menu:
* Syntax:: Syntax of the GBORN Commands
* Function:: Purpose of each of the commands
* Examples:: Usage examples of the GBORN module
Syntax of Generalized Born Solvation commands
[SYNTAX: GBORn commands]
GBORn { P1 <real> P2 <real> P3 <real> P4 <real> P5 <real>
[LAMBda <real>] [EPSILON <real>] [CUTAlpha <real>] [WEIGht] }
{ CLEAr }
Parameters of the Generalized Born Model
P1-P6 The parameters P1, P2, ..., P5 specify the particular parameters
controlling the calculation of the effective Born radius for
a particular configuration of the biomolecule as determined from
the sum over all atoms:
Alpha(i) = [-CCELEC/(2*Lambda*R(i)) + P1*CCELEC/(2*R(i)^2
+ Sum {P2*V(j)/rij^4} + Sum {P3*V(j)/rij^4}
bonds angles
+ Sum {P4*V(j)*Cij/rij^4}]^(-1)*(-CCELEC)/2
non-bonded
with Cij = 1 when (rij^2)/(R(i)+R(j))^2 > 1/P5
and Cij = 1/4(1-cos[P5*PI*(rij^2)/(R(i)+R(j))^2])^2 otherwise.
Note: R(i), V(i) correspond to the vdW radius and volume respectively,
CCELEC is the conversion from e^2/A to kcal/mol, rij is the separation
between atom i and atom j.
Lambda This is the scaling parameter for the vdW radius in the first term
of the expression above.
Note: The parameters P1-P5 and Lambda correspond to parameters for a
particular CHARMM parameter/topology set.
***The parameters P1-P5 and Lambda are required input***
EPSILON This is the value of the dielectric constant for the solvent medium.
The default value is 80.
CUTAlpha This is a maximum value for the effective Alpha for any atom
during the calculation for a particular conformation of the
biomolecule. It is necessary because in some instances the
expression above for Alpha(i) can take on negative values
of numerical problems with the expression for very buried atoms
in large globular biopolymers. The default for this value is
10^6.
WEIGht This is a flag to specify that you want the vdW radii for the
atoms to be taken from the wmain array instead of the parameter
files (from Rmin values). The default is to use the parameter
values. These values are used for the R(i) and V(i) noted above.
CLEAr Clear all arrays and logical flags used in Generalized Born
calculation.
Examples
This examples are meant to be a partial guide in setting up
an input file for PBEQ. There is one test file, tstpbeq.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 maxi 200
set ener80 = ?ENPB
scalar wmain = radius
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
END
Example(3)
----------
This example shows how to set up a membrane potential.
Note that a non-zero concentration is required for a sensible
system with a membrane potential.
PBEQ
scalar wmain = radius
scalar wmain show
SOLVE epsw 80.0 ncel 150 dcel 0.5 conc 0.150 -
Tmemb 25.0 Zmemb 0.0 epsm 2.0 vmemb 0.100
RESET
END
CHARMM Documentation / Rick_Venable@nih.gov