Reweighting Metadynamics

In my last post I described how to check convergence of metadynamics simulation. However, the 1D free energy surface showed in the previous post is not reweighted and corrected properly because I had both upper and lower wall in my metadynamics simulation. Hence we have to correct the free energy surface for the walls and we can do it using Plumed-2.4.1. It is important to mention that I ran my simulations with Plumed-1.3 and have to do some adjustments to get a 1ps sparsed COLVAR file whose header looks like following

 

#! FIELDS time d1 d2 d3 metad.rbias uwall.bias    
    1.0000      0.297346503      0.718053460      1.629577041    0.000000000     1.516799331 
    2.0000      0.312700063      0.676572204      1.619366407    0.500000000     1.017568588 
    3.0000      0.337963104      0.671517015      1.576508045    0.437842011     0.172387958 
    4.0000      0.304164201      0.674881995      1.545186877    0.981566668     0.022762353 
    5.0000      0.281478077      0.678171694      1.506772399    0.851512551     0.003159886 
    6.0000      0.343289763      0.649008214      1.564545155    1.156349301     0.086780988 
    7.0000      0.314589500      0.690125763      1.542209387    1.231884360     0.028832838 
    8.0000      0.297886729      0.662928343      1.551209688    1.714575529     0.034525435 
    9.0000      0.293216735      0.665717304      1.541974306    1.670215130     0.015825573 

 

Worth mentioning that the original COLVAR output from Plumed 1.3 looks like following

 

#! FIELDS time cv1 cv2 cv3 vbias vwall
    0.0000      0.331873775      0.679191828      1.648215294    0.000000000     2.416546345 
    0.1000      0.305706650      0.693865001      1.633054614    0.000000000     1.585586548 
    0.2000      0.289182395      0.736790895      1.625117898    0.000000000     1.509020209 
    0.3000      0.306852549      0.739897966      1.635408998    0.000000000     2.007535696 
    0.4000      0.289105147      0.741805196      1.605738878    0.000000000     0.980214059 
    0.5000      0.314303517      0.736100793      1.611924648    0.000000000     1.059433341 
    0.6000      0.316665620      0.693137527      1.657169104    0.000000000     3.068289280 
    0.7000      0.313074738      0.689189851      1.670769334    0.000000000     4.263953686 
    0.8000      0.281177700      0.716560900      1.678995728    0.000000000     5.230778694

I edited cv1, cv2, cv3 to d1, d2 and d3 and changed vbias to metad.rbias and vwall to uwall.bias just for naming convenience.

vwall in the original COLVAR defines the external wall potential as a combination of both upper and lower wall. 

Note: You have to delete the values corresponds to 0.000 ps from the final COLVAR file. 

The method is well described here.

See the following script which produces 1D reweighted free energy surface on CV1 which is d1 in the new COLVAR 1ps sparsed file. You will get two outputs fes.dat and fes_wall.dat and try to see how different the wall corrected (fes_wall.dat) and with wall free energy surfaces (fes.dat) look like. It will look something like following

fes1W50_wall.png

Figure 1. 1D reweighted free energy surface on CV1.

I used the following script

 

# Read COLVAR file
distance:    READ FILE=COLVAR_sparse  IGNORE_TIME VALUES=d1
metad:       READ FILE=COLVAR_sparse IGNORE_TIME VALUES=metad.rbias
uwall:       READ FILE=COLVAR_sparse IGNORE_TIME VALUES=uwall.bias

# Define weights
weights1: REWEIGHT_METAD TEMP=300 ARG=metad.rbias
weights2: REWEIGHT_BIAS TEMP=300 ARG=metad.rbias,uwall.bias

# Calculate histograms
HISTOGRAM ...
  ARG=distance
  GRID_MIN=0.0
  GRID_MAX=1.0
  GRID_BIN=100
  BANDWIDTH=0.002
  LOGWEIGHTS=weights1
  LABEL=hh1
... HISTOGRAM

HISTOGRAM ...
  ARG=distance
  GRID_MIN=0.0
  GRID_MAX=1.0
  GRID_BIN=100
  BANDWIDTH=0.002
  LOGWEIGHTS=weights2
  LABEL=hh2
... HISTOGRAM

# Print histograms to file
ff1: CONVERT_TO_FES GRID=hh1 TEMP=300
DUMPGRID GRID=ff1 FILE=fes.dat
DUMPGRID GRID=hh1 FILE=histo FMT=%24.16e
ff2: CONVERT_TO_FES GRID=hh2 TEMP=300
DUMPGRID GRID=ff2 FILE=fes_wall.dat
DUMPGRID GRID=hh2 FILE=histo_wall FMT=%24.16e

Ofcourse we can reweight the free energy surface on both CV1 and CV2 using the following script

# Read COLVAR file

distance1:   READ FILE=COLVAR_sparse  IGNORE_TIME VALUES=d1
distance2:   READ FILE=COLVAR_sparse  IGNORE_TIME VALUES=d2
metad:       READ FILE=COLVAR_sparse IGNORE_TIME VALUES=metad.rbias
uwall:       READ FILE=COLVAR_sparse IGNORE_TIME VALUES=uwall.bias

# Define weights

weights1: REWEIGHT_METAD TEMP=300 ARG=metad.rbias
weights2: REWEIGHT_BIAS TEMP=300 ARG=metad.rbias,uwall.bias

# Calculate histograms

HISTOGRAM ...
  ARG=distance1,distance2
  GRID_MIN=0.0,0.0
  GRID_MAX=1.0,1.0
  GRID_BIN=200,200
  BANDWIDTH=0.05,0.05
  LOGWEIGHTS=weights1
  LABEL=hh1
... HISTOGRAM

HISTOGRAM ...

  ARG=distance1,distance2
  GRID_MIN=0.0,0.0
  GRID_MAX=1.0,1.0
  GRID_BIN=200,200
  BANDWIDTH=0.05,0.05
  LOGWEIGHTS=weights2
  LABEL=hh2
... HISTOGRAM

# Print histograms to file

ff1: CONVERT_TO_FES GRID=hh1 TEMP=300
DUMPGRID GRID=ff1 FILE=fes.dat
DUMPGRID GRID=hh1 FILE=histo FMT=%24.16e
ff2: CONVERT_TO_FES GRID=hh2 TEMP=300
DUMPGRID GRID=ff2 FILE=fes_wall.dat
DUMPGRID GRID=hh2 FILE=histo_wall FMT=%24.16e

Here I plotted fes.dat again which is two dimensional free energy surface reweighted on CV1 and CV2. You can plot fes.dat and fes_wall.dat in a 2D plot to see what was the effect of wall. Here I plotted fes_wall.dat which is wall corrected 2D FES.

fes2d_wall_colornew

 Figure 2. 2D reweighted free energy surface on CV1 and CV2.

Note: Sorry for bad quality rendering (email me for better picture)

Reweighting algorithm described here

Pratyush Tiwary and Michele Parrinello. A time-independent free energy estimator for metadynamics. The Journal of Physical Chemistry B, 119(3):736–742, 2015. PMID: 25046020.

Manual of how to use reweighting metadynamics described here.

 

 

 

 

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

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

fes1W50_1d

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

conv

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.

recrossing

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.