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

Michael Lerner mglerner+gromacs at gmail.com
Thu Oct 29 23:07:05 CET 2009


On Thu, Oct 29, 2009 at 5:38 PM, Justin A. Lemkul <jalemkul at vt.edu> wrote:

>
> 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.
>
>
Oops! I actually had read some previous messages about that, but I got
confused and thought the problem was with g_energy, not g_traj. I'm not
particularly familiar with the GROMACS source code, but gmx_traj.c has this
function:

static real temp(rvec v[],real mass[],int isize,atom_id index[])
{
  real ekin2=0;
  int  i,j;

  for(i=0; i<isize; i++) {
    j = index[i];
    ekin2 += mass[j]*norm2(v[j]);
  }

  return ekin2/(3*isize*BOLTZ);
}

so it looks like g_traj is just assuming that there are 3*N degrees of
freedom. My protein topology file has 304 particles and 83 lines in the
constraints section, so I think the correction should be

295.3496*((3*304)/(3*304 - 83)) = 324.9202

for my protein vs. 322.9525 for the W's and 322.9279 for the WF's.

So, they're still a little off, but doing some block averaging and looking
at standard errors, it's well within the acceptable range.

Thanks,

-Michael


> -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
>
> ========================================
>
> _______________________________________________
> 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
>



-- 
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)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20091029/aa28330f/attachment.html>


More information about the gromacs.org_gmx-users mailing list