#!/usr/bin/perl
#
# This script reads a Gaussian log file for NBO output and constructs
# a NBO label file called <PREFIX>.46 suitable for reading by NBO
# visualizers such as xorbplot and NBOView.
#
# JPH 3rd April 2014
#

use strict;
# use warnings;
use Getopt::Long;

# Set up some vars.
#
my $help;
my $version;
my $blank;
my $lang;
my $orb_type;
my $mol_unit;
my $i;
my $logfile = '';
my $prefix = 'FILE';
my @lines;
my @items;
my @ao_array;
my @nao_array;
my @nbo_array;
my @nho_array;
my @X;
my $ao=0;
my $nao=0;
my $nbo=0;
my $nho=0;
my $num_ao;
my $num_nao;
my $num_nbo;
my $num_nho;
my $m;
my $at1;
my $at2;
my $orb_find=0;
my $inp;

#
# Now read the arguments:
#
GetOptions( 'prefix=s'  =>  \$prefix,       # --prefix=<string>
            'logfile=s' =>  \$logfile,      # --logfile=<string>
            'version'   =>  \$version,      # --version
            'help'      =>  \$help);        # --help
#
#

print "=============================================================\n";
print " NBOLABEL V 0.3: Please report bugs to jerry.hagon\@ncl.ac.uk\n";
print "=============================================================\n";
if ($version) { exit 0; }

if ($help) {
print<<EOF;

USAGE:

  nbolabel [opt] --logfile=<gaussian log file>

OPTIONS:

  --help          Print this help message.

  --version       Print version number and banner then exit.

  --prefix=<output file prefix> (default is FILE).

EXAMPLE:

  nbolabel --prefix=<prefix_string> --logfile=ammonia.log

EOF
exit(0);
}

unless ($logfile) { die "No input log file specified!" };

# Read input log file.
#
open(LOGFILE, $logfile) or die "Can't open '$logfile': $!";
@lines = <LOGFILE>;

# Set up output file.
#
if ( -e "$prefix.46" ) {
  print "$prefix.46 exists! Overwrite?(y/n): ";
  $inp = <>;
  chomp($inp);
  unless ( $inp =~ /Y|y/ ) { exit 0; }
  open (OUT, ">$prefix.46");
} else {
  open (OUT, ">$prefix.46");
}
  



print "Searching for orbitals...\n";


# Find line numbers where AO Orbital information begins;
#

for $i (0 .. $#lines) {
  if (grep { /^\s*AO  Atom No lang/ } $lines[$i]) {
    print "... AO Orbitals Found!\n";
    $orb_find=1;
    $ao=$i;
  }
}

if ($ao) {
  $num_ao=0;
  $blank=0;
  for $i ($ao+2 .. $#lines) {
    if ($lines[$i] =~/^$/) { last; }
    if ($lines[$i] =~ /^\s{10}/) {
      if ($blank == 1) { last; }
      $blank=1;
      next;
    } else {
      $num_ao++;
      $blank=0;
      @items = split(' ',$lines[$i]);
      $items[5]=~/(\d)/;
      $lang = $1.lc($items[3]);
      $ao_array[$items[0]]=sprintf("%2s%2s( %-3s)", $items[1], $items[2], $lang);
    }
  }
  @X=@ao_array;
  shift @X;
  my $str_num_ao = sprintf("   AO%10s\n",$num_ao);
  print OUT "$str_num_ao";
  $m = 7;
  print OUT " ", join("", splice @X, 0, $m), "\n" while @X;

} else {
  print "... AO Orbitals Not Found!\n";
}

# Find line numbers where NAO Orbital information begins;
#

for $i (0 .. $#lines) {
  if (grep { /NATURAL POPULATIONS:\s*Natural\s*atomic\s*orbital\s*occupancies/ } $lines[$i]) {
    print "...NAO Orbitals Found!\n";
    $orb_find=1;
    $nao=$i;
  }
}

