#!/usr/bin/perl -w
use strict;
# gelly qm helper to use MOPAC2009  
# Global Phasing Ltd 2012, NOT COPYRIGHTED! 
# author: Oliver Smart
# contact: buster-develop@GlobalPhasing.com
#
# you will need to download and obtain MOPAC2009 before using this script
# please visit
#
# http://openmopac.net/home.html
#
# 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 mopac executable is:
print "PROGRAM_TAG MOPAC2009 helper_15.03.12\n"; # will appear in final pdb
my $mopaccom="/opt/mopac/MOPAC2009.exe";
#
my $mlines;
my $line;
my $lcount=0;
my $numatoms=0;
my $gradnum=0;
my $crdin= "gelly_coords_for_qmhelper.crdin"
; # input file
my $mopacinp; # input file
my $conv_kcalmol2kjmol=4.18400; # http://mccammon.ucsd.edu/~dzhang/energy-unit-conv-table.html

print "\n";
print "gelly qm helper script to run MOPAC 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 = "PM6"; 
   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 mopac input file
my $random_number = int(rand(100000));
$mopacinp="gelly_mopac_".$$."_".$random_number.".mop";
open(MOPACINP,">".$mopacinp) || die "cannot open $mopacinp";
print "will now write out a mopac input file $mopacinp\n";
print MOPACINP "$method 1SCF GRADIENTS NOSYM CHARGE=$charge\n";
print MOPACINP "mopac input file written by gelly_qm_helper_mopac2009.pl\n";
print MOPACINP "All coordinates are Cartesian\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 MOPACINP "%d \(%s\) %f  1 %f 1 %f 1\n", $atnum, $atnam, $xcoord, $ycoord, $zcoord;
    }
}
close(CRDIN);
 
print MOPACINP "\n"; #need a couple of blank lines at the bottom
print MOPACINP "\n";
print "number of atoms converted $numatoms\n";

my $command= "$mopaccom $mopacinp ";
$command = "bash -c \"$command\"";
print "now run mopac with system command:\n   $command \n";

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

system($command);

# now check and parse output log  file
my $mopaclog= $mopacinp;
$mopaclog=~ s/\.mop$/.out/;
my ($gradx, $grady, $gradz);
if (-r $mopaclog) {
   print "Good have got output file $mopaclog\n";
   print "now try to parse out energy and gradients\n";
   open(MOPACLOG,$mopaclog) || die "cannot open $mopaclog";
   my $readgrad=0;
   while (<MOPACLOG>) {
      $line = $_;
      # energy from line like
      #          FINAL HEAT OF FORMATION =         37.28205 KCAL/MOL
      # AM1 energy look for line like
      #  Energy=   -0.084678594320 NIter=  14.
      if ($line =~/ FINAL HEAT OF FORMATION = *(.*) KCAL/) {
          my $energy=$1;
          print "read energy= $energy kcal/mol =   "; 
          $energy=$energy*$conv_kcalmol2kjmol; # convert to kJ/mol
          print                               " $energy kJ/mol\n"; 
          print "ENERGYTAG   $energy\n";
      }
      # gradients from section like
      #   PARAMETER     ATOM    TYPE            VALUE       GRADIENT
      #      1          1  C    CARTESIAN X    -0.000004    -0.026764  KCAL/ANGSTROM
      #      2          1  C    CARTESIAN Y    -0.104938    -0.019835  KCAL/ANGSTROM
      #      3          1  C    CARTESIAN Z     0.202571    -0.005600  KCAL/ANGSTROM
      elsif ($line =~/.{25}CARTESIAN X.* (\S*)  KCAL\/ANGSTROM/) {
          $gradx=$1;
      }
      elsif ($line =~/.{25}CARTESIAN Y.* (\S*)  KCAL\/ANGSTROM/) {
          $grady=$1;
      }
      elsif ($line =~/.{25}CARTESIAN Z.* (\S*)  KCAL\/ANGSTROM/) {
          $gradz=$1;
          $gradnum++;
          $gradx *= $conv_kcalmol2kjmol; # convert to kJ/molAngs
          $grady *= $conv_kcalmol2kjmol; # convert to kJ/molAngs
          $gradz *= $conv_kcalmol2kjmol; # convert to kJ/molAngs
          print "GRADIENTTAG $gradx $grady $gradz\n";
      }
   }
   close(MOPACLOG);
}
else {
   print "ERROR cannot find file $mopaclog\n";
}


if ($gradnum == $numatoms) {
   print "Success! have got required number of grads\n";
   my $clean=1;
   if ($clean) {
      $command= "rm $mopacinp";
      $command=~  s/\.mop$/.\*/;
      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 $mopaclog");
   }
  
