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

francesco oteri francesco.oteri at gmail.com
Mon Jul 16 18:08:07 CEST 2012 

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


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

More information about the gromacs.org_gmx-developers mailing list