[gmx-users] Discontinuity in first time-step velocities for diatomic molecule locally coupled to thermal drain
Berk Hess
gmx3 at hotmail.com
Tue Jul 10 20:03:06 CEST 2012
Hi,
You can't have ret_t=0 with Nose-Hoover.
We should add a check for this.
Cheers,
Berk
> Date: Thu, 28 Jun 2012 11:15:15 +0300
> From: InonShar at TAU.ac.IL
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Discontinuity in first time-step velocities for diatomic molecule locally coupled to thermal drain
>
> 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 at 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-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
More information about the gromacs.org_gmx-users
mailing list