if ($nao) {
  $num_nao=0;
  $blank=0;
  for my $i ($nao+4 .. $#lines) {
    if ($lines[$i] =~ /^\s*$|^\s*\[|^\s*[A-Z,a-z]/) {
      if ($blank == 1) { last; }
      $blank=1;
      next;
    } else {
      $num_nao++;
      $blank=0;
      @items = split(' ',$lines[$i]);
      $items[5]=~/(\d)/;
      $lang = $1.lc($items[3]);
      if ( length($items[2]) > 2 ) { $items[2] = substr( $items[2], 1, 2 ); }
      if ( length($lang) > 3 ) {
        $lang = substr( $lang, 0, 2 ).'~';
      }
      $nao_array[$items[0]]=sprintf("%2s%2s( %-3s)", $items[1], $items[2], $lang);
    }
  }
  @X=@nao_array;
  shift @X;
  my $str_num_nao = sprintf("  NAO%10s\n",$num_nao);
  print OUT "$str_num_nao";
  $m = 7;
  print OUT " ", join("", splice @X, 0, $m), "\n" while @X;
 
} else {
  print "...NAO Orbitals Not Found!\n";
}

# Find line numbers where NBO Orbital information begins;
#

for my $i (0 .. $#lines) {
  if (grep { /Molecular unit/ } $lines[$i]) {
    print "...NBO Orbitals Found!\n"; 
    $orb_find=1;
    $nbo=$i;
    last;
  }
}

if ($nbo) {
  $num_nbo=0;
  for $i ($nbo+1 .. $#lines) {
    if ($lines[$i] =~ /^\s{20}|^\s*\-\-\-|^\s*[A-Za-z]/) { next; }
    if ($lines[$i] =~ /^\s*$/) {
      if ($blank == 1) { last; }
      $blank=1;
      next;
    } else {
      $num_nbo++;
      $blank=0;
      @items = split(' ',$lines[$i]);
      $items[0]=~ s/\.$//;
      if ($lines[$i] =~ /BD \(\s*\d+\)\s*(\w+)\s+(\d+)\s*\-\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $at2 = $4;
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        if ( length($4) > 2 ) { $at2 = substr( $4, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%2s%2s\-%2s%2s ", $1, $at1, $3, $at2);
      }
      if ($lines[$i] =~ /CR \(\s*\d+\)\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $orb_type = 'cr';
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%3s%2s(%s) ", $1, $at1, $orb_type);
      }
      if ($lines[$i] =~ /LP \(\s*\d+\)\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $orb_type = 'lp';
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%3s%2s(%s) ", $1, $at1, $orb_type);
      }
      if ($lines[$i] =~ /LP\*\(\s*\d+\)\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $orb_type = 'lp*';
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%2s%2s(%s) ", $1, $at1, $orb_type);
      }
      if ($lines[$i] =~ /RY\*\(\s*\d+\)\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $orb_type = 'ry*';
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%2s%2s(%s) ", $1, $at1, $orb_type);
      }
      if ($lines[$i] =~ /BD\*\(\s*\d+\)\s*(\w+)\s+(\d+)\s*\-\s*(\w+)\s+(\d+)/) {
        $at1 = $2;
        $at2 = $4;
        if ( length($2) > 2 ) { $at1 = substr( $2, 1, 2 ); }
        if ( length($4) > 2 ) { $at2 = substr( $4, 1, 2 ); }
        $nbo_array[$items[0]]=sprintf("%2s%2s\-%2s%2s*", $1, $at1, $3, $at2);
      }
    }
  }
  @X=@nbo_array;
  shift @X;
  my $str_num_nbo = sprintf("  NBO%10s\n",$num_nbo);
  print OUT "$str_num_nbo";
  $m = 7;
  print OUT " ", join("", splice @X, 0, $m), "\n" while @X;

} else {
  print "...NBO Orbitals Not Found!\n";
}

# Find line numbers where NHO Orbital information begins;
#

for my $i (0 .. $#lines) {
  if (grep { /^\s*NHO Fock matrix/ } $lines[$i]) {
    print "...NHO Orbitals Found!\n";
    $orb_find=1;
    $nho=$i;
  }
}

if ($nho) {
  $num_nho=0;
  for $i ($nho+4 .. $#lines) {
    if ($lines[$i] =~ /^\s*$/) {
      last;
    } else {
      $num_nho++;
      $blank=0;
      @items = split(' ',$lines[$i]);
      $items[0]=~ s/\.$//;
      if ($lines[$i] =~ /\s*\d+\. (.{10})/) {
        if ( length($4) > 2 ) { $at2 = substr( $4, 1, 2 ); }
        $nho_array[$items[0]]=sprintf("%10s", $1);
      }
    }
  }
  @X=@nho_array;
  shift @X;
  my $str_num_nho = sprintf("  NHO%10s\n",$num_nho);
  print OUT "$str_num_nho";
  $m = 7;
  print OUT " ", join("", splice @X, 0, $m), "\n" while @X;

} else {
  print "...NHO Orbitals Not Found!\n";
}

unless ($orb_find) { print "No AO, NAO, NBO or NHO orbitals found!\n"; }
close(OUT);

exit 0;