Valine to Alanine Free Energy determination

The files in this directory are a simple example of setting up and running a perturbation from a valine to an alanine,
namely ACE-VAL-NME -> ACE-ALA-NME

The molecule will be surrounded by a box of TIP3P water. The solvated system is generated with xleap and is equilibrated with sander using PME for long-range electrostatics.

This tutorial is not designed to teach MD or Free Energy theory or techniques, only to instruct on the usage of the programs: xleap, sander, and gibbs.

The overview of the steps detailed in this tutorial are as follows. First the system is constructed using xleap to produce the necessary input files (paramter and coordinate files) for sander and gibbs programs. Sander is then used to move the system away from high energy configurations by minimization. The minimized system is then equilibrated with constant temperature and pressure until the system is stable. When an equilibrated system is obtained, the Free Energy of changing a valine to alanine is determined using the gibbs module. These runs are fairly short and are not supposed to actually mimic the "true results" for this perturbation.

The largest part of this tutorial is dedicated to constructing of the system and its perturbation, and pre-conditioning the system for the free-energy sampling runs. Read the GIBBS section of the manual with special attention to the introduction (sections 7.1 to 7.5), and the Discussion at the end (section 7.8).

Step 1: Generate input files in xleap

Detailed tutorial for leap beginners Condensed form of instructions

Step 2: Minimize and Equilibrate with sander

Detailed tutorial for sander beginners Condensed form of instructions

Step 3: Free Energy Sampling with gibbs

Free Energy Perturbation (FEP)

The sample tutorial run of gibbs has 10 intervals from lambda=1 to 0, and only 700 steps of sampling. These values are only meant to demonstrate the use of gibbs and smaller intervals (dlambda) and more steps are necessary. Read the manual for guidance and for references to the methods, limitations, and suggested parameters. 
% gibbs -O -i in.gibbsF -o out.gibbsF -r crd.gibbsF -p parm_pert.top -c crd.md1
Description of input file

The final output from this run shows the following statistics. Note that the free energy for forward and reverse averages are very different but should be the same absolute value. This is due to the small number of windows and the lack of convergence (too few equilibration and sampling steps in each window). 
    Current Lambda =   0.000000
    Last F.E. update: Lambda =   0.000000  Step =      11000  Method = F.E.P.     
    Accumulated "forward" quantities (Regular)       
       Lam+d_lam = 0.100000    F_energy  =    5.13896
       ELEC =     25.474    NONB =      0.656    14NB =      3.150
       14EL =    -26.101    BADH =     -0.116
    Accumulated "reverse" quantities (Regular)       
       Lam-d_lam = 0.000000    F_energy  =   -3.19955
       ELEC =    -25.400    NONB =     -2.148    14NB =     -3.026
       14EL =     25.989    BADH =      0.110
    ------------------------------------------------------------------------------
Doubling the number of windows produces: 
    Current Lambda =   0.000000
    Last F.E. update: Lambda =   0.000000  Step =      21000  Method = F.E.P.     
    Accumulated "forward" quantities (Regular)       
       Lam+d_lam = 0.050000    F_energy  =    3.88508
       ELEC =     25.316    NONB =      0.357    14NB =      3.252
       14EL =    -25.803    BADH =     -0.130
    Accumulated "reverse" quantities (Regular)       
       Lam-d_lam = 0.000000    F_energy  =   -3.17159
       ELEC =    -25.307    NONB =     -1.398    14NB =     -3.305
       14EL =     25.821    BADH =      0.122
    ------------------------------------------------------------------------------
Sampling for 9000 steps instead of 700 steps over 20 windows: 
    Current Lambda =   0.000000
    Last F.E. update: Lambda =   0.000000  Step =     210000  Method = F.E.P.     
    Accumulated "forward" quantities (Regular)       
       Lam+d_lam = 0.050000    F_energy  =    4.78230
       ELEC =     24.487    NONB =      1.175    14NB =      3.215
       14EL =    -25.175    BADH =     -0.122
    Accumulated "reverse" quantities (Regular)       
       Lam-d_lam = 0.000000    F_energy  =   -4.07408
       ELEC =    -24.706    NONB =     -2.761    14NB =     -3.247
       14EL =     25.329    BADH =      0.113
    ------------------------------------------------------------------------------
Setting isldyn=1 in the input file to select "thermodynamic integration" at the same sampling rate gives: 
    Current Lambda =   0.000000
    Last F.E. update: Lambda =   0.050000  Step =     200000  Method = T.I.       
    Accumulated "forward" quantities (Regular)       
       Lambda   =  0.050000    F_energy  =    4.16968
       Enthalpy =    0.74212   T*Entropy =   -3.42756

Further refinement at the lambda=0 endpoint is necessary. One can try dynamically modified windows, smaller delta lamda steps, and various methods of changing the vdW and electrostatic contributions in as well as which energies and derivatives are included in the free energy accumulators. This simulation only includes free energy contributions arising from the energy and entropy from perturbed groups with non-perturbed groups, and thus is mainly the solvation energy in this case.