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
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).
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.
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
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
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
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.
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory