# How to calculate the temperature for a subset of my simulation

From: Jhonatam Cordeiro Rodrigues (jcrodrig_at_aggies.ncat.edu)
Date: Mon Oct 19 2015 - 15:08:27 CDT

Dear NAMD users

I am running simulations that have different temperature couplings on
different subsets of the simulation. Please, how can I read the temperature
for a given selection of atoms?

I found this nice script written by Shailesh Pandey, but when I run it, I
get temperature readings that are certainly incorrect (over 1Bi K).

Is there something wrong with the script? or is there a script out there
that would do the trick?

########### TCL script to calculate temperature from velocity
trajectorySTART #########

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

proc getTemp { frame rigidbonds } {

set all [atomselect top all frame \$frame]
set eng {};
set esum 0.0;
set nDoF 0;
set c1 [[ atomselect top "((type 'HW' or hydrogen) or water)" ] get
index ]
set c2 [[ atomselect top "water" ] get index ]

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?
if { (\$rigidbonds == 2) && (\$id in \$c1) } {
set DoF 2
} elseif { (\$rigidbonds == 1) && (\$id in \$c2) } {
set DoF 2
}

set nDoF [expr {\$nDoF + \$DoF}]
}
\$all delete
puts "nDoF: \$nDoF, esum: \$esum"
set kB 0.00198657 ;# kcal/(mol.K)
set T [expr { \$esum * 2.0 / (\$nDoF * \$kB) }]
puts "Frame: \$frame ## Temp: \$T"
}

########### TCL script to calculate temperature from velocity trajectory
END #########

After the running the script for water I get results like this:
nDoF: 481968, esum: 638805356005.8174
Frame: 201 ## Temp: 1334370608.0432146

Sincerely
Jhonatam Cordeiro,
Ph.D. candidate
Department of Industrial & Systems Engineering
North Carolina A&T State University
1601, East Market st.
Greensboro, NC 27411
Phone: +1 336 706-9203
Fax: +1 336 334-7729
email: jcrodrig_at_aggies.nc <http://goog_74760233>at.edu

This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:21:25 CST