#! /usr/bin/perl -w
use strict;
use Getopt::Long;
use Pod::Usage;
use IO::File;


########################################################################
# adjustable parameters
my $writepdb = 0;
my $oxpbc    = 0;
my $cutoff = 1.31;
my $cell = 9.85275;
my $replicate = 0;

########################################################################
# global variables
my $natoms = 0;
my $comment = "";
my @hx = ();
my @hy = ();
my @hz = ();
my @ox = ();
my @oy = ();
my @oz = ();

my $infile;
my $outfile;
my $in =  \*STDIN;
my $out = \*STDOUT;
my $hl = 0;
my $help = 0;
my $velocity = 0;
my $first = 1;
my $hlfile = "";
my @xcelloff = (0.0);
my @ycelloff = (0.0);
my @zcelloff = (0.0);

# parse commandline arguments
Getopt::Long::Configure ("bundling");
GetOptions( 'writepdb|p'     => \$writepdb,
            'oxpbc|o'        => \$oxpbc,
            'velocity|v'     => \$velocity,
            'bondlength|b=f' => \$cutoff,
            'cell|c=f'       => \$cell,
            'hlfile|n=s'     => \$hlfile,
            'replicate|r'    => \$replicate,
            'help|h'         => \$help,)
  or pod2usage(2);

if ($replicate) {
#    @xcelloff = (-0.5 * $cell, -0.5 * $cell, -0.5 * $cell, -0.5 * $cell,
#                  0.5 * $cell,  0.5 * $cell,  0.5 * $cell,  0.5 * $cell);
#    @ycelloff = (-0.5 * $cell, -0.5 * $cell,  0.5 * $cell,  0.5 * $cell,
#                 -0.5 * $cell, -0.5 * $cell,  0.5 * $cell,  0.5 * $cell);
#    @zcelloff = (-0.5 * $cell,  0.5 * $cell, -0.5 * $cell,  0.5 * $cell,
#                 -0.5 * $cell,  0.5 * $cell, -0.5 * $cell,  0.5 * $cell);
    @xcelloff = (         0.0,          0.0,          0.0,          0.0);
    @ycelloff = (-0.5 * $cell, -0.5 * $cell,  0.5 * $cell,  0.5 * $cell);
    @zcelloff = (-0.5 * $cell,  0.5 * $cell, -0.5 * $cell,  0.5 * $cell);
}
# print man page and exit if requested.
pod2usage(-exitstatus => 0, -verbose => 2) if $help;

# open input and output file, if given.
# otherwise we will use stdout and stdin.
if (@ARGV) {
    if ($ARGV[0]) {
        $in = new IO::File $ARGV[0], "r";
        die "\ncould not open input xyz file '$ARGV[0]' for reading: $!\n\n"
          unless defined $in;
    }
    if ($ARGV[1]) {
        $out = new IO::File "$ARGV[1]", "w";
        die "\ncould not open file $ARGV[1] for writing: $!\n\n"
          unless defined $out;
    }
}

if (not "$hlfile" eq  "") {
    $hl = new IO::File "$hlfile", "w";
    die "\ncould not open file $hlfile for writing: $!\n\n"
        unless defined $hl;
    $hlfile = 1;
}

