[gmx-developers] calculating poteintial energy into do_md()

Shirts, Michael (mrs5pt) mrs5pt at eservices.virginia.edu
Mon Jul 16 21:01:58 CEST 2012 

That sounds reasonable.  I think the plan will be, then:

4.6: a good explanation of how to set up REST2 using lambda values

5.0: an implementation of REST2 (or some variant -- why does it have to be
sqrt(T_0/T_B) scaling that's best?) that is just as fast as normal replica
exchange simulations. It will probably still be build on the lambda
machinery for ease of coding, but might be toggleable by a simple mdp
switch.  If not, it will be well explained in the documentation.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821


> From: francesco oteri <francesco.oteri at gmail.com>
> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Date: Mon, 16 Jul 2012 19:50:13 +0200
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] calculating poteintial energy into do_md()
> 
> Hi Michael,
> I agree with you respect the recoding. But performance difference is of 2x.
> I started inserting the features in 4.5.5  thursday and within today or
> tomorrow
> I will end. I plan one day of debug. So totally 4-5gg...I will get back
> these
> days as performance improvement :)
> 
> Of course, implementing seriously this features will take more time both
> for coding
> and testing. But right now I am not planning to do it seriously.
> 
> Maybe later when 4.6 will be out and if it will give me problem.
> 
> 
> 2012/7/16 Berk Hess <hess at kth.se>
> 
>>  On 07/16/2012 06:08 PM, francesco oteri wrote:
>> 
>> 
>> 
>> 2012/7/16 Berk Hess <hess at kth.se>
>> 
>>>  Using lambda scaling for non-bonded interactions is going to be very
>>> slow.
>>> 
>> 
>>  If it is done with the actual implementation, you are right. But if
>> mdrun recognize rest2=yes thzn it rescales the
>> parameters at teh beginning. Then a normal MD run goes.
>> 
>> 
>>> 
>>> But this would scale the protein-protein non-bonded interactions with
>>> lambda^2.
>>> 
>> 
>> 
>>  Why? computing charge-charge interaction is:
>> 
>>  S*qi*S*qj = (S^2 )*qi*qj
>> but S=sqrt(Bi/B)
>> so what I obtain is:
>> (Bi/B)* qi*qj
>> 
>>  That is exactly what I want: rescaling by (Bi/B) protein-protein
>> interactions
>> 
>>   Ah, sorry, that works indeed.
>> 
>> Cheers,
>> 
>> Berk
>> 
>> 
>> 
>>>  Is that what you want?
>>> 
>>> Cheers,
>>> 
>>> Berk
>>> 
>>> 
>>> On 07/16/2012 05:06 PM, francesco oteri wrote:
>>> 
>>> In the paper they just say that non bonded rescaling is obtained
>>> rescaling
>>> LG epsilon by Bi/B and the protein charges by sqrt(Bi/B).
>>> 
>>>  I did simulation using the lambda approach and I qualitatively
>>> reproduced
>>> the results.
>>> 
>>>  Just to give my contribution, I think the most efficient way to
>>> implement REST2
>>> is using the lambda approach as follows:
>>> 
>>>  1) in .mdp introducing a new keyword like rest2 = yes/no
>>> 2) using the init_lambda value to rescale the parameters at the bootstrap
>>> of md
>>>    (es. in do_md or init_replica_exchange).
>>> 3) running the code as normal MD and swapping the states
>>> 
>>>  Since point 1 implies modification of grompp and point 2 requires the
>>> knowledge of
>>> the data structure,  I am just implementing point 3 while point 2 is
>>> demanded to external
>>> scripts.
>>> 
>>>  Francesco
>>> 
>>> 2012/7/16 Berk Hess <hess at kth.se>
>>> 
>>>> On 07/16/2012 04:42 PM, Shirts, Michael (mrs5pt) wrote:
>>>> 
>>>>> One question, since better REST2 is one item on the todo list for 5.0,
>>>>> and I
>>>>> want to be thinking about the possibilities.
>>>>> 
>>>>> For solute-solvent energy rescaling, how does one calculate the changes
>>>>> in
>>>>> solute-solvent energy without doing multiple PME calls to decompose the
>>>>> electrostatic energy?
>>>>> 
>>>>> I'm guessing that by since your are load multiple topologies of the
>>>>> system,
>>>>> for each of those topologies, you can explicitly write new nonbonded
>>>>> pairwise parameter terms between the water and protein.  This is a pain
>>>>> to
>>>>> do (need new ij terms for every parameter pair) but straightforward in
>>>>> the
>>>>> end to automate.  Is this how your are planning to do it?
>>>>> 
>>>>> Without PME, then it's straightforward to decompose into solute-solvent
>>>>> energy / solute-solute / solvent-solvent electrostatics, though for
>>>>> adding
>>>>> long term to Gromacs, we'd want to support PME.
>>>>> 
>>>>  No, even without PME this is not straightforward.
>>>> Using a plain cut-off's for electrostatics is horrible, so I'd say you
>>>> want to use reaction-field.
>>>> But with reaction-filed the correction terms are non pair-wise, as with
>>>> PME.
>>>> 
>>>> Cheers,
>>>> 
>>>> Berk
>>>> 
>>>>> 
>>>>> Best,
>>>>> ~~~~~~~~~~~~
>>>>> Michael Shirts
>>>>> Assistant Professor
>>>>> Department of Chemical Engineering
>>>>> University of Virginia
>>>>> michael.shirts at virginia.edu
>>>>> (434)-243-1821
>>>>> 
>>>>> 
>>>>>    From: francesco oteri <francesco.oteri at gmail.com>
>>>>>> Reply-To: Discussion list for GROMACS development <
>>>>>> gmx-developers at gromacs.org>
>>>>>> Date: Mon, 16 Jul 2012 16:13:16 +0200
>>>>>> To: Discussion list for GROMACS development <
>>>>>> gmx-developers at gromacs.org>
>>>>>> Subject: Re: [gmx-developers] calculating poteintial energy into
>>>>>> do_md()
>>>>>> 
>>>>>> Hi
>>>>>> 
>>>>>> 2012/7/16 Mark Abraham <Mark.Abraham at anu.edu.au>
>>>>>> 
>>>>>>    On 16/07/2012 10:59 PM, francesco oteri wrote:
>>>>>>> 
>>>>>>> Hi,
>>>>>>> what I am trying to do is implementing the REST2 tecnique [1]
>>>>>>> The goal is enanching the conforlational sampling splitting the
>>>>>>> system in
>>>>>>> two parts: solvent and solute.
>>>>>>> The two groups have a different APPARENT temperature because the
>>>>>>> bonded
>>>>>>> interaction within the solute atoms are rescaled by a factor Bi/B (
>>>>>>> Bi=1/Ti*Kb,  Ti = apparent temperature for the i-th replica,
>>>>>>> B=Boltzmann
>>>>>>> constant) and the  protein-water interactions are rescaled by
>>>>>>> sqrt(Bi/B).
>>>>>>> So the system temperature, regulated by the thermostat, is equal
>>>>>>> along the
>>>>>>> different replicas but the dynamic of the different replicas
>>>>>>> changes becaus of the different Boltzmann factor.
>>>>>>> To implement the tecnique, the naive approach (adopted in [2]) is:
>>>>>>> 1) Generating N-different topology (N= number of replicas).
>>>>>>> 2) Removing same check at the beginning of the Replica Exchange code,
>>>>>>> in
>>>>>>> order to run N replicas with equal run parameters (temperature,
>>>>>>> pression,
>>>>>>> ecc. ecc.)
>>>>>>> 3) Changing the acceptance ratio formula.
>>>>>>> 
>>>>>>>   This approach has been developped for gmx4.0.3 but I need to use
>>>>>>> gmx4.5.5 because it has a native implementation of the CHARMM
>>>>>>> force-field.
>>>>>>> Since the code between the two versions is different I am trying to
>>>>>>> figure
>>>>>>> out how to reinvent the wheel.
>>>>>>> 
>>>>>>>   As shown in [3] and adopted in [2] the same result can be obtained
>>>>>>> using
>>>>>>> the lamda dynamics in gromacs.In this case, only one topology has to
>>>>>>> be
>>>>>>> generated. The state A correspond to the replica at the lowest
>>>>>>> apparent
>>>>>>> temperature, while the state B represents the highest temperature
>>>>>>> replica.
>>>>>>> Intermediate replicas are generated changing the init_lambda value.
>>>>>>> 
>>>>>>>   What has been found, this approach results in slowest run.
>>>>>>> 
>>>>>>> 
>>>>>>>   So, right now I got stuck at the point 3. In fact  to find:
>>>>>>> 
>>>>>>>   delta = B( Va(Xb) + Vb(Xa) - Va(Xa) - Vb(Xb))
>>>>>>> 
>>>>>>>   I need to calculate Va(Xb) and Vb(Xa)
>>>>>>> 
>>>>>>> 
>>>>>>> Va(Xb) is determined by the potential of replica a and the
>>>>>>> coordinates of
>>>>>>> replica b, so it is most effectively evaluated on the processor that
>>>>>>> has
>>>>>>> the coordinates of replica b. If the topology is the same at the
>>>>>>> replicas a
>>>>>>> and b, then can Va(Xb) be computed by replica b by rescaling suitably
>>>>>>> the
>>>>>>> quantities that contributed to Vb(Xb), based on the knowledge that
>>>>>>> replica
>>>>>>> a is the exchange partner? Then the quantities that are exchanged
>>>>>>> before
>>>>>>> the test are Va(Xb) and Vb(Xa). Should be as fast as regular REMD.
>>>>>>> 
>>>>>>> 
>>>>>> This is my idea, indeed
>>>>>> 
>>>>>> 
>>>>>>> 
>> 
>> --
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
>> 
> 
> 
> 
> -- 
> Cordiali saluti, Dr.Oteri Francesco
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.

More information about the gromacs.org_gmx-developers mailing list