VMD-L Mailing List

From: Nuno R. L. Ferreira (nunolf_at_ci.uc.pt)
Date: Fri Feb 17 2006 - 03:50:45 CST

Hi Ana

Send your script to the mailing-list. It's easier to debugg the problem with
the tcl code ;-)

But if you want to calculate for 3 different selections the RMSD, using the
same procedure, instead of doing a foreach, call the proc 3 times with the
desired selection, something like that:

# Calculates the Root-Mean-Squared-Deviation for $selection

proc calcRmsd {{selection} {filename}} {

set sel0 [atomselect top $selection frame 0]

if {[$sel0 num] <= 0} {

error "calcRmsd: must have at least one atom in selection"

}

if {[file exists $filename]} {

error "calcRmsd: file $filename already exists!"

}

set numfrms [molinfo top get numframes]

set nframe 0

set fileId [open ./$filename w 0640]

while {$nframe<$numfrms} {

set sel1 [atomselect top $selection frame $nframe]

set tm [measure fit $sel1 $sel0]

set move_sel [atomselect top "all" frame $nframe]

$move_sel move $tm

set rmsd [measure rmsd $sel0 $sel1]

puts $fileId "$nframe \t $rmsd"

incr nframe

}

close $fileId

}

# start properties calculation
calcRmsd "protein" sel_protein.rmsd
calcRmsd "backbone" sel_backbone.rmsd
calcRmsd "protein sheet backbone" sel_prot_sheet_back.rmsd

Best,
Nuno

P.S. About the "protein sheet backbone" selection, one note. Stride will
calculate which residues are in a sheet secondary structure for frame 0, and
the RMSD will be calculated on those residues. There's no SS recalculation
during the proc.

----- Original Message -----
From: "Ana Celia Araujo Vila Verde" <avilaverde_at_engr.psu.edu>
To: <vmd-l_at_ks.uiuc.edu>
Sent: Thursday, February 16, 2006 9:49 PM
Subject: RE: vmd-l: problem with script to get RMSD

