VMD-L Mailing List

From: J T (jtibbitt_at_odu.edu)
Date: Tue Feb 24 2009 - 21:37:21 CST

VMDs 'measure dihed' command takes four atom indices as input. You
are trying to give 4 xyz coordinate sets as input instead of atom
indices. I've already written a small procedure script for
calculating the dihedral angle given 4 sets of coordinates. In your
script, you would use it like:

set angle [calc $sel11 $sel21 $sel31 $sel41]

# CALCULATE DIHEDRAL ANGLES
# n1 and n2 = normal vectors formed by atoms 1 2 3 and 2 3 4,
respectively
# n1 = (x1 y1 z1)
# n2 = (x2 y2 z2)
# angle = acos (n1 . n2)/(|n1|*|n2|)
#
# Given 3 points: (x1 y1 z1) (x2 y2 z2) (x3 y3 z3)
# Two vectors are: (x2-x1 y2-y1 z2-z1) (x3-x2 y3-y2 z3-z2)
# Normal Vector = | i j k |
# | x1-x2 y1-y2 z1-z2 |
# | x2-x3 y2-y3 z2-z3 |
#
proc calc {atom1 atom2 atom3 atom4} {
   set 12 [vecsub $atom2 $atom1]
   set 21 [vecsub $atom1 $atom2]
   set 23 [vecsub $atom3 $atom2]
   set 32 [vecsub $atom2 $atom3]
   set 34 [vecsub $atom4 $atom3]
   set n1 [veccross $12 $23]
   set n2 [veccross $23 $34]
   set cosangle [expr [vecdot $n1 $n2] / ([veclength $n1] *
[veclength $n2])]
   set angle [expr {180 * acos($cosangle) / 3.14159265358}]
   set sign [vecdot [veccross $n1 $n2] $23]
   if { $sign < 0 } {
     set angle [expr {0-$angle}]
   } elseif { $sign == 0 } {
     set angle 180.0
   }
   return $angle
}

I'm pretty sure it works well. Hope you can use it.
Jeff Tibbitt

On Feb 24, 2009, at 8:09 PM, Lee, Hwan Kyu (NIH/NHLBI) [E] wrote:

> Dear VMDer,
>
> I calculated center of mass for each residue, and then calculated
> dihedral angles for those four center of mass. I wrote tcl script
> like below, but got an error, "expected integer but got
> "11.965297699" measure dihed: bad atom index". It looks like
> $sel11 cannot be recognized as a coordinate of COM. How can I fix
> this problem?
> Thanks for your help in advance.
>
> -----------------
> set numFrames [molinfo top get numframes]
>
> for {set pp 300} {$pp < $numFrames} {incr pp} {
>
> set n 1
>
> while {$n < 35} {
>
> set n1 [expr {$n+1}]
> set n2 [expr {$n+2}]
> set n3 [expr {$n+3}]
>
> set sel1 [atomselect top "resid $n"]
> set sel2 [atomselect top "resid $n1"]
> set sel3 [atomselect top "resid $n2"]
> set sel4 [atomselect top "resid $n3"]
>
> set sel11 [measure center $sel1 weight mass]
> set sel21 [measure center $sel2 weight mass]
> set sel31 [measure center $sel3 weight mass]
> set sel41 [measure center $sel4 weight mass]
>
> set dih [measure dihed [list $sel11 $sel21 $sel31 $sel41]]
>
> $sel1 delete
> $sel2 delete
> $sel3 delete
> $sel4 delete
> $sel11 delete
> $sel21 delete
> $sel31 delete
> $sel41 delete
>
> incr n
> set result [open "dihedral.xvg" a]
> puts $result "$dih"
> close $result
> }
> }
> -----------------
> best,
> Hwankyu.
>
>