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

Shirts, Michael (mrs5pt) mrs5pt at eservices.virginia.edu
Mon Jul 16 16:42:21 CEST 2012 

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.

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.
>>> 
>> 
>> 
>> 
>>  --
>> 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