[lammps-users] Shake and energy conservation

Hi all,

I am trying to equilibrate a water droplet. The system consists of 2440 water molecules (SPC/E). I am using Shake to constrain the water molecules. An equilibration run was run 100ps (fix nve with fix temp/rescale) fixing the temperature at 220K. Following that, when I try to run a constant energy simulation (removing temp/rescale), the particle seems to have equilibrated (temp fluctuates about 220K). However the potential energy (and as a result, the total energy) seems to be all over the place.

Step Temp PotEng KinEng TotEng
        0 220.64311 -25249.309 3208.8949 -22040.415
       10 220.09195 -24686.818 3200.8791 -21485.939
       20 218.81464 -18948.086 3182.3028 -15765.783
       30 220.07978 -18136.302 3200.7022 -14935.6

Obviuosly I am doing something wrong. I am attaching the input script. If required, I can email the data file.

Any help would be hugely appreciated.

in.waterFAtestNVE (1.03 KB)

SHAKE adds fictitious forces so I wouldn't expect great
energy conservation. You're also using a Coulomb
cutoff of 10 Ang, and no long-range Ewald or PPPM, which
is not what SPC/E is parameterized for, so that would also
add to errors. However if you've run a long NVT and
equilibrated to one T, then run NVE, I think the system
should stay (on average) at the same T.