SCHEDULED DOWNTIME: GitLab's virtual host must be rebooted to apply security updates. The reboot will occur on Friday, 2025-04-25, between 7:00 and 8:00 PM.
There was an error fetching the commit references. Please try again later.
where -replex 1000 indicates that every 1000 molecular-dynamics steps
where -replex 2000 indicates that every 2000 molecular-dynamics steps
all replicas are randomly paired (e.g. 0-2 and 1-3) and exchanges are attempted
all replicas are randomly paired (e.g. 0-2 and 1-3) and exchanges are attempted
between each pair (as printed in the GROMACS *.log files).
between each pair (as printed in the GROMACS *.log files).
The same simulation can be run using WELLTEMPERED metadynamics.
\subsection bxcon Convergence of the Simulations
In the resources for this tutorial you can find the results for a 40ns long Well-Tempered Bias Exchange simulation. First of all we
can try to assess the convergence of the simulations by looking at the profiles. In the "convergence" folder there is a script that
calculates the free energy from the HILLS.0 file at incresing simulation lengths (i.e. every more 0.8 ns of simulation). The
scripts also generate two measures of the evolution of the profiles in time:
1. time-diff.HILLS.0: where it is stored the average deviation between two successive profiles
2. KL.HILLS.0: where it is stored the average deviation between profiles correctly weigheted for the free energy of the profiles themselves (Symmetrized Kullback-Lieber divergence)
From both plots one can deduce that after 8 ns the profiles don't change significantly thus suggesting that averaging over the range 8-40ns should result in
a accurate profile (we will test this using metagui).
\subsection metagui Bias-Exchange Analysis with METAGUI
In principle Bias-Exchange Metadynamics can give as a results only N 1D free energy profiles. But the information
contained in all the replicas can be used to recover multidensional free energy surfaces in >=N dimensions. A simple
way to perform this analysis is to use METAGUI. METAGUI performs the following operations:
1. Clusterizes the trajectories on a multidimensional GRID defined by at least the biased coordinates.
2. By using the 1D free energy profiles and the clustering assigns a free energy to the cluster using a WHAM procedure.
3. Lets the user visualize the clusters.
4. Approximates the kinetics among clusters.
METAGUI (Biarnes et. al) is a plugin for VMD that implements the approch developed by Marinelli et. al 2009. It can be
downloaded from the PLUMED website.
In order for the colvar and hills file to be compatible with METAGUI their header must be formatted as following:
COLVAR.#:
\verbatim
#! FIELDS time cv1 cv2 cv3 cv4
#! ACTIVE 1 1 A
#! ..
...
\endverbatim
NOTE:
1. the COLVAR.# files should contain ALL the collective variables employed (all those biased in at least one replica plus those additionaly analysed). They MUST be named cv1 ... cvN.
2. the COLVAR.# files must be synchronised with the trajectories, this means that for each frame in the trajectory at time t there must be a line in each colvar at time t and viceversa. The best option is usually to analyse the trajectories a posteriori using plumed driver.
3. a keyword #! ACTIVE NBIASEDCV BIASEDCV LABEL is needed, where NBIASEDCV is the number of biased cv in that replica (not overall), BIASEDCV is the index of the biased cv in that replica (i.e. 1 for the first replica and so on); LABEL is a letter that identify the replica (usually is simply A for the first, B for the second and so on) this is usufull if two replicas are biasing the same collective variable:
\verbatim
COLVAR.0:
#! FIELDS time cv1 cv2 cv3
#! ACTIVE 1 1 A
#! ..
...
COLVAR.1:
#! FIELDS time cv1 cv2 cv3
#! ACTIVE 1 2 B
#! ..
...
COLVAR.2:
#! FIELDS time cv1 cv2 cv3
#! ACTIVE 1 2 C
#! ..
...
COLVAR.3:
#! FIELDS time cv1 cv2 cv3
#! ACTIVE 0
#! ..
...
\endverbatim
In the above case Replica 0 biases cv1; replicas 1 and 2 biases cv2 while replica 3 is a neutral (unbiased) replica. cv3 is unbiased in all the replicas.
The ACTIVE keyword must be the FIRST LINE in the HILLS.# files:
HILLS.#:
\verbatim
#! ACTIVE 1 1 A
#! FIELDS time cv1 sigma_cv1 height biasf
#! ..
...
\endverbatim
The above notes hold for the HILLS files as well. In the folder metagui the script check_for_metagui.sh checks if the header of your file is compatible
with METAGUI, but remember that this is not enough! Synchronisation of COLVAR and trajctory files is also needed. HILLS files can be written with a
different frequency but times must be consistent.
\verbatim
./check_for_metagui.sh ../COLVAR.0
\endverbatim
will tell you that the ACTIVE keyword is missing, you need to modify all the header BEFORE proceeding with the tutorial!!
In the metagui folder there is a metagui.input file:
\verbatim
WHAM_EXE wham_bemeta.x
BASINS_EXE kinetic_basins.x
KT 2.4900
HILLS_FILE ../HILLS.0
HILLS_FILE ../HILLS.1
HILLS_FILE ../HILLS.2
HILLS_FILE ../HILLS.3
GRO_FILE VIL.pdb
COLVAR_FILE ../COLVAR.0 ../traj0.xtc "psi-1"
COLVAR_FILE ../COLVAR.1 ../traj1.xtc "phi-2"
COLVAR_FILE ../COLVAR.2 ../traj2.xtc "psi-2"
COLVAR_FILE ../COLVAR.3 ../traj3.xtc "phi-3"
TRAJ_SKIP 10
CVGRID 1 -3.1415 3.1415 15 PERIODIC
CVGRID 2 -3.1415 3.1415 15 PERIODIC
CVGRID 3 -3.1415 3.1415 15 PERIODIC
CVGRID 4 -3.1415 3.1415 15 PERIODIC
ACTIVE 4 1 2 3 4
T_CLUSTER 0.
T_FILL 8000.
DELTA 4
GCORR 1
TR_N_EXP 5
\endverbatim
where are defined the temperature in energy units, the place where to find COLVAR, HILLS and trajectory files. A reference gro or pdb file
is needed to load the trajectories. The definition of the ranges and the number of bins for the available collective variables.
Now let's start with the analysis:
1. run VMD and load metagui
2. in metagui load the metagui.input file \ref belfast-8-mg1-fig
3. In the left section of the interface "load all" the trajectories
4. Find the Microstates
In order to visualise the microstate it is convenient to align all the structures using the VMD RMSD Trajectory tool that can be found
in Extensions->Analysis.
One or more microstates can be visualised by selecting them and clicking show.
You can sort the microstate using the column name tabs, for example by clicking on size the microstates will be ordered from the larger
to the smaller. If you look at the largest one it is possible to observe that by using the four selected collective variables the backbone
conformation of the peptide is well defined while the sidechains can populate different rotameric states.
The equilibrium time in the analysis panel should be such that by averaging over the two halves of the remind of the simulation the profiles
are the same (i.e the profile averaged between Teq and Teq+(Ttot-Teq)/2 should be the same of that averaged from Teq+(Ttot-Teq)/2 and Ttot).
By clicking on COMPUTE FREE ENERGIES, the program will first generate the 1D free energy profiles from the HILLS files and then run
the WHAM analysis on the microstates. Once the analysis is done it is possible to visually check the convergence of the 1D profiles one by
one by clicking on the K bottons next to the HILLS.# files. The BLUE and the RED profiles are the two profiles just defined, while the GREEN
is the average of the two. Now it is possible for example to sort the microstates as a function of the free energy and save them by dumping
the structures for further analysis.
\anchor belfast-8-mg1-fig
\image html belfast-8-mg1.png "METAGUI Interface after loading metagui.input file"
If you look in the metagui folder you will see a lot of files, some of them can be very usefull:
metagui/MICROSTATES: is the content of the microstates list table
metagui/WHAM_RUN/VG_HILLS.#: are the opposite of the free energies calculated from the hills files
metagui/WHAM_RUN/*.gnp: are gnuplot input files to plot the VG_HILLS.# files (i.e. gnuplot -> load "convergence..")
metagui/WHAM_RUN/FES: is the result of the WHAM, for each cluster there is its free energy and the error estimate from WHAM
\verbatim
gnuplot> plot [0:40]'FES' u 2:3
\endverbatim
plots the microstate error in the free energy estimate as a function of the microstates free energy.
\Reference
This tutorial is freely inspired to the work of Biarnes et al.
More materials can be found in
1. Marinelli, F., Pietrucci, F., Laio, A. & Piana, S. A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations. PLoS Comput. Biol. 5, e1000452 (2009).
2. Biarnés, X., Pietrucci, F., Marinelli, F. & Laio, A. METAGUI. A VMD interface for analyzing metadynamics and molecular dynamics simulations. Comput. Phys. Commun. 183, 203–211 (2012).
3. Baftizadeh, F., Cossio, P., Pietrucci, F. & Laio, A. Protein folding and ligand-enzyme binding from bias-exchange metadynamics simulations. Current Physical Chemistry 2, 79–91 (2012).