>
> Hi John,
>
> Thanks for replying. I don't think I was very clear on my previous email.
> Let me see if I can explain my problem better:
> The procedure has a foreach loop, looping through the elements in variable
$selection.
>
> When $selection is set to "sheet protein backbone" (lets call this mode
A), then in the first loop it will calculate, for all frames, the RMSD of
all the atoms belonging to sheets relative to those same atoms in frame; it
will output those RMSD values into a file. In the second loop it does the
same thing but for all the atoms in the protein, and will output to a
different file. In the third time it does the same thing for all the atoms
belonging to the backbone.
>
> If I set selection to just "protein" (lets call this mode B), the foreach
loop will only take place once, and will calculate the RMSD of all the atoms
belonging to the protein relative to those same atoms in frame 0 and it
will output those into a file.
>
> The problem is that, even though I didn't make any other changes in the
tcl script, the RMSD for the protein atoms calculated using mode A is
different from the RMSD that I get using mode B, even though I'm doing
exactly the same thing in both! It's like information from the first loop
somehow passes into the second loop, even though I unset all the variables
at the end of each loop.
>
> I hope I explained my problem better.
>
> Thanks,
>
> Ana
>
>
>
> _________________________________
> Ana Célia Araújo Vila Verde
> Penn State University
> Department of Chemical Engineering
> Fenske Laboratory
> University Park, PA 16802
> USA
>
> Phone: +(1) (814) 863-2879
> Fax: +(1) (814) 865-7846
> avilaverde_at_engr.psu.edu
> _________________________________
>
>
> -----Original Message-----
> From: John Stone [mailto:johns_at_ks.uiuc.edu]
> Sent: Thursday, February 16, 2006 2:41 PM
> To: Ana Celia Araujo Vila Verde
> Cc: vmd-l_at_ks.uiuc.edu
> Subject: Re: vmd-l: problem with script to get RMSD
>
>
> Ana,
> It would be helpful to see exactly how you've changed the scripts so
there's
> no confusion on what you're actually comparing. You didn't say anything
> about how different the RMSD values are, but they are bound to be somewhat
> different just because you have different atoms in each selection, and
> different numbers of atoms, so you'll definitely get different results,
> it's just a question of how different... One thing that may help is to
> use the same selections you're using for the RMSD calculations as
> graphical representations to highlight the particular selections you're
> interested in. Then, in order to see more clearly, you may want to use
> the Trajectory tab "Draw Multiple Frames" option to show the timesteps
> that you're comparing superimposed on each other. In order to see them
> more easily, you can then change the coloring method to "Timestep".
> Give that a try and see if it helps make more sense of the numbers you're
> getting from the RMSD calculation.
>
> John
>
> On Thu, Feb 16, 2006 at 02:31:28PM -0500, Ana Celia Araujo Vila Verde
wrote:
> > Dear all,
> >
> >
> >
> > Based on the script to obtain the RMSD of one molecule relative to the
same molecule at time 0, which can be found on the VMD manual, I built a
very similar one to calculate the RMSD just for certain selections of the
protein (see below signature). My problem is that if I do
> >
> >
> >
> > selection "sheet protein backbone", the RMSD for sheet is OK but for
protein and for backbone is wrong (I think).
> >
> >
> >
> > If I take the exact same procedure and simply change to
> >
> >
> >
> > selection " protein "
> >
> >
> >
> > and then to
> >
> >
> >
> > selection "backbone ",
> >
> >
> >
> > I get different values for their respective RMSD! This leads me to think
that the RMSD I get for protein and for backbone when I use """ selection
"sheet protein backbone" """ is incorrect.
> >
> >
> >
> > Could anyone shine some light on this? Clearly I'm doing something
wrong, but I looked at the code for hours and could not find it.
> >
> >
> >
> > Thanks Ana
> >
> >
> >
> > _________________________________
> >
> > Ana Célia Araújo Vila Verde
> >
> > Penn State University
> >
> > Department of Chemical Engineering
> >
> > Fenske Laboratory
> > University Park, PA 16802
> >
> > USA
> >
> >
> >
> > Phone: +(1) (814) 863-2879
> > Fax: +(1) (814) 865-7846
> >
> > avilaverde_at_engr.psu.edu
> >
> > _________________________________
> >
> >
> >
> > To execute I type vmd -dispdev text -e nameOfFile.tcl on the unix
command line.
> >
> >
> >
> > proc print_rmsd0_through_time {{mol top}} {
> >
> >
> >
> > ##################################
> >
> > set selection "sheet protein backbone" ;# the RMSD calculation
is done for the elements of this list
> >
> > #################################
> >
> > set refFrame 0 ;# the RMDS is calculated against this reference
> >
> > set incrm 1; # the comparison is done every incrm number of
frames
> >
> > set durTimestep 0.000002; # duration of the timestep in the
simulation, in nanoseconds
> >
> > set stepsFrame 500; # number of timesteps in each frame in *dcd
fil
> >
> >
#--------------------------------------------------------------------------- 

--
> >
> >
> >
> >         foreach sel $selection {
> >
> >                 set outFile try21RMSD_${sel}_refFrame${refFrame}.txt; #
result goes to this file
> >
> >                 set out [open $outFile w]
> >
> >                 puts $out "computing RMSD for all sheets using frame
$refFrame as reference"
> >
> >                 puts $out "time (ns)        RMSD (angstrom)"
> >
> >
> >
> >                 set reference [atomselect $mol "protein and $sel" frame
$refFrame]
> >
> >                 # the frame being compared
> >
> >                 set compare [atomselect $mol "protein and $sel"]
> >
> >
> >
> >                 set num_steps [molinfo $mol get numframes]
> >
> > # Does comparison every incrm 30 frames
> >
> >                 for {set frme $refFrame} {$frme <= $num_steps} {incr
frme $incrm} {
> >
> >                         # get the correct frame
> >
> >                         $compare frame $frme
> >
> >
> >
> >                         # compute the transformation
> >
> >                         set trans_mat [measure fit $compare $reference]
> >
> >                         # do the alignment
> >
> >                         $compare move $trans_mat
> >
> >                         # compute the RMSD
> >
> >                         set rmsd [measure rmsd $compare $reference]
> >
> >                         # print the RMSD
> >
> > # puts "RMSD of $frame is $rmsd"
> >
> >                         set time [expr {$frme*$stepsFrame*$durTimestep}]
> >
> >                         puts $out "$time        $rmsd"
> >
> >                 }
> >
> >         close $out
> >
> >         $compare delete
> >
> >         unset trans_mat
> >
> >         unset rmsd
> >
> >         unset reference
> >
> >         unset num_steps
> >
> >         }
> >
> > }
> >
> > # Call procedure
> >
> > print_rmsd0_through_time
> >
> > #
> >
> > quit
> >
>
> -- 
> NIH Resource for Macromolecular Modeling and Bioinformatics
> Beckman Institute for Advanced Science and Technology
> University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
> Email: johns_at_ks.uiuc.edu                 Phone: 217-244-3349
>   WWW: http://www.ks.uiuc.edu/~johns/      Fax: 217-244-6078
>
>
>
> ---
> avast! Antivirus: Inbound message clean.
> Virus Database (VPS): 0607-2, 16-02-2006
> Tested on: 17-02-2006 9:18:05
> avast! - copyright (c) 1988-2005 ALWIL Software.
> http://www.avast.com
>
>
>
---
avast! Antivirus: Outbound message clean.
Virus Database (VPS): 0607-2, 16-02-2006
Tested on: 17-02-2006 9:50:47
avast! - copyright (c) 1988-2005 ALWIL Software.
http://www.avast.com