VMD-L Mailing List

From: PARMINDER MANKOO (pmankoo1_at_jhu.edu)
Date: Tue Jun 26 2007 - 15:46:02 CDT

Hello,

 I am trying to write a script to compute the angle bween two domains in my protein for a time trajectory. The two domains have resid (591 to 746) and (747 to 961).

 One way is to compute the center of coordinate for these two domains and use a residue in the middle of two domains as the third point. Following script is what I have written. Is this okay? And how can I do it over time trajectory (the dcd files, I am reading from CHARMM?).

 Thanks, preeti

set DC [atomselect top "resid 591 to 746"]
set AC [atomselect top "resid 747 to 961"]
 set HC [atomselect top "resid 747"]
> >
 puts $cen1 "DC center [measure center $DC]"
 puts $cen2 "AC center [measure center $AC]"
 puts $cen3 "HC center [measure center $HC]"
> >
> >
 proc angle {DC HC AC} {
    # cosine of the angle between three points
     # cos = ( v1 * v2 ) / |v1| * |v2|, where v1 and v2 are vectors
  
     # get Pi 3.14159265358979323846
    global M_PI
> >
     # setup vectors hd and ha
     set hd [vecsub $DC $HC]
     set ha [vecsub $AC $HC]
     set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] * [veclength $ha])]
> >
     # convert cosine to angle in degrees
     return [expr acos($cosine)*(180.0/$M_PI)]
> >
 }