Aim:

The aim of this tutorial is to run a Temperature Accelerated Molecular Dynamics(TAMD)/ driven-Adiabatic Free Energy Dynamics (d-AFED)  simulation using the PLUMED-AMBER interface. We will set up a simple simulation of alanine dipeptide in vacuum, analyze the output, estimate free energies and its convergence from the simulation.

Theory:

Lets consider a system with N number of particles and the potential energy defined as U(r). Let q(r) be defined as collective variable (CV). Assuming that the system is maintained at temperature T. We wish to compute free energy which is given by 

f is an arbitrary constant. The probability that CV has a value s is given by

where Z is the configuration partition function for the system, given by

It is clear from equation (2) that if the  temperature T of the system is modified by some factor, the sampling of the CV space can be enhanced. The TAMD/d-AFED  is one of the method that enhances the sampling by increasing the CV temperature. 

This method is proposed by J.B. Abrams and M. E. Tuckerman in 2008 where the sampling of collective variables (CVs) space  is accelerated by modifying the T factor such that the temperature of CV is kept high. It is not done in coordinate space directly due to some limitations instead it is done in extended space (non-physical) where the auxiliary variables are introduced which coupled with the physical/real CVs through a harmonic potential. To maintain the  adiabaticity among the real CVs and auxiliary variables, the mass of the auxiliary variable are taken to be very high in comparison to the real CVs. The harmonic constant is taken to be high enough such that the motion of auxiliary variable should be exactly followed  by the real variables. The free energy at T (normal temperature) is given by the following reweighing equation

Setup and simulation:

The alanine dipeptide system will be used in vacuum as a model system. It is modelled using the ff14SB all-atom AMBER force field and the simulations were carried out using PLUMED-AMBER interface. The Collective variables used for this simulation are phi and psi torsion angles as shown in the figure.

We need amber and plumed input files to run the simulations. The files are given as follows

The Label phi and psi represents the two torsional angels as shown in the figure and the TORSION ATOMS shows their serial numbers for their respective atoms involved in the torsion angle. KAPPA is the harmonic constant used to connect the real and auxiliary variables and TAU defines the mass of these auxiliary variables to maintain the adiabaticity. The TEMP keyword tells the temperature at which the CVs need to be heated. The plumed generates the output files named as mentioned , COLVAR and COLVAR_cc files. The COLVAR_cc file prints the information of real variables and the COLVAR file prints the information about the auxiliary variables. The TEMP file prints the velocity of the auxiliary variables to check the temperature of the auxiliary variables. 

Once the topology file, input coordinates file, amber input file, plumed input like necessary files are ready, the sander _plumed.sh file is used to submit the job

qsub sander_plumed.sh

Analysis of results:

The simulation is ran at high temperature but our aim is to calculate the free energy at normal temperature. The reweighing code assists us to reweigh it and you to make a room temperature two dimensional free energy landscape. follow these command to reweigh the data as follows:

cd reweighting/

create a link for the COLVAR file.

cp ../COLVAR .

We need to edit the COLVAR file because it contains some extra lines from PLUMED simulation.

sed -i  '1,6d'  COLVAR

The above command will delete the initial 6 lines from the COLVAR file. Now run the reweighting code by the following command. 

./tamd_reweighting.x

This run will generate the free_energy.dat file where the data for free energy is printed. You can use the plotting script named as fes_script.gnp to plot the free energy landscape.

type gnuplot in the terminal and enter

load 'fes_script.gnp'

The two-dimensional free energy surface (FES) should look like the following if enough sampling is done. 

Calculating one dimensional projection and checking convergence:

This folder contains a script which checks the convergence of free energy along phi and psi CVs. Use the following command.

sh convergence.sh

The one dimensional projection is plotted from two dimensional FES as a function of time to check the convergence. follow the following commands

cd convergence_check/

View the one-dimensional projections along phi and psi using the following commands. Copy the scripts using the following commands

cp ../phi_script.gnp .

cp ../psi_script.gnp .

load gnuplot in the terminal followed by the following commands.

load 'phi_script.gnp'

quit

load 'psi_script.gnp'

The free energy profile is converged after 10ns of simulation. The sharp peaks along phi shows that the higher free energy region is not sampled properly at 1500K. 

Try Yourself:

Is it more efficient to use a temperature of 3000K for the auxiliary variables? 

How the KAPPA values  were chosen? What about taking higher values of KAPPA?

How the values of TAUs are chosen? How would it affect the results?