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

From: Shailesh Pandey (shaileshp51_at_gmail.com)
Date: Fri Jul 11 2014 - 06:14:59 CDT

Thank you very much Norman, for your detailed answer and time. It worked
for me. I have changed your script little bit, update one is as below
Hopefully it might be useful for someone else willing to do something
similar.

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

# 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 #########

On Fri, Jul 11, 2014 at 1:21 PM, Norman Geist <
norman.geist_at_uni-greifswald.de> wrote:

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

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:20:59 CST