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

[Mark Abraham]

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