VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Wed Aug 16 2006 - 19:36:28 CDT

On Wed, 16 Aug 2006, Leonardo Sepulveda Durán wrote:

LS> Hello everyone.

hello leonardo,

LS> I am modifying the rmsd-fulltrottle.tcl script to calculate the rmsf
LS> and asociated B-factors of a trayectory. I have a question about Tcl
LS> variables. The following sentence
LS>
LS> set sel_resid [[atomselect top "protein and alpha"] get resid]

first of all, this is actually not good VMD scripting style, as the
atomselect command will (implicitely) create a new tcl function
and if you don't 'delete' them you produce a memory leak, see:
http://biocore.ks.uiuc.edu/biocore/biofs/VMD%20(Public)/tips/tclscripts.html

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

it is a list and not an array. in tcl there is a difference between
those.

LS>
LS> >Main< (Escritorio) 13 % set sel_resid [[atomselect top "protein and
LS> alpha"] get resid]
LS>
LS> 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
LS> 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
LS> 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317
LS>
LS> I was trying to access to the first value of this structure, in this
LS> case 268. If i use expr $sel_resid(1) the program says that
LS> sel_resid isn't an array. Anyway, when I run the script there are

a) tcl counts from 0, so the first element has the index 0, and
b) since you have a list, you address the element with
[lindex $sel_resid 0 ] but note, that
c) a 'get' on a selection usually returns a list of lists, since you
can request multiple properties in one go.

hope that helps,
   axel.

p.s.: i don't get, why you have to use arrays in your loob below.
regular variables should do...

LS> other errors that I dont understand. As I am new in Tcl programing, It
LS> is reasonable that i am making one or more mistakes, nevertheless I
LS> can't grasp it yet. Following is the subroutine with coments, If I
LS> comment all the lines to calculate B_factors or average rmsfs, it
LS> works for ca atoms. I hope someone could help me with this
LS>
LS> # sel_resid is a vector with all the residue numbers
LS> set sel_resid [[atomselect top "protein and alpha"] get resid]
LS>
LS> proc rmsf_residue_over_time {{mol top} res} {
LS>
LS> # Initialize vector for average rmsf of noh atoms
LS> foreach r $res {
LS> set rmsf_ave($r) 0
LS> #set rmsf_ca($r) 0
LS> #set bf_ca($r) 0
LS> }
LS>
LS> # Loop through all residues to calculate rmsfs and put information
LS> # into a file
LS> set outfile [open rmsf.dat w]
LS>
LS> foreach r $res {
LS> # select alfa carbon
LS> set ca_atom [atomselect $mol "protein and resid $r and alpha"]
LS>
LS> # select not hidrogen atoms
LS> set noh_atoms [atomselect $mol "protein and resid $r and noh"]
LS>
LS> # This expresion takes the selection $noh_atoms, which is a list.
LS> # It automatically loops over it, suming iteratively each
LS> # rmsd to rmsf_ave($r), wich is then the sum of all residue atoms rmsfs
LS> set rmsf_ave($r) [expr $rmsf_ave($r) + [measure rmsf $noh_atoms]
LS> ] # <- posible error
LS> set rmsf_ave($r) [expr $rmsf_ave($r)/[llenght $noh_atoms] ]
LS> # <- posible error
LS>
LS> # calculate the rmsf for alpha carbon
LS> set rmsf_ca($r) [measure rmsf $ca_atom]
LS>
LS> # calculate the B-factor of ca atoms from rmsf_ca using the formula
LS> # BF($ca_atom) = 8/3*(pi^2)*($rmsf_ca($r))^2
LS> set M_PI
LS> set bf_ca($r) [expr
LS> 8.0/3.0*$M_PI*$M_PI*$rmsf_ca($r)*$rmsf_ca($r)] # <- posible error
LS>
LS> # Prints all the stuff to the outputfile
LS> # set first $res(1) <- it seems there is an error here
LS> # if ($r == $first){
LS> # puts $outfile [format "%s\t%s\t%s\t%s" "resnum" "rmsf_ca"
LS> "Bf_ca" "rmsf_ave"]
LS> #}
LS> puts $outfile [format "%d\t%.3f\t%.3f\t%.3f" $r $rmsf_ca($r)
LS> $bf_ca($r) $rmsf_ave($r)]
LS> #puts $outfile [format "%d\t%.3f" $r $rmsf_ca($r) ]
LS>
LS> # delete lists
LS> $ca_atom delete
LS> $noh_atom delete
LS> }
LS> close $outfile
LS> }
LS>
LS> rmsf_residue_over_time top $sel_resid
LS>
LS>
LS> Leonardo Sepálveda Durán
LS> Universidad de Chile
LS> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.