#!/usr/bin/perl -w
use strict;
use File::Temp qw/ tempdir /;
use List::Util qw[min];
# gelly qm helper using GAUSSIAN 03 
# Global Phasing Ltd 2011, NOT COPYRIGHTED! 
# author: Oliver Smart
# contact: buster-develop@GlobalPhasing.com
#
# You are encouraged to use this as a template for running other QM packages
# as a helper. We would be happy to help buster-develop@GlobalPhasing.com 
#
# local information about where gaussian executable is:
my $gaussianver="g03";
print "PROGRAM_TAG GAUSSIAN $gaussianver helper_27_July_12\n"; # will appear in final pdb

# create scratch directory
my $SCRDIR = tempdir();
my $gaussiancom="export GAUSS_SCRDIR=$SCRDIR; g03";
print "Gaussian command=$gaussiancom\n";
#
# get number of processors on machine
my $nproc=`grep -c processor /proc/cpuinfo`; chomp($nproc);
# use up to 4 processors but not more than actual number
my $procuse=min(4,$nproc);
my $link0="\%NProc=$procuse"; # use up to 4 processors
my $mlines;
my $line;
my $lcount=0;
my $numatoms=0;
my $gradnum=0;
my $crdin= "gelly_coords_for_qmhelper.crdin"
; # input file
my $gaussianinp="gelly_gaussian.inp"; # output file
my $conv_hartee2kjmol=2625.5; # http://mccammon.ucsd.edu/~dzhang/energy-unit-conv-table.html
my $conv_bohr2angs=0.529177249; # http://www.unitconversion.org/length/bohr-radius-to-angstroms-conversion.html

print "\n";
print "gelly qm helper script to run GAUSSIAN for energy/gradient calculation\n";
print "\n";
print "helper script location $0\n";

# deal with the command line args to this script
# there will be three (always)
# 1. integer charge  
# 2. integer multiplicity
# 3. string method
my $numArgs = $#ARGV + 1;
($numArgs > 2) || die "ERROR this helper script should be supplied with 3 command line args but only have $numArgs\nERROR ";
my $charge=$ARGV[0];
my $multip=$ARGV[1];
my $method=$ARGV[2];
print "picked up charge=$charge multip=$multip method=$method from command line\n";
if ($method eq "AUTO")  {
   $method = "HF/6-31G(d)"; # small basis ab initio
   print "method is set to AUTO: selecting $method\n";
   }
print "METHOD_TAG $method\n "; # for final pdb

# open file $crdin in current directory
open(CRDIN,$crdin) || die "ERROR cannot open file: $crdin to read atoms/coords from\nERROR ";
print "have opened crdin format coordinate file $crdin\n";

# open gaussian input file
my $random_number = int(rand(100000));
$gaussianinp="gelly_gaussian_".$$."_".$random_number.".inp";
open(GAUSSIANINP,">".$gaussianinp) || die "cannot open $gaussianinp";
print "will now write out a gaussian input file $gaussianinp\n";
print GAUSSIANINP "$link0\n" if ($link0);
print GAUSSIANINP "# $method FORCE\n";
print GAUSSIANINP "\n";
print GAUSSIANINP "gaussian input file written by gelly_qm_helper_gaussian.pl\n";
print GAUSSIANINP "\n";
print GAUSSIANINP "$charge $multip\n";
while (<CRDIN>) {
    $line = $_;
    $lcount++;
    if ($lcount>8)  { # skip header
        # have a coord line
        $numatoms++;
        #print "debug $line";
        my $atnam=substr($line,10,4);
        my $atnum=substr($line,16,3);
        my $xcoord=substr($line,22,27);
        my $ycoord=substr($line,52,27);
        my $zcoord=substr($line,82,27);
        #print "debug atnam=\"$atnam\" numatoms=\"$numatoms\" \"$xcoord\"  \"$ycoord\"  \"$zcoord\"\n";
        # write the atomic number and cartesian coords for each atom
        printf GAUSSIANINP "%d %f %f %f\n", $atnum, $xcoord, $ycoord, $zcoord;
    }
}
close(CRDIN);
 
