<div dir="ltr">Hi gromacs developers,<div>after the mail exchange we had, I succeeded in the implementation and testing of REST2.</div><div style>I adopted the following strategy:</div><div style>1) From a normal .top file, N .top are generated. In each one the forcefield are resscaled according to the REST2 method</div>
<div style>2) The replica exchange is run with a modified version of mdrun. The modifications are:</div><div style> -- The temperature of the different replicas are not cecked </div><div style> -- As suggested by Mark Abraham, the energetic test Va(Xb) + Vb(Xa) - Va(Xa) - Vb(Xb) is evaluated in two steps Initially, replica A evaluates Vb(Xa), Va(Xa) and replica B Va(Xb),Vb(Xb), then the result is sent to the master of the replicas in order to decide the different exchange</div>
<div style><br></div><div style><br></div><div style>As far as I understood, implementing REST2 is in the roadmap of gromacs5, so if the community agrees and nobody is already in charge, I can send a patch for gromacs4.5.5 and, possibly for gromacs4.6, as well as the script to obtain the different topologies on the website as user contribution. Later I can work to implement REST2 in gromacs5 in a less dirty way.</div>
<div style><br></div><div style>Let me know what you think,</div><div style>Francesco</div><div style><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">2012/7/17 francesco oteri <span dir="ltr"><<a href="mailto:francesco.oteri@gmail.com" target="_blank">francesco.oteri@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br><br><div class="gmail_quote">2012/7/16 Shirts, Michael (mrs5pt) <span dir="ltr"><<a href="mailto:mrs5pt@eservices.virginia.edu" target="_blank">mrs5pt@eservices.virginia.edu</a>></span><div class="im">
<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
That sounds reasonable. I think the plan will be, then:<br>
<br>
4.6: a good explanation of how to set up REST2 using lambda values<br></blockquote><div><br></div></div><div>I agree, but in my opinion you should say something regarding the performance drop</div><div class="im"><div> </div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
5.0: an implementation of REST2 (or some variant -- why does it have to be<br>
sqrt(T_0/T_B) scaling that's best?)</blockquote><div> </div></div><div>I don't know. You shoul ask the authors :) I guess it is the most convenient because when T_B=T_0</div><div>sqrt(T_0/T_B) = 1 so you have an unmodified Hemiltonian, so you can use this replica for data analysis</div>
<div><div class="h5">
<div><br></div><div><br></div><div> that is just as fast as normal replica</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
exchange simulations. It will probably still be build on the lambda<br>
machinery for ease of coding, but might be toggleable by a simple mdp<br>
switch. If not, it will be well explained in the documentation.<br>
<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>
> 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>
</div>> Date: Mon, 16 Jul 2012 19:50:13 +0200<br>
<div><div>> 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 Michael,<br>
> I agree with you respect the recoding. But performance difference is of 2x.<br>
> I started inserting the features in 4.5.5 thursday and within today or<br>
> tomorrow<br>
> I will end. I plan one day of debug. So totally 4-5gg...I will get back<br>
> these<br>
> days as performance improvement :)<br>
><br>
> Of course, implementing seriously this features will take more time both<br>
> for coding<br>
> and testing. But right now I am not planning to do it seriously.<br>
><br>
> Maybe later when 4.6 will be out and if it will give me problem.<br>
><br>
><br>
> 2012/7/16 Berk Hess <<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>><br>
><br>
>> On 07/16/2012 06:08 PM, francesco oteri wrote:<br>
>><br>
>><br>
>><br>
>> 2012/7/16 Berk Hess <<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>><br>
>><br>
>>> Using lambda scaling for non-bonded interactions is going to be very<br>
>>> slow.<br>
>>><br>
>><br>
>> If it is done with the actual implementation, you are right. But if<br>
>> mdrun recognize rest2=yes thzn it rescales the<br>
>> parameters at teh beginning. Then a normal MD run goes.<br>
>><br>
>><br>
>>><br>
>>> But this would scale the protein-protein non-bonded interactions with<br>
>>> lambda^2.<br>
>>><br>
>><br>
>><br>
>> Why? computing charge-charge interaction is:<br>
>><br>
>> S*qi*S*qj = (S^2 )*qi*qj<br>
>> but S=sqrt(Bi/B)<br>
>> so what I obtain is:<br>
>> (Bi/B)* qi*qj<br>
>><br>
>> That is exactly what I want: rescaling by (Bi/B) protein-protein<br>
>> interactions<br>
>><br>
>> Ah, sorry, that works indeed.<br>
>><br>
>> Cheers,<br>
>><br>
>> Berk<br>
>><br>
>><br>
>><br>
>>> Is that what you want?<br>
>>><br>
>>> Cheers,<br>
>>><br>
>>> Berk<br>
>>><br>
>>><br>
>>> On 07/16/2012 05:06 PM, francesco oteri wrote:<br>
>>><br>
>>> In the paper they just say that non bonded rescaling is obtained<br>
>>> rescaling<br>
>>> LG epsilon by Bi/B and the protein charges by sqrt(Bi/B).<br>
>>><br>
>>> I did simulation using the lambda approach and I qualitatively<br>
>>> reproduced<br>
>>> the results.<br>
>>><br>
>>> Just to give my contribution, I think the most efficient way to<br>
>>> implement REST2<br>
>>> is using the lambda approach as follows:<br>
>>><br>
>>> 1) in .mdp introducing a new keyword like rest2 = yes/no<br>
>>> 2) using the init_lambda value to rescale the parameters at the bootstrap<br>
>>> of md<br>
>>> (es. in do_md or init_replica_exchange).<br>
>>> 3) running the code as normal MD and swapping the states<br>
>>><br>
>>> Since point 1 implies modification of grompp and point 2 requires the<br>
>>> knowledge of<br>
>>> the data structure, I am just implementing point 3 while point 2 is<br>
>>> demanded to external<br>
>>> scripts.<br>
>>><br>
>>> Francesco<br>
>>><br>
>>> 2012/7/16 Berk Hess <<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>><br>
>>><br>
>>>> On 07/16/2012 04:42 PM, Shirts, Michael (mrs5pt) wrote:<br>
>>>><br>
>>>>> One question, since better REST2 is one item on the todo list for 5.0,<br>
>>>>> and I<br>
>>>>> want to be thinking about the possibilities.<br>
>>>>><br>
>>>>> For solute-solvent energy rescaling, how does one calculate the changes<br>
>>>>> 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<br>
>>>>> 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<br>
>>>>> to<br>
>>>>> do (need new ij terms for every parameter pair) but straightforward in<br>
>>>>> 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<br>
>>>>> adding<br>
>>>>> long term to Gromacs, we'd want to support PME.<br>
>>>>><br>
>>>> No, even without PME this is not straightforward.<br>
>>>> Using a plain cut-off's for electrostatics is horrible, so I'd say you<br>
>>>> want to use reaction-field.<br>
>>>> But with reaction-filed the correction terms are non pair-wise, as with<br>
>>>> PME.<br>
>>>><br>
>>>> Cheers,<br>
>>>><br>
>>>> Berk<br>
>>>><br>
>>>>><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>
>>>>>> 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 16:13:16 +0200<br>
>>>>>> To: Discussion list for GROMACS development <<br>
>>>>>> <a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a>><br>
>>>>>> Subject: Re: [gmx-developers] calculating poteintial energy into<br>
>>>>>> 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>
>>>>>> 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<br>
>>>>>>> system in<br>
>>>>>>> two parts: solvent and solute.<br>
>>>>>>> The two groups have a different APPARENT temperature because the<br>
>>>>>>> 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,<br>
>>>>>>> B=Boltzmann<br>
>>>>>>> constant) and the protein-water interactions are rescaled by<br>
>>>>>>> sqrt(Bi/B).<br>
>>>>>>> So the system temperature, regulated by the thermostat, is equal<br>
>>>>>>> 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,<br>
>>>>>>> in<br>
>>>>>>> order to run N replicas with equal run parameters (temperature,<br>
>>>>>>> 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<br>
>>>>>>> force-field.<br>
>>>>>>> Since the code between the two versions is different I am trying to<br>
>>>>>>> figure<br>
>>>>>>> out how to reinvent the wheel.<br>
>>>>>>><br>
>>>>>>> As shown in [3] and adopted in [2] the same result can be obtained<br>
>>>>>>> using<br>
>>>>>>> the lamda dynamics in gromacs.In this case, only one topology has to<br>
>>>>>>> be<br>
>>>>>>> generated. The state A correspond to the replica at the lowest<br>
>>>>>>> apparent<br>
>>>>>>> temperature, while the state B represents the highest temperature<br>
>>>>>>> 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<br>
>>>>>>> coordinates of<br>
>>>>>>> replica b, so it is most effectively evaluated on the processor that<br>
>>>>>>> has<br>
>>>>>>> the coordinates of replica b. If the topology is the same at the<br>
>>>>>>> replicas a<br>
>>>>>>> and b, then can Va(Xb) be computed by replica b by rescaling suitably<br>
>>>>>>> the<br>
>>>>>>> quantities that contributed to Vb(Xb), based on the knowledge that<br>
>>>>>>> replica<br>
>>>>>>> a is the exchange partner? Then the quantities that are exchanged<br>
>>>>>>> before<br>
>>>>>>> the test are Va(Xb) and Vb(Xa). Should be as fast as regular REMD.<br>
>>>>>>><br>
>>>>>>><br>
>>>>>> This is my idea, indeed<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/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" target="_blank">gmx-developers-request@gromacs.org</a>.<br>
>><br>
><br>
><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/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" target="_blank">gmx-developers-request@gromacs.org</a>.<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/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" target="_blank">gmx-developers-request@gromacs.org</a>.<br>
</div></div></blockquote></div></div></div><div class="HOEnZb"><div class="h5"><br><br clear="all"><div><br></div>-- <br>Cordiali saluti, Dr.Oteri Francesco<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br>Cordiali saluti, Dr.Oteri Francesco
</div>