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

Wade wadelv at foxmail.com
Sat Oct 18 16:34:39 CEST 2014

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.

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?. 
Dose it means that the conserved En. in ener file is H-tilde?


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)?

With best wishes,


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


> Is that right? I really need your help.
> Any suggestion is valuable.
> Wade?
> --‍

More information about the gromacs.org_gmx-users mailing list