#!/usr/local/bin/perl -w

# dm ("Distance Matrix/Measurement") v. 7 Apr 1996
# Dave Schweisguth <dcs@proton.chem.yale.edu>

### Preliminaries

require 5.001;			    # Perl 5.001 required, 5.001m recommended
use strict;			    # Require optional-but-desirable practices

### Parameters

# Environment

(my $whatami = $0) =~ s|.*/||;	    # `basename $0`

# Configuration

my $res_diff	= 1000000;	    # i.e. no limit
my $d_diff	= 5;

# Initialization (don't change these)

my($atom_type, $res_type, $chain, $res_num, $x, $y, $z, $seg, $tag,
    @tag, @res_num, @x, @y, @z, $i, $j, $d);

### Arguments and error-checking

my($arg, $sign, $first, $rest);
while (@ARGV and ($sign, $first, $rest) = ($ARGV[0] =~ /^([\-+])(.)(.*)/)) {
    if ($sign eq '+' && $first !~ /[\0]/) { # -/+ switches (none at the moment)
	&usage("$whatami: $sign$first is not an option.\n");
    }
    if ($first =~ /[dr]/) {	# Switches with arguments
        shift;
    	$arg = $rest ne '' ? $rest : @ARGV ? shift :
            &usage("$whatami: -$first requires an argument.\n");
    } elsif ($rest eq '') {
	shift;
    } else {
    	$ARGV[0] = "$sign$rest";
    }
    if ($first eq 'd')	  { $d_diff = $arg; }
    elsif ($first eq 'r') { $res_diff = $arg; }
    elsif ($first eq 'u') { &usage(0); }
    else		  { &usage("$whatami: $sign$first is not an option.\n"); }
}

sub usage {
    my($message) = $_[0];

    warn $message if $message;
    warn <<EOP;
Usage: $whatami [-d #] [-r #] PDB-file
-d #	Print distances less than # (default $d_diff)
-r #	Print distances between residue numbers differing by less than #
	    (default $res_diff, i.e. no limit)
-u	This message
EOP
    exit !! $message;
}

### Do it

while(<>) {
    next unless ($atom_type, $res_type, $chain, $res_num, $x, $y, $z, $seg) =
	/^(?:ATOM  |HETATM).{5} (.{5})(.{3}) (.)(.{4}).   (.{8})(.{8})(.{8}).{6}.{6} .{3}  (.{4})/;
	# Column 17 ("alternate location") appended to columns 13-16
	#     (atom type), mostly to compensate for bad Insight PDB
	# Columns 77-80 (element and charge fields) ignored
    foreach $i ($atom_type, $res_type, $chain, $res_num, $x, $y, $z, $seg) {
	$i =~ s/^\s*(\S*)\s*/$1/;
    }
    foreach $i ($atom_type, $res_type, $res_num, $x, $y, $z) {
	# Chain and segment IDs not always present, so not required
        warn "$whatami: Bad ATOM record on $ARGV line $.!\n" if $i eq '';
    }
    $tag = $seg ne '' ? $seg : $chain;
    next unless $atom_type =~ /H/;
    push(@tag, $tag . ($tag ? ' ' : '') . "$res_type $res_num $atom_type");
    push(@res_num, $res_num);
    push(@x, $x);
    push(@y, $y);
    push(@z, $z);
}

foreach $i (0 .. $#tag) {
    foreach $j (($i + 1) .. $#tag) {
	next if abs($res_num[$i] - $res_num[$j]) > $res_diff;
        $d = sqrt(
	    ($x[$i] - $x[$j]) ** 2 +
	    ($y[$i] - $y[$j]) ** 2 +
	    ($z[$i] - $z[$j]) ** 2
	);
	printf("$tag[$i]\t$tag[$j]\t%8.3f\n", $d) unless $d > $d_diff;
    }
}