[gmx-users] Umbrella Sampling

Discussion:

[gmx-users] Umbrella Sampling - Justin tutorial

Steven Neumann

2011-11-16 09:31:32 UTC

Hi GMX Users,

I am doing Justin tutorial of Umbrella sampling. I have just finished
continous pulling of chainA from the reference chainB. I have some
questions. I looked at the trajectory of pulling and it has began with
dissociating residue 27Alanine from the ChainB following 26, 25, 24...1. My
question is why? As you apply pulling with the constant force to the COM of
the whole chain why does it start with terminal residue following then one
by one? Why not the middle one or any other?

The second thing I would like to extract starting configurations from from
my pulling. Till frame 189 the COM varies from 0.49 to 0.56 - makes sense
as the ChainA is still within ChainB. I would like to use configurations:

0 - 0.50
50 - 0.52
100 - 0.51
150 - 0.51
200 - 0.62
250 - 2.21
....
500 - 5.48

My question is: Do I have to use exactly the same e.g. 0.2 nm spacing or
this configuration above is ok? Can the spacing in nm vary?
And the last thing - is it required to use frames till 189 as the COM
varies in this area?

Thank you!

Steven
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111116/2362d418/attachment.html>

Justin A. Lemkul

2011-11-16 14:31:56 UTC

Permalink

By pulling on the COM of any molecule, the forces are then redistributed to the
other atoms. In the case given in the tutorial, the region you see dissociate
first is held together by relatively weak interactions. If you pull directly on
a specific residue or atom, that residue will dissociate first in all
likelihood. For further discussion pertinent to this specific system, please
refer to our paper linked from the tutorial. You're observing what we observed,
so there is no problem and I am happy that the behavior is reproducible for
others :)

Post by Steven Neumann
The second thing I would like to extract starting configurations from
from my pulling. Till frame 189 the COM varies from 0.49 to 0.56 - makes
sense as the ChainA is still within ChainB. I would like to use
0 - 0.50
50 - 0.52
100 - 0.51
150 - 0.51
200 - 0.62
250 - 2.21
....
500 - 5.48
My question is: Do I have to use exactly the same e.g. 0.2 nm spacing or
this configuration above is ok? Can the spacing in nm vary?

I'm not clear on what you mean here. In constructing a reaction coordinate, you
need to sample at finite intervals along the given path. 0.2 nm is a common
spacing used for many systems, and if you want to reproduce the tutorial's
results, you should use that spacing. Otherwise, I can't guarantee what you
will see. The peptide chains remain bound for a long time, until the force
applied by the harmonic spring overcomes the intermolecular interactions. This
was useful in our study (again, see the paper for why).

Post by Steven Neumann
And the last thing - is it required to use frames till 189 as the COM
varies in this area?

You only need one to represent the center of that particular window, since the
COM distance isn't changing much here.

-Justin
--
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================

Steven Neumann

2011-11-16 14:58:27 UTC

Permalink

Thank you Justin, now I get what you mean!
As I understood I should pick just one frame untill 189 ps (when
dissociation occured) - does it matter which one I will choose on the final
binding free energy?
Then spacing should be 0.2 nm. Right? But how is it possible to do spacing
like this from such results of summary_distances.dat:

189 0.5769072

190 0.5753590

191 0.5711223

192 0.5719844

193 0.5856807

194 0.5886649

195 0.5958101

.......

211 0.7111782

212 0.7322708

213 0.7675126

...

378 4.2101626

379 4.1993918

380 4.2223649

......

500 5.48839

Do you mean the approximate valueof 0.2 nm spacing? As it would be
difficult with these numbers. Can you e.g. do spacing of app. 0.1 nm until
e.g. 3 nm and then increase it to 0.2 nm (I think as in your paper) until
5nm?

Thank you,

Steven

Post by Justin A. Lemkul

By pulling on the COM of any molecule, the forces are then redistributed
to the other atoms. In the case given in the tutorial, the region you see
dissociate first is held together by relatively weak interactions. If you
pull directly on a specific residue or atom, that residue will dissociate
first in all likelihood. For further discussion pertinent to this specific
system, please refer to our paper linked from the tutorial. You're
observing what we observed, so there is no problem and I am happy that the
behavior is reproducible for others :)
The second thing I would like to extract starting configurations from from

