[gmx-users] Temperature coupling problem: reference T is different to given T in mdp file

[Andrey Frolov]

 Dear gmx-users,

1) I am running simulation of two paracetamol molecules dissolved in supercritical CO2 in NVT ensemble.  The problem is that the resulting temperature of the system is larger than the temperature specified in mdp file by ~ 20 K (i set ref_t = 410 K in mdp file).  Berendsen and v-rescale thermostats  behave similarly. The output from g_energy:
"
Statistics over 250001 steps [ 0.0000 through 500.0000 ps ], 1 data sets
All statistics are over 2501 points
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
T-System     430.342    1.1 14.1796 6.30923 (K)
"


2) When i run the same simulation but set 3 different subsystems (all CO2 molecules, the first paracetamol molecule, the second paracetamol molecule)  coupled to different bathes with the same temperature, the resulting temperatures do not match the specified temperature = 410 K. The output from g_energy: 
"
Statistics over 250001 steps [ 0.0000 through 500.0000 ps ], 4 data sets
All statistics are over 2501 points
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Temperature    428.82            1.2     14.356 -3.39055 (K)
T-CO2              427.432          1.2     14.6848 -5.24 (K)
T-PARA            452.434          3.7     91.2873 20.656 (K)
T-PAR2            454.211          10      93.7436 37.8375 (K) "

I would appreciate very much if someone gives me a hint what could be wrong.
Thank you.

Andrey
----
Andrey I. Frolov, PhD
Postdoctoral fellow
Institute of Solution Chemistry
Russian Academy of Sciences


Some information about system:
There are 200 CO2 molecules, cubic box: 4 x 4 x 4 nm3. I set "coulombtype = Cut-off" just for testing purposes, PME gives similar results. 

Some mdp control options (for the 2) case):
"
; RUN CONTROL PARAMETERS = 
integrator = md
dt = 0.002
nsteps = 250000
comm-mode = Linear
nstcomm = 1
; NEIGHBORSEARCHING PARAMETERS = 
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.1
; OPTIONS FOR ELECTROSTATICS AND VDW = 
coulombtype = Cut-off
rcoulomb = 1.1
vdw_type = Shift 
rvdw = 1.0
; OPTIONS FOR WEAK COUPLING ALGORITHMS = 
tcoupl = v-rescale
tc-grps = CO2 PARA PAR2
tau_t = 1.0 1.0 1.0
ref_t = 410 410 410
; GENERATE VELOCITIES FOR STARTUP RUN = 
gen_vel = yes
gen_temp = 410 
gen_seed = 473529
; OPTIONS FOR BONDS = 
constraints = hbonds
constraint_algorithm = lincs
unconstrained_start = no
shake_tol = 0.00001
lincs_order = 4
lincs_warnangle = 30
morse = no
lincs_iter = 2 "