# AW: Regarding Temperature calculation from velocity trajectory of NAMD-2.9

From: Norman Geist (norman.geist_at_uni-greifswald.de)
Date: Fri Jul 11 2014 - 02:51:44 CDT

Seriously? A contrained hydrogen does only have 2 not 3 DOF. And for TIP3P water each atom has only 2DOF.

Therefore:

#rigidbonds none=0; water=1; all=2

proc getTemp { frame rigidbonds } {
set all [atomselect top all frame \$frame]
set eng {}; set esum 0; set nDoF 0;
foreach m [\$all get mass] v [\$all get {x y z}] id [\$all get index] {
set e [expr {0.5 * \$m * [vecdot \$v \$v]}]
lappend eng \$e
set esum [expr {\$esum + \$e}]

set all [atomselect top all frame \$frame]

set DoF 3.

#SHAKE?

set c1 [atomselect top "index \${id} and ((type 'HW' or hydrogen) or water)" ]

set c2 [atomselect top "index \${id} and water" ]

if { \$rigidbonds == 2 && [\$c1 num] > 0 } {

set DoF 2.

} elseif { \$rigidbonds == 1 && [\$c2 num] > 0 } {

set DoF 2.

}

\$c1 delete

\$c2 delete

incr nDOF \$DoF
}

set T [expr { \$esum * 2.0 / (\$nDoF * 0.00198657) }]
puts "frame: \$frame ## Temp: \$T"
}

Von: Shailesh Pandey [mailto:shaileshp51_at_gmail.com]
Gesendet: Donnerstag, 10. Juli 2014 16:30
An: Norman Geist
Betreff: Re: namd-l: Regarding Temperature calculation from velocity trajectory of NAMD-2.9

Yes, I am using rigidBonds water, how reduction in D.F is to be taken into account?

On Thu, Jul 10, 2014 at 3:38 PM, Norman Geist <norman.geist_at_uni-greifswald.de> wrote:

You didn’t consider the reduced degrees of freedom for constrained hydrogen due SHAKE aka rigidbonds ;)

Norman Geist.

Von: owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] Im Auftrag von Shailesh Pandey
Gesendet: Donnerstag, 10. Juli 2014 11:30
An: namd-l_at_ks.uiuc.edu
Betreff: namd-l: Regarding Temperature calculation from velocity trajectory of NAMD-2.9

Dear NAMD users,

I have performed a MD simulation of a protein+ligand system with explicit solvent and PBC with NPT ensemble.

Now as one of the tests of equilibration, I want to look at temperature profile of Solvated system, protein+ligand, water, and center of mass of complex. As discussed in

Gallo, M. T., Grant, B. J., Teodoro, M. L., Melton, J., Cieplak, P., Phillips, G. N., & Stec, B. (2009). Novel procedure for thermal equilibration in molecular dynamics simulation. Molecular Simulation, 35(June 2014), 349–357. doi:10.1080/08927020802647272

And, we can calculate temperatures of different configurations if we have corresponding velocity trajectory.

I am using following tcl code to calculate temperature after loading psf file and velocity trajectory as molecule data

proc getTemp { frame } {
set all [atomselect top all frame \$frame]
set eng {}; set esum 0;
foreach m [\$all get mass] v [\$all get {x y z}] {
set e [expr {0.5 * \$m * [vecdot \$v \$v]}]
lappend eng \$e
set esum [expr {\$esum + \$e}]
}
set T [expr { \$esum * 2.0 / (3.0 * [\$all num] * 0.00198657) }]
puts "frame: \$frame ## Temp: \$T"
}

I am following it from http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node14.html

It produces temperature in range of 220-230 for solvated system, but log file shows temperature fluctuating around 310, as 310K was the simulation temperature.

I am unable to find where I am doing wrong?

Any help in this regard is appreciated.

Thank you.

_____

Diese E-Mail ist frei von Viren und Malware, denn der avast! Antivirus <http://www.avast.com/> Schutz ist aktiv.

```---
Diese E-Mail ist frei von Viren und Malware, denn der avast! Antivirus Schutz ist aktiv.
http://www.avast.com
```

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:35 CST