VMD-L Mailing List

From: Qasim Pars (qasimpars_at_gmail.com)
Date: Wed Oct 12 2016 - 07:42:28 CDT

Dear users,

I want to calculate the helix (residues 90-115) orientation angle between
the first/reference and final frame of a protein. To do this I aligned the
reference structure to the z-axis... For this aim I have two pdb files and
the below script. The script gives me large angle (178 degree) for helix,
which is not compatible with literature. Could you please tell me what I am
doing wrong in the script?

mol load pdb [lindex $argv 0].pdb
mol load pdb [lindex $argv 1].pdb

proc angle { a b } {
    set amag [veclength $a]
    set bmag [veclength $b]
    set dotprod [vecdot $a $b]
    return [expr 57.2958 * acos($dotprod / ($amag * $bmag))]
}

package require Orient
namespace import Orient::orient

set sel [atomselect 0 "chain A and name CA"]
set I [draw principalaxes $sel]
set A [orient $sel [lindex $I 2] {0 0 1}]
[atomselect 0 all] move $A
set I [draw principalaxes $sel]

#fit the final structure on the residues 2-50 of reference structure
set ref [atomselect 0 "name CA and resid 2 to 50 and protein and chain A"]
set to_move [atomselect 1 "name CA and resid 2 to 50 and protein and chain
A"]
set trans_mat [measure fit $to_move $ref]
[atomselect 1 all] move $trans_mat

set sel_h1 [atomselect 0 "name CA and resid 90 to 115 and chain A"]
set I_hel_1 [draw principalaxes $sel_h1]

set sel_h2 [atomselect 1 "name CA and resid 90 to 115 and chain A"]
set I_hel_2 [draw principalaxes $sel_h2]

set dot1 [vecdot [lindex $I_hel_1 2] [lindex $I 2]]
set hel1 [vecsub [lindex $I_hel_1 2] [vecscale $dot1 [lindex $I 2]]]

set dot2 [vecdot [lindex $I_hel_2 2] [lindex $I 2]]
set hel2 [vecsub [lindex $I_hel_2 2] [vecscale $dot2 [lindex $I 2]]]

set rot [angle $hel1 $hel2]
puts "Orientation angle is $rot"

quit 

-- 
Qasim Pars