We will play with a toy system, alanine dipeptide simulated in vacuo using the AMBER99SB force field (see Fig. \ref trieste-4-ala-fig).
This rather simple molecule is useful to benchmark data analysis and free-energy methods.
This system is a nice example because it presents two metastable states separated by a high free-energy barrier.
It is conventional use to characterize the two states in terms of Ramachandran dihedral angles, which are denoted with \f$ \Phi \f$ and \f$ \Psi \f$ in Fig. \ref trieste-4-transition-fig .
\anchor trieste-4-ala-fig
\image html belfast-2-ala.png "The molecule of the day: alanine dipeptide."
\anchor trieste-4-transition-fig
\image html belfast-2-transition.png "Two metastable states of alanine dipeptide are characterized by their Ramachandran dihedral angles."
\section trieste-4-ex-1 Exercize 1: my first metadynamics calculation
In this excercise we will setup and perform a well-tempered metadynamics run using the backbone dihedral \f$ \phi \f$
as collective variable. During the calculation, we will also monitor the behavior of the other backbone dihedral \f$ \psi \f$.
Here you can find a sample `plumed.dat` file that you can use as a template.
Whenever you see an highlighted \highlight{FILL} string, this is a string that you should replace.
\plumedfile
# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C
phi: TORSION ATOMS=__FILL__
# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N
psi: TORSION ATOMS=__FILL__
# Activate well-tempered metadynamics in phi
metad: __FILL__ ARG=phi...
# Deposit a Gaussian every 500 time steps, with initial height equal
# to 1.2 kJoule/mol, biasfactor equal to 10.0
PACE=500 HEIGHT=1.2 BIASFACTOR=10.0
# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run
SIGMA=__FILL__
# Gaussians will be written to file and also stored on grid
FILE=HILLS GRID_MIN=-pi GRID_MAX=pi
...
# Print the collective variables and the value of the bias potential on COLVAR file
PRINT ARG=phi,psi,__FILL__ FILE=COLVAR STRIDE=10
\endplumedfile
Once your `plumed.dat` file is complete, you can run a XX-ns long metadynamics simulations with the following command
\verbatim
> gmx mdrun -s topol.tpr -plumed plumed.dat
\endverbatim
During the metadynamics simulation, PLUMED will create two files, named COLVAR and HILLS.
The COLVAR file contains all the information specified by the PRINT command, in this case
the value of the CVs every 10 steps of simulation, along with the current value of the metadynamics bias potential.
We can use `gnuplot` to visualize the behavior of the CV during the simulation, as reported in the COLVAR file:
\verbatim
gnuplot> p "COLVAR" u 1:2
\endverbatim
\anchor trieste-4-phi-fig
\image html munster-metad-phi.png "Time evolution of the CV phi during the first 2 ns of a metadynamics simulation of alanine dipeptide in vacuum."
By inspecting Figure \ref trieste-4-phi-fig, we can see that the system is initialized in one of the two metastable
states of alanine dipeptide. After a while (t=0.1 ns), the system is pushed
by the metadynamics bias potential to visit the other local minimum. As the simulation continues,
the bias potential fills the underlying free-energy landscape, and the system is able to diffuse in the
entire phase space.
The HILLS file contains a list of the Gaussians deposited along the simulation.
If we give a look at the header of this file, we can find relevant information about its content:
\verbatim
#! FIELDS time phi psi sigma_phi sigma_psi height biasf
#! SET multivariate false
#! SET min_phi -pi
#! SET max_phi pi
#! SET min_psi -pi
#! SET max_psi pi
\endverbatim
The line starting with FIELDS tells us what is displayed in the various columns of the HILLS file:
the time of the simulation, the value of phi and psi, the width of the Gaussian in phi and psi,
the height of the Gaussian, and the biasfactor.
We can use the HILLS file to visualize the decrease of the Gaussian height during the simulation,
according to the well-tempered recipe:
\anchor trieste-4-phihills-fig
\image html munster-metad-phihills.png "Time evolution of the Gaussian height."
If we look carefully at the scale of the y-axis, we will notice that in the beginning the value
of the Gaussian height is higher than the initial height specified in the input file, which should be 1.2 kJoule/mol.
In fact, this column reports the height of the Gaussian rescaled by the pre-factor that
in well-tempered metadynamics relates the bias potential to the free energy.
One can estimate the free energy as a function of the metadynamics CVs directly from the metadynamics
bias potential. In order to do so, the utility \ref sum_hills should be used to sum the Gaussians
deposited during the simulation and stored in the HILLS file.
To calculate the free energy as a function of phi, it is sufficient to use the following command line:
\verbatim
plumed sum_hills --hills HILLS
\endverbatim
The command above generates a file called fes.dat in which the free-energy surface as function
of phi is calculated on a regular grid. One can modify the default name for the free energy file,
as well as the boundaries and bin size of the grid, by using the following options of \ref sum_hills :
\verbatim
--outfile - specify the outputfile for sumhills
--min - the lower bounds for the grid
--max - the upper bounds for the grid
--bin - the number of bins for the grid
--spacing - grid spacing, alternative to the number of bins
\endverbatim
The result should look like this:
\anchor trieste-4-metad-phifes-fig
\image html munster-metad-phifes.png "Estimate of the free energy as a function of the dihedral phi from a 10ns-long well-tempered metadynamics simulation."
To assess the convergence of a metadynamics simulation, one can calculate the estimate of the free energy as a function
of simulation time. At convergence, the reconstructed profiles should be similar.
The option --stride should be used to give an estimate of the free energy every N Gaussians deposited, and
the option --mintozero can be used to align the profiles by setting the global minimum to zero.
