[gmx-users] inter-helical angles

[Vincent Lemaitre]

Hi Sadhna,

The alternative I found was to write a script to do the job for me. Here 
is the code and some explanations.

Hope it helps,

Vincent

#!/usr/bin/perl

# A script by Vincent Lemaitre called cross_angle.pl which parse the.xvg 
file
# generated by the following serie of commands:
#
#   # 1. [g_bundle] 	Calculate angle between helix 1 & 2
#
#   echo 0 1 | g_bundle -nice 10 -f name.xtc -s name.tpr \
#                       -n HxHy.ndx -dt 50 -b 0 -e 10000 -na 2 -oa 
H1H2.pdb\
#                       >& g_bundle.H1H2.0.out

#   # Note the HxHy.ndx is an index file containing CA of the residues 
belonging to the helices.
#
#   # 2. rm, Remove some unuseful files
#
#   rm bun_len.xvg bun_dist.xvg bun_z.xvg bun_tilt.vxg bun_tiltl.xvg 
bun_tiltr.xvg
#
#   # 3. [g_traj]	Get the coordinates of the vectors describing the 
helices
#   echo 0 | g_traj -nice 10 -f H1H2.pdb -s H1H2.pdb -ox H1H2.xvg
#
# Where the file generated by the g_traj with the option -ox contains 
coordinates of
# the vector describing the helices as a function of time.
#
# This script is intended to compute the inter-helical angle as a 
function of time.
#
# The output contains a header and two columns of number: time (ps), 
cross angle between 2 helices.
# The script allows also to do some further averaging if n is set to a 
different value than 1.
# In this case, make sure the total number of time step is a multiple of 
n !
#
# The input file, a.pdo file, contains 19 columns of numbers:
# column 1: time
# Columns 2-4: x,y and z coodinates of top point of helix 1 axe.
# Columns 5-7: x,y and z coodinates of middle point of helix 1 axe.
# Columns 8-10: x,y and z coodinates of bottom point of helix 1 axe.
# Columns 11-13: x,y and z coodinates of top point of helix 2 axe.
# Columns 14-16 x,y and z coodinates of middle point of helix 2 axe.
# Columns 17-19: x,y and z coodinates of bottom point of helix 2 axe.

# The script need .pdo file as an input. It also need to set-up several 
parameters:
# 1. The name of the input file, $file_name: data.xvg.
# 2. The root name for output file.
# 3. The variable $n, which control the number of point output and the 
number of
#    steps over which averaging is done (1st momentum).
#

# Some Perl statements
use strict;
use warnings;
# Program invoque an asin function described in the following package:
use Math::Trig;

##########################################################################
# Section of the script to set-up some parameters used for parsing and
# averaging data in .xvg file
##########################################################################
# Set name of the file to open
# Compulsory parameter to be put in the input line, 1st argument.
#my $name = '/home/vincent/Projects/H1H2.xvg';
my $file_name = $ARGV[0];
#
##########################################################################
# Name for output file
# Compulsory parameter to enter in the input line, 2nd Argument.
#my $name = '/home/vincent/Projects/H1H2_angle.xvg';
my $name = $ARGV[1];
#
##########################################################################
# Set variable n , the number of step per point in the output file.
# Average data and output them in ouput file every n steps.
# Set the per default value of $n
my $n = 1;
# $n can be set using input line, 3rd argument, optional
if (exists $ARGV[2] ) {
$n = $ARGV[2];
}
#
##########################################################################
#

# Start Script
my $scriptname = 'cross_angle.pl';

# Initialize variables.
my $time;

# Coordinates of Helix 1 axe
my $H1_top_x;
my $H1_top_y;
my $H1_top_z;

my $H1_mid_x;
my $H1_mid_y;
my $H1_mid_z;

my $H1_bottom_x;
my $H1_bottom_y;
my $H1_bottom_z;

# Coordinates of Helix 2 axe
my $H2_top_x;
my $H2_top_y;
my $H2_top_z;

my $H2_mid_x;
my $H2_mid_y;
my $H2_mid_z;

my $H2_bottom_x;
my $H2_bottom_y;
my $H2_bottom_z;

# Sum $time and all the coordinates for n steps of calculation, then
# average them by dividing the sum by n.
# sumforce and sumz gather the result of the sum over these n steps.
my $sumtime = 0;

my $sumH1_top_x = 0;
my $sumH1_top_y = 0;
my $sumH1_top_z = 0;

