[gmx-developers] calculating poteintial energy into do_md()
Shirts, Michael (mrs5pt)
mrs5pt at eservices.virginia.edu
Tue Feb 26 18:47:57 CET 2013
I don't know that we want to add any functionality into 4.6 at this point --
that's not what minor versions are for. In terms of getting REST2 into 5.0,
I'm happy to make sure we have the conversation so that REST2 easily fits
into the general formalism for 5.0. It SHOULD be very simple to put in if
we organize things properly.
The best way to do this, Francesco, is for you to outline a list of what
would need to be changed/added as a feature for 5.0 at
http://redmine.gromacs.org/, and we can start having the conversation there.
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>
> Date: Tue, 26 Feb 2013 18:41:19 +0100
> To: "michael.shirts at virginia.edu" <michael.shirts at virginia.edu>, Discussion
> list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] calculating poteintial energy into do_md()
>
> 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
More information about the gromacs.org_gmx-developers
mailing list