[gmx-users] Question about conserved energy in MTTK

[Shun Sakuraba]

Dear list,

I am trying to use MTTK barostat in GROMACS 4.5.5.
After analyzing the result for a while, I found that the conserved energy (not total energy) of MTTK is drifting during the simulation.
The .xvg, .edr files are uploaded at [1] and [2]. It is drifting with a constant ratio of ca. -185 kJ/mol/ps.

I cannot believe this is an expected behavior, so could anyone point out where I am wrong in my simulation setup? I found similar report at [3] but seems it was when 4.5 was in pre-release stage.

Thanks in advance for your help!

* Simulation detail
The system consists of 1000 SPC-E water molecules, and the time step is set to 0.5 fs, just in case the long timestep harms the conserevation (c.f. [3]). The interaction energy is set to switching version, just in case, too. Changing these parameter does not seem to improve the conservation.
The double precision version of GROMACS is used (single precision version also has the same problem).
The system has been pre-equilibrated with Berendsen pressure coupling simulation with the same pressure and temperature.

[1] https://www.dropbox.com/s/16bgfhvavqcw1f8/mttk05.xvg
[2] https://www.dropbox.com/s/tg8butw39rmz8qk/mttk05.edr
[3] http://lists.gromacs.org/pipermail/gmx-developers/2010-January/003979.html

== .mdp file contents follow

integrator = md-vv
define =

dt = 0.0005
nsteps = 1000000 ; 500 ps

coulombtype = PME-Switch
vdwtype = Switch
pbc = xyz

rlist = 1.2
rcoulomb = 1.0
rcoulomb_switch = 0.9
rvdw = 1.0
rvdw_switch = 0.9
nstlist = 1

tinit = 0
tcoupl = nose-hoover
tc_grps = System
tau_t = 0.5
ref_t = 300.0
nsttcouple = 1

pcoupl = MTTK
pcoupltype = isotropic
compressibility = 4.5e-5
ref_p = 1.01325
tau_p = 0.5
refcoord_scaling = no
nstpcouple = 1

constraints = hbonds
constraint_algorithm = LINCS

nstxtcout = 100
nstlog = 100
nstenergy = 100
nstvout = 0
nstxout = 1000

-- 
Shun SAKURABA, Ph.D.
Postdoc @ Molecular Modeling & Simulation Group, Japan Atomic Energy Agency