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

francesco oteri francesco.oteri at gmail.com
Mon Jul 16 19:50:13 CEST 2012 

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120716/812217d6/attachment.html>

More information about the gromacs.org_gmx-developers mailing list