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

Berk Hess hess at kth.se
Mon Jul 16 16:46:56 CEST 2012 

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
>>> [2]http://pubs.acs.org/doi/abs/10.1021/ct100493v
>>> [3] 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
>>>>>>> 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.
>>>>
>>>
>>>
>>>   --         * 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
>>> 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.

More information about the gromacs.org_gmx-developers mailing list