Fwd: help to improve the script to predict mean square displacement

From: jampani srinivas (jampanis_at_gmail.com)
Date: Fri Mar 12 2010 - 14:51:55 CST

---------- Forwarded message ----------
From: jampani srinivas <jampanis_at_gmail.com>
Date: Thu, Mar 11, 2010 at 8:18 PM
Subject: help to improve the script to predict mean square displacement
To: vmd-l_at_ks.uiuc.edu

Dear VMD users,

I have written a small script to predict the mean square displacement in
Tcl. The script is working fine but it is taking huge computational time, it
took around three hours to finish 200 frames with 6000 water molecules. Can
anybody help me to improve the script in such way that it
takes reasonably less time to calculate MSD.

Here is my script

proc MSD {{file ""} {mol top}} {
        # open file if given
        set file_given 0
        if {$file != ""} \
        {
                set fileId [open $file w]
                set file_given 1
        }
set num_steps [molinfo top get numframes]
#set frozen [[atomselect top "beta 1.0" frame 0] get index]
set mobile1 [atomselect top "beta 0"]
set mobile2 [atomselect top "beta 0"]
  for {set i 0} {$i < $num_steps} {incr i } {
    $mobile1 frame $i
    set coord1 [$mobile1 get {x y z}]
    set list2 []
                         for {set j [expr $i+1]} {$j < $num_steps} {incr j}
{
                         set z 0
                         $mobile2 frame $j
                         set coord2 [$mobile2 get {x y z}]
                              for {set l 0} {$l <= [$mobile1 num]} {incr l}
{
                                     set dist1 [vecdist [lindex $coord1 $l]
[lindex $coord2 $l]]
                                     set dist2 [expr $dist1*$dist1]
                                     set z [expr $dist2 + $z]
                                     }
                           set conv [expr $z/[$mobile1 num]]
                           lappend list2 "$conv"
                          }
     for {set k 0} {$k < $i} {incr k} {
     lappend list2 [expr 1.0-1.0]
     }
           set final $list2
                  if {$k == 0} {
                  set avrfinal $final
                      } else {
                          set avrfinal [vecadd $avrfinal $final]
                      }
  }
             # print to file
                    if {$file_given == 1} \
                          {
                                for {set m 0} {$m < [llength $avrfinal]}
{incr m} {
                                set last [expr [expr 1.0/[expr [llength
$avrfinal] - $m]]*[lindex $avrfinal $m]]
                                puts $fileId "$last"
                                }
                           }
              #close file
                       if {$file_given == 1} \
                          {
                            close $fileId
                          }
}

Thanks in advance
Srinivas.

-- 
*********************************************
J. Srinivasa Rao
Post-doctoral Research Associate
C/o Prof. Luis R Cruz Cruz
Computational Biophysics Group
Department of Physics
Drexel University
3141 Chestnut St
Philadelphia, PA 19104
Ph:  Off:     215-895-1989
       Mob:  704-706-4191
Web:http://jsrao.web.officelive.com/default.aspx
**********************************************
-- 
*********************************************
J. Srinivasa Rao
Post-doctoral Research Associate
C/o Prof. Luis R Cruz Cruz
Computational Biophysics Group
Department of Physics
Drexel University
3141 Chestnut St
Philadelphia, PA 19104
Ph:  Off:     215-895-1989
       Mob:  704-706-4191
Web:http://jsrao.web.officelive.com/default.aspx
**********************************************

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:53 CST