<br><br><div class="gmail_quote">2012/7/16 Berk Hess <span dir="ltr"><<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
<div>Using lambda scaling for non-bonded
interactions is going to be very slow.<br></div></div></blockquote><div><br></div><div>If it is done with the actual implementation, you are right. But if mdrun recognize rest2=yes thzn it rescales the</div><div>parameters at teh beginning. Then a normal MD run goes.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div text="#000000" bgcolor="#FFFFFF"><div>
<br>
But this would scale the protein-protein non-bonded interactions
with lambda^2.<br></div></div></blockquote><div><br></div><div><br></div><div>Why? computing charge-charge interaction is:</div><div><br></div><div>S*qi*S*qj = (S^2 )*qi*qj</div><div>but S=sqrt(Bi/B) </div><div>so what I obtain is:</div>
<div>(Bi/B)* qi*qj</div><div><br></div><div>That is exactly what I want: rescaling by (Bi/B) protein-protein interactions</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF"><div>
Is that what you want?<br>
<br>
Cheers,<br>
<br>
Berk<div><div class="h5"><br>
<br>
On 07/16/2012 05:06 PM, francesco oteri wrote:<br>
</div></div></div><div><div class="h5">
<blockquote type="cite">
In the paper they just say that non bonded rescaling is obtained
rescaling
<div>LG epsilon by Bi/B and the protein charges by sqrt(Bi/B).</div>
<div><br>
</div>
<div>I did simulation using the lambda approach and I
qualitatively reproduced </div>
<div>the results.</div>
<div><br>
</div>
<div>Just to give my contribution, I think the most efficient way
to implement REST2</div>
<div>is using the lambda approach as follows:</div>
<div><br>
</div>
<div>1) in .mdp introducing a new keyword like rest2 = yes/no</div>
<div>2) using the init_lambda value to rescale the parameters at
the bootstrap of md </div>
<div> (es. in do_md or init_replica_exchange). </div>
<div>3) running the code as normal MD and swapping the states </div>
<div><br>
</div>
<div>Since point 1 implies modification of grompp and point 2
requires the knowledge of </div>
<div>the data structure, I am just implementing point 3 while
point 2 is demanded to external </div>
<div>scripts.</div>
<div>
<br>
</div>
<div>Francesco</div>
<div>
<div><br>
<div class="gmail_quote">2012/7/16 Berk Hess <span dir="ltr"><<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>On 07/16/2012 04:42 PM, Shirts, Michael
(mrs5pt) wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
One question, since better REST2 is one item on the
todo list for 5.0, and I<br>
want to be thinking about the possibilities.<br>
<br>
For solute-solvent energy rescaling, how does one
calculate the changes in<br>
solute-solvent energy without doing multiple PME calls
to decompose the<br>
electrostatic energy?<br>
<br>
I'm guessing that by since your are load multiple
topologies of the system,<br>
for each of those topologies, you can explicitly write
new nonbonded<br>
pairwise parameter terms between the water and
protein. This is a pain to<br>
do (need new ij terms for every parameter pair) but
straightforward in the<br>
end to automate. Is this how your are planning to do
it?<br>
<br>
Without PME, then it's straightforward to decompose
into solute-solvent<br>
energy / solute-solute / solvent-solvent
electrostatics, though for adding<br>
long term to Gromacs, we'd want to support PME.<br>
</blockquote>
</div>
No, even without PME this is not straightforward.<br>
Using a plain cut-off's for electrostatics is horrible, so
I'd say you want to use reaction-field.<br>
But with reaction-filed the correction terms are non
pair-wise, as with PME.<br>
<br>
Cheers,<br>
<br>
Berk<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div>
<br>
Best,<br>
~~~~~~~~~~~~<br>
Michael Shirts<br>
Assistant Professor<br>
Department of Chemical Engineering<br>
University of Virginia<br>
<a href="mailto:michael.shirts@virginia.edu" target="_blank">michael.shirts@virginia.edu</a><br>
(434)-243-1821<br>
<br>
<br>
</div>
</div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div>
From: francesco oteri <<a href="mailto:francesco.oteri@gmail.com" target="_blank">francesco.oteri@gmail.com</a>><br>
Reply-To: Discussion list for GROMACS development
<<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a>><br>
Date: Mon, 16 Jul 2012 16:13:16 +0200<br>
To: Discussion list for GROMACS development <<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a>><br>
Subject: Re: [gmx-developers] calculating
poteintial energy into do_md()<br>
<br>
Hi<br>
<br>
2012/7/16 Mark Abraham <<a href="mailto:Mark.Abraham@anu.edu.au" target="_blank">Mark.Abraham@anu.edu.au</a>><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On 16/07/2012 10:59 PM, francesco oteri wrote:<br>
<br>
Hi,<br>
what I am trying to do is implementing the REST2
tecnique [1]<br>
The goal is enanching the conforlational
sampling splitting the system in<br>
two parts: solvent and solute.<br>
The two groups have a different APPARENT
temperature because the bonded<br>
interaction within the solute atoms are rescaled
by a factor Bi/B (<br>
Bi=1/Ti*Kb, Ti = apparent temperature for the
i-th replica, B=Boltzmann<br>
constant) and the protein-water interactions
are rescaled by sqrt(Bi/B).<br>
So the system temperature, regulated by the
thermostat, is equal along the<br>
different replicas but the dynamic of the
different replicas<br>
changes becaus of the different Boltzmann
factor.<br>
To implement the tecnique, the naive approach
(adopted in [2]) is:<br>
1) Generating N-different topology (N= number of
replicas).<br>
2) Removing same check at the beginning of the
Replica Exchange code, in<br>
order to run N replicas with equal run
parameters (temperature, pression,<br>
ecc. ecc.)<br>
3) Changing the acceptance ratio formula.<br>
<br>
This approach has been developped for gmx4.0.3
but I need to use<br>
gmx4.5.5 because it has a native implementation
of the CHARMM force-field.<br>
Since the code between the two versions is
different I am trying to figure<br>
out how to reinvent the wheel.<br>
<br>
As shown in [3] and adopted in [2] the same
result can be obtained using<br>
the lamda dynamics in gromacs.In this case, only
one topology has to be<br>
generated. The state A correspond to the replica
at the lowest apparent<br>
temperature, while the state B represents the
highest temperature replica.<br>
Intermediate replicas are generated changing the
init_lambda value.<br>
<br>
What has been found, this approach results in
slowest run.<br>
<br>
<br>
So, right now I got stuck at the point 3. In
fact to find:<br>
<br>
delta = B( Va(Xb) + Vb(Xa) - Va(Xa) - Vb(Xb))<br>
<br>
I need to calculate Va(Xb) and Vb(Xa)<br>
<br>
<br>
Va(Xb) is determined by the potential of replica
a and the coordinates of<br>
replica b, so it is most effectively evaluated
on the processor that has<br>
the coordinates of replica b. If the topology is
the same at the replicas a<br>
and b, then can Va(Xb) be computed by replica b
by rescaling suitably the<br>
quantities that contributed to Vb(Xb), based on
the knowledge that replica<br>
a is the exchange partner? Then the quantities
that are exchanged before<br>
the test are Va(Xb) and Vb(Xa). Should be as
fast as regular REMD.<br>
<br>
</blockquote>
<br>
This is my idea, indeed<br>
<br>
</div>
</div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div>
<br>
Mark<br>
<br>
<br>
<br>
[1] <a href="http://pubs.acs.org/doi/abs/10.1021/jp204407d" target="_blank">http://pubs.acs.org/doi/abs/10.1021/jp204407d</a><br>
[2]<a href="http://pubs.acs.org/doi/abs/10.1021/ct100493v" target="_blank">http://pubs.acs.org/doi/abs/10.1021/ct100493v</a><br>
[3] <a href="http://onlinelibrary.wiley.com/doi/10.1002/jcc.21703/abstract" target="_blank">http://onlinelibrary.wiley.com/doi/10.1002/jcc.21703/abstract</a><br>
<br>
<br>
<br>
<br>
<br>
2012/7/16 Berk Hess <<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On 7/16/12 14:18 , Shirts, Michael (mrs5pt)
wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Actually I did it but I the performance
decrease a lot (2x).<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
In particular I run the same run with
free_energy = yes and<br>
free_energy =<br>
no and in the second case the run is
2times faster.<br>
I found in the mailing list that this
performance drop is known.<br>
<br>
</blockquote>
What are the particular parameters you are
trying to change? I wrote the<br>
new Hamiltonian Exchange code, so I'm
interested in making sure it's as<br>
useful as possible. There might be some
cases where this performance<br>
drop<br>
can be avoided, and It's likely the PME
parameters, since doing<br>
repeated<br>
that's the only thing that really makes free
energy slow.<br>
<br>
UNLESS of course you are changing ALL of the
atoms, which which case all<br>
of<br>
them will go through the free energy code,
and it certainly will be<br>
slower.<br>
There are some possible ways to fix this,
but that will not be until 5.0.<br>
Best,<br>
~~~~~~~~~~~~<br>
Michael Shirts<br>
Assistant Professor<br>
Department of Chemical Engineering<br>
University of Virginia<br>
<a href="mailto:michael.shirts@virginia.edu" target="_blank">michael.shirts@virginia.edu</a><br>
(434)-243-1821<br>
<br>
<br>
From: francesco oteri <<a href="mailto:francesco.oteri@gmail.com" target="_blank">francesco.oteri@gmail.com</a>><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Reply-To: Discussion list for GROMACS
development <<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a>><br>
Date: Mon, 16 Jul 2012 14:08:49 +0200<br>
To: Discussion list for GROMACS
development <<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
Subject: Re: [gmx-developers] calculating
poteintial energy into do_md()<br>
<br>
Dear Berk<br>
<br>
2012/7/16 Berk Hess <<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>><br>
<br>
On 7/16/12 13:05 , francesco oteri
wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Dear gromacs developers,<br>
I am trying to implement in gromacs
4.5.5 a particular Hemiltonian<br>
Replica<br>
Exchange tecnique.<br>
Right now, at every exchange attempt in
do_md, I figured out how to<br>
access at the potential<br>
energy of replica A (=configuration A at
temperature A) and B<br>
(=configuration B at temperature B) and
so on.<br>
<br>
In my case, each replica has the same
temperature, but there is a<br>
different Hemiltonian equation for every
replica.<br>
The different Hemiltonian are obtained
simply changing the force field<br>
parameters in the input topology so I
dont<br>
need to modify anything in gromacs.<br>
<br>
But at every exchange attempt I have
to test if the configuration B<br>
can<br>
exist in the state A so I need to
calculate<br>
its potential energy using the force
field data of replica A.<br>
<br>
I found that function sum_epot
calculates the potential energy but I<br>
suspect that it uses values calculated
in<br>
do_force since sum_epot is called by
do_force_lowlevel in turn called<br>
by<br>
do_force.<br>
<br>
<br>
<br>
So my question is, should I call
do_force in replica A with<br>
coordinates<br>
from replica B reach my goal?<br>
<br>
Yes.<br>
Are you changing bonded and/or
non-bonded parameters?<br>
Some non-bonded parameters might be
preprocessed, so you might need<br>
reprocess those<br>
and them reprocesses again to get back
to the original state.<br>
<br>
<br>
Which parameters need such a
preprocessing?<br>
</blockquote>
<br>
<br>
Note that 4.6 will have proper
Hamiltonian replica exchange, but that<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
will<br>
use the lambda coupling<br>
parameter approach. If you need to do
something similar, it might be<br>
much<br>
simpler to use this code.<br>
<br>
<br>
Actually I did it but I the
performance decrease a lot (2x).<br>
</blockquote>
In particular I run the same run with
free_energy = yes and<br>
free_energy =<br>
no and in the second case the run is
2times faster.<br>
I found in the mailing list that this
performance drop is known.<br>
<br>
</blockquote>
I think the performance drop is due to
the, temporarily, missing sse<br>
</blockquote>
non-bonded kernels.<br>
They should be back somewhere this week.<br>
<br>
Cheers,<br>
<br>
Berk<br>
<br>
</blockquote>
</div>
</div>
</blockquote>
</blockquote>
</blockquote>
</blockquote>
</div>
</div>
</div>
</blockquote>
<br>
</div></div></div>
<br>--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</a><br>
Please don't post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br></blockquote></div><br><br clear="all"><div><br></div>-- <br>Cordiali saluti, Dr.Oteri Francesco<br>