VMD-L Mailing List

From: Ashar Malik (asharjm_at_gmail.com)
Date: Wed Sep 20 2023 - 04:13:25 CDT

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