[gmx-users] Conserved energy ("Conserved En.") in NVT simulation (Wade)

Johnny Lu johnny.lu128 at gmail.com
Sun Oct 19 01:35:23 CEST 2014 

>From Bussi paper: "If we use this analysis and go to the limit of an
infinitesimal time step, we find [Equation 15] where the last two terms
come from the integration of Eq. 7 along the trajectory. Note that a
similar integration along the path is present in H_{Nosé}. However, in our
scheme a stochastic integration is also necessary. In the continuum limit
the changes in energy induced by the rescaling compensate exactly the
fluctuations in H. For a finite time step this compensation is only
approximate and the conservation of tilde{H} provides a measure on the
accuracy of the integration. This accuracy has to be interpreted in the
sense of the ability of generating configurations representative of the
ensemble. The physical meaning of Eq. 15 is that the fluxes of energy
between the system and the thermostat are exactly balanced."

Does that mean even with a finite timestep, if the conformations in the
trajectory are all correctly sampled from the intended ensemble without any
biases, tilde H should still be constant?

It is a bit unfortunate that they didn't say how the drift rate of tilde H
related to deviation of distribution of conformations from the intended
ensemble.

On Sat, Oct 18, 2014 at 7:02 PM, Johnny Lu <johnny.lu128 at gmail.com> wrote:

> The bussi paper is at [
> marge.uochb.cas.cz/~roesel/teach/Clanky/paper5_Bussi_v_rescale.pdf
> <http://marge.uochb.cas.cz/%7Eroesel/teach/Clanky/paper5_Bussi_v_rescale.pdf>
> ].
>
> I seem to vaguely remember someone wrote if tilde H drop linearly, it
> means the energy is conserved. But I'm not sure now.
>
> On Sat, Oct 18, 2014 at 6:58 PM, Johnny Lu <johnny.lu128 at gmail.com> wrote:
>
>> Figure 2 of the Bussi paper has a diagram with constant tilde H for NVT
>> simulation with 2fs time step, and another one with rising tilde H for NVT
>> simulation with 40fs timestep, and say: "for tilde H in the NVT
>> calculation. We notice that, while in the NVE calculation the system
>> explodes at some point, the NVT calculation is always stable, thanks to the
>> thermostat. However, in spite of the stability, the drift in tilde H
>> indicates that thesampling is inaccurate under these conditions."
>>
>> So, tilde H should be constant even in a NVT simulation, and I don't
>> exactly agree with the old mail-list post at [
>> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-August/083330.html
>> ]
>>
>> On Sat, Oct 18, 2014 at 6:51 PM, Johnny Lu <johnny.lu128 at gmail.com>
>> wrote:
>>
>>> Does this old maillist post say tilde H is conserved energy? [
>>> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-August/083330.html
>>> ]
>>>
>>> On Sat, Oct 18, 2014 at 12:46 PM, Mark Abraham <mark.j.abraham at gmail.com
>>> > wrote:
>>>
>>>> On Sat, Oct 18, 2014 at 4:32 PM, Wade <wadelv at foxmail.com> wrote:
>>>>
>>>> > Dear Mark,
>>>> >     Thank your very much for your reply. However, I still confused
>>>> because
>>>> > results showed that the conserved En. was time dependent (with small
>>>> > fluctuation); but, the total energy in the ener file has small drift
>>>> with
>>>> > high fluctuations.
>>>> >
>>>>
>>>> "Both these things vary" doesn't mean anything. Numbers, or it didn't
>>>> happen ;-)
>>>>
>>>>
>>>> > Further, I rechecked the code  (gmx-4.5.5) ?in these days, and I
>>>> noticed a
>>>> > part of code might be useful ("enerd->term[F_ECONSERVED] =
>>>> > enerd->term[F_ETOT] +
>>>> compute_conserved_from_auxiliary(ir,state,&MassQ);").
>>>> > It will be actived when the leap-frog integrator is employed.
>>>> > To assure the the conserved En. is either H or H-tilde, I changed the
>>>> code
>>>> > as this: "enerd->term[F_ECONSERVED] = enerd->term[F_ETOT];". Result
>>>> shown
>>>> > that the "total energy" and "conserved En." are same now in the
>>>> output of
>>>> > ener file?.
>>>> >
>>>>
>>>> Great. A = A...
>>>>
>>>>
>>>> > Dose it means that the conserved En. in ener file is H-tilde?
>>>> >
>>>>
>>>> ...which doesn't say what B is... The question is whether there's an
>>>> intrinsic time dependence on the value returned by
>>>> compute_conserved_from_auxiliary. You can only do that by removing the
>>>> *other* term.
>>>>
>>>>
>>>> >
>>>> > ?
>>>> >
>>>> > If yes, why does the conserved En. drifted with a constant rate in a
>>>> 10 ns
>>>> > ?simulation for a pure water box (216 tip4p water)?
>>>> >
>>>> >
>>>> Because it's characteristic of the drift of your setup or of the
>>>> implementation of what is reported as the conserved quantity. I
>>>> suggested a
>>>> diagnostic procedure. Did you try it? If not, why not?
>>>>
>>>> Mark
>>>>
>>>> With best wishes,
>>>> >
>>>> >
>>>> > Wade
>>>> >
>>>> >
>>>> > (part of the code please see below, I hope I have got the key point).
>>>> > ?
>>>> >
>>>> >
>>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>>> >         /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
>>>> > /* sum up the foreign energy and dhdl terms */
>>>> > sum_dhdl(enerd,state->lambda,ir);         /* use the directly
>>>> determined
>>>> > last velocity, not actually the averaged half steps */        if
>>>> (bTrotter
>>>> > && ir->eI==eiVV)         {            enerd->term[F_EKIN] = last_ekin;
>>>> >   }        enerd->term[F_ETOT] = enerd->term[F_EPOT] +
>>>> > enerd->term[F_EKIN];                if (bVV)        {
>>>> > enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
>>>> > saved_conserved_quantity;        }        else         {
>>>> > enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
>>>> > compute_conserved_from_auxiliary(ir,state,&MassQ);?
>>>> >
>>>> >
>>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> >
>>>> > >That could well be true. I did look carefully at the code earlier
>>>> this
>>>> > >year, and decided it was correct and wrote output of sensible things.
>>>> > >However I don't remember whether the conserved quantity reported is
>>>> H or
>>>> > >H-tilde, and don't have time now to go back. But IIRC because one of
>>>> them
>>>> > >has a time dependence, you can run double-precision GROMACS in a
>>>> > >conservative NVE setup and observe no drift, then add this
>>>> thermostat, and
>>>> > >from that run you will know which of H and H-tilde is reported as the
>>>> > >conserved quantity.
>>>> >
>>>> > >Mark
>>>> >
>>>> > > Is that right? I really need your help.
>>>> > > Any suggestion is valuable.
>>>> > >
>>>> > > Wade?
>>>> > > --‍
>>>> > --
>>>> > Gromacs Users mailing list
>>>> >
>>>> > * Please search the archive at
>>>> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
>>>> > send a mail to gmx-users-request at gromacs.org.
>>>> >
>>>> --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>>
>>>
>>>
>>
>

More information about the gromacs.org_gmx-users mailing list