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

francesco oteri francesco.oteri at gmail.com
Mon Jul 16 14:59:32 CEST 2012 

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)

[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<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/48229a9a/attachment.html>

More information about the gromacs.org_gmx-developers mailing list