VMD-L Mailing List

From: Mark Cunningham (cunningham_at_utpa.edu)
Date: Thu Jul 07 2011 - 09:22:27 CDT

Sara:

In trying to convert from radians to degrees, 180/3.14159 is only correct to
5 decimal places, where double precision arithmetic has approximately 12
decimal places of precision. This may seem unimportant but can lead to
hard-to-trace bugs in more complex calculations. Better practice is to use
the value of pi to machine precision, which can be obtained from calling the
arc cosine function with the argument -1. Also, tcl is an interpreted language,
meaning that every time you go through the loop tcl has to parse the string
180/3.14159 and decide what to do with it. Compiled languages like Fortran
and even C++ will recognize this to be a constant value and will just compute
it once. Pulling this outside the loop will improve your efficiency. As Axel
pointed out in his comments, you need to be careful with what you do inside
the loops. Recomputing the conversion factor 180/pi each time just makes
the code inefficient. Trying to reopen the output file will cause different errors
on different machines: some will simply rewind the output file and start over,
others may complain and die, others may just die. Repeating the atom
selections (without performing a $sel delete) will allocate more and more memory
to the selections, which can also result in your program crashing.

The %5i and %7.2f are format specifiers. Tcl uses a C-like syntax. They mean
an integer with a field width of 5 spaces and a floating point number with a field width of
7 spaces and 2 places after the decimal. This will produce lines like
    5 90.00 80.00 -80.00 120.00 130.00 -45.00
in your output file. You can, of course, change these to suit your needs. This sort of
structure to your output file will make it simpler to be read by some other program for plotting.
The available format specifiers can be found in the tcl documentation.

Mark
________________________________
From: Sara baretller [sarabiocomputation_at_gmail.com]
Sent: Thursday, July 07, 2011 8:08 AM
To: Mark Cunningham
Cc: vmd-l_at_ks.uiuc.edu
Subject: Re: vmd-l: Re: Problem with for loop.

Thank you so much for you help, i really appreciate your time and effort.

I am sorry but i am going to ask you what did you mean by set conv [expr 180/acos(-1.0)] # Using just the digits of pi that you remember is poor practice
also what is the value "%5i" and "%7.2f " are they just examples...thank you again

Sara

On Wed, Jul 6, 2011 at 1:34 PM, Mark Cunningham <cunningham_at_utpa.edu<mailto:cunningham_at_utpa.edu>> wrote:
 Sara:

If the values are not changing, it is undoubtably due to the fact that you are not
incrementing the frame. If i is your frame index, you need a
$sel1 frame $i

As a cautionary note, atomselect is computationally expensive. It would be best
to not recompute them. You can use the tcl array construct:

 for {set j 0} {$j < 6} {incr j} {
    set sel1($j) [atomselect $mol "residue $j and resname LEU and name BBc"]
    set sel2($j) [atomselect $mol "residue $j and resname GLU and name BBc"]
}

Inside your loop over the frames, you can access the stored selections:
set vec2 {0 1 0} # This isn't changing so don't put it inside the loop
set conv [expr 180/acos(-1.0)] # Using just the digits of pi that you remember is poor practice
for {set i 0} {$i < $Nframes} {incr i} {
   puts -nonewline $outfile [format "%5i" $i]

   for {set j 0} {$j < 6}{incr j} {
       $sel1($j) frame $i
       $sel2($j) frame $i
       set coord1 [lindex [$sel1($j) get { x y z }] 0]
       set coord2 [lindex [$sel2($j) get { x y z }] 0]
       set vec1 [vecnorm[vecsub $coord1 $coord2 ]] #vecnorm will not fail if the vector length is 0 so
                                                                          #there's no particular reason to do the checking you do
       set angle [expr $conv*acos([vecdot $vec1 $vec2])]
       puts -nonewline [format "%7.2f " $angle] # this will write all 6 angles on one line
   }
   puts $outfile " " # this will force a newline
}

Good luck

Mark
________________________________
From: owner-vmd-l_at_ks.uiuc.edu<mailto:owner-vmd-l_at_ks.uiuc.edu> [owner-vmd-l_at_ks.uiuc.edu<mailto:owner-vmd-l_at_ks.uiuc.edu>] on behalf of Sara baretller [sarabiocomputation_at_gmail.com<mailto:sarabiocomputation_at_gmail.com>]
Sent: Wednesday, July 06, 2011 10:04 AM
To: vmd-l_at_ks.uiuc.edu<mailto:vmd-l_at_ks.uiuc.edu>
Subject: vmd-l: Re: Problem with for loop.

Hi all

please if you have an idea or suggestion , i would really appreciated it

