[gmx-users] Anomaly in Temperature Behavior
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Tue Jan 22 21:12:29 CET 2013
Dear All,
We are doing a simulation on a heterogeneous three phase system (gas,
liquid, solid) surrounded by two silica surfaces which are kept frozen.
To set up the simulation at 300 K, we first perform a 0.5 ns NVT simulation
which runs without a problem. During this run we freeze the solid phase in
all three dimensions. We then take the final configuration of this run,
unfreeze the solid phase, and subject it to a long NVE simulation to study
the dissociation of the solid phase (this is similar to melting of ice)
which is an endothermic process.
However, the initial temperature of the NVE simulation is around 185 K
while the end temperature of NVT is ~300K. Also, the temperature keeps
increasing which is contrary to the physical fact that the dissociation
process is endothermic.
(I have run the same system with DL_POLY 2.20 and observed gradual
temperature decrease as the result of solid decomposition.)
I appreciate being advised on what might be wrong with my system?
Here are the mdp files as well as sample temperature profiles:
For NVT:
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Run control
integrator = md ; Leap-frog algorithm
tinit = 0 ; starting time [ps]
dt = 0.001 ; time step [ps]
nsteps = 500000 ; number of steps
nstcomm = 100 ; frequency for center of mass motion
removal [steps]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Output control
nstxout = 0 ; frequency to write coordinates to
output trajectory file [steps]
nstvout = 0 ; frequency to write velocities to
output trajectory file [steps]
nstfout = 0 ; frequency to write forces to output
trajectory file [steps]
nstlog = 500 ; frequency to write energies to log
file [steps]
nstenergy = 500 ; frequency to write energies to energy
file [steps]
nstxtcout = 0 ; frequency to write coordinates to xtc
trajectory [steps]
xtc-precision = 1000 ; precision to write to xtc trajectory
[real]
xtc_grps = HYDW HYDG SOL
energygrps = HYDW HYDG SOL SiO2 SiOH
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Neighborsearching and short-range nonbonded interactions
nstlist = 10 ; frequency to update the neighbor list
[steps]
ns_type = grid ; (grid / simple) search for
neighboring list
pbc = xyz ; priodic boundary conditions (xyz / no
/ xy)
rlist = 1.0 ; cut-off distance for the short-range
neighbor list [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Electrostatics
coulombtype = PME
rcoulomb = 1.0 ; distance for the Coulomb cut-off [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; van der Waals
vdw-type = switch ; (the LJ normal out at rvdw_switch to
reach zero at rvdw)
rvdw-switch = 0.8 ; where to strat switching the LJ
potential [nm]
rvdw = 0.9 ; cut-off distance for vdw potenrial
[nm]
DispCorr = EnerPres; (Apply long range dispersion
corrections for Energy and Pressure)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; EWALD/PME/PPPM parameters
fourierspacing = 0.12 ; maximum grid spacing for the FFT when
using PPPM or PME [nm]
pme_order = 6 ; interpolation order for PME
ewald_rtol = 1e-06 ; relative strength of the
Ewald-shifted direct potential at rcoulomb
epsilon_surface = 0 ; dipole correction to the Ewald
summation in 3d
optimize_fft = no ; optimal FFT plan for the grid at
startup (yes / no)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Temperature coupling
tcoupl = nose-hoover
tc_grps = system
tau_t = 1.0 ; time constatn for coupling [ps], 1
for each group
ref_t = 300 ; refernce temperature [K], 1 for each
group
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Pressure coupling
Pcoupl = no
;Pcoupl = Parrinello-Rahman
tau_p = 0.5 ; time constant for coupling
compressibility = 4.5e-05
ref_p = 100.0 ; reference pressure for coupling
[bar]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Velocity generation
gen_vel = yes ; generating velocities according to
maxwell distribution
gen_temp = 300 ; [K]
gen_seed = -1 ; initialize random generator based on
the process ID number [integer]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Energy group exclusions
energygrp_excl = SiOH SiOH SiO2 SiO2 SiOH SiO2
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Bonds
constraints = all-angles
constraint-algorithm = shake
;continuation = yes ; do not apply constraints to the
start configuration
shake-tol = 1e-5
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Non-equilibrium MD
freezegrps = HYDW SiOH SiO2
freezedim = Y Y Y Y Y Y Y Y Y
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
For NVE:
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Run control
integrator = md ; Leap-frog algorithm
tinit = 0 ; starting time [ps]
dt = 0.001 ; time step [ps]
nsteps = 1000000 ; number of steps
nstcomm = 100 ; frequency for center of mass motion
removal [steps]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Output control
nstxout = 0 ; frequency to write coordinates to
output trajectory file [steps]
nstvout = 0 ; frequency to write velocities to
output trajectory file [steps]
nstfout = 0 ; frequency to write forces to output
trajectory file [steps]
nstlog = 500 ; frequency to write energies to log
file [steps]
nstenergy = 500 ; frequency to write energies to energy
file [steps]
nstxtcout = 1000 ; frequency to write coordinates to xtc
trajectory [steps]
xtc-precision = 1000 ; precision to write to xtc trajectory
[real]
xtc_grps = HYDW HYDG SOL
energygrps = HYDW HYDG SOL SiO2 SiOH
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Neighborsearching and short-range nonbonded interactions
nstlist = 10 ; frequency to update the neighbor list
[steps]
ns_type = grid ; (grid / simple) search for
neighboring list
pbc = xyz ; priodic boundary conditions (xyz / no
/ xy)
rlist = 1.0 ; cut-off distance for the short-range
neighbor list [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Electrostatics
coulombtype = PME
rcoulomb = 1.0 ; distance for the Coulomb cut-off [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; van der Waals
vdw-type = switch ; (the LJ normal out at rvdw_switch to
reach zero at rvdw)
rvdw-switch = 0.8 ; where to strat switching the LJ
potential [nm]
rvdw = 0.9 ; cut-off distance for vdw potenrial
[nm]
DispCorr = EnerPres; (Apply long range dispersion
corrections for Energy and Pressure)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; EWALD/PME/PPPM parameters
fourierspacing = 0.12 ; maximum grid spacing for the FFT when
using PPPM or PME [nm]
pme_order = 6 ; interpolation order for PME
ewald_rtol = 1e-06 ; relative strength of the
Ewald-shifted direct potential at rcoulomb
epsilon_surface = 0 ; dipole correction to the Ewald
summation in 3d
optimize_fft = no ; optimal FFT plan for the grid at
startup (yes / no)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Temperature coupling
tcoupl = no
;tcoupl = nose-hoover
tc_grps = system
tau_t = 1.0 ; time constatn for coupling [ps], 1
for each group
ref_t = 300 ; refernce temperature [K], 1 for each
group
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Pressure coupling
Pcoupl = no
;Pcoupl = Parrinello-Rahman
tau_p = 0.5 ; time constant for coupling
compressibility = 4.5e-05
ref_p = 100.0 ; reference pressure for coupling
[bar]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Velocity generation
gen_vel = no ; generating velocities according to
maxwell distribution
gen_temp = 300 ; [K]
gen_seed = -1 ; initialize random generator based on
the process ID number [integer]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Energy group exclusions
energygrp_excl = SiOH SiOH SiO2 SiO2 SiOH SiO2
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Bonds
constraints = all-angles
constraint-algorithm = shake
;continuation = yes ; do not apply constraints to the
start configuration
shake-tol = 1e-5
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Non-equilibrium MD
freezegrps = SiOH SiO2
freezedim = Y Y Y Y Y Y
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
NVT ending T sample points:
Time (ps) Temperature (K) 495.5 296.6523 496 301.1944 496.5 301.899 497
298.3572 497.5 302.4087 498 302.5511 498.5 299.8555 499 299.1132 499.5
298.5053 500 302.8503
NVE starting T sample points:
Time (ps) Temperature (K) 0 184.2413 0.5 253.9046 1 256.6712 1.5
253.8046 2 256.9058 2.5 259.6746 3 257.2886 3.5 257.8648 4 260.7298
4.5 257.4706 5 259.5725 5.5 257.4381 6 261.8065 6.5 261.5226 7 258.4074
7.5 257.2928 8 257.778 8.5 260.1319 9 258.1589 9.5 261.3586 10 261.9669
This brings the following question to my mind too:
Does gromacs calculate temperature based on the kinetic energy content of
the unfrozen particles? If this is the case, it makes sense if I am getting
lower T at the beginning of my NVE run, because after starting NVE the
kinetic energy content (coming from NVT with more frozen particles) is
distributed over larger number of free (unfrozen) molecules and hence the T
is lower.
--
S. Alireza Bagherzadeh
PhD Candidate
Dept. of Chem. & Bio. Eng.
University of BC
2360 East Mall
Vancouver BC V6T1Z3
More information about the gromacs.org_gmx-users
mailing list