[gmx-users] Different temperatures for different groups, even with Nose-Hoover

Justin A. Lemkul jalemkul at vt.edu
Thu Oct 29 22:38:48 CET 2009



Michael Lerner wrote:
> Hi,
> 
> I'm using GROMACS 4.0.5 to do a MARTINI simulation of a CG protein in a 
> CG water box, looking at the diffusion constant of the protein among 
> other things. I'm using the Nose-Hoover thermostat and Parrinello-Rahman 
> barostat with tc-grps = System. In CHARMM, a protein-in-water simulation 
> with an extended ensemble will show equal temperatures for the protein 
> and the water. However, I'm finding that, with ref_t = 323, I get a SOL 
> (both W and WF particles) temperature of 323 and a Protein temperature 
> of 295. Is this what I should expect?
> 
> The system has 62.5k particles, the box is 20 x 20 x 20, and the 
> relevant .mdp parameters I've used are:
> 
> tcoupl                   = Nose-Hoover
> tc-grps                  = System
> tau_t                    = 2
> ref_t                    = 323
> Pcoupl                   = Parrinello-Rahman
> Pcoupltype               = isotropic
> tau_p                    = 5
> compressibility          = 5e-5
> ref_p                    = 1.0
> 
> These parameters seem to lead to stable systems for bilayer simulations. 
> I've included the whole .mdp file at the end in case I've forgotten to 
> mention something relevant.
> 
> I calculated temperatures with
> 
> g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5
> 

I see you are using "constraints = none," but are there any constraints defined 
in the topology?  If so, as noted in g_traj -h:

"Option -ot plots the temperature of each group, provided velocities are
present in the trajectory file. No corrections are made for constrained
degrees of freedom! This implies -com."

So you will get an inherently incorrect answer when using g_traj, if there are 
any constraints.

-Justin

> Thanks,
> 
> -Michael
> 
> ------ begin md2.mdp ------
> 
> ; RUN CONTROL PARAMETERS =
> integrator               = md
> ; start time and timestep in ps =
> tinit                    = 0.0
> dt                       = 0.020
> nsteps                   = 25500000
> ; mode for center of mass motion removal
> comm-mode                = linear
> ; number of steps for center of mass motion removal
> nstcomm                  = 1
> ; group(s) for center of mass motion removal
> comm-grps                = System
> 
> 
> ; OUTPUT CONTROL OPTIONS =
> ; Output frequency for coords (x), velocities (v) and forces (f) =
> nstxout                  = 50000
> nstvout                  = 50000
> nstfout                  = 0
> ; Output frequency for energies to log file and energy file =
> nstlog                   = 2500
> nstenergy                = 2500
> ; Output frequency and precision for xtc file =
> nstxtcout                = 2500
> xtc_precision            = 100
> ; This selects the subset of atoms for the xtc file. You can =
> ; select multiple groups. By default all atoms will be written. =
> xtc-grps                 =
> ; Selection of energy groups =
> energygrps               = System Protein W WF
> 
> ; NEIGHBORSEARCHING PARAMETERS =
> ; nblist update frequency =
> nstlist                  = 5
> ; ns algorithm (simple or grid) =
> ns_type                  = grid
> ; Periodic boundary conditions: xyz or none =
> pbc                      = xyz
> ; nblist cut-off         =
> rlist                    = 1.2
> 
> ; OPTIONS FOR ELECTROSTATICS AND VDW =
> ; Method for doing electrostatics =
> coulombtype              = Shift
> rcoulomb_switch          = 0.0
> rcoulomb                 = 1.2
> ; Dielectric constant (DC) for cut-off or DC of reaction field =
> epsilon_r                = 15
> ; Method for doing Van der Waals =
> vdw_type                 = Shift
> ; cut-off lengths        =
> rvdw_switch              = 0.9
> rvdw                     = 1.2
> ; Apply long range dispersion corrections for Energy and Pressure =
> DispCorr                 = No
> ; Spacing for the PME/PPPM FFT grid =
> fourierspacing           = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used =
> fourier_nx               = 10
> fourier_ny               = 10
> fourier_nz               = 10
> ; EWALD/PME/PPPM parameters =
> pme_order                = 4
> ewald_rtol               = 1e-05
> epsilon_surface          = 0
> optimize_fft             = no
> 
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS =
> ; Temperature coupling   =
> tcoupl                   = Nose-Hoover
> ; Groups to couple separately =
> tc-grps                  = System
> ; Time constant (ps) and reference temperature (K) =
> tau_t                    = 2
> ref_t                    = 323
> ; Pressure coupling      =
> Pcoupl                   = Parrinello-Rahman
> Pcoupltype               = isotropic
> ; Time constant (ps), compressibility (1/bar) and reference P (bar) =
> tau_p                    = 5
> compressibility          = 5e-5
> ref_p                    = 1.0
> 
> ; SIMULATED ANNEALING CONTROL =
> annealing                = no
> 
> ; GENERATE VELOCITIES FOR STARTUP RUN =
> continuation      = yes
> gen_vel                  = no
> ;gen_temp                 = 323
> ;gen_seed                 = 473529
> 
> ; OPTIONS FOR BONDS     =
> constraints              = none
> ; Type of constraint algorithm =
> constraint_algorithm     = Lincs
> ; Relative tolerance of shake =
> shake_tol                = 0.0001
> ; Highest order in the expansion of the constraint coupling matrix =
> lincs_order              = 4
> ; Lincs will write a warning to the stderr if in one step a bond =
> ; rotates over more degrees than =
> lincs_warnangle          = 30
> ; Convert harmonic bonds to morse potentials =
> morse                    = no
> 
> ; NMR refinement stuff  =
> ; Distance restraints type: No, Simple or Ensemble =
> disre                    = No
> ; Force weighting of pairs in one distance restraint: Equal or 
> Conservative =
> disre_weighting          = Equal
> ; Use sqrt of the time averaged times the instantaneous violation =
> disre_mixed              = no
> disre_fc                 = 1000
> disre_tau                = 1.25
> ; Output frequency for pair distances to energy file =
> nstdisreout              = 100
> ------- end md2.mdp -------
> 
> 
> 
> -- 
> Michael Lerner, Ph.D.
> IRTA Postdoctoral Fellow
> Laboratory of Computational Biology NIH/NHLBI
> 5635 Fishers Lane, Room T909, MSC 9314
> Rockville, MD 20852 (UPS/FedEx/Reality)
> Bethesda MD 20892-9314 (USPS)
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================



More information about the gromacs.org_gmx-users mailing list