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 class="im">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 class="h5">
<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 class="h5">
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 class="h5">
<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/<u></u>10.1021/jp204407d</a><br>
[2]<a href="http://pubs.acs.org/doi/abs/10.1021/ct100493v" target="_blank">http://pubs.acs.org/doi/<u></u>abs/10.1021/ct100493v</a><br>
[3] <a href="http://onlinelibrary.wiley.com/doi/10.1002/jcc.21703/abstract" target="_blank">http://onlinelibrary.wiley.<u></u>com/doi/10.1002/jcc.21703/<u></u>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>
Since I dont need to change lambda along the simulation, I prefer using<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
a<br>
modified input topology and running a normal REMD changing something in<br>
repl_ex.c and md.c.<br>
<br>
Of course, any suggestions to improve performance are wellcome :)<br>
<br>
<br>
<br>
<br>
Cheers,<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Berk<br>
<br>
<br>
Thanks in advance,<br>
<br>
Francesco<br>
<br>
<br>
<br>
<br>
<br>
<br>
--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>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" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
<br>
<br>
</blockquote>
--<br>
Cordiali saluti, Dr.Oteri Francesco<br>
--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>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" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
<br>
</blockquote></blockquote>
--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>developers</a><br>
Please don't post (un)subscribe requests to the list. Use the www<br>
interface or send it to <a href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
<br>
</blockquote>
<br>
<br></div></div>
-- * This assumes close to cubic cells, which the code tries to make.<div class="im"><br>
<br>
Cordiali saluti, Dr.Oteri Francesco<br>
<br>
<br>
<br>
<br>
<br>
--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>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" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
<br>
</div></blockquote><div class="im">
<br>
<br>
-- <br>
Cordiali saluti, Dr.Oteri Francesco<br>
-- <br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>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" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
</div></blockquote></blockquote><div class="HOEnZb"><div class="h5">
<br>
<br>
-- <br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-<u></u>developers</a><br>
Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@<u></u>gromacs.org</a>.<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br>Cordiali saluti, Dr.Oteri Francesco<br>
</div></div>