[gmx-developers] energies from previous step used for hamiltonian replica exchange evaluation?

[Berk Hess]

Hi,

I think Floris is right. The exchange is done based on all quantities at 
step n. But what we exchange is actually the state at step n+1 (and v at 
n+0.5 with leap-frog). This means that the system has evolved with the 
wrong Hamiltonian and/or at the wrong state point for one MD step. In 
practice this will be unnoticable, since with e.g. REMD one MD step at a 
few K different temperature will have negligible effect.
Fixing this is non-trivial, since we need the update in case we need the 
pressure or the kinetic energy. One solution is to make a backup of the 
state at step n and then redo the update after exchange. In the current 
code this is hard to achieve, but when we modularize the main MD loop 
this should be easier.

Cheers,

Berk

On 06/11/2014 12:16 AM, Shirts, Michael R. (mrs5pt) wrote:
> Hi, Floris-
>
> I will check this when I can (Travelling at a conference), but the
> bookkeeping in gromacs is a somewhat convoluted (on the list of things to
> clean up) -- a lot of quantities are saved and not stored to file or used
> immediately for various annoying reasons. So I THINK that everything is
> properly accounted for.
>
> I will take an extra look, however, when I get a chance (unless Berk gets
> to it first).  I've certainly compared it to non-replica exchange, and
> gotten statistical identical results -- which does not guarantee it's 100%
> right, however. Note that md-vv should to the bookkeeping correctly as
> well.
>
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> (434)-243-1821
>
>
>
> On 6/10/14, 3:47 PM, "Floris Buelens" <floris_buelens at yahoo.com> wrote:
>
>> Right - all the required quantities come from a single call of do_force
>> before update. But the way I see it (please do correct me), as soon as
>> update() is executed, all the calculated energy quantities are those from
>> the previous set of coordinates, and no longer valid as the basis for a
>> monte carlo move for the new state of the system.
>>
>>
>>
>>
>>
>>
>> On Tuesday, 10 June 2014, 15:37, Berk Hess <hess at kth.se> wrote:
>> Hi,
>>
>> Gromacs doesn't do Hamiltonian replica exchange by re-calculating
>> potentials, does it?
>> All potentials and lambda dependencies are calculated using the
>> coordinates before update, so I don't think there is any issue at all.
>>
>> Cheers,
>>
>> Berk
>>
>>
>>
>>
>> On 06/10/2014 02:56 PM, Floris Buelens wrote:
>>> Hi,
>>>
>>> While trying to evaluate Hamiltonian replica exchange probabilities
>>> 'offline', e.g. not at runtime but from checkpoint files, I think I
>>> might have run into an inconsistency in how they're evaluated in
>>> Gromacs. My worry is that exchange probabilities are calculated at the
>>> very end of the MD step loop, after update() has been called. This would
>>> mean that the exchange criterion is evaluated using energies from step
>>> x, but coordinates from step x+1.
>>>
>>> If I've got this right, the strictly correct way would be to evaluate
>>> replica exchanges before update() is called, and then to recalculate the
>>> energies and forces where the replica exchange attempt was successful
>>> (since they were calculated according to a different U(x) and are
>>> strictly no longer valid) - and only then to call update(). Or
>>> equivalently to implement hybrid MC/MD, where x MD steps are followed by
>>> a single MC step, where the exchange move is either accepted or rejected
>>> based on a fresh energy evaluation.
>>>
>>> Am I missing something?
>>>
>>> I had a look at replica exchange probabilities between successive steps
>>> and predictably the correlation is very high, but some of the deviations
>>> are not so small (some sample data pasted in below).
>>>
>>> thanks,
>>>
>>> Floris
>>>
>>>
>>>    p(n)  p(n+1)
>>> 0.0080  0.0102
>>> 1.0000  1.0000
>>> 0.0137  0.0142
>>> 0.7724  0.6141
>>> 0.0418  0.0619
>>> 0.0731  0.0859
>>> 0.2771  0.2374
>>> 0.0092  0.0130
>>> 0.3064  0.2375
>>> 0.0804  0.0770
>>> 0.0349  0.0325
>>> 1.0000  1.0000
>>> 0.0506  0.0718
>>> 0.1385  0.2349
>>> 0.0003  0.0004
>>> 0.9125  0.8968
>>> 0.0000  0.0000
>>> 0.7520  0.7820
>>> 0.0116  0.0124
>>> 1.0000  0.8525
>>> 0.2255  0.2747
>>> 0.0112  0.0130
>>> 0.0241  0.0230
>>> 1.0000  1.0000
>>> 0.0015  0.0010
>>> 0.0018  0.0027
>>> 0.1365  0.2044
>>> 1.0000  1.0000
>>> 0.0104  0.0062
>>> 0.0238  0.0420
>>> 0.4723  0.5447
>> -- 
>> Gromacs Developers mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
>> or send a mail to gmx-developers-request at gromacs.org.