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

[Johnny Lu]

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