Post by Steven Neumann
my pulling. Till frame 189 the COM varies from 0.49 to 0.56 - makes sense
0 - 0.50
50 - 0.52
100 - 0.51
150 - 0.51
200 - 0.62
250 - 2.21
....
500 - 5.48
My question is: Do I have to use exactly the same e.g. 0.2 nm spacing or
this configuration above is ok? Can the spacing in nm vary?

I'm not clear on what you mean here. In constructing a reaction
coordinate, you need to sample at finite intervals along the given path.
0.2 nm is a common spacing used for many systems, and if you want to
reproduce the tutorial's results, you should use that spacing. Otherwise,
I can't guarantee what you will see. The peptide chains remain bound for a
long time, until the force applied by the harmonic spring overcomes the
intermolecular interactions. This was useful in our study (again, see the
paper for why).
And the last thing - is it required to use frames till 189 as the COM

Post by Steven Neumann
varies in this area?

You only need one to represent the center of that particular window, since
the COM distance isn't changing much here.
-Justin
--
==============================**==========
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
==============================**==========
--
gmx-users mailing list gmx-users at gromacs.org
http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
Please search the archive at http://www.gromacs.org/**
Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111116/724944f9/attachment.html>

Justin A. Lemkul

2011-11-16 15:06:13 UTC

Permalink

Post by Steven Neumann
Thank you Justin, now I get what you mean!
As I understood I should pick just one frame untill 189 ps (when
dissociation occured) - does it matter which one I will choose on the
final binding free energy?

The free energy minimum should (theoretically) be when the peptides are still
interacting, i.e. at time zero. Use this frame and base the spacing off of it.

Post by Steven Neumann
Then spacing should be 0.2 nm. Right? But how is it possible to do
189 0.5769072
190 0.5753590
191 0.5711223
192 0.5719844
193 0.5856807
194 0.5886649
195 0.5958101
.......
211 0.7111782
212 0.7322708
213 0.7675126
...
378 4.2101626
379 4.1993918
380 4.2223649
......
500 5.48839
Do you mean the approximate valueof 0.2 nm spacing? As it would be
difficult with these numbers. Can you e.g. do spacing of app. 0.1 nm
until e.g. 3 nm and then increase it to 0.2 nm (I think as in your
paper) until 5nm?

Steven Neumann

2011-11-16 15:17:23 UTC

Permalink

Thank you Justin! That helped a lot!

Steven

Post by Justin A. Lemkul

The free energy minimum should (theoretically) be when the peptides are
still interacting, i.e. at time zero. Use this frame and base the spacing
off of it.
Then spacing should be 0.2 nm. Right? But how is it possible to do spacing

Post by Steven Neumann
189 0.5769072
190 0.5753590
191 0.5711223
192 0.5719844
193 0.5856807
194 0.5886649
195 0.5958101
.......
211 0.7111782
212 0.7322708
213 0.7675126
...
378 4.2101626
379 4.1993918
380 4.2223649
......
500 5.48839
Do you mean the approximate valueof 0.2 nm spacing? As it would be
difficult with these numbers. Can you e.g. do spacing of app. 0.1 nm until
e.g. 3 nm and then increase it to 0.2 nm (I think as in your paper) until
5nm?

Choose approximate values.
My protocol of uneven spacing was necessary to refine the energy minimum
for the suitable data for the paper. I kept the tutorial simple; you
should still achieve a reasonable result. Uneven spacing is uncommon in
the literature, so I did not want to give anyone the impression that it was
a standard protocol, leaving it as an exercise for the reader to understand
the paper and the motivation for uneven spacing.
-Justin
--
==============================**==========
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
==============================**==========
--
gmx-users mailing list gmx-users at gromacs.org
http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
Please search the archive at http://www.gromacs.org/**
Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111116/3b6629ba/attachment.html>

Steven Neumann

2011-11-17 11:55:31 UTC

Permalink

Hi Justin,

