Finding free energy differences between conformation clusters on plumed

[billy.willi...@gmail.com]

Hi All,

• I ran a well-tempered metadynamics simulation a cyclic peptide [say cyclo(ABCDEF) ].
• I biased it in terms of three internal distances, d7, d8 and d9, where d7 B->E, d8 was A->D and d9 was C->F from the example above.
• This simulation ran to convergence until the Gaussian height reached zero and the metadynamics bias did not change, and neither did the FES for 40ns.  The simulation ran for 100ns all up.
• I have an FES for d7, d8 and d9.
• I have clustred all conformers across the entire metadynamics simulation into 6 distinct groups.  This clustering was based on the c-alpha coordinates and how they varied across the metaD simulation.
• I found that cluster 6 was the lowest energy conformation based on measurements of d7, d8 and d9 in those clusters, in conjunction with the FES of d7, d8 and d9.
  

My question is, how do I get the difference in free energy between my 6 clusters?  If I had a 4 dimensional free energy surface ( dG, d7, d8 and d9 ) then I would be able to locate each of my clusters on that surface.  Can plumed do this?

Best regards,

Billy

[Alejandro Gil-Ley]

Hi

To find the reweighted population of each cluster and then calculate the free-energy difference, you can follow the procedure of tutorial http://plumed.github.io/doc-v2.2/user-doc/html/belfast-4.html, section "Reweighting the results". Just, instead of using a plumed file similar to the one used in the example (for umbrella sampling) you can modified your plumed.dat like this:

from

phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
a: METAD BIASFACTOR=15 ARG=phi,psi PACE=500 HEIGHT=0.4 SIGMA=0.25,0.25 .... FILE=HILLS
PRINT STRIDE=100 ARG=phi,psi FILE=COLVAR

to

phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
a: METAD RESTART=YES BIASFACTOR=15 ARG=phi,psi PACE=100000000 HEIGHT=0 SIGMA=0.25,0.25 .... FILE=HILLS
PRINT STRIDE=100 ARG=a.bias FILE=BIAS

once you have the bias potential on the BIAS file it would possible to reweight the simulation with a simple command line like the one used in the tutorial (e.g. your final file could have for each time a column with a binary variable that indicates if the frame belong to cluster #1, and second column to cluster #2, ..., final column would be the bias potential):

> grep -v \# CLUSTER_BIAS | awk '{w=exp($NF/2.5); if($2==1)a+=w; else b+=w;}END{print a/(a+b)}'

Best

A

[Alejandro Gil-Ley]

I forgot, you have to use plumed driver to print the BIAS file:

plumed driver --plumed new_plumed.dat --mf_xtc traj.xtc 

(if a distance is calculated between the center of masses of two group of atoms, you will need to created a mcfile (DUMPMASSCHARGE) with the atom masses and complete the previous command line with --mc mcfile)

[Andrea Spitaleri]

Hi,

have look to our work on very similar subject:

http://onlinelibrary.wiley.com/doi/10.1002/chem.201501196/abstract

HTH

and

--
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...@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/ac58fb67-78a4-49e6-8b7d-52a41362c32a%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

"finche' andate in bici e fischia il vento
non sara' mai una brutta giornata"
----------------------------------------------------------------
http://sites.google.com/site/andreaspitaleri/
----------------------------------------------------------------
http://it.linkedin.com/in/andreaspitaleri
----------------------------------------------------------------

[Pratyush Tiwary]

Dear Billy

Firstly, I would request you to dispel the notion that "simulation ran to convergence until the Gaussian height reached zero and the metadynamics bias did not change". Gaussian height reaching zero has absolutely nothing to do with metadynamics converging - neither of these two implies the other. Your FES might be converged, but bias still non-zero, or your bias could be close to 0, but FES not converged. Be careful, and check whether your 1-d profiles (or 2-d profiles) as function of any of the three biased CVs you have has actually converged.

Secondly, if the FES has actually converged,  then you are fine using the umbrella sampling type reweighting using the final bias suggested earlier in this thread by Alejandro . If it has not converged (which won't entirely be surprising given you are biasing 3 CVs), you should do a slightly more careful reweighting, as discussed previously in this mailing list for example here: https://groups.google.com/forum/#!searchin/plumed-users/bogdan%7Csort:relevance/plumed-users/pTHUKYShPro/Vy52N2UZBgAJ

- Pratyush

--

Pratyush Tiwary

Post-doctoral Scientist, Columbia University

Department of Chemistry

sites.google.com/site/pratyushtiwary/

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/ac58fb67-78a4-49e6-8b7d-52a41362c32a%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

"finche' andate in bici e fischia il vento
non sara' mai una brutta giornata"
----------------------------------------------------------------
http://sites.google.com/site/andreaspitaleri/
----------------------------------------------------------------
http://it.linkedin.com/in/andreaspitaleri
----------------------------------------------------------------

--

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/CAA0kTsDbqEfBat%3DydfjatoSJJF1_mQbbqb%2BVjSx%2BYke%3D%2B9W8ww%40mail.gmail.com.

[Alejandro Gil-Ley]

Thanks Pratyush for the correction, I didn't specify one should check the convergence to use that kind of reweighting

Best

A

To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@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/ac58fb67-78a4-49e6-8b7d-52a41362c32a%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

"finche' andate in bici e fischia il vento
non sara' mai una brutta giornata"
----------------------------------------------------------------
http://sites.google.com/site/andreaspitaleri/
----------------------------------------------------------------
http://it.linkedin.com/in/andreaspitaleri
----------------------------------------------------------------

--
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...@googlegroups.com.

To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

[billy.willi...@gmail.com]

Thank you for all the advice everyone! :) 

Just a quick question though.  If I were interested in just the relative free energies between the 6 clusters could I simply find the differences in the bias potential required to hold d7, d8 and d9 in place for each cluster?

The lower the energy required to contrain a conformer into position, the more energetically favourable that conformer must be.  If I set the lowest energy as the reference point, could I not just find how less favourable every other conformer is based on the biasing potential of the 5 separate umbrella sampling simulations?

Best regards,

Billy