my $sumH1_bottom_x = 0;
my $sumH1_bottom_y = 0;
my $sumH1_bottom_z = 0;

my $sumH2_top_x = 0;
my $sumH2_top_y = 0;
my $sumH2_top_z = 0;

my $sumH2_bottom_x = 0;
my $sumH2_bottom_y = 0;
my $sumH2_bottom_z = 0;

# Input. Array where .xvg input is stored
my @kinkrecord;
# Input. Counter
my $c = 0;

# counter used to generate the output file.
my $line = 0;
# Array where output is stored
my @angle;

# Wellcome message
print STDERR "\nScript $scriptname.";
print STDERR "\nParse .xvg file generated by g_traj using the pdb file 
output by g_bundle -oa .\n";
print STDERR "Output file with cross angle.\n Check comments of this 
script for more information\n\n";
print STDERR "Parameters are given to the script throught the command 
line.\nNo need to edit the script.\n\n";
# Open input file
open (INPUT, $file_name) || die "\nCan't open Input File 
$file_name.\nUsage: ./$scriptname INPUT.xvg OUTPUT ($n - optional, 
perdefault $n = 1)\n\n";
my @xvgrecord = <INPUT>;

# Parsing .xvg data file

# Getting rid of the header of the .xvg file.

foreach my $xvgline (@xvgrecord) {

   # Discard blank line
   if ($xvgline =~ /^\s*$/) {
      next;

   # Discard comment line
   }  elsif($xvgline =~ /^\s*#/) {
      next;

   # Discard header line
   }  elsif($xvgline =~ /^@/) {
      next;

   # Discard other comment line
   }  elsif($xvgline =~ /^&/) {
      next;

   # Keep line, and add to array of line
   }  else {
      # Remove space at the beginning of the line
      $xvgline =~ s/^\s*//;

      $kinkrecord[$c] = $xvgline;
      $c = $c +1;
   }
}

