[gmx-users] Energy conservation problem in NVE simulation for a system with frozen atoms

[Yi Ren]

Hi all,

I am simulating a system with two slabs and water molecules in between.
The slabs were frozen during simulations and I've exclude the interactions
between slab atoms. After energy minimization, I run nvt at 300 K for 10
ns, followed by annealing the system to 250 K in 1 ns. Then another nvt
was run at 250 K for 6 ns. So far, everything was normal. Then I switched
to nve simulation by turning off the thermostat. However, the energy of
the system kept decreasing (temperature was reaching 0 K). I was using
Group cut-off for all the simulations. I've checked the energy between
different groups. The energies between slab and slab, both coulombic and
LJ are 0, which make sense because I've exclude the interactions. The
energies between water and water, on the other hand, decreased to a huge
extend. I am not sure why this happened. If I remove the slabs, only
simulate water molecules, then it works just fine.

Please let me know if you have ever encountered such problem. Any
suggestions would be greatly appreciated.

I am using Gromacs/5.1.4 and the parameters in .mdp file are listed as
follows:
Input Parameters:
   integrator                     = md
   tinit                          = 0
   dt                             = 0.002
   nsteps                         = 25000000
   init-step                      = 0
   simulation-part                = 1
   comm-mode                      = Linear
   nstcomm                        = 100
   bd-fric                        = 0
   ld-seed                        = 603782470
   emtol                          = 10
   emstep                         = 0.01
   niter                          = 20
   fcstep                         = 0
   nstcgsteep                     = 1000
   nbfgscorr                      = 10
   rtpi                           = 0.05
   nstxout                        = 250000
   nstvout                        = 250000
   nstfout                        = 0
   nstlog                         = 250000
   nstcalcenergy                  = 100
   nstenergy                      = 250000
   nstxout-compressed             = 0
   compressed-x-precision         = 1000
   cutoff-scheme                  = Group
   nstlist                        = 10
   ns-type                        = Grid
   pbc                            = xyz
   periodic-molecules             = FALSE
   verlet-buffer-tolerance        = 0.005
   rlist                          = 1
   rlistlong                      = 1
   nstcalclr                      = 0
   coulombtype                    = PME
   coulomb-modifier               = None
   rcoulomb-switch                = 0
   rcoulomb                       = 1
   epsilon-r                      = 1
   epsilon-rf                     = inf
   vdw-type                       = Cut-off
   vdw-modifier                   = None
   rvdw-switch                    = 0
   rvdw                           = 1
   DispCorr                       = EnerPres
   table-extension                = 1
   fourierspacing                 = 0.16
   fourier-nx                     = 36
   fourier-ny                     = 36
   fourier-nz                     = 48
   pme-order                      = 4
   ewald-rtol                     = 1e-05
   ewald-rtol-lj                  = 0.001
   lj-pme-comb-rule               = Geometric
   ewald-geometry                 = 0
   epsilon-surface                = 0
   implicit-solvent               = No
   gb-algorithm                   = Still
   nstgbradii                     = 1
   rgbradii                       = 1
   gb-epsilon-solvent             = 80
   gb-saltconc                    = 0
   gb-obc-alpha                   = 1
   gb-obc-beta                    = 0.8
   gb-obc-gamma                   = 4.85
   gb-dielectric-offset           = 0.009
   sa-algorithm                   = Ace-approximation
   sa-surface-tension             = 2.05016
   tcoupl                         = No
   nsttcouple                     = -1
   nh-chain-length                = 0
   print-nose-hoover-chain-variables = FALSE
   pcoupl                         = No
   pcoupltype                     = Isotropic
   nstpcouple                     = -1
   tau-p                          = 1
   compressibility (3x3):
      compressibility[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      compressibility[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      compressibility[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   ref-p (3x3):
      ref-p[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   refcoord-scaling               = No
   posres-com (3):
      posres-com[0]= 0.00000e+00
      posres-com[1]= 0.00000e+00
      posres-com[2]= 0.00000e+00
   posres-comB (3):
      posres-comB[0]= 0.00000e+00
      posres-comB[1]= 0.00000e+00
      posres-comB[2]= 0.00000e+00
   QMMM                           = FALSE
   QMconstraints                  = 0
   QMMMscheme                     = 0
   MMChargeScaleFactor            = 1
qm-opts:
   ngQM                           = 0
   constraint-algorithm           = Lincs
   continuation                   = TRUE
   Shake-SOR                      = FALSE
   shake-tol                      = 0.0001
   lincs-order                    = 6
   lincs-iter                     = 1
   lincs-warnangle                = 30
   nwall                          = 0
   wall-type                      = 9-3
   wall-r-linpot                  = -1
   wall-atomtype[0]               = -1
   wall-atomtype[1]               = -1
   wall-density[0]                = 0
   wall-density[1]                = 0
   wall-ewald-zfac                = 3
   pull                           = FALSE
   rotation                       = FALSE
   interactiveMD                  = FALSE
   disre                          = No
   disre-weighting                = Conservative
   disre-mixed                    = FALSE
   dr-fc                          = 1000
   dr-tau                         = 0
   nstdisreout                    = 100
   orire-fc                       = 0
   orire-tau                      = 0
   nstorireout                    = 100
   free-energy                    = no
   cos-acceleration               = 0
   deform (3x3):
      deform[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   simulated-tempering            = FALSE
   E-x:
      n = 0
   E-xt:
      n = 0
   E-y:
      n = 0
   E-yt:
      n = 0
   E-z:
      n = 0
   E-zt:
      n = 0
   swapcoords                     = no
   adress                         = FALSE
   userint1                       = 0
   userint2                       = 0
   userint3                       = 0
   userint4                       = 0
   userreal1                      = 0
   userreal2                      = 0
   userreal3                      = 0
   userreal4                      = 0
grpopts:
   nrdf:       24285
   ref-t:           0
   tau-t:           0
annealing:          No
annealing-npoints:           0
   acc:            0           0           0
   nfreeze:           Y           Y           Y           N           N
       N
   energygrp-flags[  0]: 1 0
   energygrp-flags[  1]: 0 0


Thanks again and best regards,
Yi