[gmx-users] Issues with energy when running in the microcanonical ensemble

[Mark Abraham]

On Wed, Jun 12, 2013 at 11:26 AM, C.M.Sampson <cms1g11 at soton.ac.uk> wrote:

> Dear all,
>
> I'm using Gromacs 4.5.5 with the amber ff99SB forcefield to run NVE
> simulations on ethane in 725 TIP3P water molecules.  I used antechamber
> to parametrise my ethane.
>
> I have 2 issues;
>
> Firstly the energy conservation of my simulations seems to get worse as
> I lower my timestep, here is an example .mdp file that I'm using:
>
>
> integrator      = md-vv
> gen_vel         = no
> nsteps          = 1000000
> dt              = 0.001
> nstlist         = 10
> rlist           = 1.3
> coulombtype     = pme
> rcoulomb        = 1.0
> vdw-type        = Shift
> rvdw            = 1.0
> tcoupl          = no
> pcoupl          = no
> nstenergy       = 100
>
> I found a link describing that it could be an issue with the precision
> and that Gromacs can be installed with double precision:
>
> http://www.gromacs.org/Documentation/Terminology/Conserving_Energy
>
> but is this the only way to ensure energy conservation?
>

Many issues with md-vv have been fixed since 4.5.5. Please update and try
again.


>
> My second issue is with the calculation of the kinetic energy, I have
> 2183 atoms in my simulation box which means the kinetic energy should be
> (roughly)
>
> 0.5 * 3*2183-3(degrees of freedom) * 0.001987204118(kB) * 300(T) *
> 4.184(conversion from kcal/mol to kJ/mol) = 8163.97
>
> however, when I look in my simulation the kinetic energy Gromacs outputs
> is 5498.20
>
>            Step           Time         Lambda
>               0        0.00000        0.00000
>
>    Energies (kJ/mol)
>        Bond          Angle Ryckaert-Bell.          LJ-14     Coulomb-14
> 5.95743e+00    1.96046e+01    4.71247e-01    5.73520e-01    3.92661e+00
>     LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
> 5.12842e+03   -3.20604e+04   -2.67111e+03   -2.95726e+04    5.49820e+03
>     Total Energy    Temperature Pressure (bar)
>    -2.40744e+04    3.02575e+02   -1.82251e+03
>
> By running in NVE I'm not using a thermostat and can't find any reason
> in the manual to scale the kinetic energy.  Is there an a reason I'm
> missing as to why my kinetic energy seems lower than it should be?
>

I couldn't follow your equation above, but the most likely issue would be
the use of some algorithm that reduces the number of degrees of freedom,
e.g. rigid water/bonds.

Mark