#!/usr/bin/perl

## Calculates the interaction energies between every residue pair in a range
##
## Author: Gadi Oron, gadi@proteologics.com
##
use PDL;
use strict;

my $rng1 = shift;
my $rng2 = shift;

STDERR->autoflush;

die "Usage energy_matrix.pl <RANGE1> <RANGE2>" unless defined $rng2;

my ($seg1,$start1,$end1);
my ($seg2,$start2,$end2);

## Parse ranges
die "Wrong range format! Should be SEG:<st>-<end>" unless($rng1 =~ /(.+):(\d+)-(\d+)/ and 
							  ($seg1,$start1,$end1) = ($1, $2, $3) and 
							  $rng2 =~ /(.+):(\d+)-(\d+)/ and 
							  ($seg2,$start2,$end2) = ($1, $2, $3));

for my $i ($start1..$end1) {
    for my $j ($start2..$end2) {
	my ($mean,$rms,$median,$min,$max) = calc_pair($seg1,$i,$seg2,$j);
	print "$seg1:$i\t$seg2:$j\t$mean\t$rms\n";
    }
}


sub calc_pair($$$$)
{
    my $seg1 = shift;
    my $res1 = shift;
    my $seg2 = shift;
    my $res2 = shift;

    open(PAIR,"| vmd -dispdev none > /dev/null") or die "Cannot open vmd pipe";
    
    print PAIR <<___; 
mol load pdb complete.pdb
set all [atomselect top all]
set grp1 [atomselect top "segname $seg1 and resid $res1"]
set grp2 [atomselect top "segname $seg2 and resid $res2"]
\$all set beta 0
\$grp1 set beta 1
\$grp2 set beta 2
\$all writepdb pair_tmp.pdb
quit
___

close(PAIR);

open(NAMD, "./namd2.exe pairs.namd +nice 10 | ")
    or die "Cannot run NAMD";
    
## alocate vectors
    my @timesteps = ();
    my @elects = ();
    my @vdws = ();
    my @totals = ();
## Read the NAMD output
    print STDERR "Parsing NAMD output\n";
    while(<NAMD>) {
	chomp;
	## Look for energy lines
#	print STDERR "$_\n";
	if(/^ENERGY: +(\d+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+) +(\S+)/) {
	    my ($ts,$elect,$vdw,$total) = ($1,$6,$7,$11);
	    push @timesteps, $ts;
	    push @elects , $elect;
	    push @vdws , $vdw;
	    push @totals, $total;
	    print STDERR "$ts...";
	}
    }
    print STDERR "\n\n";
    close(NAMD);

    print STDERR "\n";
    my $totpdl = pdl(@totals);
    
## Stats are ($mean,$rms,$median,$min,$max) = 
    my @stats = stats($totpdl);

    return @stats;
}