VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Wed Jul 06 2011 - 10:26:10 CDT

On Wed, Jul 6, 2011 at 10:38 AM, Sara baretller
<sarabiocomputation_at_gmail.com> wrote:
>
> I have been having a problem with this script ,  I will really appreciate
> your help. so the script does get me the angle between two vectors but then
> it will not go to the next loop. it stops in the first. by giving the angle
> between the first vector and the vector {0 1 0 }.
>
>
>
> Script:
>
>  for {set j 0} {$j < 6} {incr j} {
>
>                 set outfile [open "out.dat" "w"]

you open the file _inside_ a loop. i don't think
that is intended.

that will reopen the file in every iteration,
_or_ if you are on an operating system like
windows, will may even refuse to open it for
writing a second time.

>
>                 set mol [molinfo top]

this line can also be hoisted outside the loop.

>                 set sel1 [atomselect $mol "residue $j and resname LEU and
> name BBc"]
>
>                 set sel2 [atomselect $mol "residue $j and resname GLU and
> name BBc"]

>                 $sel1 frame $j
>                 $sel2 frame $j

hmmm.... you have j as variable for residue number
what does this have to do with the frame number?

you most likely want to write a nested loop, where the
"outside" loop would loop over another variable, say 'i'
and then over the number of frames.

>                 $sel1 update
>                 $sel2 update

you selection text is so that the selected atoms won't change over
the trajectory. no need to update the selection again. better use
the frame key word to the atomselect keyword right away.

(even faster would be to make a pair of selection for each of the
'j' residues and store them in an array for the loop over the trajectory
frames and then delete them only at the end. this way you avoid the
time consuming selection update process entirely).

>                 puts "residue $j of 6 peptides"
>                 puts "$sel1 $sel2"
>
>                 set nLEU [llength [$sel1 list]]
>                 set nGLU [llength [$sel2 list]]
>
>                 set coord1  [$sel1 get { x y z }]
>                 set coord2  [$sel2 get { x y z }]
>
>
>                 set vec1 [vecsub [lindex $coord1 0] [lindex $coord2 0]]
>
>                 set  vec2 {0 1 0}
>
>                 set  vec1l [veclength $vec1]
>                 set  vec2l [veclength $vec2]
>
>                 if {$vec1l != 0} {
>                 set vec1 [vecnorm $vec1]
>                 }
>                 if {$vec2l != 0} {
>                 set vec2 [vecnorm $vec2]
>                 }
>
>                 set dotprod [vecdot $vec1 $vec2]
>
>                 set angle [expr acos($dotprod)]
>
>                 set angle [expr ($angle * 180/3.14159)]
>                 puts "Final angle : $angle"
>                 puts $angle

you are at the end of a nested loop in which you
generated atom selections with the atomselect command,
but didn't "delete" them. this will cause a memory leak
and depending on the number of trajectory frames you may
run out of memory and/or VMD will crash or even crash
your machine. this is probably the most common mistake
VMD users make when writing analysis scripts. have a closer
look at the VMD user's guide section on writing trajectory
analysis scripts and also search the mailing list archives
for previous discussions on the subject of memory leaks.
there are plenty and there is a lot to learn from them.

cheers,
   axel.

>                 }
>
> close $outfile
>
> tcl councel:
>
> angle 1.5284741404344042
> Final angle : 87.57519131337723
> 87.57519131337723
>
>
> Thank you
>
> Sara
>
> 

-- 
Dr. Axel Kohlmeyer
akohlmey_at_gmail.com  http://goo.gl/1wk0
Institute for Computational Molecular Science
Temple University, Philadelphia PA, USA.