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:

mds.in

plumed.dat

prmtop

ala_dp.rst

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