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

[Justin Lemkul]

On 1/22/13 4:58 AM, Andrey Frolov wrote:
>   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:
> "

This approach is unsound.  There are insufficient degrees of freedom in single 
molecules to justify distinct thermostats.

> 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.

Considering all data points will include early frames that are not necessarily 
anywhere close to being equilibrated.  You could try block averaging to see if 
the results converge to the correct result over the latter part of the 
simulation, or at least see if the system is approaching the target temperature.

> 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.
>

It would be better to diagnose the system using PME.  Even if the results are 
"similar," rounding errors are compounded when using plain cutoffs, which render 
the results inherently unreliable.

-Justin

> 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 "
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================