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
PRINT W_STRIDE 50 HILLS RESTART W_STRIDE 1000 HEIGHT 0.5 WELLTEMPERED SIMTEMP 298 BIASFACTOR 10 DISTANCE LIST 1143 1221 SIGMA 0.05 DISTANCE LIST 1143 520 SIGMA 0.05 DISTANCE LIST 513 1153 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 ENDMETA
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 PRINT ARG=d1,d2 STRIDE=1 FILE=COLVAR
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 done
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.
