VMD-L Mailing List

From: Leonardo Sepulveda Durán (leonardosepulveda_at_gmail.com)
Date: Wed Aug 16 2006 - 16:41:46 CDT

Hello everyone.

I am modifying the rmsd-fulltrottle.tcl script to calculate the rmsf
and asociated B-factors of a trayectory. I have a question about Tcl
variables. The following sentence

set sel_resid [[atomselect top "protein and alpha"] get resid]

makes sel_resid some sort of list or array, because if I use that
command in the Tk console in vmd it gives the following output:

>Main< (Escritorio) 13 % set sel_resid [[atomselect top "protein and
alpha"] get resid]

268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317

I was trying to access to the first value of this structure, in this
case 268. If i use expr $sel_resid(1) the program says that
sel_resid isn't an array. Anyway, when I run the script there are
other errors that I dont understand. As I am new in Tcl programing, It
is reasonable that i am making one or more mistakes, nevertheless I
can't grasp it yet. Following is the subroutine with coments, If I
comment all the lines to calculate B_factors or average rmsfs, it
works for ca atoms. I hope someone could help me with this

# sel_resid is a vector with all the residue numbers
set sel_resid [[atomselect top "protein and alpha"] get resid]

proc rmsf_residue_over_time {{mol top} res} {

   # Initialize vector for average rmsf of noh atoms
   foreach r $res {
      set rmsf_ave($r) 0
      #set rmsf_ca($r) 0
      #set bf_ca($r) 0
   }

   # Loop through all residues to calculate rmsfs and put information
   # into a file
   set outfile [open rmsf.dat w]

   foreach r $res {
      # select alfa carbon
      set ca_atom [atomselect $mol "protein and resid $r and alpha"]

      # select not hidrogen atoms
      set noh_atoms [atomselect $mol "protein and resid $r and noh"]

      # This expresion takes the selection $noh_atoms, which is a list.
      # It automatically loops over it, suming iteratively each
      # rmsd to rmsf_ave($r), wich is then the sum of all residue atoms rmsfs
      set rmsf_ave($r) [expr $rmsf_ave($r) + [measure rmsf $noh_atoms]
] # <- posible error
      set rmsf_ave($r) [expr $rmsf_ave($r)/[llenght $noh_atoms] ]
  # <- posible error

      # calculate the rmsf for alpha carbon
      set rmsf_ca($r) [measure rmsf $ca_atom]

      # calculate the B-factor of ca atoms from rmsf_ca using the formula
      # BF($ca_atom) = 8/3*(pi^2)*($rmsf_ca($r))^2
      set M_PI
      set bf_ca($r) [expr
8.0/3.0*$M_PI*$M_PI*$rmsf_ca($r)*$rmsf_ca($r)] # <- posible error

      # Prints all the stuff to the outputfile
      # set first $res(1) <- it seems there is an error here
      # if ($r == $first){
      # puts $outfile [format "%s\t%s\t%s\t%s" "resnum" "rmsf_ca"
"Bf_ca" "rmsf_ave"]
      #}
      puts $outfile [format "%d\t%.3f\t%.3f\t%.3f" $r $rmsf_ca($r)
$bf_ca($r) $rmsf_ave($r)]
      #puts $outfile [format "%d\t%.3f" $r $rmsf_ca($r) ]

      # delete lists
      $ca_atom delete
      $noh_atom delete
   }
   close $outfile
}

rmsf_residue_over_time top $sel_resid

Leonardo Sepálveda Durán
Universidad de Chile