From: Ana Celia Araujo Vila Verde (avilaverde_at_engr.psu.edu)
Date: Fri Feb 17 2006 - 08:41:20 CST

Hi Nuno,

 

Thanks for the suggestion. I'll try that and see if the error disappears.

 

About the source code: it was in the first message I sent, below my signature, but here it is again:

 

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

 

 

 

 

_________________________________

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: Nuno R. L. Ferreira [ <mailto:nunolf_at_ci.uc.pt> mailto:nunolf_at_ci.uc.pt]

Sent: Friday, February 17, 2006 4:51 AM

To: vmd mailing list; Ana Celia Araujo Vila Verde

Subject: Re: vmd-l: problem with script to get RMSD

 

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 <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 <http://www.avast.com/>