VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Sun Sep 24 2006 - 10:40:06 CDT

On Sun, 24 Sep 2006, Ivan Degtyarenko wrote:

ID>
ID> Dear All, to anyone interested. Since that utilities measuring angles and
ID> dihedrals are still missing from VMD, here is a small dirty Tcl code which
ID> takes coordinates of three points as input and returns angle between them
ID> in degrees. It is not a state of art, but works and might be useful at
ID> some point. Corrections/improvements are welcome.

well, you could make more use of vmd provided primitives.
e.g. with:

proc angle {D H A} {
 # cosine of the angle between three points
 # get Pi
 global M_PI

 # setup vectors hd and ha
 set hd [vecsub $D $H]
 set ha [vecsub $A $H]
 set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] * [veclength $ha])]

 # convert cosine to angle in degrees
 return [expr acos($cosine)*(180.0/$M_PI)]
}

axel.

ID>
ID> YT, Ivan Degtyarenko
ID>
ID>
ID> # compute angle between three arbitrary points in space
ID> # (initially has been written for hydrogen bond angles,
ID> # notation has been kept: D-H...A)
ID> # usage: angle {x1 y1 z1} {x2 y2 z2} {x3 y3 z3}
ID>
ID> proc angle {D H A} {
ID> # cosine of the angle between three points
ID> # cos = ( v1 * v2 ) / |v1| * |v2|, where v1 and v2 are vectors
ID>
ID> # get Pi 3.14159265358979323846
ID> global M_PI
ID>
ID> # initialize arrays
ID> set d() 0; set h() 0; set a() 0; set hd() 0; set ha() 0;
ID>
ID> # split coordinates
ID> set d(x) [lindex $D 0]; set d(y) [lindex $D 1]; set d(z) [lindex $D 2];
ID> set h(x) [lindex $H 0]; set h(y) [lindex $H 1]; set h(z) [lindex $H 2];
ID> set a(x) [lindex $A 0]; set a(y) [lindex $A 1]; set a(z) [lindex $A 2];
ID>
ID> # setup vectors hd and ha
ID> set hd(x) [expr $d(x) - $h(x)];
ID> set hd(y) [expr $d(y) - $h(y)];
ID> set hd(z) [expr $d(z) - $h(z)];
ID>
ID> set ha(x) [expr $a(x) - $h(x)];
ID> set ha(y) [expr $a(y) - $h(y)];
ID> set ha(z) [expr $a(z) - $h(z)];
ID>
ID> # compute cosine
ID> set cosine [expr \
ID> ($hd(x)*$ha(x) + $hd(y)*$ha(y) + $hd(z)*$ha(z)) / \
ID> (sqrt(pow($hd(x),2) + pow($hd(y),2) + pow($hd(z),2)) * \
ID> sqrt(pow($ha(x),2) + pow($ha(y),2) + pow($ha(z),2)))];
ID>
ID> # convert cosine to angle in degrees
ID> set angle [expr acos($cosine)*(180.0/$M_PI)]
ID>
ID> return $angle;
ID> }
ID> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.