VMD-L Mailing List

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Wed Oct 12 2016 - 11:28:40 CDT

Hi Qasim,

Its not one script, and it might need some edits, since I couldn't test
it (for instance, looking back, I see that the colvar definition is
incomplete, since I never told the colvars module to actually USE the
orientationAngle). Basically you'd copy the colvars {} piece into its
own file, and the next two pieces can be inserted into a fresh .tcl
file, which you can source from the tkconsole. From there, it *should*
just work on a trajectory, so long as frame 0 is the reference you want.

-Josh

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

On 10/12/2016 10:19 AM, Qasim Pars wrote:
> Hi Joshua,
>
> Thanks for your reply. Could you please attach the script? I couldn't figure out how to run the script? I think I need only a trajectory file (e.g. xtc file)? The reference file is not necessary, right?
>
> Thanks in advance
>
>
>
>> On 12 Oct 2016, at 18:36, "Vermaas, Joshua" <Joshua.Vermaas_at_nrel.gov> wrote:
>>
>> 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
>>
>