print GAUSSIANINP "\n"; #need a couple of blank lines at the bottom
print GAUSSIANINP "\n";
print "number of atoms converted $numatoms\n";

my $command= "$gaussiancom < $gaussianinp >& gelly_gaussian.log";
$command = "bash -c \"$command\"";
print "now run gaussian with system command:\n   $command \n";

# DEBUG test helper without gaussian by simply copying emailed file to gelly_gaussian.log
# DEBUG file localation as $method
####$command = "cp $method  gelly_gaussian.log";
####print "DEBUG  test without gaussian instead use emailed file:\n   $command \n";

system($command);

# now check and parse output log  file
my $gaussianlog= "gelly_gaussian.log";
my $skipfour = -1; # to skip 4 lines 
if (-r $gaussianlog) {
   print "Good have got output file $gaussianlog\n";
   print "now try to parse out energy and gradients\n";
   open(GAUSSIANLOG,$gaussianlog) || die "cannot open $gaussianlog";
   my $readgrad=0;
   while (<GAUSSIANLOG>) {
      $line = $_;
      # AM1 energy look for line like
      #  Energy=   -0.084678594320 NIter=  14.
      if ($line =~/ Energy= *(.*) NIter=/) {
          my $energy=$1;
          print "AM1 type energy= $energy hartrees =   "; 
          $energy=$energy*$conv_hartee2kjmol; # convert to kJ/mol
          print                               " $energy kJ/mol\n"; 
          print "ENERGYTAG   $energy\n";
      }
      # HF energy look for line like
      #  SCF Done:  E(RHF) =  -115.030043902     A.U. after   11 cycles
      elsif ($line =~/ SCF Done:  E\(.HF\) = *(.*) *A\.U\./) {
          my $energy=$1;
          print "RHF type energy= $energy hartrees =   "; 
          $energy=$energy*$conv_hartee2kjmol; # convert to kJ/mol
          print                               " $energy kJ/mol\n"; 
          print "ENERGYTAG   $energy\n";
      }
      elsif ($line =~/ \*\*\*\*\* Axes restored to original set \*\*\*\*\*/) {
          #print "debug got to Axes... $line";
          $skipfour = 4; # need to skip next 4 lines
      }
      elsif ($skipfour > 0) {
          #print "debug skipping:    $line";
          $skipfour--;
      }
      elsif ( ($skipfour == 0) && ($gradnum < $numatoms) ) {
          $gradnum++;
          #print "debug grad #$gradnum from: $line";
          # gradient from lines like
          #                       123456789012345
          #    1          6           0.007083447    0.012979385   -0.010013744
          #    2          8          -0.024807510    0.013305065    0.008396735
          #1234567890123456789012345678901234567890123456789012345678901234567890
          #         1         2         3         4         5         6 
          my $xgrad=substr($line,23,15);
          my $ygrad=substr($line,38,15);
          my $zgrad=substr($line,53,15);
          #print "debug parsed to '$xgrad' '$ygrad' '$zgrad'\n";
          $xgrad=$xgrad*$conv_hartee2kjmol/$conv_bohr2angs; # convert to kJ/mol/Angs
          $ygrad=$ygrad*$conv_hartee2kjmol/$conv_bohr2angs; # convert to kJ/mol/Angs
          $zgrad=$zgrad*$conv_hartee2kjmol/$conv_bohr2angs; # convert to kJ/mol/Angs
          # Gaussian reports force not gradient, change signs
          $xgrad = -$xgrad; $ygrad = -$ygrad; $zgrad = -$zgrad;
          print "GRADIENTTAG $xgrad $ygrad $zgrad\n";
      }
      
   }
   close(GAUSSIANLOG);
}
else {
   print "ERROR cannot find file $gaussianlog\n";
}


if ($gradnum == $numatoms) {
   print "Success! have got required number of grads\n";
   $command= "rm $gaussianinp; rm -rf $SCRDIR";
   print("clean up with system command:\n   $command\n");
   system($command);
   }
else {
   print "ERROR job failed try to cat the log file\n";
   system("cat gelly_gaussian.log");
   }