I am sorry for so many questions but I do not understand something.
First we run the simulation of pulling Chain A from ChainB with constant
force (pull_k1=1000) and constant velocity of pulling (pull_rate1=0.01) We
extract windows as we discussed and then run simulations with those
configurations as a starting point. I saw the trajectory of one of these
simulations and it looks like normal MD simulation. My question is: Why do
we have in mdp file still pull code as we do not pull protein chain any
more? Pull rate is set to zero but force is still applied... why? Is this
code just used to extract pullf.xvg and pullx.xvg which does not change too
much?
I would appreciate the explanation as without undesratnding the basics its
not good to do any simulation like this.

Thank you

Steven

Are my questions to trivial or noones knows? Please, help!

Post by Steven Neumann
Hi GMX Users,
I am doing Justin tutorial of Umbrella sampling. I have just finished
continous pulling of chainA from the reference chainB. I have some
questions. I looked at the trajectory of pulling and it has began with
dissociating residue 27Alanine from the ChainB following 26, 25, 24...1. My
question is why? As you apply pulling with the constant force to the COM of
the whole chain why does it start with terminal residue following then one
by one? Why not the middle one or any other?
The second thing I would like to extract starting configurations from from
my pulling. Till frame 189 the COM varies from 0.49 to 0.56 - makes sense
0 - 0.50
50 - 0.52
100 - 0.51
150 - 0.51
200 - 0.62
250 - 2.21
....
500 - 5.48
My question is: Do I have to use exactly the same e.g. 0.2 nm spacing or
this configuration above is ok? Can the spacing in nm vary?
And the last thing - is it required to use frames till 189 as the COM
varies in this area?
Thank you!
Steven

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111117/8b9ec4a6/attachment.html>

Justin A. Lemkul

2011-11-17 13:25:46 UTC

Permalink

Post by Steven Neumann
Hi Justin,
I am sorry for so many questions but I do not understand something.
First we run the simulation of pulling Chain A from ChainB with constant
force (pull_k1=1000) and constant velocity of pulling (pull_rate1=0.01)
We extract windows as we discussed and then run simulations with those
configurations as a starting point. I saw the trajectory of one of these
simulations and it looks like normal MD simulation. My question is: Why
do we have in mdp file still pull code as we do not pull protein chain
any more? Pull rate is set to zero but force is still applied... why? Is
this code just used to extract pullf.xvg and pullx.xvg which does not
change too much?
I would appreciate the explanation as without undesratnding the basics
its not good to do any simulation like this.

Steven Neumann

2011-11-17 13:48:04 UTC

Permalink

Thank you Justin, I will read it for sure and come back to my simulations!

Steven

Post by Justin A. Lemkul

Post by Steven Neumann
Hi Justin,
I am sorry for so many questions but I do not understand something.
First we run the simulation of pulling Chain A from ChainB with constant
force (pull_k1=1000) and constant velocity of pulling (pull_rate1=0.01) We
extract windows as we discussed and then run simulations with those
configurations as a starting point. I saw the trajectory of one of these
simulations and it looks like normal MD simulation. My question is: Why do
we have in mdp file still pull code as we do not pull protein chain any
more? Pull rate is set to zero but force is still applied... why? Is this
code just used to extract pullf.xvg and pullx.xvg which does not change too
much?
I would appreciate the explanation as without undesratnding the basics
its not good to do any simulation like this.

Have you read the WHAM paper, or the one specific for g_wham, or perhaps
papers about umbrella sampling in general? You should start there before
diving in to doing the simulations.
The harmonic force applied in the SMD and US simulations is simply a
biasing force to make something happen. With a non-zero pull_rate, motion
in a particular direction is forced to happen. With a pull_rate of zero,
the COM distance is simply restricted - the biasing force maintains the COM
distance between the two defined species, while allowing it to oscillate
according to a harmonic force defined by pull_k1. Thus you establish
sampling overlap between neighboring windows along the reaction coordinate.
-Justin
--
==============================**==========
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
==============================**==========
--
gmx-users mailing list gmx-users at gromacs.org
http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
Please search the archive at http://www.gromacs.org/**
Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111117/53692dbe/attachment.html>