CHARMM c30b1 gbmv.doc

File: GBMV ]-[ Node: Top
Up: (commands.doc) -=- Next: Syntax\n

                 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 (

* 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\n


    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 is slower than the other GB methods in CHARMM.
However, GBMV can be made quite fast through the use of multiple time steps.

     GBMV is one the most accurate GB methods in the literature as of 2002.
For PARAM22, GBMV is the best (only) choice as of version 29. However, for 
PARAM19, GBORN and ACE may suffice.

Papers: M. S. Lee, F. R. Salsbury, Jr., and C. L. Brooks III. J. Chem. Phys.,
         116, 10606 (2002).
        M. S. Lee, M. Feig, F. R. Salsbury, Jr., and C. L. Brooks III. 
         J. Comp. Chem., submitted (2002)
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\n

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>} 
       CORR <int> 
        { ESHIFT <real> SHIFT <real> TT <real>   (CORR = 0) }             
        { SHIFT <real> SLOPE <real>              (CORR = 1) }

[method II: slower but more accurate (recommended)]

       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) }             
        { SHIFT <real> SLOPE <real>              (CORR = 1) } 

[ modify alpha update frequency ]

     { UPDATE <int> }

[ grid-based method ]

            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]


File: GBMV ]-[ Node: Function
Up: Top -=- Previous: Syntax -=- Next: Examples\n

          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)

CORR      Coloumb field correction method: 1 for R^7 method (default), 
          0 for R^5 method. The R^7 method using SHIFT/SLOPE. The R^5
          method used SHIFT/ESHIFT.

TT        Multiplicative factor for correction term (CORR = 0 only). 

SHIFt     The shifting factor of Alpha(i). CORR=0 or 1. MUST be set! 

ESHIft    Energy shifting factor of the self-polarization 
          energies: 1/Alpha(i). CORR=0 only. (Default 0.0)

SLOPE     Multiplicative factor of the Alpha(i). CORR=1 only. (Default 1)

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
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. For single points and
          minimizations, value of 1 is recommended. For dynamics,
          values of 5-10 may be used with Nose thermostatting (NOSE). 
          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

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 -20 to 
          -50 are recommended to GBMV II. Value of -20 is necessary 
          for stable dynamics in 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.

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)

          Parameters specific to GBMV I:

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


                 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.

          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\n

                  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.

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:


  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

  <Known Compatible with>
   - 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)
   - EWALD 
   - multiple dielectric
   - QUANTUM* (single energy with original charges is ok)
   - FLUCQ
   - GRID

CHARMM .doc Homepage

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute