<br><br><div class="gmail_quote">On Thu, Oct 29, 2009 at 5:38 PM, Justin A. Lemkul <span dir="ltr">&lt;<a href="mailto:jalemkul@vt.edu">jalemkul@vt.edu</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div class="im"><br></div>
I see you are using &quot;constraints = none,&quot; but are there any constraints defined in the topology?  If so, as noted in g_traj -h:<br>
<br>
&quot;Option -ot plots the temperature of each group, provided velocities are<br>
present in the trajectory file. No corrections are made for constrained<br>
degrees of freedom! This implies -com.&quot;<br>
<br>
So you will get an inherently incorrect answer when using g_traj, if there are any constraints.<br>
<br></blockquote><div><br>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&#39;m not particularly familiar with the GROMACS source code, but gmx_traj.c has this function:<br>
<br>static real temp(rvec v[],real mass[],int isize,atom_id index[])<br>{<br>  real ekin2=0;<br>  int  i,j;<br><br>  for(i=0; i&lt;isize; i++) {<br>    j = index[i];<br>    ekin2 += mass[j]*norm2(v[j]);<br>  }<br><br>  return ekin2/(3*isize*BOLTZ);<br>
}<br><br>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 <br><br>295.3496*((3*304)/(3*304 - 83)) = 324.9202<br>
<br>for my protein vs. 322.9525 for the W&#39;s and 322.9279 for the WF&#39;s.<br><br>So, they&#39;re still a little off, but doing some block averaging and looking at standard errors, it&#39;s well within the acceptable range.<br>
<br>Thanks,<br><br>-Michael<br> </div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
-Justin<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div><div></div><div class="h5">
Thanks,<br>
<br>
-Michael<br>
<br>
------ begin md2.mdp ------<br>
<br>
; RUN CONTROL PARAMETERS =<br>
integrator               = md<br>
; start time and timestep in ps =<br>
tinit                    = 0.0<br>
dt                       = 0.020<br>
nsteps                   = 25500000<br>
; mode for center of mass motion removal<br>
comm-mode                = linear<br>
; number of steps for center of mass motion removal<br>
nstcomm                  = 1<br>
; group(s) for center of mass motion removal<br>
comm-grps                = System<br>
<br>
<br>
; OUTPUT CONTROL OPTIONS =<br>
; Output frequency for coords (x), velocities (v) and forces (f) =<br>
nstxout                  = 50000<br>
nstvout                  = 50000<br>
nstfout                  = 0<br>
; Output frequency for energies to log file and energy file =<br>
nstlog                   = 2500<br>
nstenergy                = 2500<br>
; Output frequency and precision for xtc file =<br>
nstxtcout                = 2500<br>
xtc_precision            = 100<br>
; This selects the subset of atoms for the xtc file. You can =<br>
; select multiple groups. By default all atoms will be written. =<br>
xtc-grps                 =<br>
; Selection of energy groups =<br>
energygrps               = System Protein W WF<br>
<br>
; NEIGHBORSEARCHING PARAMETERS =<br>
; nblist update frequency =<br>
nstlist                  = 5<br>
; ns algorithm (simple or grid) =<br>
ns_type                  = grid<br>
; Periodic boundary conditions: xyz or none =<br>
pbc                      = xyz<br>
; nblist cut-off         =<br>
rlist                    = 1.2<br>
<br>
; OPTIONS FOR ELECTROSTATICS AND VDW =<br>
; Method for doing electrostatics =<br>
coulombtype              = Shift<br>
rcoulomb_switch          = 0.0<br>
rcoulomb                 = 1.2<br>
; Dielectric constant (DC) for cut-off or DC of reaction field =<br>
epsilon_r                = 15<br>
; Method for doing Van der Waals =<br>
vdw_type                 = Shift<br>
; cut-off lengths        =<br>
rvdw_switch              = 0.9<br>
rvdw                     = 1.2<br>
; Apply long range dispersion corrections for Energy and Pressure =<br>
DispCorr                 = No<br>
; Spacing for the PME/PPPM FFT grid =<br>
fourierspacing           = 0.12<br>
; FFT grid size, when a value is 0 fourierspacing will be used =<br>
fourier_nx               = 10<br>
fourier_ny               = 10<br>
fourier_nz               = 10<br>
; EWALD/PME/PPPM parameters =<br>
pme_order                = 4<br>
ewald_rtol               = 1e-05<br>
epsilon_surface          = 0<br>
optimize_fft             = no<br>
<br>
; OPTIONS FOR WEAK COUPLING ALGORITHMS =<br>
; Temperature coupling   =<br>
tcoupl                   = Nose-Hoover<br>
; Groups to couple separately =<br>
tc-grps                  = System<br>
; Time constant (ps) and reference temperature (K) =<br>
tau_t                    = 2<br>
ref_t                    = 323<br>
; Pressure coupling      =<br>
Pcoupl                   = Parrinello-Rahman<br>
Pcoupltype               = isotropic<br>
; Time constant (ps), compressibility (1/bar) and reference P (bar) =<br>
tau_p                    = 5<br>
compressibility          = 5e-5<br>
ref_p                    = 1.0<br>
<br>
; SIMULATED ANNEALING CONTROL =<br>
annealing                = no<br>
<br>
; GENERATE VELOCITIES FOR STARTUP RUN =<br>
continuation      = yes<br>
gen_vel                  = no<br>
;gen_temp                 = 323<br>
;gen_seed                 = 473529<br>
<br>
; OPTIONS FOR BONDS     =<br>
constraints              = none<br>
; Type of constraint algorithm =<br>
constraint_algorithm     = Lincs<br>
; Relative tolerance of shake =<br>
shake_tol                = 0.0001<br>
; Highest order in the expansion of the constraint coupling matrix =<br>
lincs_order              = 4<br>
; Lincs will write a warning to the stderr if in one step a bond =<br>
; rotates over more degrees than =<br>
lincs_warnangle          = 30<br>
; Convert harmonic bonds to morse potentials =<br>
morse                    = no<br>
<br>
; NMR refinement stuff  =<br>
; Distance restraints type: No, Simple or Ensemble =<br>
disre                    = No<br>
; Force weighting of pairs in one distance restraint: Equal or Conservative =<br>
disre_weighting          = Equal<br>
; Use sqrt of the time averaged times the instantaneous violation =<br>
disre_mixed              = no<br>
disre_fc                 = 1000<br>
disre_tau                = 1.25<br>
; Output frequency for pair distances to energy file =<br>
nstdisreout              = 100<br>
------- end md2.mdp -------<br>
<br>
<br>
<br>
-- <br>
Michael Lerner, Ph.D.<br>
IRTA Postdoctoral Fellow<br>
Laboratory of Computational Biology NIH/NHLBI<br>
5635 Fishers Lane, Room T909, MSC 9314<br>
Rockville, MD 20852 (UPS/FedEx/Reality)<br>
Bethesda MD 20892-9314 (USPS)<br>
<br>
<br></div></div>
------------------------------------------------------------------------<div class="im"><br>
<br>
_______________________________________________<br>
gmx-users mailing list    <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/search" target="_blank">http://www.gromacs.org/search</a> before posting!<br>
Please don&#39;t post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can&#39;t post? Read <a href="http://www.gromacs.org/mailing_lists/users.php" target="_blank">http://www.gromacs.org/mailing_lists/users.php</a><br>
</div></blockquote>
<br>
-- <br>
========================================<br>
<br>
Justin A. Lemkul<br>
Ph.D. Candidate<br>
ICTAS Doctoral Scholar<br>
Department of Biochemistry<br>
Virginia Tech<br>
Blacksburg, VA<br>
jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> | (540) 231-9080<br>
<a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin</a><br>
<br>
========================================<div><div></div><div class="h5"><br>
_______________________________________________<br>
gmx-users mailing list    <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/search" target="_blank">http://www.gromacs.org/search</a> before posting!<br>
Please don&#39;t post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can&#39;t post? Read <a href="http://www.gromacs.org/mailing_lists/users.php" target="_blank">http://www.gromacs.org/mailing_lists/users.php</a><br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Michael Lerner, Ph.D.<br>IRTA Postdoctoral Fellow<br>Laboratory of Computational Biology NIH/NHLBI<br>5635 Fishers Lane, Room T909, MSC 9314<br>Rockville, MD 20852 (UPS/FedEx/Reality)<br>
Bethesda MD 20892-9314 (USPS)<br>