VMD-L Mailing List

From: Castro Martinez, Camila (camila.castro_at_tu-dortmund.de)
Date: Wed Sep 20 2023 - 02:53:52 CDT

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--_000_092405a8be5346ea8ba47fcc50853f16tudortmundde_--