Good morning.
While doing some verification, I tried to simulate a harmonically bound
diatomic molecule where one of the atoms is coupled to a thermal reservoir
held at zero Kelvin (a thermal "drain" or "sink"). The initial conditions
were such that the atom coupled to the bath starts from rest and the other
has some initial velocity (set in conf.gro). I notice that while the
initial velocities are read in properly, the velocities in subsequent
time-steps are far smaller, and continue to fluctuate among these smaller
values (by small I mean a difference of two orders of magnitude, when the
initial velocity was 1 nm/ps). This is also visible in a sudden drop in
the kinetic energy and temperature immediately after the first time step.
The problem is that I've set COMM=None.Is there a separate mechanism which
subtracts velocities after the beginning of the first time-step? This does
not happen when there is no coupling to the thermal bath.
Many thanks,
--
Inon Sharony
J+N+W+N% ShR+W+N+J+
972-3-6407634
Please consider your environmental responsibility before printing this e-mail.
Further information:
( conf.gro , topol.itp , traj.trr , energy.xvg , topol.tpr )
conf.gro:
Molecule starting from Minimal Potential Energy Configuration with
non-zero initial velocity
2
1SAU S 1 3.150 3.150 3.270 0.000 0.000 0.000
1SAU S 2 3.150 3.150 3.030 0.000 0.000 1.000
6.37511 6.37511 6.37511
___________________________________________________________________________
grompp.mdp:
integrator = md-vv-avek
dt = 0.001
nsteps = 999
;
;------------------------------------------------------------------------------------------------
nstlog = 1
nstcalcenergy = 1
;------------------------------------------------------------------------------------------------
;
pbc = no ; no periodic boundary conditions (only one
molecule)
comm_mode = None
nstcomm = 1
;
;------------------------------------------------------------------------------------------------
; 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 1
;
;-------------------------------------------------------------------------------------------------
; neighbor list search
nstype = simple ; particle (as opposed to domain)
decomposition
coulombtype = cut-off
;
;-------------------------------------------------------------------------------------------------
;
; Langevin dynamics temperature coupling
;
ld_seed = 1 ; (1993) [integer]
tcoupl = v-rescale
nsttcouple = 1
;
tc-grps = 0 1
tau_t = 1 0 ; mass/gamma [a.m.u. nanosecond]
ref_t = 0 0 ; reference (bath) temperature
;-------------------------------------------------------------------------------------------------
;
; Velocity generation
gen_vel = no ; (no)
gen_temp = 0 ; (300) [K]
gen_seed = 1 ; (-1) = from PID
___________________________________________________________________________
topol.itp:
[ moleculetype ]
; Name nrexcl
SAU 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass
1 S 1 SAU S 1 0.000 32.0600
2 S 1 SAU S 2 0.000 32.0600
[ bonds ]
; ai aj fu c0, c1, ...
1 2 1 0.238 680000.0 ; harmonic
___________________________________________________________________________
traj.trr:
traj.trr frame 0:
natoms= 2 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 (2x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.27000e+00}
x[ 1]={ 3.15000e+00, 3.15000e+00, 3.03000e+00}
v (2x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, -2.12102e-02}
v[ 1]={ 0.00000e+00, 0.00000e+00, 1.02121e+00}
f (2x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, -1.36000e+03}
f[ 1]={ 0.00000e+00, 0.00000e+00, 1.36000e+03}
traj.trr frame 1:
natoms= 2 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 (2x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.26996e+00}
x[ 1]={ 3.15000e+00, 3.15000e+00, 3.03002e+00}
v (2x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, -6.29244e-02}
v[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
f (2x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, -1.31673e+03}
f[ 1]={ 0.00000e+00, 0.00000e+00, 1.31673e+03}
traj.trr frame 2:
natoms= 2 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 (2x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.26987e+00}
x[ 1]={ 3.15000e+00, 3.15000e+00, 3.03004e+00}
v (2x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, -1.02810e-01}
v[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
f (2x3):
f[ 0]={ 0.00000e+00, 0.00000e+00, -1.24604e+03}
f[ 1]={ 0.00000e+00, 0.00000e+00, 1.24604e+03}
.
.
.
___________________________________________________________________________
energy.xvg:
# This file was created Tue Jun 26 22:00:38 2012
# by the following command:
# g_energy_d -s -dp -quiet
#
# g_energy_d is part of G R O M A C S:
#
# Go Rough, Oppose Many Angry Chinese Serial killers
#
@ title "Gromacs Energies"
@ xaxis label "Time (ps)"
@ yaxis label "(kJ/mol), (K)"
@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 "Bond"
@ s1 legend "Potential"
@ s2 legend "Kinetic En."
@ s3 legend "Total Energy"
@ s4 legend "Conserved En."
@ s5 legend "Temperature"
0.000000 1.360000000000 1.360000000000 8.033028696195
9.393028696195 9.393028696195 322.048544262812
0.001000 1.274838872373 1.274838872373 0.077195406901
1.352034279274 17.389274589248 3.094806374580
0.002000 1.141621384530 1.141621384530 0.181863160294
1.323484544824 17.374836589202 7.290994249209
0.003000 0.971972108898 0.971972108898 0.325139897720
1.297112006618 17.361628208451 13.035037555885
0.004000 0.780311882060 0.780311882060 0.493462222612
1.273774104672 17.350136510882 19.783172256830
0.005000 0.582583355004 0.582583355004 0.671485683506
1.254069038510 17.340700175901 26.920230842567
0.006000 0.394878246021 0.394878246021 0.843413133674
1.238291379694 17.333484424611 33.812897001156
0.007000 0.232083734770 0.232083734770 0.994330449068
1.226414183837 17.328471077977 39.863255286285
0.008000 0.106661122703 0.106661122703 1.111437369236
1.218098491939 17.325464232219 44.558136207263
0.009000 0.027656180079 0.027656180079 1.185072612382
1.212728792461 17.324110715250 47.510213656278
0.010000 0.000018869457 0.000018869457 1.209452000918
1.209470870375 17.323933267090 48.487596768555
0.011000 0.024282388081 0.024282388081 1.183064285501
1.207346673582 17.324373367123 47.429698725646
0.012000 0.096620173650 0.096620173650 1.108699349473
1.205319523123 17.324839892590 44.448367487097
.
.
.
___________________________________________________________________________
topol.tpr:
inputrec:
integrator = md-vv-avek
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 = V-rescale
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 3
ref_t: 0 0
tau_t: 1 0
anneal: No No
ann_npoints: 0 0
acc: 0 0 0
nfreeze: N N N
energygrp_flags[ 0]: 0 0
energygrp_flags[ 1]: 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 = 2
lambda = 0.000000e+00
topology:
name="Free Sulfur Molecule"
#atoms = 2
molblock (0):
moltype = 0 "SAU"
#molecules = 1
#atoms_mol = 2
#posres_xA = 0
#posres_xB = 0
ffparams:
atnr=1
ntypes=2
functype[0]=LJ_SR, c6= 9.98400640e-03, c12= 1.30754560e-05
functype[1]=BONDS, b0A= 2.38000e-01, cbA= 6.80000e+05, b0B=
2.38000e-01, cbB= 6.80000e+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="SAU"
atoms:
atom (2):
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]={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 (2):
atom[0]={name="S"}
atom[1]={name="S"}
type (2):
type[0]={name="S",nameB="S"}
type[1]={name="S",nameB="S"}
residue (1):
residue[0]={name="SAU", nr=1, ic=' '}
cgs:
nr=2
cgs[0]={0..0}
cgs[1]={1..1}
excls:
nr=2
nra=4
excls[0][0..1]={0, 1}
excls[1][2..3]={0, 1}
Bond:
nr: 3
iatoms:
0 type=1 (BONDS) 0 1
grp[T-Coupling ] nr=2, name=[ 0 1]
grp[Energy Mon. ] nr=2, name=[ 0 1]
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 (6):
grpname[0]={name="System"}
grpname[1]={name="Other"}
grpname[2]={name="SAU"}
grpname[3]={name="0"}
grpname[4]={name="1"}
grpname[5]={name="rest"}
groups T-Cou Energ Accel Freez User1 User2 VCM XTC Or. R
QMMM
allocated 2 2 0 0 0 0 0 0
0 0
groupnr[ 0] = 0 0 0 0 0 0 0 0
0 0
groupnr[ 1] = 1 1 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 (2x3):
x[ 0]={ 3.15000e+00, 3.15000e+00, 3.27000e+00}
x[ 1]={ 3.15000e+00, 3.15000e+00, 3.03000e+00}
v (2x3):
v[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
v[ 1]={ 0.00000e+00, 0.00000e+00, 1.00000e+00}
Group statistics
T-Coupling : 1 1 (total 2 atoms)
Energy Mon. : 1 1 (total 2 atoms)
Acceleration: 2 (total 2 atoms)
Freeze : 2 (total 2 atoms)
User1 : 2 (total 2 atoms)
User2 : 2 (total 2 atoms)
VCM : 2 (total 2 atoms)
XTC : 2 (total 2 atoms)
Or. Res. Fit: 2 (total 2 atoms)
QMMM : 2 (total 2 atoms)
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* 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
