#!/usr/bin/perl

use lib "/home/lmt09/SOFTWARE/taopackage/";
use ESPT::Glog 0.07;
use File::Basename;

# When adding new options you need to update:
# 1. SYNOPSIS manpage
# 2. Argument check (default set to zero and if check section)
# 3. OPTIONS manpage
# 4. Input error if no input
# 5. main while loop and data output
# 6. call to plotopt subroutine
# 7. usage subroutine

=head1 NAME

GREGAP - Gaussian Regular Expression Graphical Analysis in Perl 

=head1 SYNOPSIS

B<gregap>  [ B<-d> F<Gaussian_log_file> ] [ B<-e> F<Gaussian_log_file> ] [ B<-f> F<Gaussian_log_file> ] [ B<-h> ] 

=head1 DESCRIPTION

This program uses GNUPlot to allow graphical monitoring of Gaussian jobs.

=cut

### Main Program ###
our $version = "1.0";

#######################################################################
# GREGAP                                                              #
# CREATED 05/03/11                                                    #
# LAST MODIFIED 15/03/11                                              #
# LEE THOMPSON                                                        #
# This Program takes the forces from each geometry optimization step  #
# and shows it in graphical format for debugging purpose. Future      #
# intention to extend this further with options.                      #
#######################################################################


# check for arguments 
usage() if ( $#ARGV < 0 );

help() if $ARGV[0] eq "-h"; 
help() if $ARGV[0] eq "--help";

$haveInput = 0;
$isLogFile = 0;
$DoForces = 0;
$DoDisp = 0;
$DoEnerg = 0;

for (my $i=0; $i<=$#ARGV; $i++) { 
	if ($ARGV[$i] eq "-f"){
		$infile = $ARGV[$i+1];
		$haveInput = 1;
		$DoForces = 1;
	}
	if ($ARGV[$i] eq "-d"){
		$infile = $ARGV[$i+1];
		$haveInput = 1;
		$DoDisp = 1;
	}
        if ($ARGV[$i] eq "-e"){
		$infile = $ARGV[$i+1];
                $haveInput = 1;
                $DoEnerg = 1;
	} 
}

=head1 OPTIONS

Command line option specifications are processed from left to right and may 
be specified more than once. If conflicting options are specified, later ones take precedence.

=over 16

=item B<-d>

Display RMS and maximum displacement against iteration number.

=item B<-e> 

Display Energy in Hartree against iteration number.

=item B<-f>	

Display RMS and maximum force against iteration number.

=item B<-h>

=item B<--help> 	

Print full GREGAP documentation via perldoc. Cannot be used with other options.

=back

=cut

print "\nGREGAP $version : Gaussian RegEx Graphical Analysis \n\n";

if ($haveInput == 0){
	if (@ARGV[0] =~ /^-[def]/){
		print "Gaussian log file missing.\n";
		die "Exit.\n $!\n";
	}elsif (@ARGV[0] =~ /^-/){
		print "### Invalid option ###\n";
		usage();
	}else{
		print "### No option specified ###\n";
		usage();
	}
}

open(INFILE,$infile) or die "Could not read $infile\n!$\n";

$[ = 1;			# set array base to 1

open(DATA,'>plot.dat') || die "Cannot open plot data file!\n";

#my $Glog = ESPT::Glog->new();
#$Glog->analyze($infile);
#$route = $Glog->get(ROUTE);
#print "$route\n";
#$ftype = $Glog->get(JOBTYPE);
#print "$ftype\n";
#@keywords[$i++] = $Glog->get(KEYWORDS);
#foreach (@keywords){
#	print "$_\n";
#}
#$fcomp = $Glog->get(COMPLETE);
#print "$fcomp\n";
#@eelec = $Glog->get(EELEC);
#foreach (@eelec){
#	print "$_\n";
#}


while(<INFILE>){
    next if /^$/;

    if ($isLogFile == 0){
	if( /^\s+Entering\s+Gaussian\s+System/ ){
		$isLogFile = 1;
		print "Input file $infile is a Gaussian log file.\n";
    	}else{
		print "Input file $infile is not a Gaussian log file.\n";
		die "Exit.\n\n";
    	}
    }
    if (/Maximum Force\s+([\d\.]+)/) {
    	push @MF, $1;
    }
    if (/RMS     Force\s+([\d\.]+)/) {
    	push @RF, $1;
    }
    if (/Maximum Displacement\s+([\d\.]+)/){
	push @MD, $1;
    }
    if (/RMS     Displacement\s+([\d\.]+)/){
        push @RD, $1;
    }
    if (/^\s+SCF\s+Done:\s+.*\s+=\s+(-*\d+\.\d+)/){
	push @EEl, $1;
    }
}

print DATA "# ITN\tRMS Forc\tMax Forc\tRMS Disp\tMax Disp\tSCF Ener\n";
for ($i = 1; $i <= scalar @MF; $i++) {
	printf DATA ("%d\t%f\t%f\t%f\t%f\t%f\n", $i, $RF[$i], $MF[$i], $RD[$i], $MD[$i], $EEl[$i]);
}

close DATA;
close INFILE;

if ($DoForces == 1){
	&plotopt('plot.dat', 'Iteration', 'Force', 'CONVERGENCE', '', '', 'lines','2|3', 'RMS Force|Max Force');
}
if ($DoDisp == 1){
	&plotopt('plot.dat', 'Iteration', 'Displacement', 'CONVERGENCE', '', '', 'lines', '4|5', 'RMS Displacement|Max Displacement');
}
if ($DoEnerg == 1){
	&plotopt('plot.dat', 'Iteration', 'Energy (Har)', 'ENERGY', '', '', 'lines', '6', '');
}

    #######################################################################
    # SUBROUTINE PLOTOPT                                                  #
    # CREATED 10/05/12                                                    #
    # LAST MODIFIED 10/05/12                                              #
    # LEE THOMPSON                                                        #
    # This subroutine takes data from plot.dat, constructs a gnuplot file #
    # plot.plt and displays the gnuplot graphic. Series and columns input #
    # use "|" as a delimiter. First column is plotted as x data and each  #
    # subsequent column corresponds to y data for each series. Extend by  #
    # inserting option for surface plot.                                  #
    #######################################################################

sub plotopt {
    local($input, $xlabel, $ylabel, $title, $xrange, $yrange, $style, $column,$series) = @_;

    open FHDL, ">plot.plt" or die "Cannot open plot input\n";

    @names = split(/\|/, $series);
    @columns = split(/\|/, $column);
    if (scalar @columns > 1) {
	$command = sprintf("plot \"%s\" using 1:%s title \"%s\"\n", $input,$columns[1], $names[1]);
	for ($i = 2; $i <= scalar @columns; $i++) {
	    $command = sprintf("%sreplot \"%s\" using 1:%s title \"%s\"\n",$command, $input, $columns[$i], $names[$i]);
	}
    }
    elsif (scalar @names == $column) {
	$command = sprintf("plot \"%s\" using 1 title \"%s\"\n", $input,$names[1]);
	for ($i = 2; $i <= $column; $i++) {
	    $command = sprintf("%sreplot \"%s\" using %s title \"%s\"\n",$command, $input, $i, $names[$i]);
	}
    }
    else {
	$command= sprintf("plot \"%s\" using 1:%s title \"%s\"\n", $input,$column, $series);
    }

    printf FHDL "set xlabel \"%s\"\n" . "set ylabel \"%s\"\n" . "set xrange [%s]\n" . "set yrange [%s]\n" . "set grid\n" .  "set data style %s\n" . "set title \"%s\"\n" . '%s' . 'pause -1', $xlabel, $ylabel, $xrange, $yrange, $style, $title, $command;
    close FHDL;

    system('gnuplot plot.plt');
}


    #######################################################################
    # SUBROUTINE HELP                                                     #
    # CREATED 15/05/12                                                    #
    # LAST MODIFIED 10/05/12                                              #
    # LEE THOMPSON                                                        #
    # Full help documentation                                             #
    #######################################################################

sub help {
	system("perldoc /home/lmt09/PHD_Y3/PROGRAMS/PERL/LOGANALYSIS/VERSION_1/gregap_1-0.pl");
	exit;
}


    #######################################################################
    # SUBROUTINE USAGE                                                    #
    # CREATED 15/05/12                                                    #
    # LAST MODIFIED 10/05/12                                              #
    # LEE THOMPSON                                                        #
    # Brief usage documentation                                           #
    #######################################################################

sub usage {
	print "\nGREGAP $version : Gaussian Regular Expression Graphical Analysis in Perl \n";
	print "\nUsage: gregap [options] <Gaussian_logfile>\n";
	print "\t-d\t\t\tdisplay RMS and Maximum displacement\n";
	print "\t-e\t\t\tdisplay energy per iteration\n";
	print "\t-f\t\t\tdisplay RMS and Maximum force\n";
	print "\t-h\t\t\tprint full documentation\n"; 
	print "\n";
	exit;
}

=head1 EXAMPLES

=over

=item gregap

Called without any parameters, GREGAP will display usage information.
If B<-h> or B<--help> is passed, then the full GREGAP documentation is displayed via perldoc.

=item gregap -f foo.log

GREGAP reads foo.log and creates plot of Max Force and RMS Force.

=back

=head1 NOTES

GREGAP

=head1 VERSION

1.0

=head1 AUTHOR

Lee Thompson, E<lt>l.thompson09@imperial.ac.ukE<gt>

=head1 COPYRIGHT

Copyright (c) 2012 Lee Thompson

This is free software; modification and redistribution is made under the same terms as Perl itself.

=cut