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

[Johnny Lu]

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.
>>>
>>
>>
>