[gmx-users] Discontinuity in first time-step velocities for diatomic molecule locally coupled to thermal drain

Inon Sharony InonShar at TAU.ac.IL
Thu Jun 28 10:15:15 CEST 2012


   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)



More information about the gromacs.org_gmx-users mailing list