Convergence of Metadynamics simulation

Convergence of metadynamics or any enhanced sampling calculations plays an important role in deciding how much you can trust the outcomes from metadynamics simulations. here is a step by step tutorial and codes on how to check convergence in metadynamics simulations.

I ran two metadynamics simulations with DISTANCE CVs using Plumed 1.3 patched with Gromacs 4.6.2. As an output I got HILLS file whose first part looks like following

 2.000 0.312700063 0.676572204 0.050000001 0.050000001 0.555555556 10.000 
 4.000 0.304164201 0.674881995 0.050000001 0.050000001 0.543420527 10.000 
 6.000 0.343289763 0.649008214 0.050000001 0.050000001 0.539083481 10.000 
 8.000 0.297886729 0.662928343 0.050000001 0.050000001 0.525467237 10.000 
 10.000 0.291971087 0.575445116 0.050000001 0.050000001 0.547548082 10.000 
 12.000 0.343111038 0.602572322 0.050000001 0.050000001 0.531554010 10.000 
 14.000 0.285310268 0.667471945 0.050000001 0.050000001 0.513465246 10.000 
 16.000 0.319890589 0.655008733 0.050000001 0.050000001 0.497615549 10.000 
 18.000 0.389135718 0.529440880 0.050000001 0.050000001 0.550858286 10.000 
 20.000 0.286446780 0.555136204 0.050000001 0.050000001 0.533876313 10.000

DISTANCE CVs used in the calculation

DISTANCE LIST 1143 1221 SIGMA 0.05 
DISTANCE LIST 1143 520 SIGMA 0.05 

UWALL CV 1 LIMIT 0.45 KAPPA 5000.0 EXP 4.0
LWALL CV 1 LIMIT 0.25 KAPPA 5000.0 EXP 4.0

UWALL CV 2 LIMIT 0.65 KAPPA 5000.0 EXP 4.0
LWALL CV 2 LIMIT 0.25 KAPPA 5000.0 EXP 4.0

UWALL CV 3 LIMIT 1.50 KAPPA 5000.0 EXP 4.0


I also post processed all my trajectories to generate a file named trajfit.xtc using the following script

#Combining the .xtc files
/home/per/lun/gromacs-4.6.2/build/src/tools/trjcat_mpi -f traj.part00*.xtc -o trajall.xtc
#PBC fixed, echo 1 signifies only protein as output, 1was the indx number
echo 1 | /home/per/lun/gromacs-4.6.2/build/src/tools/trjconv_mpi -f trajall.xtc -pbc nojump -o trajnoj.xtc -s md.tpr 
# Aligning on CA and protein as output
echo 3 1 | /home/per/lun/gromacs-4.6.2/build/src/tools/trjconv_mpi -s md.tpr -f trajnoj.xtc -fit rot+trans -o trajfit.xtc

I used plumed driver on trajfit.xtc file to generate a COLVAR file with 1 ps time step. The plumed.dat file used with the driver looks like following

d1: DISTANCE ATOMS=1143,1221
d2: DISTANCE ATOMS=1143,520


I issued the driver command on the trajfit.xtc file

plumed driver --mf_xtc trajfit.xtc --plumed plumed.dat

I got a COLVAR file as output and removed the header as well as removed the value corresponds to 0.00 ps. The next step was to match the end part of COLVAR and HILLS files.  

Now I dumped fes files corresponds to each ns using the following command

~/plumed-1.3/utilities/WT_reweight/reweight -colvar COLVAR -hills HILLS -ncv 2 -nvar 2 -stride 2 -fes 1 -temp 300 -welltemp -timeout 1000 -kjoule

So the one dimensional reweighted free energy surface looks like following


Figure 1. 1D FES reweighted on CV1.

The FES above shows two basins, basin 1: averaged over 0.27-0.33 nm and basin 2: averaged over 0.62-0.74 nm.

We created a script which is named calcfediff

#Usage: calcfediff fes.dat a b c d
#Calculates difference between the two regions a-b and c-d.

#This is the exponentially averaged energy
grep . $1 | awk 'BEGIN{s1='$2'; s2='$3'; s3='$4'; s4='$5'; kt=2.5}{e=$2; c=exp(-e/kt); if ($1>s1 && $1<=s2) {q1+=c} if($1>s3 && $1<s4) {q2+=c} }END{print -kt*log(q1/q2)}'
#grep . $1 | awk 'BEGIN{s1='$2'; s2='$3'; s3='$4'; s4='$5'; kt=2.5}{e=$2; c=exp(-e/kt); if ($1>s1 && $1<=s2) {q1+=c} if($1>s3 && $1<s4) {q2+=c} }END{print -kt*log(q1/q2), q1, q2}'

We incorporated calcfediff inside a .sh file by following away

shopt -s nullglob
rm conv2.dat
for a in fes.dat.? fes.dat.?? fes.dat.???; do
echo -n ${a#fes.dat.} " " >>conv2.dat
calcfediff $a 0.27 0.33 0.62 0.74 >>conv2.dat

This gives a plot of Free Energy as a function of time which looks something like following


Figure 2. Free Energy as a function of time

We can say from Figure 2 that free energy value is converged to a certain value after ~300 ns. 

Another important part to see if the simulation is converged is to see if there is enough recrossings in CVs. Figure 3 shows the fluctuation of CVs during our metadynamics simulations. As our CVs are hydrogen bond distances we can definitely say from the fluctuation that we have enough alternative fluctuations between CV1 and CV2.


Figure 3. Fluctuation of CV1 and CV2 in our metadynamics simulation.


However the free energy profiles showed in Figure 1 is not reweighted and corrected for the wall bias. We need to reweight it properly and apply the correction for the applied wall potential. Keep an eye on the upcoming post where I will explain how to reweight and correct for wall potential.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s