[gmx-developers] User defined reaction coordinate for umbrella sampling

[Sheeba Jem]

Hi all,

I am trying to evaluate the free energy of peptide folding using umbrella
sampling and I need some help in implementing the reaction coordinate in
gromacs.  The reaction coordinate I want to use is the average C-alpha
distance between the i and i+4 th residues. The umbrella potential would
look like:

psi = 1/2 * k (avg - avg0)^2

avg = (r1+r2....rN)/N

r = CA_i - CA_i+4

avg0 = umbrella position

So if I run five umbrella simulations fixing the average calpha distance at
say 0.6, 0.7, 0.8, 0.9 and 1.0 nm then I should be able to use g_wham to get
the PMF as a function of the average calpha distance of the peptide. I am
studying a short alpha helical peptide, therefore when the peptide is
helical the average calpha distance should be close to 0.6 nm.

As far as I understand it is not possible to define 'average values' as pull
groups and what I want to be able to do is: say when I fix the window at 0.7
nm, at each step during the simulation, I need to calculate the average
Calpha distance of the peptide and calculate the force (and hence the new
set of coordinates) due to the deviation from 0.7 and split the force to
each Calpha atom either equally or depending on the Calpha distance the atom
is involved in.  Is it possible to modify the pull.c code to apply such a
reaction coordinate?

I understand there are other free energy methods like replica exchange etc
to study the thermodynamics of folding and this is not an ideal reaction
coordinate for folding but once I get a hold of what I can do with the pull
code I can try better, more complicated reaction coordinates.

The gromacs version I am using is 4.0.5.

I appreciate any help.

Thanks

Sheeba
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20110722/662fdf8c/attachment.html>