Home » Metadynamics on Fly

Metadynamics on Fly

metadynamics_page

Metadynamics is an enhanced sampling technique which was introduced by Laio and Parrinello in 2002 and often used to construct multidimensional free energy landscape of certain system and to understand rare biological events which is now become in a fashion in biophysical chemistry. Since its inception several reviews and tutorials were published in order to provide both theoretical and hands on experience. Some of these tutorials and reviews are mentioned below:

1. Theory of Metadynamics http://people.sissa.it/~laio/Research/Res_metadynamics.php

2. Plumed, an on the fly package which can be integrated with some widely used MD codes and will be discussed here. It contains theory of Metadynamics, several inputs as well as some examples on how to execute metadynamics with an integrated MD code http://plumed.github.io/doc-v2.1/user-doc/html/index.html

3. A nice review on Metadynamics is here http://onlinelibrary.wiley.com/doi/10.1002/wcms.31/pdf , though in literature you can find many recent reviews and application. You can keep a track of all those research articles and reviews by typing “metadynamics” as search keyword in PubMed and in any other database such as ACS, ScienceDirect etc.

If you are following the Plumed Manual and tutorial which can be accessed here http://plumed.github.io/doc-v2.1/user-doc/html/index.html then we will directly start with some real life examples which we were testing. This tutorial will help us to remind our old work and will provide some beginners some ideas regarding parameters and how it works.

General advice for parameter setting in metadynamics

Sigma (if constant): use approximately half of the standard deviation of the CV in an unbiased MD simulation

Sigma (if adaptive=DIFF):  use a carefully selected number of steps for collecting fluctuations. With ADAPTIVE=DIFF we were using SIGMA=20, as it seems to be quite well in our test case but you can try different sigma. See the tutorial for the usage and example.

Sigma (if adaptive=GEOM): Unfortunately we are not lucky with ADAPTIVE=GEOM. May be we need to test a bit more. The parameter specifying the length-scale for the characterization of the landscape is suggested to be 0.005 nm in the manual.

Speed of adding bias (i.e. hill height divided by pace):

  • rough metadynamics: 2-4 kJ/mol per picosecond (very rough, not well-tempered)
  • production metadynamics: 1 kJ/mol per picosecond or below, and using well-tempered metadynamics (but sometimes a faster speed is required even in production)

Bias factor for well-tempered metadynamics: should be chosen to explore the desired region of the energy landscape, based on preliminary calculations or sometimes educated guesses. Rule of thumb: A biasfactor of 10 makes the hills become significantly smaller at a bias level of ~10kT, so in practice it explores well up to ~15 kT (very rough estimates here).

Pace for adding hills: 0.2 to 1 picosecond   (rough runs use smaller), i.e. 100-500 steps.

Hill height: should be adjusted to give the desired speed (see above). The greater the height it will explore regions quickly but its prone to error. But for rough calculation you can try greater height such as 2 kJ/mol, HEIGHT=2.0

In our test case we are using the following CVs. Note: These are atom numbers not residue numbers.

Definition of CVs:

dist1: distance between 424:N:Lys28 and 488:O:Arg31

dist2: distance between 445:O:Lys28 and 465:N:Arg31

dist3: distance between 391:N:Val26 and 528:O:Val33

dist4: distance between 406:O:Val26 and 513:N:Val33

dist5: distance between 367:N:Ala24 and 563:O:Phe35

dist6: distance between 376:O:Ala24 and 544:N:Phe35

#example script for a 1D-Metadynamics with dist1 as a CV

DISTANCE ATOMS=424,488 LABEL=d1
 METAD ARG=d1 SIGMA=0.05 HEIGHT=0.5 PACE=100 LABEL=metad
 PRINT ARG=d1,metad.bias STRIDE=100 FILE=COLVAR

You will get files named as COLVAR and HILLS. Open the COLVAR and HILL files to see the fluctuation in d1 and d2 during Metadynamics. You can increase the HEIGHT or PACE to increase the pace of metadynamics so that it can quickly explore certain regions.