# Turn the Scalar into an array of time lines
# be sure to start from 0 to $#kinkrecord. Array element starts at 0!
for my $i (0 .. $#kinkrecord ) {
   ($time, $H1_top_x, $H1_top_y, $H1_top_z, $H1_mid_x, $H1_mid_y, 
$H1_mid_z, $H1_bottom_x, $H1_bottom_y, $H1_bottom_z, $H2_top_x, 
$H2_top_y, $H2_top_z, $H2_mid_x, $H2_mid_y, $H2_mid_z, $H2_bottom_x, 
$H2_bottom_y, $H2_bottom_z) = split /\s*\s/, $kinkrecord[$i];

   # Summing data for later averaging
   $sumtime = $sumtime + $time;

   $sumH1_top_x = $sumH1_top_x + $H1_top_x;
   $sumH1_top_y = $sumH1_top_y + $H1_top_y;
   $sumH1_top_z = $sumH1_top_z + $H1_top_z;

   $sumH1_bottom_x = $sumH1_bottom_x + $H1_bottom_x;
   $sumH1_bottom_y = $sumH1_bottom_y + $H1_bottom_y;
   $sumH1_bottom_z = $sumH1_bottom_z + $H1_bottom_z;

   $sumH2_top_x = $sumH2_top_x + $H2_top_x;
   $sumH2_top_y = $sumH2_top_y + $H2_top_y;
   $sumH2_top_z = $sumH2_top_z + $H2_top_z;

   $sumH2_bottom_x = $sumH2_bottom_x + $H2_bottom_x;
   $sumH2_bottom_y = $sumH2_bottom_y + $H2_bottom_y;
   $sumH2_bottom_z = $sumH2_bottom_z + $H2_bottom_z;

   # Every n step, average and store the result in a new array
    if ( (int(($i+1)/$n) - (($i+1)/$n))  == 0 ) {
     # averaged time and kink angle over n steps
     my $time_aver = $sumtime/$n;

     my $H1_top_x_aver = $sumH1_top_x/$n;
     my $H1_top_y_aver = $sumH1_top_y/$n;
     my $H1_top_z_aver = $sumH1_top_z/$n;

     my $H1_bottom_x_aver = $sumH1_bottom_x/$n;
     my $H1_bottom_y_aver = $sumH1_bottom_y/$n;
     my $H1_bottom_z_aver = $sumH1_bottom_z/$n;

     my $H2_top_x_aver = $sumH2_top_x/$n;
     my $H2_top_y_aver = $sumH2_top_y/$n;
     my $H2_top_z_aver = $sumH2_top_z/$n;

     my $H2_bottom_x_aver = $sumH2_bottom_x/$n;
     my $H2_bottom_y_aver = $sumH2_bottom_y/$n;
     my $H2_bottom_z_aver = $sumH2_bottom_z/$n;

# Start Real Calculation !
     # Get coordinate of average axe, helix 1
     my $V1_x = $H1_bottom_x_aver - $H1_top_x_aver;
     my $V1_y = $H1_bottom_y_aver - $H1_top_y_aver;
     my $V1_z = $H1_bottom_z_aver - $H1_top_z_aver;

     # Get coordinate of average axe, helix 2
     my $V2_x = $H2_bottom_x_aver - $H2_top_x_aver;
     my $V2_y = $H2_bottom_y_aver - $H2_top_y_aver;
     my $V2_z = $H2_bottom_z_aver - $H2_top_z_aver;

     # Calculate norm of vector V1 and V2 (helix1 and 2 axe)
     my $norm1 = sqrt( ($V1_x * $V1_x) + ($V1_y * $V1_y) + ($V1_z * 
$V1_z));
     my $norm2 = sqrt( ($V2_x * $V2_x) + ($V2_y * $V2_y) + ($V2_z * 
$V2_z));

     # Normalize Coordinate of vector V1 and V2
     my $V1_x_n = $V1_x / $norm1;
     my $V1_y_n = $V1_y / $norm1;
     my $V1_z_n = $V1_z / $norm1;

     my $V2_x_n = $V2_x / $norm2;
     my $V2_y_n = $V2_y / $norm2;
     my $V2_z_n = $V2_z / $norm2;

     # Compute dot product of the normalized vector 1 and 2
     my $dot = ( $V1_x_n * $V2_x_n ) + ( $V1_y_n * $V2_y_n ) + 
( $V1_z_n * $V2_z_n );

     # Get the angle between Helix 1 and 2, in rad !
     my $ihangle_rad = asin($dot);
     # Convert angle from rad to degree
     my $ihangle = ( 360 / ( 2*3.1428571 )) * $ihangle_rad;

     # Store information in array
     $angle[$line] = sprintf("%6.1f %6.3f", $time_aver, $ihangle);

     # Reset $sumstep, $sumf and $sumz before new averaging loop
     $sumtime = 0;

     $sumH1_top_x = 0;
     $sumH1_top_y = 0;
     $sumH1_top_z = 0;

     $sumH1_bottom_x = 0;
     $sumH1_bottom_y = 0;
     $sumH1_bottom_z = 0;


     $sumH2_top_x = 0;
     $sumH2_top_y = 0;
     $sumH2_top_z = 0;

     $sumH2_bottom_x = 0;
     $sumH2_bottom_y = 0;
     $sumH2_bottom_z = 0;

     $line = $line + 1;
     }
   }

# Print the results in an .xvg file format

# Print Ouput File: cross angle = f(t)
open (OUTFILE, "> $name");
# Print Header of file
print OUTFILE "# This file was created by $scriptname\n";
print OUTFILE "# which is a perl script written by Vincent Lemaitre.";
print OUTFILE "\n\@    title \"Inter-Helices Angle\"\n";
print OUTFILE "\@    xaxis  label \"Time (ps)\"\n";
print OUTFILE "\@    yaxis  label \"Inter-helices angle [degree]\"\n";
print OUTFILE "\@TYPE xy\n";
print OUTFILE "\@ view 0.15, 0.15, 0.75, 0.85\n";
print OUTFILE "\@ legend on\n";
print OUTFILE "\@ legend box on\n";
print OUTFILE "\@ legend loctype view\n";
print OUTFILE "\@ legend 0.78, 0.8\n";
print OUTFILE "\@ legend length 2\n";
print OUTFILE "\@ s0 legend \"Inter-helices angle\"\n";

# Print averaged data
for my $j ( 0 .. $#angle ) {
     print OUTFILE "$angle[$j]\n";
     }

print OUTFILE "\n";
# End message
print STDERR "\nOutput file $name\.xvg with $line data points.\n";
print STDERR "\nHave fun !\n";

exit;



On Mercredi, mars 12, 2003, at 05:06 , sadhna wrote:

> hi all,
>         Can gromacs be used to calculate inter-helical angles between 
> the
> two helices? Could anyone suggest an alternative wroking program for the
> same.
>
> regards
>
>
> Sadhna Joshi
> Research Scholar
> Dept of Chemical Engg
> Indian Institute of Technology,Powai
> Mumbai-400076
> India
>
>
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
>