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

francesco oteri francesco.oteri at gmail.com
Mon Jul 16 17:06:21 CEST 2012 

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
>>>
>>>
>>>> Mark
>>>>
>>>>
>>>>
>>>>   [1]  http://pubs.acs.org/doi/abs/**10.1021/jp204407d<http://pubs.acs.org/doi/abs/10.1021/jp204407d>
>>>> [2]http://pubs.acs.org/doi/**abs/10.1021/ct100493v<http://pubs.acs.org/doi/abs/10.1021/ct100493v>
>>>> [3] http://onlinelibrary.wiley.**com/doi/10.1002/jcc.21703/**abstract<http://onlinelibrary.wiley.com/doi/10.1002/jcc.21703/abstract>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>   2012/7/16 Berk Hess <hess at kth.se>
>>>>
>>>>    On 7/16/12 14:18 , Shirts, Michael (mrs5pt) wrote:
>>>>>
>>>>>    Actually I did it but I the performance decrease a lot (2x).
>>>>>>
>>>>>>> In particular I run the same run with free_energy  = yes and
>>>>>>> free_energy =
>>>>>>> no and in the second case the run is 2times faster.
>>>>>>> I found in the mailing list that this performance drop is known.
>>>>>>>
>>>>>>>  What are the particular parameters you are trying to change?  I
>>>>>> wrote the
>>>>>> new Hamiltonian Exchange code, so I'm interested in making sure it's
>>>>>> as
>>>>>> useful as possible.  There might be some cases where this performance
>>>>>> drop
>>>>>> can be avoided, and   It's likely the PME parameters, since doing
>>>>>> repeated
>>>>>> that's the only thing that really makes free energy slow.
>>>>>>
>>>>>> UNLESS of course you are changing ALL of the atoms, which which case
>>>>>> all
>>>>>> of
>>>>>> them will go through the free energy code, and it certainly will be
>>>>>> slower.
>>>>>> There are some possible ways to fix this, but that will not be until
>>>>>> 5.0.
>>>>>>    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 14:08:49 +0200
>>>>>>> To: Discussion list for GROMACS development <
>>>>>>> gmx-developers at gromacs.org
>>>>>>> Subject: Re: [gmx-developers] calculating poteintial energy into
>>>>>>> do_md()
>>>>>>>
>>>>>>> Dear Berk
>>>>>>>
>>>>>>> 2012/7/16 Berk Hess <hess at kth.se>
>>>>>>>
>>>>>>>     On 7/16/12 13:05 , francesco oteri wrote:
>>>>>>>
>>>>>>>>    Dear gromacs developers,
>>>>>>>> I am trying to implement in gromacs 4.5.5 a particular Hemiltonian
>>>>>>>> Replica
>>>>>>>> Exchange tecnique.
>>>>>>>> Right now, at every exchange attempt in do_md, I figured out  how to
>>>>>>>> access at the potential
>>>>>>>> energy of replica A (=configuration A at temperature A) and B
>>>>>>>> (=configuration B at temperature B) and so on.
>>>>>>>>
>>>>>>>> In my case, each replica has the same temperature, but there is a
>>>>>>>> different Hemiltonian equation for every replica.
>>>>>>>> The different Hemiltonian are obtained simply changing the force
>>>>>>>> field
>>>>>>>> parameters in the input topology so I dont
>>>>>>>> need to modify anything in gromacs.
>>>>>>>>
>>>>>>>>    But at every exchange attempt I have to test if the
>>>>>>>> configuration B
>>>>>>>> can
>>>>>>>> exist in the state A so I need to calculate
>>>>>>>> its potential energy using the force field data of replica A.
>>>>>>>>
>>>>>>>>    I found that function sum_epot calculates the potential energy
>>>>>>>> but I
>>>>>>>> suspect that it uses values calculated in
>>>>>>>> do_force since sum_epot is called by do_force_lowlevel in turn
>>>>>>>> called
>>>>>>>> by
>>>>>>>> do_force.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>    So my question is, should I call do_force in replica A with
>>>>>>>> coordinates
>>>>>>>> from replica B reach my goal?
>>>>>>>>
>>>>>>>> Yes.
>>>>>>>> Are you changing bonded and/or non-bonded parameters?
>>>>>>>> Some non-bonded parameters might be preprocessed, so you might need
>>>>>>>> reprocess those
>>>>>>>> and them reprocesses again to get back to the original state.
>>>>>>>>
>>>>>>>>
>>>>>>>>   Which parameters need such a preprocessing?
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>   Note that 4.6 will have proper Hamiltonian replica exchange, but
>>>>>>> that
>>>>>>>
>>>>>>>> will
>>>>>>>> use the lambda coupling
>>>>>>>> parameter approach. If you need to do something similar, it might be
>>>>>>>> much
>>>>>>>> simpler to use this code.
>>>>>>>>
>>>>>>>>
>>>>>>>>   Actually I did it but I the performance decrease a lot (2x).
>>>>>>>>
>>>>>>> In particular I run the same run with free_energy  = yes and
>>>>>>> free_energy =
>>>>>>> no and in the second case the run is 2times faster.
>>>>>>> I found in the mailing list that this performance drop is known.
>>>>>>>
>>>>>>>     I think the performance drop is due to the, temporarily, missing
>>>>>> sse
>>>>>>
>>>>> non-bonded kernels.
>>>>> They should be back somewhere this week.
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Berk
>>>>>
>>>>>   Since I dont need to change lambda along the simulation, I prefer
>>>>> using
>>>>>
>>>>>> a
>>>>>>> modified input topology and running a normal REMD changing something
>>>>>>> in
>>>>>>> repl_ex.c and md.c.
>>>>>>>
>>>>>>> Of course, any suggestions to improve performance are wellcome :)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Cheers,
>>>>>>>
>>>>>>>  Berk
>>>>>>>>
>>>>>>>>
>>>>>>>>    Thanks in advance,
>>>>>>>>
>>>>>>>>    Francesco
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> gmx-developers mailing list
>>>>>>>> gmx-developers at gromacs.org
>>>>>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-**developers<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@**gromacs.org<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<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@**gromacs.org<gmx-developers-request at gromacs.org>
>>>>>>> .
>>>>>>>
>>>>>>>  --
>>>>> gmx-developers mailing list
>>>>> gmx-developers at gromacs.org
>>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-**developers<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@**gromacs.org<gmx-developers-request at gromacs.org>
>>>>> .
>>>>>
>>>>>
>>>>
>>>>   --         * This assumes close to cubic cells, which the code tries
>>>> to make.
>>>>
>>>>
>>>> Cordiali saluti, Dr.Oteri Francesco
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> gmx-developers mailing list
>>>> gmx-developers at gromacs.org
>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-**developers<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@**gromacs.org<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<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@**gromacs.org<gmx-developers-request at gromacs.org>
>>> .
>>>
>>
>
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-**developers<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@**gromacs.org<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/672c3f64/attachment.html>

More information about the gromacs.org_gmx-developers mailing list