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 ]