#!/usr/bin/perl
#
# Extract time-averaged HBond per residue data from a trajectory
# Original 2006 w/ Thomas Chun
# re-written 4/2010 
#
# usage: reshb.pl contact.out
#
# data must be grepped from contact.out 
# ($ charmm < contact.inp > contact.out)
#  ..contact.out has pertinent details:
#
#     I-atom                 J-atom          <occupancy>     <time>  # events
# ---------------------------------------------------------------------------
#      A15  ALA  1    HN   -  A15  ALA  1    O      0.03        5.8      26
#                A15  ALA  1    HN       Summary:   0.06        6.7      55
#
#  we want only the data lines, not the summary lines.
# The first part of this program will grep the charmm output (contact.out)
####
# shell> grep " -  A15  ALA" contact.tmp > contact.dat
####
# I don't want to use Scanf any more;  will fix this with a 
# 'split()' function...
#
use String::Scanf; # imports sscanf()
open(INFILE, $ARGV[0]) or die "Cannot open $ARGV[0]\nYou must specify a file to fix!!\n\n";
$filen = $ARGV[0];
substr($filen, -4, 4) = ""; # drop that suffix
open(OUTFILE, ">$filen-HBres.dat") or die "couldn't open $filen-HBres-dat for writing...\n";
flock(OUTFILE, 2);
#print OUTFILE "#\n#RESNO   0    1    2      310    Alpha  Pi     six    seven  Total\n#\n";
my $x = 0;
my $y = 0;
my $resno = 0;
my $span = 0;
my $infile = 0;
my @count = 0;
# Initialize array (do we need to initialize?)
for $x (0 .. 20) { # this should be number of residues;
#  need to extract this information from the file name (?) or something
#  ..maybe from the data list?  Max $resnoi?
    $sum[$x] = 0;
    for $y (0 .. 7) { # five cols of data?
	$count[$x][$y] = 0;
    }
}
$totalresidues = 0;
foreach $line (<INFILE>)
{
    # some sort of conditional here....
    ($junk1,$residi,$resnoi,$atomi,$dash1,$junk2,$residj,$resnoj,$atomj,$occ,$lifetime,$numhb) = split(' ', $line);
    print "$span\n";
    if ($resnoi > $totalresidues)
    {
	$totalresidues = $resnoi;
    }
    $span=abs($resnoi - $resnoj);
    $count[$resnoi][$span] += $occ;
    $count[$resnoj][$span] += $occ;
    $sum[$resnoi] += $occ;
    $sum[$resnoj] += $occ;
}
print "$totalresidues\n";
format OUTFILE =
#RESNO 0     1     2      310    Alpha  Pi     six    seven  Total\n#\n";
  @>>> @>>>> @>>>> @>>>>  @>>>>  @>>>>  @>>>>  @>>>>  @>>>>  @>>>>
$x   $count[$x][0]  $count[$x][1]  $count[$x][2]   $count[$x][3]   $count[$x][4]   $count[$x][5]   $count[$x][6]   $count[$x][7]    $sum[$x]
.
for $x (1 .. $totalresidues) { # this should be residue number
    write (OUTFILE);
}

#test to dump array
#print "$count[33][3]\n"; # "33" means t = 66.00
#for $time ( @count ) {
#    print OUTFILE " $time @$time\n";
#}

close(OUTFILE);

