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

francesco oteri francesco.oteri at gmail.com
Mon Jul 16 16:13:16 CEST 2012 

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

More information about the gromacs.org_gmx-developers mailing list