#! /usr/bin/perl -w
# $Id: out2xyz.pl,v 1.2 2004/05/18 06:38:47 akohlmey Exp $
# small script to extract coordinates from cpmd output file and 
# convert them into a multi step xyz file.
#
# an updated version may be found at
# http://www.theochem.ruhr-uni-bochum.de/go/cpmd-vmd.html
#
# (c) 2004 by axel.kohlmeyer@theochem.ruhr-uni-bochum.de

use strict;
use Getopt::Long;
use Pod::Usage;
use IO::File;

# initialize variables
my @x = ();
my @y = ();
my @z = ();
my @a = ();
my ($atoms, $atmcount, $step, $nfi) = (0,0,0,0);
my $bohr = 0.529177249;

# variables for commandline parsing.
my $in = \*STDIN;
my $out = \*STDOUT;
my $help = 0;
my $first = 1;

# parse commandline arguments
Getopt::Long::Configure ("bundling");
GetOptions( 'input|i=s'  => sub { $in  = new IO::File "$_[1]", "r";
                                   die "\ncould not open file $_[1] for reading: $!\n\n"
                                       unless defined $in;},,
            'output|o=s' => sub { $out = new IO::File "$_[1]", "w";
                                   die "\ncould not open file $_[1] for writing: $!\n\n"
                                       unless defined $out;},
            'help|h'      => \$help,)
  or pod2usage(2);

# print man page and exit if requested.
pod2usage(-exitstatus => 0, -verbose => 2) if $help;

# parse the the reference geometry to get the number of atoms.
while (<$in>) {
    last if (/\Q*********************** ATOMS *************************/);
}
<$in>;
$atoms=0;
while (<$in>) {
    last if (/\Q*********************************************************/);
    ++$atoms;
}


# read and convert the trajectory file
$atmcount = 0;
my $nbeads   = 0;
my $oldnfi   = 0;
my $savedline = "";

while (<$in>) {
    # coordinate format without forces
    if (/.*\*                      ATOMIC COORDINATES                      \*.*/) {
        ++$nfi;
        $_ = <$in>;
        print "$atoms\n";
        print "step: $nfi\n";
        for (my $i=0; $i < $atoms; ++$i) {
            $_ = <$in>;
            @_ = split;
            printf "%4s %12.5f %12.5f %12.5f\n", $_[1], $_[2] * $bohr,
                                                 $_[3] * $bohr, $_[4] * $bohr;
        }
    }
    # coordinate format with forces
    if (/.*ATOM          COORDINATES            GRADIENTS \(-FORCES\).*/) {
        ++$nfi;
        print "$atoms\n";
        print "step: $nfi\n";
        for (my $i=0; $i < $atoms; ++$i) {
            $_ = <$in>;
            @_ = split;
            printf "%4s %12.5f %12.5f %12.5f\n", $_[1], $_[2] * $bohr,
                                                 $_[3] * $bohr, $_[4] * $bohr;
        }
    }
}

exit 0;


# embedded documentation.

__END__

=head1 NAME

out2xyz.pl - convert cpmd output files into multistep XYZ format.

=head1 SYNOPSIS

out2xyz.pl [options] 

Options:

 --input=<name>, -i <name>   set input file to '<name>'.  defaults to stdin.
 --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 atom coordinates printed in the
output of a cpmd run into a multistep xyz format movie file. 
The output has been tested with VMD. 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
reads from stdin (i.e. your keyboard).

=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<-h, --help>

Prints this manual page and exits.

=back

=head1 KNOWN BUGS

=head1 VERSION

$Id: out2xyz.pl,v 1.2 2004/05/18 06:38:47 akohlmey Exp $

=cut

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