while (<$in>) {

    # read xyz header
    chomp;
    $natoms = int($_);
    $comment = <$in>;

    for (my $i = 0; $i < $natoms; ++$i) {
        $_ = <$in>;
        chomp $_;

        # store hydrogen positions after translation back to unit cell
        if (/^\s*H\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*).*$/) {
            push @hx, &fixpbc($1);
            push @hy, &fixpbc($2);
            push @hz, &fixpbc($3);
        }

        # store oxygen positions after translation back to unit cell
        if (/^\s*O\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*).*$/) {
	    if ($oxpbc) {
             	push @ox, &fixpbc($1);
            	push @oy, &fixpbc($2);
            	push @oz, &fixpbc($3);
	    } else {	
            	push @ox, $1;
            	push @oy, $2;
            	push @oz, $3;
	    }
        }
    }
    $hlfile = 1;

    # account for replicas
    $natoms *= $#xcelloff + 1;
    # print output header
    if ($writepdb) {
	print $out "HEADER    $comment";
	print $out "REMARK    natoms = $natoms\n";
	printf $out "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n",
	  $cell, $cell, $cell, 90.0, 90.0, 90.0;
    } else {
        print $out "$natoms\n$comment";
    }
    # bond analysis for all oxygens
    my $atmcount = 0;
    for (my $i = 0; $i < $#ox + 1; ++$i) {

        my $nbond;
        my ($x,$y,$z,$dist) = (0,0,0,0);

        $nbond = 0;
        for (my $j = 0; $j < $#hx + 1; ++$j) {

            # translate hydrogen to minimum image position
            # with respect to the current oxygen.
            $x = &fixdist($ox[$i], $hx[$j]);
            $y = &fixdist($oy[$i], $hy[$j]);
            $z = &fixdist($oz[$i], $hz[$j]);

            # if the distance is below the bond cutoff
            # keep the hydrogen translation so that we
            # only replace complete molecules with their images.
            # count the bond so that we can mark oxygens
            # with an unusual number of bonds.
            $dist = sqrt(  ($ox[$i]-$x)**2
                         + ($oy[$i]-$y)**2
                         + ($oz[$i]-$z)**2);
            if ($dist < $cutoff) {
                $hx[$j] = $x;
                $hy[$j] = $y;
                $hz[$j] = $z;
                ++ $nbond
            }
        }

        for (my $r = 0; $r <= $#xcelloff; ++$r) {
            if ($writepdb) {
                ++$atmcount;
                if ($nbond == 2) {
                    printf $out "HETATM% 5d  OW  WAT % 5d    % 8.3f% 8.3f% 8.3f" .
                      "% 6.2f% 6.2f    O   \n",  $atmcount, $i,
                        $ox[$i] + $xcelloff[$r], $oy[$i] + $ycelloff[$r],
                          $oz[$i] + $zcelloff[$r], $nbond, 0.0;
                } else {
                    printf $out "HETATM% 5d  NW  WAT % 5d    % 8.3f% 8.3f% 8.3f" .
                      "% 6.2f% 6.2f    N   \n",  $atmcount, $i,
                        $ox[$i] + $xcelloff[$r], $oy[$i] + $ycelloff[$r],
                          $oz[$i] + $zcelloff[$r], $nbond, 0.0;
                }
            } else {
                if ($nbond == 2) {
                    printf $out "O  %12.6f%12.6f%12.6f\n", $ox[$i] + $xcelloff[$r],
                      $oy[$i] + $ycelloff[$r], $oz[$i] + $zcelloff[$r];
                } else {
                    printf $out "N  %12.6f%12.6f%12.6f\n", $ox[$i] + $xcelloff[$r],
                      $oy[$i] + $ycelloff[$r], $oz[$i] + $zcelloff[$r];
                    if ($hl && ($hlfile == 1)) {
                        printf $hl "1\nFirst atom with nbonds!=2. idx=$i\n"
                          . "O  %12.6f%12.6f%12.6f\n", $ox[$i], $oy[$i], $oz[$i];
                        $hlfile=0;
                    }
                }
            }
        }

    }

    for (my $r = 0; $r <= $#xcelloff; ++$r) {
        for (my $i = 0; $i <= $#hx; ++$i) {
            if ($writepdb) {
                ++$atmcount;
                printf $out "HETATM% 5d  HW  WAT % 5d    % 8.3f% 8.3f% 8.3f" .
                  "% 6.2f% 6.2f    H   \n",  $atmcount, $i,
                    $hx[$i] + $xcelloff[$r], $hy[$i] + $ycelloff[$r],
                      $hz[$i] + $zcelloff[$r], 1.0, 0.0;
            } else {
                printf $out "H  %12.6f%12.6f%12.6f\n", $hx[$i] + $xcelloff[$r],
                  $hy[$i] + $ycelloff[$r], $hz[$i] + $zcelloff[$r];
            }
        }
    }

    if ($writepdb) {
	print $out "TER\nEND\n";
    }
    @hx = ();
    @hy = ();
    @hz = ();
    @ox = ();
    @oy = ();
    @oz = ();
}

exit 0;

sub fixpbc {
    my ($coord) = @_;
    while ($coord*1.0 > 0.5*$cell) {
        $coord -= $cell;
    } 
    while ($coord*1.0 < -0.5*$cell) {
        $coord += $cell;
    }
    return $coord;
}

sub fixdist {
    my ($a, $b) = @_;

    if ($a - $b > 0.5*$cell) {
        $b += $cell;
    } elsif ($a - $b < -0.5*$cell) {
        $b -= $cell;
    }
    return $b;
}

# embedded documentation.

__END__

=head1 NAME

highlight.pl - rearrange an xyz movie file to visualize proton transfer.

=head1 SYNOPSIS

highlight.pl [options] [<input file> [<output file>]]

Options:

 --input=<name>, -i <name>    set input file to '<name>'.
                              defaults to TRAJECTORY.
 --geometry=<name>, -g <name> set name of reference geometry file
                              (in xyz format) defaults to GEOMETRY.xyz.
 --output=<name>, -o <name>  set output file to '<name>'. defaults to stdout.
 --help, -h                  print full manual page.

=head1 DESCRIPTION

B<traj2xyz.pl> tries to convert a cpmd style TRAJECTORY file into
a multistep xyz format movie file. The output has been tested with
molden. Other programs to view the resulting xyz file may or may not work.

=head1 OPTIONS

=over 8

=item B<-i, --input>

Set the input file name. If you do not set this option the program 
tries to read a file named 'TRAJECTORY' in the current directory.

=item B<-o, --output>

Set the output file name. If you do not set this option the program
writes to the standard output (i.e. your terminal).

=item B<-g, --geometry>

Set the reference geometry file name. If you do not set this option, 
the program  tries to read a file named 'GEOMETRY.xyz' in the current 
directory. This file B<has> to be in xyz-format.

=item B<-h, --help>

Prints this manual page and exits.

=back

=head1 KNOWN BUGS

no coordinate boxshift postprocessing when using PBC.

=head1 VERSION

$Id: highlight.pl,v 1.1 2003/08/18 17:28:07 akohlmey Exp $

=cut

# Local Variables:
# mode: cperl
# cperl-indent-level: 4
# End:
