umbrella sampling and metadynamics

[qingh...@gmail.com]

Hello,

I have done both umbrella sampling and well-tempered metadynamics on the same reaction coordinates of the same system using Gromacs-4.6.7 interfaced with Plumed-2.1-hrex. The reaction coordinates was defined as the distance of the center of mass of two loops (only backbone atoms were considered)

For the US, the reaction coordinates were ranging from 5.0 to 16.5 angstrom, and the bin size was 0.5 angstrom. For each window, 1 ns equilibration and 10 ns sampling were perfromed. Here is the plumed file for US:

    UNITS LENGTH=A TIME=ps ENERGY=kcal/mol

    chainA: GROUP ATOMS=1-4129
    WHOLEMOLECULES ENTITY0=chainA

    com1: COM ATOMS=2866,2868,2887,2889,2891,2894,2896,2898,2910,2912,2914,2917
    com2: COM ATOMS=3525,3527,3530,3532,3534,3541
     dist: DISTANCE ATOMS=com1,com2
     restraint-dist: RESTRAINT ARG=dist KAPPA=40 AT=XXX
     PRINT STRIDE=25 ARG=dist,restraint-dist.bias FILE=dist.XXX.pr

200-ns WT metadynamics simulation was perfromed on the same reaction coordinates, and here is the plumed file:

UNITS LENGTH=A TIME=ps ENERGY=kcal/mol
  chainA: GROUP ATOMS=1-4129
 WHOLEMOLECULES ENTITY0=chainA

  com1: COM ATOMS=2866,2868,2887,2889,2891,2894,2896,2898,2910,2912,2914,2917
  com2: COM ATOMS=3525,3527,3530,3532,3534,3541

dist: DISTANCE ATOMS=com1,com2
metad: METAD ARG=dist PACE=500 HEIGHT=0.1 SIGMA=0.2 FILE=HILLS BIASFACTOR=3.0 TEMP=300.0
PRINT STRIDE=50 ARG=dist,metad.bias FILE=COLVAR

Then I got the two free energy profiles, I can understand the profile extract from wham:

There is an activation barrier of around 17 kcal/mol, and the transition state is at ~11 angstrom.

This plots were generated directly by "plumed sum_hills --hills HILLS". According to the tutorial, the plots are the free energy surfaces, then the global minima locates at

~9 angstrom, which contradicts with the US profile (~6.0 angstrom) though there is also local minima on the metaD profile. And the energy barrier from the close state (6.0 A) to the open state (~12 A) is <3 kcal/mol, which is much smaller than the one of US profile (~17 kcal/mol).

For the metaD profile, some of our group member thought that the plots presented here is not the free energy surface, and it need to be flipped up. But I don't agree with that.

Could you please help me understand the differences of the two profiles, especially the metaD profile? Thanks very much!

All the best,

Qinghua

[Giovanni Bussi]

Hi,

notice that if the barrier is really 12 kcal/mol your bias factor is probably too small.

In any case, I don't think your WHAM simulation is converged. As a check, you can repeat it initializing all the structures from the right-most restraint (at 16.5), assuming that now you initialized them from the left-most one (at 5). In a recent work we found a huge discrepancy, clearly showing that our WHAM results were not converged:

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-16-S9-S6

Notice that metadynamics will try to push the system back and forth, so if the CV is not sufficient to converge the simulation in the simulated time you will clearly see it from the time series (or from the bias time series, as in your case). Detecting problems in umbrella sampling is a bit more tricky.

My guess it that this single CV is not sufficient in your system.

Giovanni

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users+unsubscribe@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/51fbd5c4-fda2-4aeb-837d-570205aaa362%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[qingh...@gmail.com]

Dear Prof. Bussi,

Thanks very much for your explanation and suggestions.

As we are still at the exploration stage, I started with a small bias factor, I will try a bigger one later.

Yeah, I have read the paper you sent. Actually, I did the US simulations following the protocols of your paper. I have already done several US simulations,
backward and forward. As you said, the differences of the free energy profiles are big.

Yeah, one CV could be not enough to describe the loop motion, and I am also thinking of adding some more CVs, using replica exchange or multiple walker.
Furthermore, I am also thinking of using path metadynamics to describe the loop motion. As I am still learning metadynamics, I will try these different variations.

Thanks again!

All the best,
Qinghua

To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.

[liyib...@gmail.com]

Hell,

I am  also a newer to plumed by using MetaD. By selecting the distance  as CV I got similar energy profile as yours.I would like to know why the energy becomes Zero at lowest distance and highest distance respectively?

Best wishes

Liyi

在 2017年3月14日星期二 UTC+9上午6:29:58，qingh...@gmail.com写道：