I solved the first problem , and the script worked and gave me the angle between the six peptides and the vector { 0 1 0} so i added another loop for the frames. for {set i 0} {$i < 8} {incr i} { , 8 as an example.

 it gave me the same results over and over again for the frames.

here the loop is only till frame 8 as an example:

the repeat for all the angle but i highlighted two anglesas an example that keep repeating between from the first frame to the other frame. angle 117.72161220856759 between (peptides 0)'s vector and the vector { 0 1 0} in frame 0. is the same as angle 117.72161220856759 between peptides 1 vector and the vector { 0 1 0} in frame 1 and so on.

 i alos hilighted : 55.63259488419982 for the height frame

Tcl Counsel results:

residue 0 of 6 peptides
 vec1l vec2l : 19.99249894859944 1.0
Final angle : 117.72161220856759
117.72161220856759
0
residue 1 of 6 peptides
 vec1l vec2l : 18.69491870304378 1.0
Final angle : 135.78875014576658
135.78875014576658
0
residue 2 of 6 peptides
 vec1l vec2l : 27.989461118438296 1.0
Final angle : 55.63259488419982
55.63259488419982
0
residue 3 of 6 peptides
 vec1l vec2l : 27.10166056426399 1.0
Final angle : 93.17285877063786
93.17285877063786
0
residue 4 of 6 peptides
 vec1l vec2l : 29.5839470751503 1.0
Final angle : 145.5654673648764
145.5654673648764
0
residue 5 of 6 peptides
 vec1l vec2l : 25.998842961169046 1.0
Final angle : 87.57519131337723
87.57519131337723
0
residue 0 of 6 peptides
 vec1l vec2l : 19.99249894859944 1.0
Final angle : 117.72161220856759
117.72161220856759
1
residue 1 of 6 peptides
 vec1l vec2l : 18.69491870304378 1.0
Final angle : 135.78875014576658
135.78875014576658
1
residue 2 of 6 peptides
 vec1l vec2l : 27.989461118438296 1.0
Final angle : 55.63259488419982
55.63259488419982
1
residue 3 of 6 peptides
 vec1l vec2l : 27.10166056426399 1.0
Final angle : 93.17285877063786
93.17285877063786
1
residue 4 of 6 peptides
 vec1l vec2l : 29.5839470751503 1.0
Final angle : 145.5654673648764
145.5654673648764
1
residue 5 of 6 peptides
 vec1l vec2l : 25.998842961169046 1.0
Final angle : 87.57519131337723
87.57519131337723
1
residue 0 of 6 peptides
 vec1l vec2l : 19.99249894859944 1.0
Final angle : 117.72161220856759
117.72161220856759
2
residue 1 of 6 peptides
 vec1l vec2l : 18.69491870304378 1.0
Final angle : 135.78875014576658
135.78875014576658
2
residue 2 of 6 peptides
 vec1l vec2l : 27.989461118438296 1.0
Final angle : 55.63259488419982
55.63259488419982
2
residue 3 of 6 peptides
 vec1l vec2l : 27.10166056426399 1.0
Final angle : 93.17285877063786
93.17285877063786
2
residue 4 of 6 peptides
 vec1l vec2l : 29.5839470751503 1.0
Final angle : 145.5654673648764
145.5654673648764
2
residue 5 of 6 peptides
 vec1l vec2l : 25.998842961169046 1.0
Final angle : 87.57519131337723
87.57519131337723
2
residue 0 of 6 peptides
 vec1l vec2l : 19.99249894859944 1.0
Final angle : 117.72161220856759
117.72161220856759
3
residue 1 of 6 peptides
 vec1l vec2l : 18.69491870304378 1.0
Final angle : 135.78875014576658
135.78875014576658
3
residue 2 of 6 peptides
 vec1l vec2l : 27.989461118438296 1.0
Final angle : 55.63259488419982
55.63259488419982
3
residue 3 of 6 peptides
 vec1l vec2l : 27.10166056426399 1.0
Final angle : 93.17285877063786
93.17285877063786
3
residue 4 of 6 peptides
 vec1l vec2l : 29.5839470751503 1.0
Final angle : 145.5654673648764
145.5654673648764
3
residue 5 of 6 peptides
 vec1l vec2l : 25.998842961169046 1.0
Final angle : 87.57519131337723
87.57519131337723
3
residue 0 of 6 peptides
 vec1l vec2l : 19.99249894859944 1.0
Final angle : 117.72161220856759
117.72161220856759
4
residue 1 of 6 peptides
 vec1l vec2l : 18.69491870304378 1.0
Final angle : 135.78875014576658
135.78875014576658
4
residue 2 of 6 peptides
 vec1l vec2l : 27.989461118438296 1.0
Final angle : 55.63259488419982
55.63259488419982
4
residue 3 of 6 peptides
 vec1l vec2l : 27.10166056426399 1.0
Final angle : 93.17285877063786
93.17285877063786
4
residue 4 of 6 peptides
 vec1l vec2l : 29.5839470751503 1.0
Final angle : 145.5654673648764
145.5654673648764
4
residue 5 of 6 peptides
 vec1l vec2l : 25.998842961169046 1.0
Final angle : 87.57519131337723
87.57519131337723
4

On Wed, Jul 6, 2011 at 10:38 AM, Sara baretller <sarabiocomputation_at_gmail.com<mailto: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"]

                set mol [molinfo top]

                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

                $sel1 update
                $sel2 update

                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
                }

close $outfile

tcl councel:

angle 1.5284741404344042
Final angle : 87.57519131337723
87.57519131337723

Thank you

Sara