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

francesco oteri francesco.oteri at gmail.com
Tue Feb 26 18:41:19 CET 2013 

Hi gromacs developers,
after the mail exchange we had, I succeeded in the implementation and
testing of REST2.
I adopted the following strategy:
1) From a normal .top file, N .top are generated. In each one the
forcefield are resscaled according to the REST2 method
2) The replica exchange is run with a modified version of mdrun. The
modifications are:
    -- The temperature of the different replicas are not cecked
    -- 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


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.

Let me know what you think,
Francesco



2012/7/17 francesco oteri <francesco.oteri at gmail.com>

>
>
> 2012/7/16 Shirts, Michael (mrs5pt) <mrs5pt at eservices.virginia.edu>
>
> That sounds reasonable.  I think the plan will be, then:
>>
>> 4.6: a good explanation of how to set up REST2 using lambda values
>>
>
> I agree, but in my opinion you should say something regarding the
> performance drop
>
>
>>
>> 5.0: an implementation of REST2 (or some variant -- why does it have to be
>> sqrt(T_0/T_B) scaling that's best?)
>
>
> I don't know. You shoul ask the authors :) I guess it is the most
> convenient because when T_B=T_0
> sqrt(T_0/T_B)  = 1 so you have an unmodified Hemiltonian, so you can use
> this replica for data analysis
>
>
>  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.
>>
>> --
>> 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
>



-- 
Cordiali saluti, Dr.Oteri Francesco
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20130226/babbcad8/attachment.html>

More information about the gromacs.org_gmx-developers mailing list