VMD-L Mailing List

From: Vermaas, Josh (vermaasj_at_msu.edu)
Date: Wed Sep 20 2023 - 06:54:23 CDT

Definitely in radians. acos -1 gives you pi, not 180, but I think that is accounted for by "conv". I think this is coming from your definitions for p1 and p2. As written, p1 is the vector from residue 1 to 2, and p2 is the vector from residue 2 to 3. That will usually be an angle > 90 degrees, and so the cross product won't be as well behaved as if p1 and p2 are made by vectors that share an origin. To be honest, as written, the answer should be 0, not 70, since residue 1 and 2 are used to define the vector AND the plane, so both 10 and 70 are clearly wrong here.

-Josh

On 9/20/23 05:13, Ashar Malik wrote:
Not sure I recall this correctly, but isn’t the angle supposed to be handled in radians? Might be wrong but woke a look.



On Wed, 20 Sep 2023 at 6:37 PM, Castro Martinez, Camila <camila.castro_at_tu-dortmund.de<mailto:camila.castro_at_tu-dortmund.de>> wrote:
Dear all


I'm looking to calculate the angle between a vector and a plane over a trajectory using tcl script. The code I've been using is the following:


set outfile [open output.dat w]
set nf [molinfo top get numframes]

#for the vector
set v1 [atomselect top "residue 1"]
set v2 [atomselect top "residue 2"]

#for the plane
set n1 [atomselect top "residue 1"]
set n2 [atomselect top "residue 2"]
set n3 [atomselect top "residue 3"]

for { set i 0 } { $i < $nf } { incr i } {

  $v1 frame $i
  $v2 frame $i
  set com1 [measure center $v1]
  set com2 [measure center $v2]
  set m2 [vecsub $com2 $com1]
  set m2mag [veclength $m2]

  $n1 frame $i
  $n2 frame $i
  $n3 frame $i
  set com3 [measure center $n1]
  set com4 [measure center $n2]
  set com5 [measure center $n3]
  set p1 [vecsub $com4 $com3]
  set p2 [vecsub $com5 $com4]
  set cn [veccross $p1 $p2]
  set cnmag [veclength $cn]

  set dotp [vecdot $m2 $cn]
  set theta [expr $conv * acos($dotp / ($m2mag * $cnmag))]
  if { $theta > 90.0 } {
    set theta [expr 180.0 - $theta]
  }

  puts $outfile "[expr ($i + 1)] $theta"
}

close $outfile

This commands actually worked, but the results are not what I expected. When I checked in Chimera for example the angle should be 70°, and the angle calculated with this script was 10°. This is clearly not correct. Where have I gone wrong in the code? or am I missing something? Thanks in advance


--
Josh Vermaas
vermaasj_at_msu.edu<mailto:vermaasj_at_msu.edu>
Assistant Professor, Plant Research Laboratory and Biochemistry and Molecular Biology
Michigan State University
vermaaslab.github.io