VMD-L Mailing List

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Wed Oct 12 2016 - 10:36:21 CDT

Hi Qasim,

Have you checked the arrows drawn graphically? My recollection is that while orient will correctly give the principle axes, it isn't consistent about its reference frame, so while one arrow representing the primary principle axis might be drawn "up", it might be drawn "down" for a related structure, even if the structures are pretty similar. This would mean that you could get unlucky and have your arrows drawn in opposite directions, meaning that the angle between would be 180-178, which is probably pretty reasonable.

A more robust method might be to use the collective variables module built into VMD. OrientationAngle does exactly what you want to do, but is slightly more involved to setup:

First, you'd need a colvars configuration file, with contents similar to this:

colvar {
    name orientation
    atoms {
        atomNameResidueRange CA 90-115
    }
    refPositionsFile reference.pdb
    refPositionsCol B
    refPositionsColValue 1
}

This tells the colvars module that you want to track the orientation of residues 90-115 after fitting the geometry to minimize the RMSD for the piece specified in reference.pdb
Next, you'd need to make a reference.pdb file, which determines what gets aligned during the fitting process:

set all [atomselect top "all" frame 0]
set ref [atomselect top "name CA and resid 2 to 50 and chain A"]
$all set beta 0
$ref set beta 1
$all writepdb reference.pdb

The final step would be to calculate the colvar at each frame of your trajectory:

cv molid top
cv configfile whateveryounamedyourconfigfile.conf
for { set f 0 } { $f < [molinfo top get numframes] } { incr f } {
    cv frame $f
    cv update
    cv printframe
}

That should result in an answer that makes much more sense.

-Josh

On 10/12/2016 09:11 AM, Qasim Pars wrote:
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