[gmx-users] Nose-Hoover does not couple
Inon Sharony
InonShar at TAU.ac.IL
Tue Jul 3 16:38:37 CEST 2012
I'm performing a verification of the energy dissipation of a single atom,
thermally coupled to a heat drain (or sink). In other words, the Langevin
equation of motion is d/dt v(t) = - gamma * v(t).
If using the Stochastic Dynamics integrator I do indeed get dissipation,
however if using the Nose-Hoover thermostat with either the leap-frog or
velocity Verlet integrators the energy does not seem to dissipate. I get
the same thing for two harmonically-bound atoms. Could anyone please take
a look and tell me why this is?
Also, the "Conserved Energy" I get is "not a number". I don't know what
this means, either.
energy.xvg:
@ title "Gromacs Energies"
@ xaxis label "Time (ps)"
@ yaxis label "(kJ/mol), (K), (bar)"
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "LJ (SR)"
@ s1 legend "Kinetic En."
@ s2 legend "Total Energy"
@ s3 legend "Conserved En."
@ s4 legend "Temperature"
@ s5 legend "Pressure"
0.000000 0.000000000000 16.030000000000
16.030000000000 nan 1285.303055615475 0.000000000000
0.001000 0.000000000000 16.030000000000
16.030000000000 nan 1285.303055615475 0.000000000000
0.002000 0.000000000000 16.030000000000
16.030000000000 nan 1285.303055615475 0.000000000000
.
.
.
0.998000 0.000000000000 16.030000000000
16.030000000000 nan 1285.303055615475 0.000000000000
0.999000 0.000000000000 16.030000000000
16.030000000000 nan 1285.303055615475 0.000000000000
=============================================
Supporting files:
traj.trr:
traj-tss.trr frame 0:
natoms= 1 step= 0 time=0.0000000e+00 lambda=
0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.15000e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
f (1x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
traj-tss.trr frame 1:
natoms= 1 step= 1 time=1.0000000e-03 lambda=
0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.15100e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
f (1x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
traj-tss.trr frame 2:
natoms= 1 step= 2 time=2.0000000e-03 lambda=
0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.15200e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
f (1x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
.
.
.
traj-tss.trr frame 998:
natoms= 1 step= 998 time=9.9800000e-01 lambda=
0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 4.14800e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
f (1x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
traj-tss.trr frame 999:
natoms= 1 step= 999 time=9.9900000e-01 lambda=
0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 4.14900e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
f (1x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
=============================================
conf.gro:
Single Sulfur Atom
1
1ATM S 1 3.150 3.150 3.150 0.000 0.000 1.000
6.37511 6.37511 6.37511
=============================================
topol.itp:
[ moleculetype ]
; Name nrexcl
ATM 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass
1 S 1 ATM S 1 0.000 32.0600
=============================================
grompp.mdp:
integrator = md
dt = 0.001
nsteps = 999
;
;------------------------------------------------------------------------------------------------
nstlog = 1
nstcalcenergy = 1
;------------------------------------------------------------------------------------------------
;
pbc = no ; no periodic boundary conditions (only one
molecule)
comm_mode = None ;"Angular" = remove center of mass
translatory AND rotational motion
nstcomm = 1 ; number of steps between removal of center
of mass motion
;
;------------------------------------------------------------------------------------------------
; number of steps between writing out of different phase data and
energetics
nstxout = 1
nstvout = 1
nstfout = 1
nstenergy = 1
; groups to write energetics of
energygrps = 0
;
;-------------------------------------------------------------------------------------------------
; neighbor list search
nstype = simple ; particle (as opposed to domain)
decompositon
coulombtype = cut-off
;
;-------------------------------------------------------------------------------------------------
;
; Langevin dynamics temperature coupling
;
ld_seed = 1 ; (1993) [integer]
tcoupl = nose-hoover
nsttcouple = 1
;
tc-grps = 0
tau_t = 1 ; mass/gamma [a.m.u. nanosecond]
ref_t = 0 ; reference (bath) temperature
;-------------------------------------------------------------------------------------------------
;
; Velocity generation
gen_vel = no ; (no)
gen_temp = 0 ; (300) [K]
gen_seed = 1 ; (-1) = from PID
=============================================
topol-tss.tpr:
inputrec:
integrator = md
nsteps = 999
init_step = 0
ns_type = Simple
nstlist = 10
ndelta = 2
nstcomm = 0
comm_mode = None
nstlog = 1
nstxout = 1
nstvout = 1
nstfout = 1
nstcalcenergy = 1
nstenergy = 1
nstxtcout = 0
init_t = 0
delta_t = 0.001
xtcprec = 1000
nkx = 0
nky = 0
nkz = 0
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 0
epsilon_surface = 0
optimize_fft = FALSE
ePBC = no
bPeriodicMols = FALSE
bContinuation = FALSE
bShakeSOR = FALSE
etc = Nose-Hoover
nsttcouple = 1
epc = No
epctype = Isotropic
nstpcouple = -1
tau_p = 1
ref_p (3x3):
ref_p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref_p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref_p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress (3x3):
compress[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
refcoord_scaling = No
posres_com (3):
posres_com[0]= 0.00000e+00
posres_com[1]= 0.00000e+00
posres_com[2]= 0.00000e+00
posres_comB (3):
posres_comB[0]= 0.00000e+00
posres_comB[1]= 0.00000e+00
posres_comB[2]= 0.00000e+00
andersen_seed = 815131
rlist = 1
rlistlong = 1
rtpi = 0.05
coulombtype = Cut-off
rcoulomb_switch = 0
rcoulomb = 1
vdwtype = Cut-off
rvdw_switch = 0
rvdw = 1
epsilon_r = 1
epsilon_rf = 1
tabext = 1
implicit_solvent = No
gb_algorithm = Still
gb_epsilon_solvent = 80
nstgbradii = 1
rgbradii = 1
gb_saltconc = 0
gb_obc_alpha = 1
gb_obc_beta = 0.8
gb_obc_gamma = 4.85
gb_dielectric_offset = 0.009
sa_algorithm = Ace-approximation
sa_surface_tension = 2.05016
DispCorr = No
free_energy = no
init_lambda = 0
delta_lambda = 0
n_foreign_lambda = 0
sc_alpha = 0
sc_power = 0
sc_sigma = 0.3
sc_sigma_min = 0.3
nstdhdl = 10
separate_dhdl_file = yes
dhdl_derivatives = yes
dh_hist_size = 0
dh_hist_spacing = 0.1
nwall = 0
wall_type = 9-3
wall_atomtype[0] = -1
wall_atomtype[1] = -1
wall_density[0] = 0
wall_density[1] = 0
wall_ewald_zfac = 3
pull = no
disre = No
disre_weighting = Conservative
disre_mixed = FALSE
dr_fc = 1000
dr_tau = 0
nstdisreout = 100
orires_fc = 0
orires_tau = 0
nstorireout = 100
dihre-fc = 1000
em_stepsize = 0.01
em_tol = 10
niter = 20
fc_stepsize = 0
nstcgsteep = 1000
nbfgscorr = 10
ConstAlg = Lincs
shake_tol = 0.0001
lincs_order = 4
lincs_warnangle = 30
lincs_iter = 1
bd_fric = 0
ld_seed = 1
cos_accel = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
grpopts:
nrdf: 3
ref_t: 0
tau_t: 1
anneal: No
ann_npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp_flags[ 0]: 0
efield-x:
n = 0
efield-xt:
n = 0
efield-y:
n = 0
efield-yt:
n = 0
efield-z:
n = 0
efield-zt:
n = 0
bQMMM = FALSE
QMconstraints = 0
QMMMscheme = 0
scalefactor = 1
qm_opts:
ngQM = 0
header:
bIr = present
bBox = present
bTop = present
bX = present
bV = present
bF = not present
natoms = 1
lambda = 0.000000e+00
topology:
name="Single Sulfur Atom"
#atoms = 1
molblock (0):
moltype = 0 "ATM"
#molecules = 1
#atoms_mol = 1
#posres_xA = 0
#posres_xB = 0
ffparams:
atnr=1
ntypes=1
functype[0]=LJ_SR, c6= 9.98400640e-03, c12= 1.30754560e-05
reppow = 12
fudgeQQ = 1
cmap
atomtypes:
atomtype[ 0]={radius=-1.00000e+00, volume=-1.00000e+00,
gb_radius=-1.00000e+00, surftens=-1.00000e+00, atomnumber= 16,
S_hct=-1.00000e+00)}
moltype (0):
name="ATM"
atoms:
atom (1):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m=
3.20600e+01, q= 0.00000e+00, mB= 3.20600e+01, qB= 0.00000e+00, resind=
0, atomnumber= 16}
atom (1):
atom[0]={name="S"}
type (1):
type[0]={name="S",nameB="S"}
residue (1):
residue[0]={name="ATM", nr=1, ic=' '}
cgs:
nr=1
cgs[0]={0..0}
excls:
nr=1
nra=1
excls[0][0..0]={0}
grp[T-Coupling ] nr=1, name=[ 0]
grp[Energy Mon. ] nr=1, name=[ 0]
grp[Acceleration] nr=1, name=[ rest]
grp[Freeze ] nr=1, name=[ rest]
grp[User1 ] nr=1, name=[ rest]
grp[User2 ] nr=1, name=[ rest]
grp[VCM ] nr=1, name=[ rest]
grp[XTC ] nr=1, name=[ rest]
grp[Or. Res. Fit] nr=1, name=[ rest]
grp[QMMM ] nr=1, name=[ rest]
grpname (5):
grpname[0]={name="System"}
grpname[1]={name="Other"}
grpname[2]={name="ATM"}
grpname[3]={name="0"}
grpname[4]={name="rest"}
groups T-Cou Energ Accel Freez User1 User2 VCM XTC Or. R
QMMM
allocated 0 0 0 0 0 0 0 0
0 0
groupnr[ *] = 0 0 0 0 0 0 0 0
0 0
box (3x3):
box[ 0]={ 6.37511e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.37511e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.37511e+00}
box_rel (3x3):
box_rel[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
box_rel[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
box_rel[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
boxv (3x3):
boxv[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
boxv[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
boxv[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
pres_prev (3x3):
pres_prev[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
pres_prev[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
pres_prev[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
svir_prev (3x3):
svir_prev[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
svir_prev[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
svir_prev[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
fvir_prev (3x3):
fvir_prev[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
fvir_prev[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
fvir_prev[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
nosehoover_xi: not available
x (1x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.15000e+00}
v (1x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
Group statistics
T-Coupling : 1 (total 1 atoms)
Energy Mon. : 1 (total 1 atoms)
Acceleration: 1 (total 1 atoms)
Freeze : 1 (total 1 atoms)
User1 : 1 (total 1 atoms)
User2 : 1 (total 1 atoms)
VCM : 1 (total 1 atoms)
XTC : 1 (total 1 atoms)
Or. Res. Fit: 1 (total 1 atoms)
QMMM : 1 (total 1 atoms)
--
Inon Sharony
J+N+W+N% ShR+W+N+J+
972-3-6407634
Please consider your environmental responsibility before printing this e-mail.
More information about the gromacs.org_gmx-users
mailing list