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

Mark Abraham mark.j.abraham at gmail.com
Tue Feb 26 18:59:22 CET 2013 

Agreed. People are most welcome to create whatever they wish, but no new
functionality will go in the official 4.5 or 4.6 branches.

And while I'm here, any new functionality for 5.0 will have to be prepared
to deal with the fact that there will be plenty of API changes going on
while we transition to C++. If any new functionality causes problems and
its author is not around or unable to help fix it, then that functionality
won't make 5.0 either. That might well apply to features from the 4.x era
and existing new stuff in master branch too! Our resources are very
stretched, and we can't promise to maintain everything that's ever been
done. The good news is that version control means nothing is ever truly
gone, if someone comes along who wants to put work into it. :-)

Cheers,

Mark

On Tue, Feb 26, 2013 at 6:47 PM, Shirts, Michael (mrs5pt) <
mrs5pt at eservices.virginia.edu> wrote:

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

More information about the gromacs.org_gmx-developers mailing list