[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