VMD-L Mailing List
From: Ana Celia Araujo Vila Verde (avilaverde_at_engr.psu.edu)
Date: Fri Feb 17 2006 - 08:41:20 CST
- Next message: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Previous message: Axel Kohlmeyer: "Re: vmd 1.8.b does not start up"
- Maybe in reply to: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Next in thread: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
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
_________________________________
-----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
> _________________________________
>
>
> -----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/>
- Next message: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Previous message: Axel Kohlmeyer: "Re: vmd 1.8.b does not start up"
- Maybe in reply to: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Next in thread: Ana Celia Araujo Vila Verde: "RE: problem with script to get RMSD"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
