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
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
