VMD-L Mailing List

From: Sara baretller (sarabiocomputation_at_gmail.com)
Date: Thu Jul 07 2011 - 16:13:03 CDT

Hi Mark;

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

Sara

On Thu, Jul 7, 2011 at 10:22 AM, Mark Cunningham <cunningham_at_utpa.edu>wrote:

> 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>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 [owner-vmd-l_at_ks.uiuc.edu] on behalf of
>> Sara baretller [sarabiocomputation_at_gmail.com]
>> *Sent:* Wednesday, July 06, 2011 10:04 AM
>> *To:* 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> 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
>>>
>>>
>>
>