#!/usr/bin/perl -w
use strict;
print "PROGRAM_TAG GAMESS (US) helper_17.05.11\n"; # will appear in final pdb
# gelly qm helper using GAMESS (US);
# 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 gamess is:
my $gamess_rungms="/home/osmart/gamess/gamess/rungms";
my $gamess_scr="/home/osmart/scr";
my $gamess_verno="01"; # argument $2 of rungms VERNO "revision number of the executable created by 'lked' step"
my $gamess_ncpus=8;    # argument $3 of rungms NCPUS "number of compute processes to be run"
#
my $scftyp="scftyp=rhf"; 
my $mlines;
my $line;
my $lcount=0;
my $numatoms=0;
my $gradnum=0;
my $crdin= "gelly_coords_for_qmhelper.crdin"
; # input file
my $gamesinp="gelly_games.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 GAMESS (US) for energy/gradient calculation\n";
print "You must request, download, compile/link, setup and test GAMESS separately!\n";
print "GAMESS web site: http://www.msg.ameslab.gov/gamess/\n";
print "\n";
print "helper script location $0\n";
(-x $gamess_rungms) || die "ERROR cannot find gamess rungms script $gamess_rungms\nYou must localize $0 script";
(-d $gamess_scr) ||    die "ERROR cannot find gamess scratch directory $gamess_scr\nYou must localize $0 script";

# 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";



# method stuff
$method = "HF/3-21G(d)_pcm" if ($method eq "AUTO"); # default AUTO method

# strip off any "_pcm" from method
my $pcm = 0;
if ($method =~ /(.*)_pcm/) {
  $method = $1;
  print "method includes _pcm so QM will use PCM (polarizable continuum model) for water.\n";
  $pcm = 1;
}

if ($method eq "AM1") {
  $mlines=" \$basis  gbasis=AM1  \$end\n"; 
  $gamess_ncpus = 1; # only one cpu allowed for semi-empirical
}
elsif ($method eq "PM3") {
  $mlines=" \$basis  gbasis=PM3  \$end\n"; 
  $gamess_ncpus = 1; # only one cpu allowed for semi-empirical
}
elsif ($method eq "MNDO") {
  $mlines=" \$basis  gbasis=MNDO  \$end\n"; 
  $gamess_ncpus = 1; # only one cpu allowed for semi-empirical
}
elsif ($method eq "HF/3-21G") {
  $mlines=" \$basis  gbasis=n21 ngauss=3  \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
elsif ($method eq "HF/3-21G(d)") {
  $mlines=" \$basis  gbasis=n21 ngauss=3 ndfunc=1 \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
elsif ($method eq "HF/6-31G(d)") {
  $mlines=" \$basis  gbasis=n31 ngauss=6 ndfunc=1  \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
elsif ($method eq "HF/6-31+G(2d,p)") {
  $mlines=" \$basis  gbasis=n31 ngauss=6 ndfunc=2 npfunc=1 diffsp=.t. \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
# small basis DFT - not a great idea
#elsif ($method eq "B3LYP/3-21G(d)") {
#  $scftyp="dfttyp=b3lyp"; 
#  $mlines=" \$basis  gbasis=n21 ngauss=3 ndfunc=1 \$end\n";
#  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
#}
elsif ($method eq "B3LYP/6-31G(d)") {
  $scftyp="dfttyp=b3lyp"; 
  $mlines=" \$basis  gbasis=n31 ngauss=6 ndfunc=1 \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
elsif ($method eq "B3LYP/6-31+G(2d,p)") {
  $scftyp="dfttyp=b3lyp"; 
  $mlines=" \$basis  gbasis=n31 ngauss=6 ndfunc=2 npfunc=1 diffsp=.t. \$end\n";
  $mlines=$mlines." \$scf dirscf=.t. \$end\n";
}
else {
  die "ERROR unrecognized method=$method\nERROR ";
}

if ($pcm) {
  $mlines=$mlines." \$pcm solvnt=water \$end\n";
  $method = $method. " with PCM water ";
}


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 "will use $gamess_ncpus cpus for calculation\n";
print "have opened crdin format coordinate file $crdin\n";

# open output file
my $random_number = int(rand(100000));
$gamesinp="gelly_gamess_".$$."_".$random_number;
open(GAMESSINP,">".$gamesinp.".inp") || die "cannot open $gamesinp.inp";
print "will now write out a gamess input file $gamesinp.inp\n";
print GAMESSINP "! gamess control file writen by $0\n";
print GAMESSINP " \$contrl $scftyp  ICHARG=$charge MULT=$multip runtyp=gradient  \$end\n";
print GAMESSINP " \$system mwords=55 timlim=99999.0   \$end\n";
print GAMESSINP $mlines;
print GAMESSINP " \$data \n";
print GAMESSINP " \$gelly helper\n";
print GAMESSINP " C1\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";
        printf GAMESSINP " %s %4.1f %f %f %f\n", $atnam, $atnum, $xcoord, $ycoord, $zcoord;
    }
}
close(CRDIN);
 
print GAMESSINP " \$end\n";
print "number of atoms converted $numatoms\n";

my $command= "$gamess_rungms $gamesinp $gamess_verno $gamess_ncpus   >& gelly_games.log";
$command = "bash -c \"$command\"";
print "now run gamess with system command:\n   $command \n";
system($command);

# now check to see whether the job has produced output .dat file
my $thisgamessdat= $gamess_scr."/".$gamesinp.".dat";
if (-r $thisgamessdat) {
   print "Good have got output file $thisgamessdat\n";
   print "now try to parse out energy and gradients\n";
   open(GAMESSOUT,$thisgamessdat) || die "cannot open $thisgamessdat";
   my $readgrad=0;
   while (<GAMESSOUT>) {
      $line = $_;
      if ($line =~/E=  (.*)  GMAX=   .*  GRMS=  .*/) {
          my $energy=$1;
          print "energy= $energy hartrees =   "; 
          $energy=$energy*$conv_hartee2kjmol; # convert to kJ/mol
          print                               " $energy kJ/mol\n"; 
          print "ENERGYTAG   $energy\n";
          $readgrad=1;
      }
      elsif ($line =~/^ \$END/) {
          $readgrad=0;
      }
      elsif ($readgrad) {
          $gradnum++;
          # gradient line 
          #print "read gradient line: $line";
          my $xgrad=substr($line,18,17);
          my $ygrad=substr($line,38,17);
          my $zgrad=substr($line,58,17);
          $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
          print "GRADIENTTAG $xgrad $ygrad $zgrad\n";
      }
   }
   close(GAMESSOUT);
}
else {
   print "ERROR cannot find file $thisgamessdat\n";
}

if ($gradnum == $numatoms) {
   print "Success! have got required number of grads\n";
   $command= "rm $gamesinp.inp $thisgamessdat";
   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_games.log");
   }
  
