[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,
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?
> --
More information about the gromacs.org_gmx-users
mailing list