Metadynamics
Introduction
The aim of this tutorial is to run a Metadynamics simulation using the PLUMED-AMBER interface. We will set up a simple simulation of alanine dipeptide in vacuum and reconstruct the free energy surface using the two dihedral angles phi and psi
Alanine dipeptide system
Alanine dipeptide system is 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 below
Metadynamics
Setup and running a simulation: In order to run this simulation we need to prepare the required inputs for AMBER and PLUMED both separately.
The following files will be required for running the simulation:
01-sander_plumed.sh
mds.in
ALanine dipeptide
&cntrl
imin=0, ntx=1,
irest=0, nstlim=5000000,
dt=0.001,
ntf=2, ntc=2,
tempi=300.0, temp0=300.0,
ntpr=100, ntwx=100,
cut=999.0,
ntp=0, ntt=3,
gamma_ln=2.0,
ig=-1,
ntb=0, igb=0,
plumed=1,
plumedfile='plumed.dat'
&end
plumed.dat
phi: TORSION ATOMS=5,7,9,15
psi:TORSION ATOM=7,9,15,17
metad: METAD ARG=phi,psi PACE=500 HEIGHT=2.4 SIGMA=0.25,0.25 FILE=HILLS
PRINT STRIDE=100 ARG=phi,psi FILE=COLVAR
Here, the labels phi and psi define the torsion φ and ψ as shown in the figure of alanine dipeptide. The serial numbers defines the atoms involved in making the torsion angles. The label metad has the keywords METAD is used to start the metadynamics, ARG defines the CVs that will be used , PACE is used for depositing the two-dimensional Gaussian at every 500 steps with height equal to 2.4 kJ/mol and width 0.5 rad.
The output will be COLVAR and HILLS files. The information that these files contain can be seen from the first line printed on each of the files.
Choosing Parameters to run the simulation for Metadynamics:
- How to choose height of the Gaussian(HEIGHT)?
- How to choose width of the Gaussian(SIGMA)?
- How to choose the frequency at which the Gaussians would be deposited(PACE)?
To start the simulation we need input file from AMBER mds.in, the toplogy file prmtop, restart file ala_dp.rst and plumed file plumed.dat. Once the necessary files are ready, the script 01-sander_plumed.sh will be used for running the simulation.
Calculating free energy landscape along CVs
Once the simulation is over, the free energy as a function of φ and ψ can be be generated by using sum_hills tool of plumed. The HILLS file contains information about the bias. From there,one can reconstruct free energy surface by summing up the Gaussians deposited at each value of CVs. For this, the following can be used
plumed sum_hills --hills HILLS --mintozero
This will create an output file fes.dat from which the free energy can be plotted in gnuplot by loading the script plot_fes.gnp
Free energy convergence
The free energy can be plotted at different simulation time in order to check the convergence. For this we need to use the following command of plumed
plumed sum_hills --hills HILLS --idw phi --kt 2.4 --stride 2500 --mintozero
This will plot the free energy as a function of phi and print it in a file after every 2500 gaussian deposition. Thus we can plot the free energy at different time in gnuplot by using the following script 03-plot_convergence.gpi
Well-Tempered Metadynamics
In well-tempered metadynamics, the bias is added with progressive decrease in the height of the Gaussian that has already been added at the
particular value of the CV. In order to run WT-MTD simulation, we need to add the keyword BIASFACTOR and TEMP in the input file of plumed in the line METAD. The input should look something like this
phi: TORSION ATOMS=5,7,9,15
psi:TORSION ATOM=7,9,15,17
metad: METAD ARG=phi,psi PACE=500 HEIGHT=2.4 SIGMA=0.5,0.5 BIASFACTOR=6 TEMP=300.0
PRINT STRIDE=100 ARG=phi,psi FILE=COLVAR
After the simulation is complete we can generate the two dimensional free energy surface by the same command used for the metadynamics:
Free energy convergence
The free energy can be plotted at different simulation time in order to check the convergence. For this we need to use the following command of plumed
plumed sum_hills --hills HILLS --idw phi --kt 2.4 --stride 2500 --mintozero
This will plot the free energy as a function of phi and print it in a file after every 2500 gaussian deposition. Thus we can plot the free energy at different time in gnuplot by using the following script 03-plot_convergence.gnp