If you plot the COLVAR file using GNUPLOT you will see something like this

image1

This will give you an idea about the fluctuation of hydrogen bond between the selected atoms during Metadynamics simulation. If we construct the free energy surface for this using plumed sum_hills then we can see something like the following (run the following command)

plumed sum_hills --hills HILLS --outfile fes --stride 500 (#you can change the stride to 100, plot and see the difference)

we can plot all fes files using GNUPLOT. See the free energy landscape

image2

Monitor the growing of the bias during simulation by plotting metad.bias against time (inside COLVAR file)

image3

Oh!! Yes here is the script on how to run Metadynamics on a cluster which already has Gromacs and Plumed patched with Gromacs

#!/bin/bash

#SBATCH -t 48:00:00
 #SBATCH -J test2
 #SBATCH -A SNIC2015-1-135
 #SBATCH -N 1
 #SBATCH -n 18

# Load the module for gromacs.
 module add gromacs/5.0.4
 if [ -e state.cpt  ]; then
 extra='-cpi state.cpt'
 fi
 srun  mdrun_plumed_mpi_d -ntomp 1 -s gpwtest.tpr $extra -noappend -maxh 48 -plumed plumed.dat

Note: the plumed.dat files contains the inputs which we discussed above

Now we are trying the same 1D-Metadynamics wit dist2 as a CV (collective variable)

DISTANCE ATOMS=445,465 LABEL=d2
 METAD ARG=d2, SIGMA=0.05 HEIGHT=0.5 PACE=100 LABEL=metad
 PRINT ARG=d2,metad.bias STRIDE=100 FILE=COLVAR

Now see the change in d2 vs time (inside COLVAR file)

image5

You can monitor change in metad.bias as we described earlier. But lets see the free energy landscape

image4In this point we will try ADAPTIVE=DIFF (Diffusion) keyword with d1 in case of 1D- metadynamics. Lets see the input

DISTANCE ATOMS=424,488 LABEL=d1
 METAD ARG=d1 SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF
 PRINT ARG=d1,metad.bias STRIDE=100 FILE=COLVAR

Now lets plot the d1 in the COLVAR file

image6Now we will see how the corresponding FES looks like (with stride 500)

image7Well lets try the metadynamics on d1 CV but with adaptive GEOM (geometric)

DISTANCE ATOMS=424,488 LABEL=d1
 METAD ARG=d1 SIGMA=0.005 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=GEOM
 PRINT ARG=d1,metad.bias STRIDE=100 FILE=COLVAR

Ok, its time to see the variation in d1 as well as the free energy landscape (FES with stride 500)

#time vs d1 (inside the COLVAR file)

image8Lets see the FES file

image9

I hope you the get the key of setting up the input. Now you can bias all 6 CVs using ADAPTIVE=DIFF by using the following input

DISTANCE ATOMS=424,488 LABEL=d1
 DISTANCE ATOMS=445,465 LABEL=d2
 DISTANCE ATOMS=367,563 LABEL=d3
 DISTANCE ATOMS=376,544 LABEL=d4
 DISTANCE ATOMS=391,528 LABEL=d5
 DISTANCE ATOMS=406,513 LABEL=d6

METAD ARG=d1,d2,d3,d4,d5,d6 SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF
 PRINT ARG=d1,d2,d3,d4,d5,d6,metad.bias STRIDE=100  FILE=COLVAR

Now we will shift to 2D-metadynamics by taking in account both d1 and d2 as a CV. An advantage of 2D-metadynamics is it will give a more realistic understanding of how some rare events occurred. But that discussion we will continue later. Now lets see how the input looks like. We will start with ADAPTIVE=DIFF on 2D-metadynamics. Lets start

this is how plumed.dat looks like

DISTANCE ATOMS=424,488 LABEL=d1
 DISTANCE ATOMS=445,465 LABEL=d2

METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF

PRINT ARG=d1,d2,metad.bias STRIDE=100  FILE=COLVAR

Now type head COLVAR to see the title of the each columns created inside COLVAR file and we can start plotting

Time vs d1

time_d1

Time vs d2

time_d2

Deposition of bias

bias_depos

Sigma d1_d1 (inside the COLVAR file)

sigma_d1_d1

Sigma d2_d2

sigma_d2_d2

Sigma d1_d2

sigma_d1_d2

Now we are integrating out the FES at a particular temperature (see the sum_hills page here for the keywords http://plumed.github.io/doc-v2.0/user-doc/html/sum_hills.html and also http://plumed.github.io/doc-v2.1/user-doc/html/kt.html )

#projection d1 -kt 2.494

plumed sum_hills --hills HILLS --outfile fes --stride 500 --bin 500,500 --idw d1 --kt 2.494

projection_d1You can see that taking in account d1 CV is started to visiting another conformational spce as evident from the another basin in the FES plot.

# projection d2 -kt 2.494

plumed sum_hills --hills HILLS --outfile fesint --stride 500 --bin 500,500 --idw d2 --kt 2.494

projection_d2

Lets plot a 2D free energy landscape by taking in account Fes4.dat and Fes10.dat just as an example using GNUPLOT. You can try with the folowing set of commands in GNUPLOT to carry out this operation on Fes1-10.dat files.

~/abpfs/Gpw/Gpwtest/Meta2d_diff $ fg
 gnuplot

gnuplot>
 gnuplot> sp 'fes2.dat' u 1:2:3
 gnuplot> set multiplot layout 4,2
 multiplot> sp 'fes2.dat' u 1:2:3
 multiplot> sp 'fes3.dat' u 1:2:3
 multiplot> sp 'fes4.dat' u 1:2:3
 multiplot> sp 'fes5.dat' u 1:2:3
 multiplot> sp 'fes6.dat' u 1:2:3
 multiplot> sp 'fes7.dat' u 1:2:3
 multiplot> sp 'fes8.dat' u 1:2:3
 multiplot> sp 'fes9.dat' u 1:2:3
 multiplot> unset multiplot
 gnuplot>

Fes4.dat looks something like this

fes4_pcaWhile the Fes10.dat looks like the following

fes10_pca

 Now we will shift our view to switching functions and path CV. Lets see the original tutorials here

http://plumed.github.io/doc-v2.0/user-doc/html/switchingfunction.html

Description of contact map is here http://plumed.github.io/doc-v2.1/user-doc/html/_c_o_n_t_a_c_t_m_a_p.html

How to define path CV on the space of contact map http://plumed.github.io/doc-v2.1/user-doc/html/_f_u_n_c_p_a_t_h_m_s_d.html

Firstly lets make histogram of all 6 CVs from d1 to d6 during simulation time. It will look something like the following

histogram_distanceThen we will plot histogram and distribution in each plot for hist2, hist3…hist6 using GNUPLOT using a command something like the following

p 'hist2.dat' u ($1*0.1):(0.003*$2) w l, (1-((x-d0)/r0)**n)/(1-((x-d0)/r0)**m), d0=0.29, n=6,m=12,r0=0.05
 set xrange [0.25:0.65]

Lets see hist2 how it looks like by using the above mentioned script

histo2

Now this plot will show you how to adjust the RATIONAL_R0 and D_0. If we project the peak of the hist2 on the x axis it will come something around 0.28 which is the value of D_0 whereas the end part of the green line will end somewhere near 0.70 (though it seems to be not in the plot range but you can adjust it by setting the xrange [0.25:0.80]) hence you can set the RATIONAL_R0 to 0.7. We are using NN=3 MM=8 as default (see the links provided above for more details).

Now lets plot hist3 and you can see D_0 is .3 whereas the RATIONAL_R0 is similar to that of previous one that means 0.007.

histo3You can set RATIONAL_R0 and D_0 for all 6 histograms and now we can vi plumed.dat to specify path CV based on PATHMSD (pathRMSD) for all 6 CVs using contact map. See the script below.

DISTANCE ATOMS=424,488 LABEL=d1
 DISTANCE ATOMS=445,465 LABEL=d2
 DISTANCE ATOMS=391,528 LABEL=d3
 DISTANCE ATOMS=406,513 LABEL=d4
 DISTANCE ATOMS=367,563 LABEL=d5
 DISTANCE ATOMS=376,544 LABEL=d6

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.0
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.0
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.0
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.0
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.0
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.0
 LABEL=cm1
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.1
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.1
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.1
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.1
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.1
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.1
 LABEL=cm2
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.2
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.2
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.2
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.2
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.2
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.2
 LABEL=cm3
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.3
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.3
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.3
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.3
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.3
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.3
 LABEL=cm4
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.4
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.4
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.4
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.4
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.4
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.4
 LABEL=cm5
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.5
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.5
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.5
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.5
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.5
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.5
 LABEL=cm6
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.6
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.6
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.6
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.6
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.6
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.6
 LABEL=cm7
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.7
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.7
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.7
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.7
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.7
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.7
 LABEL=cm8
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.8
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.8
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.8
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.8
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.8
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.8
 LABEL=cm9
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=0.9
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=0.9
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=0.9
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=0.9
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=0.9
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=0.9
 LABEL=cm10
 CMDIST
 ... CONTACTMAP

CONTACTMAP ...
 ATOMS1=424,488 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8} REFERENCE1=1.0
 ATOMS2=445,465 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=3 MM=8} REFERENCE2=1.0
 ATOMS3=391,528 SWITCH3={RATIONAL R_0=0.08 D_0=0.29 NN=3 MM=8} REFERENCE3=1.0
 ATOMS4=406,513 SWITCH4={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE4=1.0
 ATOMS5=367,563 SWITCH5={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE5=1.0
 ATOMS6=376,544 SWITCH6={RATIONAL R_0=0.07 D_0=0.29 NN=3 MM=8} REFERENCE6=1.0
 LABEL=cm11
 CMDIST
 ... CONTACTMAP

p1: FUNCPATHMSD ARG=cm1,cm2,cm3,cm4,cm5,cm6,cm7,cm8,cm9,cm10,cm11 LAMBDA=39
 METAD ARG=p1.s,p1.z SIGMA=20 HEIGHT=0.5 PACE=1000 ADAPTIVE=DIFF LABEL=metad
 PRINT ARG=d1,d2,d3,d4,d5,d6,cm1,cm2,cm3,cm4,cm5,cm6,cm7,cm8,cm9,cm10,cm11,p1.s,p1.z,metad.bias STRIDE=1 FILE=COLVAR FMT=%8.4f

We will see the results after a bit but for now lets see how a general input taking in account the co-ordination or contact map looks like (based on dist1 and dist2 as a CV)

#corodination
 DISTANCE ATOMS=424,488 LABEL=d1
 DISTANCE ATOMS=445,465 LABEL=d2
 COORDINATION GROUPA=424 GROUPB=488 LABEL=c1 SWITCH={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8}
 COORDINATION GROUPA=445 GROUPB=465 LABEL=c2 SWITCH={RATIONAL R_0=0.07 D_0=0.31 NN=6 MM=8}
 METAD ARG=d1,d2,c1,c2 SIGMA=0.35,0.35 HEIGHT=0.5 PACE=1000 LABEL=metad
 PRINT ARG=d1,d2,c1,c2,metad.bias STRIDE=100 FILE=COLVAR

#contactmap
 DISTANCE ATOMS=424,488 LABEL=d1
 DISTANCE ATOMS=445,465 LABEL=d2
 ATOMS1=424,488 LABEL=cm1 SWITCH1={RATIONAL R_0=0.07 D_0=0.28 NN=3 MM=8}
 ATOMS2=445,465 LABEL=cm2 SWITCH2={RATIONAL R_0=0.07 D_0=0.31 NN=6 MM=8}
 METAD ARG=cm1,cm2 SIGMA=0.35,0.35 HEIGHT=0.5 PACE=1000 LABEL=metad
 PRINT ARG=d1,d2,cm1,cm2,metad.bias STRIDE=100 FILE=COLVAR

Now we will see how centre of mass command works with Plumed

The COM commands creates a centre of mass for a group of atoms and stored them in a virtual atom which can be used to computed distances from centre of mass and subsequently metadynamics

#centre of mass among atoms 1-7
 COM ATOMS=1-7         LABEL=c1
 #centre of mass among atoms 15 and 20
 COM ATOMS=15,20       LABEL=c2
 #distance between these two centre of masses
 DISTANCE ATOMS=c1,c2  LABEL=d1
 #metadynamics with ADAPTIVE=DIFF
 METAD ARG=d1 SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF
 PRINT ARG=d1,metad.bias STRIDE=100  FILE=COLVAR

Lets see how the ANGLE keyword works for (make sure you specify atom numbers not the residue number)

ANGLE ATOMS=50,25,124 LABEL=a
 ANGLE ATOMS=50,25,124,149 LABEL=b
 METAD ARG=a,b SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF
 PRINT ARG=a,b,metad.bias STRIDE=100  FILE=COLVAR

 Now we will see how the script of TORSION atom will look like

TORSION ATOMS=1,2,3,4 LABEL=t
 METAD ARG=t SIGMA=20 HEIGHT=0.5 PACE=100 LABEL=metad ADAPTIVE=DIFF
 PRINT ARG=t,metad.bias STRIDE=100  FILE=COLVAR

To monitor the perform Metadynamics with ALPHABETA keyword http://plumed.github.io/doc-v2.1/user-doc/html/_a_l_p_h_a_b_e_t_a.html . here we combine ALPHABETA keyword with path CV with FUNCPATHMSD.

#ATOMS1 the contributing atoms to create a torsional angle. REFERENCE1 in first part is the reference (starting) value for the torsion angle (in radian) and the second part is the torsion value of the fluctuation. For example in case of ATOMS1 the REFERENCE1 in the first part is the initial value (-2.43 radian) for torsion angle and the later value (-0.69 radian) is for the fluctuated value.

ALPHABETA ...
 ATOMS1=1062,1064,1066,1074 REFERENCE1=-2.43
 ATOMS2=1074,1076,1078,1081 REFERENCE2=1.22
 ATOMS3=1081,1083,1085,1095 REFERENCE3=-0.69
 ATOMS4=1179,1181,1183,1200 REFERENCE4=-2.44
 ATOMS5=1200,1202,1204,1216 REFERENCE5=-1.05
 ATOMS6=1216,1218,1220,1227 REFERENCE6=-1.40
 ATOMS7=1227,1229,1231,1234 REFERENCE7=2.44
 ATOMS8=1234,1236,1238,1248 REFERENCE8=-2.61
 ATOMS9=1248,1250,1252,1264 REFERENCE9=-2.44
 ATOMS10=1275,1277,1279,1282 REFERENCE10=-3.14
 ATOMS11=1282,1284,1286,1302 REFERENCE11=-2.44
 LABEL=ab1
 ... ALPHABETA

ALPHABETA ...
 ATOMS1=1062,1064,1066,1074 REFERENCE1=-0.69
 ATOMS2=1074,1076,1078,1081 REFERENCE2=3.14
 ATOMS3=1081,1083,1085,1095 REFERENCE3=3.14
 ATOMS4=1179,1181,1183,1200 REFERENCE4=-1.05
 ATOMS5=1200,1202,1204,1216 REFERENCE5=1.05
 ATOMS6=1216,1218,1220,1227 REFERENCE6=-2.27
 ATOMS7=1227,1229,1231,1234 REFERENCE7=-3.14
 ATOMS8=1234,1236,1238,1248 REFERENCE8=-1.39
 ATOMS9=1248,1250,1252,1264 REFERENCE9=-1.22
 ATOMS10=1275,1277,1279,1282 REFERENCE10=3.14
 ATOMS11=1282,1284,1286,1302 REFERENCE11=-1.57
 LABEL=ab2
 ... ALPHABETA

p1: FUNCPATHMSD ARG=ab1,ab2 LAMBDA=0.351
 METAD ARG=p1.s,p1.z SIGMA=20 HEIGHT=0.5 PACE=1000 ADAPTIVE=DIFF LABEL=metad
 PRINT ARG=ab1,ab2,p1.s,p1.z,metad.bias STRIDE=100 FILE=COLVAR FMT=%8.4f

From manual: λ should be chosen so as to have λ*d(Xi, Xi+1) ≃ 2.3 on average.

For the example above: d(X1, X2) = d(formed, broken) = (1-0)^2 + (1-0)^2 = 2, so one should set LAMBDA=1.15

Use of COMBINE keyword in order to construct a 1D MetaD path with torsion atoms and ALPHABETA keyword

ALPHABETA ...
 ATOMS1=1062,1064,1066,1074 REFERENCE1=-2.43
 ATOMS2=1074,1076,1078,1081 REFERENCE2=1.22
 ATOMS3=1081,1083,1085,1095 REFERENCE3=-0.69
 ATOMS4=1179,1181,1183,1200 REFERENCE4=-2.44
 ATOMS5=1200,1202,1204,1216 REFERENCE5=-1.05
 ATOMS6=1216,1218,1220,1227 REFERENCE6=-1.40
 ATOMS7=1227,1229,1231,1234 REFERENCE7=2.44
 ATOMS8=1234,1236,1238,1248 REFERENCE8=-2.61
 ATOMS9=1248,1250,1252,1264 REFERENCE9=-2.44
 ATOMS10=1282,1284,1286,1302 REFERENCE10=-2.44
 LABEL=s1
 ... ALPHABETA

ALPHABETA ...
 ATOMS1=1062,1064,1066,1074 REFERENCE1=-1.2
 ATOMS2=1074,1076,1078,1081 REFERENCE2=3.14
 ATOMS3=1081,1083,1085,1095 REFERENCE3=-2.7
 ATOMS4=1179,1181,1183,1200 REFERENCE4=-1.05
 ATOMS5=1200,1202,1204,1216 REFERENCE5=1.05
 ATOMS6=1216,1218,1220,1227 REFERENCE6=-2.6
 ATOMS7=1227,1229,1231,1234 REFERENCE7=2.8
 ATOMS8=1234,1236,1238,1248 REFERENCE8=-1.39
 ATOMS9=1248,1250,1252,1264 REFERENCE9=-1.5
 ATOMS10=1282,1284,1286,1302 REFERENCE10=-1.57
 LABEL=s2
 ... ALPHABETA

t1: TORSION ATOMS=1062,1064,1066,1074
 t2: TORSION ATOMS=1074,1076,1078,1081
 t3: TORSION ATOMS=1081,1083,1085,1095
 t4: TORSION ATOMS=1179,1181,1183,1200
 t5: TORSION ATOMS=1200,1202,1204,1216
 t6: TORSION ATOMS=1216,1218,1220,1227
 t7: TORSION ATOMS=1227,1229,1231,1234
 t8: TORSION ATOMS=1234,1236,1238,1248
 t9: TORSION ATOMS=1248,1250,1252,1264
 t10: TORSION ATOMS=1282,1284,1286,1302

d1: COMBINE ARG=s1,s2 POWERS=1,0 COEFFICIENTS=-1,10 PERIODIC=NO
 d2: COMBINE ARG=s2,s1 POWERS=1,0 COEFFICIENTS=-1,10 PERIODIC=NO
 f: COMBINE ARG=d1,d2 POWERS=1,1 COEFFICIENTS=1,-1 PERIODIC=NO

METAD ARG=f SIGMA=1000 HEIGHT=0.5 PACE=500 ADAPTIVE=DIFF LABEL=metad
 PRINT ARG=s1,s2,d1,d2,f,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,metad.bias STRIDE=100 FILE=COLVAR FMT=%8.4f

Metadynamics with Principal Component Analysis

Here is a PCA based metadynamics simulation which I tried.

PRINT W_STRIDE 50
HILLS RESTART W_STRIDE 15000 HEIGHT 2.0
WELLTEMPERED SIMTEMP 298 BIASFACTOR 10
PCA FRAME ref.dat EIGENVEC ev1.dat DIFF SIGMA 0.05
UWALL CV 1 LIMIT 1.0 KAPPA 5000.0 EXP 4.0
LWALL CV 1 LIMIT -3.0 KAPPA 5000.0 EXP 4.0

ALIGN_ATOMS LIST <prot>
prot->
5 24 31 42 53 67 79 93 112 127 146 162 174 194 211 225 244 261 281 302 309 321 331 346 362 369 381 395 412 437 443 463 477 497 516 535 547 561 568 579 589 603 622 646 670 676 687 703 725 735 749 763 773 780 790 809 823 845 863 882 903 915 926 937 959 970 994 1008 1029 1044 1066 1078 1085 1099 1121 1137 1152 1169 1183 1204 1220 1231 1238 1252 1268 1279 1286 1306 1326 1337 1359 1371 1390 1406 1420 1436 1443 1457 1476 1487 1514 1520 1541 1563 1583 1602 1617 1633 1652 1664 1678 1692 1699 1719 1742 1748 1762 1783 1797 1807 1818 1832 1852 1864 1871 1890 1909 1916 1935 1942 1966 1988 2000 2019 2030 2049 2056 2067 2083 2103 2109 2128 2144 2160 2175 2194 2216 2230 2247 2261 2283 2302 2317 2331 2341 2360 2380 2394 2414 2435 2462 2468 2484 2502 2514 2536 2554 2568 2575 2595 2614 2628 2647 2654 2661 2680 2695 2710 2734 2754 2775 2790 2805 2811 2830 2844 2865 2880 2902 2921 2935 2953 2965 2984 3005 3029 3046 3065 3079 3098 3110 3120 3138 3154 3161 3175 3194 3211 3230 3245 3267 3277 3291 3302 3321 3337 3350 3361 3368 3382 3393 3403 3422 3436 3460 3466 3480 3492 3512 3531 3545 3567 3584 3603 3620 3634 3653 3665 3681 3700 3722 3746 3752 3772 3799 3805 3825 3846 3862 3876 3895 3905 3919 3933 3944 3966 3993 3999 4013 4033 4048 4068 4082 4093 4108 4122 4129 4151 4172 4186 4205 4228 4234 4249 4270 4291 4310 4327 4345 4364 4379 4391 4407 4422 4428 4435 4454 4464 4481 4500 4514 4533 4552 4559 4578 4590 4618 4624 4648 4654 4668 4688 4707 4726 4733 4754 4760 4780 4797 4821 4843 4864 4884 4898 4914 4934 4946 4967 4979 4993 5011 5022 5038 5045 5064 5074 5093 5103 5125 5147 5161  
prot<-

ENDMETA

 Steered MD:

Steered MD is often used to drag a system from its initial configuration to its final configuration. I used this technique once to drag a ligand from the active site of the protein in order to study the extend of flap opening. In my opinion, the steer MD is not an ideal way to study ligand binding and unbinding as it is highly biased and very hard to converge. Here is a script which I used to study the ligand unbinding from the protein cavity. The velocity is set at 0.0003 per 1000 steps (2 ps, as in gromacs) by this velocity it will take 20ns to reach 3 nm in terms of distance CV which I defined.

DISTANCE LIST <f1> <g1> NOPBC

#centre of mass of catalytic aspartyl residues
f1->

535 3337

f1<-

#centre of mass of the ligand
g1->

5179 5180 5181 5182 5183 5184 5185 5186 5187 5188 5189 5190 5191 5192 5193 5194 5195 5196 5197 5198 5199 5200 5201 5202 5203 5204 5205 5206 5207 5208 5209 5210 5211 5212 5213 5214 5215 5216 5217 5218 5219 5220 5221 5222 5223 5224 5225 5226 5227 5228 5229 5230 5231 5232 5233 5234 5235 5236 5237 5238 5239 5240 5241 5242 5243 5244 5245 5246 5247 5248 5249 5250 5251 5252 5253 5254 5255 5256 5257 5258 5259 5260 5261 5262 5263 5264 5265 5266 5267 5268 5269 5270 5271 5272 5273 5274 5275 5276 5277

g1<-


STEER CV 1 TO 3.0 VEL 0.0003 KAPPA 500.0

PRINT W_STRIDE 50

ALIGN_ATOMS LIST <prot>
prot->
5 24 31 42 53 67 79 93 112 127 146 162 174 194 211 225 244 261 281 302 309 321 331 346 362 369 381 395 412 437 443 463 477 497 516 535 547 561 568 579 589 603 622 646 670 676 687 703 725 735 749 763 773 780 790 809 823 845 863 882 903 915 926 937 959 970 994 1008 1029 1044 1066 1078 1085 1099 1121 1137 1152 1169 1183 1204 1220 1231 1238 1252 1268 1279 1286 1306 1326 1337 1359 1371 1390 1406 1420 1436 1443 1457 1476 1487 1514 1520 1541 1563 1583 1602 1617 1633 1652 1664 1678 1692 1699 1719 1742 1748 1762 1783 1797 1807 1818 1832 1852 1864 1871 1890 1909 1916 1935 1942 1966 1988 2000 2019 2030 2049 2056 2067 2083 2103 2109 2128 2144 2160 2175 2194 2216 2230 2247 2261 2283 2302 2317 2331 2341 2360 2380 2394 2414 2435 2462 2468 2484 2502 2514 2536 2554 2568 2575 2595 2614 2628 2647 2654 2661 2680 2695 2710 2734 2754 2775 2790 2805 2811 2830 2844 2865 2880 2902 2921 2935 2953 2965 2984 3005 3029 3046 3065 3079 3098 3110 3120 3138 3154 3161 3175 3194 3211 3230 3245 3267 3277 3291 3302 3321 3337 3350 3361 3368 3382 3393 3403 3422 3436 3460 3466 3480 3492 3512 3531 3545 3567 3584 3603 3620 3634 3653 3665 3681 3700 3722 3746 3752 3772 3799 3805 3825 3846 3862 3876 3895 3905 3919 3933 3944 3966 3993 3999 4013 4033 4048 4068 4082 4093 4108 4122 4129 4151 4172 4186 4205 4228 4234 4249 4270 4291 4310 4327 4345 4364 4379 4391 4407 4422 4428 4435 4454 4464 4481 4500 4514 4533 4552 4559 4578 4590 4618 4624 4648 4654 4668 4688 4707 4726 4733 4754 4760 4780 4797 4821 4843 4864 4884 4898 4914 4934 4946 4967 4979 4993 5011 5022 5038 5045 5064 5074 5093 5103 5125 5147 5161 5179 5180 5181 5182 5183 5184 5185 5186 5187 5188 5189 5190 5191 5192 5193 5194 5195 5196 5197 5198 5199 5200 5201 5202 5203 5204 5205 5206 5207 5208 5209 5210 5211 5212 5213 5214 5215 5216 5217 5218 5219 5220 5221 5222 5223 5224 5225 5226 5227 5228 5229 5230 5231 5232 5233 5234 5235 5236 5237 5238 5239 5240 5241 5242 5243 5244 5245 5246 5247 5248 5249 5250 5251 5252 5253 5254 5255 5256 5257 5258 5259 5260 5261 5262 5263 5264 5265 5266 5267 5268 5269 5270 5271 5272 5273 5274 5275 5276 5277
prot<-
ENDMETA


Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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 )

Google+ photo

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

Connecting to %s

%d bloggers like this: