CHARMM c24 pert.doc



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



                    Free Energy Perturbation Calculations


    The PERTurbe command allows the scaling of energy between PSFs for use in
energy analysis, comparisons, slow growth free energy simulations,
widowing free energy simulation, and for slow growth homology modelling.
This is a rather flexible implementation of free energy perturbation
that allows connectivity to change.  Also, three energy restraint
terms (harmonic, dihedral and NOE) and the SKIP command flags are subject
to change which allows a flexible way in which to compute free energy
differences between different conformations.  This code in implemented
in a non-intrusive manner and works with all minimizers and
integrators.  SHAKE can now be applied to bonds which are mutated as
well; an appropriate constraint corrections is calculated
automatically in these cases.


* Menu:

* Syntax::           Syntax of PERT Commands
* Description::      Description of PERT Commands
* Restrictions::     Restrictions in PERT Command usage
* References::       References regarding Free Energy Perturbations
* Examples::         A Sample Input Files
* Constraints::      Special considerations regarding SHAKE



File: Pert ]-[ Node: Syntax
Up: Top -=- Next: Description -=- Previous: Top


                Syntax of Free Energy Perturbation Commands

[Syntax PERT]

PERTurb   [OFF] [INBFrq int nonbond-specs] [RESEt] atom-selection
                [INBFrq  0               ]

      The PERT OFF command disables the free energy routines.  The
atom-selection indicated which atoms have changed.  This is to make
the calculation run more efficiently.  If only a small percent of the
atoms have changed, this doubles the performance.  The nonbond specs
are included to make sure the nonbond exclusion lists are properly
setup.  This then allows the connectivity to change during the
simulation.  INBFrq should not be set to zero here unless the
exclusion lists have already been setup in a previous command.

----------------------------------------------------------------------------
ENERgy   ...   [ RESET ]  [ free-energy-step-spec ]
DYNAmics ...              [ PUNIt integer         ]    
MINImize ...

  [ RESET ]     ! Resets all all accumulation data and counters.
                 (automatic for the first PERT or after a PERT OFF command)

free-energy-step-spec::=
   [PWINdow  [LAMBda ] ] [PSTArt int] [PSTOp int] [LSTArt real] [LSTOp  real] -
   [PSLOwgrowth        ]   
                 [PINCrement int]  [PEQUilibrate int]  [LAVErage]  [LINCrement]


  [PWINdow   ]  ! specifies the windowing algorithm (default)
  [PSLOwgrowth] ! specifies the slow growth algorithm
  
  [LAMBda    ]  ! specifies the lambda value for windowing methods or for energy
                  or minimization calculations.

  [PSTArt int]  ! starting dynamics step number for accumulation (default 1)
  [PSTOp  int]  ! ending dynamics step number for accumulation (default 0)

  [LSTArt real] ! specifies the starting lambda value (default 0.0)
  [LSTOp  real] ! specifies the final lambda value (default 1.0)

  [PINCrement int] ! specifies number of steps to next window (auto mode).
  [PEQUil int]  ! specifies number of steps for equilibration (auto mode).
  [LAVErage]    ! specifies that lambda = (LSTART+LSTOP)/2    (auto mode).
  [LINCrement]  ! Specifies the lambda increment between windows (auto mode).

  [ END   ]     ! Turns off the free energy code
  

The PSTArt and PSTOp values are relative to the number of dynamics steps
since PERT command was first enabled, or if a PERT RESET command is used
(regardless of whether the DYNAmics commands is run more than once or
whether the dynamic run involved the use of restart files).

By specifying the auto mode parameters (PINCrement, PEQUilibrate, LINCrement),
a new window will start at the conclusion of the current window with modified
parameters.  Also, in auto mode, the run will terminate at the end of a window
where the LSTOP value is 1.0 or 0.0.

The PUNIt option allows the free-energy-step-spec to be specified 
more than once and acts as a scheduler for a particular simulation.
The format of the PUNIt file is;

* title
*
repeat-lines-of(free-energy-step-specs)
END

The end is optional and terminates the free energy run.
Lack of an END (i.e. and end-of-file or blank lines) will put PERT into
auto mode which will continue until LSTOP becomes 1.0 or 0.0 (based on
the LINCrement value).



File: Pert ]-[ Node: Description
Up: Top -=- Previous: Syntax -=- Next: Restrictions


                      Description of PERT Commands

      The PERTurb command copies and saves the current PSF and restraint
data for harmonic, NOE and dihedral restraints to the initial (lambda=0)
saved state.  The SKIP command flags are also saved to allow linear scaling
of entire energy terms.

      The structure may then be modified or perturbed with patches,
SCALar commands, with the DELEte command, or by generating or reading
a new PSF.

The Basic mode of operation is;

....
PERTurb        ! Define the lambda=0 state.
PATCH ....     ! Define the lambda=1 state.
DYNAmics ....  ! Run MD on intermediate surface...
....

The PSF in use when dynamics or energy minimization is invoked
becomes the final (lambda=1) state.  The actual energy computed
is a linear combination of these two endpoints.

      The PATCH command may be replaced with any other command that
modifies the PSF.  Some examples which modify the PSF;

SCALAR CHARGE SET -0.55 SELE ATOM A 1 O* SHOW END ! change a charge

DELETE ATOM SELE ALL END
READ SEQUENCE ....
GENERATE ...    ! generate a new different PSF.

DELETE CONNECTIVITY ....  ! modify the PSF by changing the connectivity.

SCALAR TYPE SET 14 SELE ATOM A 1 O SHOW END ! change the vdw atom type 

OPEN ....     ! Read a new PSF
READ PSF ... 

      It is not required that the PSF be modified.  If one wants to
carry out coordinate perturbation only, it is sufficient to modify
the harmonic restraints, the NOE restraints, or the dihedral restraints.
In this way, it is possible to calculate the free energy differences
between different conformers.  (However, this option should not be
used with simulataneous change of SHAKE constraints)

      Note that with this implementation, because two PSFs are used,
that the connectivity may change.  The use of 1-4 interactions and
nonbond exclusions is fully supported.  This allows this method to be
used for examining changes that involve bond changes, such as cystine
bridge formation.

The advanced mode of operation is;

....
PERTurb
PATCH ....
DYNAmics ....
....
PERTurb
PATCH ....
DYNAmics ....
....
PERTurb
PATCH ....
DYNAmics ....
....

In this way, several changes can be affected in a single CHARMM run.
For example, the first patch might be the removal of charges, and the
second patch could correspond to a change in atom size, and the third
patch could simply consist of modifying dihedral restraints so as to
affect a conformational change.  The free energy differences and
fluctuations will be calculated for each window as well as the total
for all previous windows.


File: Pert ]-[ Node: Restrictions
Up: Top -=- Previous:Description -=- Next: References


RESTRICTIONS:

      The number of atoms in both sets must match!  If the system of
interest has different numbers of atoms, then dummy atoms must be
used.  The mapping of atoms between the first and last structure is
one to one. 

The following CHARMM features are not currently supported for use
with free energy perturbation;

INTEraction_energy

These commands will continue to work, but will only use the final
(lambda=1) structure.  Most other energy related CHARMM features
are supported.  The IMAGE facility will be supported in the near
future.

The following CHARMM energy related features cannot be modified with
the PERT command (e.g. cannot be part of what is changing, and are
only determined by the final state).

HBON - hydrogen bond energy
ST2  - ST2 energy
CIC  - internal coordinate constraint energy
CDRO - quartic droplet potential energy
USER - user supplied energy (USERLINK)
RXNF - Reaction field energy
IMNB - image van der Waal energy
IMEL - image electrostatic energy
IMHB - image hydrogen bond energy
IMST - image ST2 energy
SBOU - solvent boundary energy
UREY - Urey Bradley energy terms
XTLV - Crystal vdw terms
XTLE - Crystal electrostatics
EXTE - Extended electrostatics



File: Pert ]-[ Node: References
Up: Top -=- Previous: Restrictions -=- Next: Examples


Some References:

M Mezei and D.L. Beveridge, in Annals of the NYAS, "Free Energy
Simulations" 482 (1986)

T. P. Straatsma Ph.D. Thesis "Free Energy Evaluation by
Molecular Dynamics Simulations"

Kollman, P. A.; et al.  J. Am. Chem. Soc. 1987, 109, 1607.

Kollman, P. A.; et al.  J. Am. Chem. Soc. 1987, 109, 6283.

Kollman, P. A.; et al.  J. Chem. Phys. 1989, 91, 7831.

Bell, C. D.; Harvey, S. C.,   J. Phys. Chem. 1986, 90, 6595.

van Gunsteren, W.F. et al. in: Computer Simulation of Biomolecular
Systems: Theoretical and Experimental Applications, vol. 2, eds. van
Gunsteren W.F. and Weiner P.K. (Escom, Leiden, 1994), p. 349



File: Pert ]-[ Node: Examples
Up: Top -=- Previous: References -=- Next: Constraints


Examples:

The input file:

* A SIMPLE TEST RUN FOR PERT
*
bomlev -1
OPEN READ FILE UNIT 1 NAME ~/c22pt/toph19.mod
READ      RTF   UNIT  1
OPEN READ FILE UNIT 2 NAME ~/c22pt/param19.mod
READ      PARAMETER  UNIT  2
READ      SEQUENCE  CARD
*  FIRST SEQUENCE FOR SECOND DERIVATIVE TEST
*
    2
AMN CBX
GENERATE  A
GENERATE  B DUPLICATE A

OPEN UNIT 3 READ CARD NAME perttest.crd
READ COOR CARD UNIT 3

! modify the charge for the lambda=0 state
SCALAR CHARGE SET -0.55 SELE ATOM A 1 O* SHOW END

! minimize initial state so initial forces will be small.
MINI ABNR NSTEP 100 CTOFNB 12.0 CUTNB 14.0

PERT  ! save all PSF data for the lambda=0 state

! modify the charge again for the lambda=1 state
SCALAR CHARGE SET -0.15 SELE ATOM A 1 O* SHOW END

! carry out free energy run from first to final state

OPEN READ CARD UNIT 88 NAME perttest.punit
DYNA VERLET STRT  NSTEP 12000 TIMESTEP 0.001 -
    IPRFRQ 100 IHTFRQ 0 IEQFRQ 100 NTRFRQ 2000  -
    IUNCRD 50  ISEED 314159  -
    NPRINT 100 NSAVC 0 NSAVV 0 INBFRQ 25 IHBFRQ 0 -
    CTOFNB 12.0 CUTNB 14.0 -
    FIRSTT 300.0 FINALT 300.0 TEMINC 100.0   -
    IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 20.0 TWINDL -20.0 -
    PUNIT 88

PERT OFF
energy ! just a check at lamda=1
STOP

The punit file:

* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
  LSTART 0.0  LAMBDA 0.0  LSTOP 0.05  PSTART  12000  PSTOP  20000  PWIND
  LSTART 0.05 LAMBDA 0.1  LSTOP 0.15  PSTART  22000  PSTOP  30000  PWIND
  LSTART 0.15 LAMBDA 0.2  LSTOP 0.25  PSTART  32000  PSTOP  40000  PWIND
  LSTART 0.25 LAMBDA 0.3  LSTOP 0.35  PSTART  42000  PSTOP  50000  PWIND
  LSTART 0.35 LAMBDA 0.4  LSTOP 0.45  PSTART  52000  PSTOP  60000  PWIND
  LSTART 0.45 LAMBDA 0.5  LSTOP 0.55  PSTART  62000  PSTOP  70000  PWIND
  LSTART 0.55 LAMBDA 0.6  LSTOP 0.65  PSTART  72000  PSTOP  80000  PWIND
  LSTART 0.65 LAMBDA 0.7  LSTOP 0.75  PSTART  82000  PSTOP  90000  PWIND
  LSTART 0.75 LAMBDA 0.8  LSTOP 0.85  PSTART  92000  PSTOP 100000  PWIND
  LSTART 0.85 LAMBDA 0.9  LSTOP 0.95  PSTART 102000  PSTOP 110000  PWIND
  LSTART 0.95 LAMBDA 1.0  LSTOP 1.0   PSTART 112000  PSTOP 120000  PWIND
  END

Or equivalently using auto mode:

* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
 LSTART 0.0  LAMBDA 0.0  LSTOP 0.05  PSTART  12000  PSTOP  20000  PWIND PEQUIL 2000 PINCR 10000 LINCR 0.1


Or also equivalently as:

* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
  LSTART 0.0  LAMBDA 0.0  LSTOP 0.05  PSTART  12000  PSTOP  20000  PWIND
  LSTART 0.05 LAMBDA 0.1  LSTOP 0.15  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.15 LAMBDA 0.2  LSTOP 0.25  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.25 LAMBDA 0.3  LSTOP 0.35  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.35 LAMBDA 0.4  LSTOP 0.45  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.45 LAMBDA 0.5  LSTOP 0.55  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.55 LAMBDA 0.6  LSTOP 0.65  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.65 LAMBDA 0.7  LSTOP 0.75  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.75 LAMBDA 0.8  LSTOP 0.85  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.85 LAMBDA 0.9  LSTOP 0.95  PEQUIL   2000  PINCR  10000  PWIND
  LSTART 0.95 LAMBDA 1.0  LSTOP 1.0   PEQUIL   2000  PINCR  10000  PWIND
  END



File: Pert ]-[ Node: Constraints
Up: Top -=- Previous: Examples -=- Next: Top


	If SHAKE is applied to bond terms which are changed as the
result of an alchemical mutation a constraint correction is calculated
where required in slow-growth mode and TI in windowing mode.  The
exponential formula in windowing mode is not supported.  The user has
to beware of subtle problems regarding a possible "moment of inertia"
term that may be or may be not included in this correction (S. Boresch
& M. Karplus, to be published) In order for the constraint correction
to work properly, attention has to be given to the following points:

(1) SHAKE must not be applied to angle terms
(2) the PARA option has to be used (it's not clear, how to support
    reference coordinates in the context of an alchemical mutation)
(3) the SHAKE command has to issued after the PERT command.  (This
    is similar to setting up IMAGEs in connection with PERT).  A
    typical input will look something like

	PERT SELE ... END
	!change psf; after ALL changes have been made
	SHAKe BOND PARA
	DYNA ... ! carry out MD simulations etc.
	PERT OFF

(4) One should not mix situations where a constraint correction for
SHAKE is necessary with the use of harmonic, NOE and dihedral
restraints to calculate conformational free energy differences.

Items (2) and (3) can lead to error messages in situations where there
is actually no problem, e.g. you just want to apply SHAKe to your
solvent which is not affected by the mutation, so you specify SHAKe
before PERT and "bomb".  Nevertheless, I thought better safe than
sorry and if wanted one can override the warnings with a BOMLEV -3.
Item (4) is simply untested.

CHARMM .doc Homepage


Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory