CHARMM c33b2 gbmv.doc
File: GBMV -=- Node: Top
Up: (commands.doc) -=- Next: Syntax
Generalized Born using Molecular Volume (GBMV)
Solvation Energy and Forces Module
- and -
Surface Area
Questions and comments regarding GBMV should be directed to
Michael S. Lee c/o
Charles L. Brooks, III (brooks@scripps.edu)
* Menu:
* Description:: Description of GBMV and related commands
* Syntax:: Syntax of the GBMV Commands
* Function:: Purpose of each of the commands
* Examples:: Usage examples of the GBMV module
File: GBMV -=- Node: Description
Up: Top -=- Previous: Top -=- Next: Syntax
Background:
The GBMV module is a Generalized Born method for
mimicking the Poisson-Boltzmann (PB) electrostatic solvation energy. The PB
method for obtaining solvation energies is considered a benchmark for implicit
solvation calculations. However, the PB method is slow and the derivatives,
i.e. forces, are ill-defined unless one changes the definition of the m
olecular volume.
The Generalized Born equation, as prescribed by Still, et. al. allows
one to compute solvation energies very similar to the PB equations.
As it is an analytical expression, forces are available as well:
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
D = r^2 / (K_s * alpha * alpha )
ij ij i j
where K_s = 4 for Still's original equation, or 8 for modified equation.
The only problem is that one needs to calculate the alpha's, a.k.a.
Born radii for each atom, accurately. There are various methods available,
such as the GBORN, ACE, and PBEQ modules in CHARMM.
The GBMV method obtains the Born radii very accurately,
i.e, w/ greater than 0.99 correlation. It is available as three approaches:
1) grid-based (Most accurate)
2) analytical method I (Least accurate, fastest)
2) analytical method II (preferred for dynamics)
approach and an analytical method. The analytical method has derivatives
and thus can be used in molecular dynamics simulations. The grid-based method
has no derivatives, however, is the most accurate and still can be used in
energy ranking and Monte-Carlo methods.
When should you use GBMV?
Because the analytical and grid-based methods are quite accurate, the
parameters change very little when optimized for a particular force-field.
Hence, forcefields besides those of CHARMM can be used with GBMV with
out refitting of parameters. The GBMV method II approximates the molecular
surface directly which leads to very good agreement (<1% relative error)
with respect to electrostatic solvatione energies from standard Poisson
theory. The higher accuracy comes at a price, however, and GBMV is
slower than all of the other GB methods in CHARMM.
Papers: (1) M. S. Lee, F. R. Salsbury, Jr., and C. L. Brooks III. J. Chem. Phys.,
116, 10606 (2002).
(2) M. S. Lee, M. Feig, F. R. Salsbury, Jr., and C. L. Brooks III.
J. Comp. Chem., 24, 1348 (2003)
(3) M. Feig, W. Im, C. L. Brooks III.
J. Chem. Phys., 120, 903 (2004)
(4) S. Tanizaki, M. Feig, J. Chem. Phys., 122, 124706 (2005)
Surface Area
We have implemented a solvent accessible surface area (SASA) calculation
within the GB module. It is relatively FREE compared to the GB calculation
and 5 times faster than the exact calculation. It is accurate to within 1%
of the exact SASA and is much more accurate than the SASA module.
File: GBMV -=- Node: Syntax
Up: Top -=- Previous: Description -=- Next: Function
Syntax of Generalized Born Molecular Volume (GBMV) Solvation commands
[SYNTAX: GBMV commands]
[method I: faster but less accurate]
GBMV { P1 <real> P2 <real> LAMBda1 <real> DN <real> SHIFT <real>
WATR <real> BETA <real> EPSILON <real>
SA <real> SB <real> SCUT <real>}
[WEIGHT]
CORR <int>
{ ESHIFT <real> SHIFT <real> TT <real> (CORR = 0) }
{ SHIFT <real> SLOPE <real> (CORR = 1) }
}
[method II: slower but more accurate (recommended)]
GBMV { [GEOM] [ARITH] [WEIGHT] [FIXA]
BETA <real> EPSILON <real> DN <real> WATR <real>
LAMBDA1 <real> TOL <real> BUFR <real> MEM <int> CUTA <int>
HSX1 <real> HSX2 <real> ONX <real> OFFX <real>
ALFRQ <int> EMP <real>
P1 <real> P2 <real> P3 <real> P4 <real> P6 <real>
SA <real> SB <real> SCUT <real>
KAPPA <real>
WTYP <int> NPHI <int>
CORR <int>
{ ESHIFT <real> SHIFT <real> TT <real> (CORR = 0,2) }
{ SHIFT <real> SLOPE <real> (CORR = 1,2) }
{ A1 <real> A2 <real> A3 <real> (CORR = 2) }
GCUT <int>
RADG <int> <real ...>
FAST 1|0 SGBFRQ <int> SXD <real>
}
[HDGB method: heterogeneous dielectric / membrane model]
GBMV { { GBMV II options }
CORR 3
A1 <real> A3 <real> A4 <real> A5 <real>
UNEPS <int>
ZS <real> ZM <real> ZT <real> ST0 <real>
[HDGBRC]
}
[ modify alpha update frequency ]
{ UPDATE <int> }
[ grid-based method ]
GBMV GRID { [GEOM] [ARITH] [CONV] [WEIGHT]
EPSILON <int> DN <real> WATR <real>
P6 <real>
KAPPA <real>
WTYP <int> NPHI <int>
CORR <int>
{ SHIFT <real> SLOPE <real> (CORR = 1) }
{ ESHIFT <real> SHIFT <real> (CORR = 0) } }
}
[ free-up memory and/or start over]
GBMV CLEAr
File: GBMV -=- Node: Function
Up: Top -=- Previous: Syntax -=- Next: Examples
---------------------------------------------------------------
Parameters of the Generalized Born using Molecular Volume Model
common to all methods:
---------------------------------------------------------------
WTYP Angular integration grid type:
0 - Dodecahedron
1 - Spherical polar
2 - Lebedev (DEFAULT)
3 - Alternating octahedron/cube
NPHI Used when WTYP equals 1 or 2. When WTYP=1, it corresponds to number
of phi angles. When WTYP=2, it corresponds to size of
Lebedev grid, which can only have values of 6,26 (Default), and 38
at the present time.
CUTA Extent of radial integration points in Angstroms. (Default 20)
GCUT radial spacing of integration grid
1 - default spacing:
0.1 0.2 0.3 0.4 0.5 0.75 1.0 1.25
1.5 1.75 2.0 2.5 3.0 3.5 4.0 5.0
6.0 7.0 8.0 10.0 12.0 16.0 20.0
2 - finer spacing for small radii
0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.0
1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6
2.8 3.0 3.2 3.7 4.1 5.1 6.1 7.1
8.1 10.1 12.1 16.1 20.1
3 - custom grid, specify with RADG
RADG custom grid spacing, first argument is number of intervals
following arguments are interval limits
CORR Coloumb field correction method:
0 for R^5 method. use: SHIFT/ESHIFT/TT
alpha(i) = - 1/( r4 - TT * r5 + ESHIFT ) + SHIFT
1 for R^7 method (default) use: SHIFT/SLOPE
alpha(i) = SLOPE/( (1-1/sqrt(2)) * r4 + r7) + SHIFT
2 for R^5/R^7 use: A1/A2/A3/SHIFT/SLOPE/ESHIFT
alpha(i) = SLOPE/( A1 * r4 + A2 * r5 + A3 * r7 + ESHIFT) + SHIFT
this mode is intended for the calculation of Born radii
in different dielectric environments
3 for R^5/R^7 use: A1/A3/A4/A5/SLOPE
alpha(i) = SLOPE/( A1 * r4 + A3(i) * r7) + A4 + A5/(eps(i) + 1)
A3(i) = A3 * 3 * eps(i) / (3 * eps(i) + 2 * EPS)
this mode is intended for the implicit membrane model (see below)
where r4 is volume integral over 1/r^4, r5 is square root
of integral over 1/r^5 and r7 is integral over 1/r^7 to the
power of 1/4.
TT Multiplicative factor for correction term (CORR = 0 only).
SHIFt The shifting factor of Alpha(i). MUST be set!
ESHIft Energy shifting factor of the self-polarization
energies: 1/Alpha(i). CORR=0 or 2 only. (Default 0.0)
SLOPE Multiplicative factor of the Alpha(i). CORR=1 or 2 only. (Default 1)
A1,A2, A3 Multiplicative factors in calculation of Alpha(i).
WATR The radius of the water probe. Usually this is set to 1.4
Angstroms. If this were changed, other parameters would have
to be modified.
EPSILON This is the value of the dielectric constant for the solvent medium.
The default value is 80.
KAPPA Debye-Huckel ionic term: Units of inverse length (Angs). Default
is 0 (no salt).
GEOM Select geometric cross-term in Still equation (default).
ARITH Select arithmetic cross-term in Still equation.
P6 Exponent in exponential of Still equation. Default is 4, for
historical reasons. Value of 8 is RECOMMENDED for GEOM, 6.5 for
ARITH.
WEIGHT Use WMAIN array for radii. (Default uses vdW radii array)
CLEAr Clear all arrays and logical flags used in Generalized Born
calculation. Use command by itself.
-----------------------------------------------------------
Parameters specific to GBMV I and II:
-----------------------------------------------------------
FIXA Update alphas only if coordinates have changed more
than expected for finite differences. Useful for
static pka calculations. With FIXA keyword, finite-difference
wouldn't work correctly, hence it must be specified. Not
on by default.
ALFRQ Update frequency of Born radii. Use with caution!
One of LIMP,IMP, or EMP options must be selected. (Default 1)
LIMP Use ALFRQ*(dE/dalpha)(dalpha/dx) part of GB force every ALFRQ
steps. For ALFRQ <= 5.
EMP Decay constant of the impulse force. Default is 1.5, which
is meant for ALFRQ of 5. Generally, EMP ~= ALFRQ/4. For
ALFRQ <= 10. (Recommended option)
IMP Use (dE/dalpha)(dalpha/dx) part of GB force every ALFRQ
steps. Any ALFRQ can be used. Only meant for equilibrium
calculations.
DN The cell width of the lookup grid. Larger values make program
slower. Smaller values use up more memory. Default of 1.0 A is best
compromise between speed and memory.
BETA Smoothing factor for tailing off of volume.
Values of around -100 are fine for GBMV I. Values of -8 to
-50 are recommended to GBMV II. (Default -20)
In general for GBMV I, larger values will make the
calculation go faster but potentially introduce jumps in the
potential energy surface.
Smaller values of beta lead to more stable dynamics, but compromise
the agreement with Poisson theory. In GBMV II the choice of BETA
also affects P3. Recommended pairs of values for GBMV II are:
BETA = -20, P3 = 0.70
BETA = -12, P3 = 0.65
BETA = -10, P3 = 0.57
BETA = -8, P3 = 0.35
LAMBda The threshold value for the atomic volumes. In GBMV I, smaller
values produces shorter Born radii and wide variance w/respect
to accurate PB radii. Large values produce larger radii but smaller
variance. In GBMV II, value should be kept at 0.5.
BUFR Distance that any atom is allowed to move before lookup table
is rebuilt. Larger values lead to less lookup table update but larger
memory usage. Use 0.0 for static structure.
Values between 0.2 and 1.0 Angstrom. (Default 0.5)
MEM Percentage extra memory beyond hypothetical calculation of table
size. (Default 10)
TOL Accuracy of the switching function used to determine accuracy of the
first derivatives, i.e. forces. (Default 1e-8)
SA Surface area coefficient (KCAL/(MOL*A**2)). (Default 0.0)
SASA Energy term shows up under EXTERN/ASP.
SB Surface area constant (KCAL/MOL) (no effect on forces) (Default 0.0)
SON The startpoint for the switching function of each hard sphere.
(Default 1.2) Units in Angstroms
SOFF The endpoint for the switching function of each hard sphere.
(Default 1.5)
P1 The multiplicative factor for the exponent of the
quartic exponential atomic function:
Gamma(i) = P1 * log(lambda)/(Rad(i)^4)
Parameters specific to GBMV II:
P1,P2 Variables which affect the shape of the VSA atomic function in the
region of R to R+2.
F(x) = A^2 / (A + x^2 - R^2)^2
where
A = P1 * R + P2 (Defaults: P1 = 1.25/P2 = 0.45)
P3 Scaling factor of VSA function. Default = 0.7
P4 Scaling coefficient for correction term to Still's equation.
(set to 0.0 for now)
P5 Exponent to the Still correction term. (use default for now)
HSX1/HSX2 Start and stop of hard-sphere tail with R(vdW) as origin.
(Defaults: -0.125/0.25).
ONX/OFFX Start and stop of VSA tail. Increasing values up to 2.8 A
makes better accuracy, however slows calculation. Compromise
of 1.9/2.1 is default.
FAST Turns on fast GBMV routine.
SGBFRQ Update frequency of internal lookup list in fast GBMV mode
(Default 1). Values between 1 and 10 are recommended.
SXD Delta used in fast GBMV mode lookup buffer. (Default 0).
Recommended values between 0.1 and 0.5. Requires 'FAST 1'
-----------------------------------------------------------
Parameters specific to HDGB (CORR = 3):
-----------------------------------------------------------
UNEPS Unit number of an input file holding the delectric profile
values (Use -1 for the default profile). The format of this
input file is restricted. Comments are not allowed in a file.
The dielectric profile must be sampled in equal intervals.
The first line needs the number of sampling points n and the
sampling interval h (Angstrom). Two columns of the z coordinates
(Angstrom) and dielectric constants. The example is given
in test/data/hdgb_eps.dat.
A4,A5 Parameters in calculation of Alpha(i). A4 and A5 correspond
to the parameter D and E of Equation (15) in the reference (3)
respectively.
ZS,ZM,ZT Parameters for a switching function for the nonpolar energy.
ST0 ZS, ZM, ZT, and ST0 corresponds to Za, Zb, Zc, and C of
Equation (11) in the reference (4).
HDGBRC If this flag is specified, the radius will be corrected
upon insertion to an implicit membrane. At this point,
this option is only available for the static single point
calculations, and not for dynamics simulations. It need
to be compiled with HDGBRC option.
-----------------------------------------------------------
Parameters specific to Grid-based GBMV:
-----------------------------------------------------------
ML Number of surface points to carve out re-entrant surface
CONV Smear grid with cross-shaped blur function to improve accuracy
File: GBMV -=- Node: Examples
Up: Top -=- Previous: Function -=- Next: Top
Usage Examples and Compatibility
The examples below illustrate some of the uses of the generalized Born
Molecular Volume (GBMV) module. See c29test/gbmvtest.inp for more examples.
--------------------------------------
THERE ARE TWO REQUIREMENTS TO RUN GBMV
--------------------------------------
1) Coordinates MUST be defined for all atoms before invoking the GBMV keyword.
Otherwise, "infinite" grid is established which uses too much memory.
2) CUTOFF Parameters MUST be defined.
For non-infinite cutoffs, "switch" in nonbonded parameters is NECESSARY.
Example 1
!To perform a single-point energy calculation w/infinite cutoffs using
!GBMV I algorithm (any forcefield):
scalar wmain = radii
GBMV BETA -100 EPSILON 80 DN 1.0 WATR 1.4 TT 2.92 -
SHIFT -0.5 ESHIFT 0.0 LAMBDA1 0.1 P1 0.44 -
BUFR 0.5 Mem 20 CUTA 20 WTYP 0 -
WEIGHT ! Radii from wmain
ENERGY ctonnb 979 ctofnb 989 cutnb 999
Example 2
!To perform a single-point energy calculation w/infinite cutoffs using
!the GBMV II algorithm (any forcefield):
GBMV BETA -20 EPSILON 80 DN 1.0 watr 1.4 GEOM -
TOL 1e-8 BUFR 0.5 Mem 10 CUTA 20 HSX1 -0.125 HSX2 0.25 -
ALFRQ 1 EMP 1.5 P4 0.0 P6 8.0 P3 0.70 ONX 1.9 OFFX 2.1 -
WTYP 2 NPHI 38 SHIFT -0.102 SLOPE 0.9085 CORR 1
ENERGY ctonnb 979 ctofnb 989 cutnb 999
GBMV CLEAR ! Clear GB arrays
Example 3
!To perform a minimization w/ cutoffs using
!the GBMV II algorithm:
GBMV BETA -20 EPSILON 80 DN 1.0 watr 1.4 GEOM -
TOL 1e-8 BUFR 0.5 Mem 10 CUTA 20 HSX1 -0.125 HSX2 0.25 -
ALFRQ 1 EMP 1.5 P4 0.0 P6 8.0 P3 0.70 ONX 1.9 OFFX 2.1 -
WTYP 2 NPHI 38 SHIFT -0.102 SLOPE 0.9085 CORR 1
MINI sd nstep 100 ctonnb 12 ctofnb 14 cutnb 16 vswitch switch
!Then run some dynamics with multiple time step alpha update:
GBMV UPDATE 5
DYNAMICS vver start timestep 0.001 nstep 1000 nprint 100 iprfrq 100 -
firstt 298 finalt 298 ichecw 0 iasors 0 iasvel 1 isvfrq 1000 -
ntrfrq 100 - !GB is not rotationally invariant due to finite int. grid
tstruc 298 - !not a fixed water molecule
NOSE RSTN TREF 298.0 QREF 10 NCYC 10 ! UPDATE 5 is not energy-conserving
Example 4
!Minimal parameter specification (Analytical Method II):
GBMV P6 8.0 SHIFT -0.102 SLOPE 0.9085
ENERGY ctonnb 979 ctofnb 989 cutnb 999
Example 5
!Grid-based GBMV:
GBMV GRID EPSILON 80 DN 0.2 watr 1.4 GEOM P6 8.0 -
WTYP 0 NPHI 10 SHIFT -0.007998 SLOPE 0.9026 CORR 1 CONV
ENERGY ctonnb 979 ctofnb 989 cutnb 999
Example 6
! HDGB DPPC membrane
! If you want the default DPPC profile used in the reference (4),
! comment out the open file statement and set UNEPS to -1.
! The input file will be closed automatically, so you don't need
! the explicit close statement.
open unit 1 name eps.dat read form
GBMV A1 0.3255 A3 1.085 A4 -0.14 A5 -0.15 -
UNEPS 1 -
ZS 0.5 ZM 9.2 ZT 25 ST0 0.32
ENERGY ctonnb 979 ctofnb 989 cutnb 999
----------------------------------------------------------------------
<Known Compatible with>
- PARALLEL
- CONS FIX
- INTE
- PHMD
- VIBRAN (finite difference second derivatives)
- MMFF (WEIGHT keyword must be used)
<Known Incompatible with (so far)>
- VIBRAN (no analytic second derivatives)
- BLOCK (hence not compat. w/ PERT/PIMPLEM/PERTURB/REPLICA)
- IMAGE/CRYSTAL
- EWALD
- multiple dielectric
- QUANTUM* (single energy with original charges is ok)
- FLUCQ
- GAMESS
- GENBORN
- GRID
- PRESSURE
- SBOUND
<PREFX keywords required for compilation>
- GENBORN
- GBMV
- GBMVFAST
- HDGB
- HDGBRC (optional)
Up: (commands.doc) -=- Next: Syntax
Up: Top -=- Previous: Top -=- Next: Syntax
Up: Top -=- Previous: Description -=- Next: Function
Up: Top -=- Previous: Syntax -=- Next: Examples
Up: Top -=- Previous: Function -=- Next: Top