#!/usr/bin/perl -s

our ($h,$help,$orca);

if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print "
  Perl script to transform hessian from PRIRODA/GAMESS/ORCA 
  to \$HESS group for GAMESS or to hess-file for ORCA (with option -orca).
  
         $program hessian_file > gamess.dat
         $program hessian_file >> gamess.inp
         $program -orca hessian_file > file.hess
  
  hessian_file may be (auto detected):
  - PRIRODA output after task=hessian
  - GAMESS output after runtyp=hessian
  - GAMESS punch (last \$HESS group)
  - ORCA hess-file 
  
  For ORCA, \$hessian block is extracted from all hessian_file, 
  \$atoms block only from GAMESS and PRIRODA outs.
  \n";
  exit;
}

our %ISOTOP = qw(
H  1.007947    He 4.0026022   Li 6.9412      Be 9.0121823   B  10.8117
C  12.01078    N  14.00672    O  15.99943    F  18.99840325 Ne 20.17976
Na 22.9897702  Mg 24.30506    Al 26.9815382  Si 28.08553    P  30.9737612
S  32.0655     Cl 35.4532     Ar 39.9481     K  39.09831    Ca 40.0784
Sc 44.9559108  Ti 47.8671     V  50.94151    Cr 51.99616    Mn 54.9380499
Fe 55.8452     Co 58.9332009  Ni 58.69342    Cu 63.5463     Zn 65.4094
Ga 69.7231     Ge 72.641      As 74.921602   Se 78.963      Br 79.9041
Kr 83.7982     Rb 85.46783    Sr 87.621      Y  88.905852   Zr 91.2242
Nb 92.906382   Mo 95.942      Tc 98          Ru 101.072     Rh 102.905502
Pd 106.421     Ag 107.86822   Cd 112.4118    In 114.8183    Sn 118.7107    
Sb 121.7601    Te 127.603     I  126.904473  Xe 131.2936    Cs 132.905452  
Ba 137.3277    La 138.90552   Ce 140.1161    Pr 140.907652  Nd 144.243     
Pm 145         Sm 150.363     Eu 151.9641    Gd 157.253     Tb 158.925342  
Dy 162.5001    Ho 164.930322  Er 167.2593    Tm 168.934212  Yb 173.043     
Lu 174.9671    Hf 178.492     Ta 180.94791   W  183.841     Re 186.2071    
Os 190.233     Ir 192.2173    Pt 195.0782    Au 196.966552  Hg 200.592     
Tl 204.38332   Pb 207.21      Bi 208.980382  Po 209         At 210         
Rn 222         Fr 223         Ra 226         Ac 227         Th 232.03811   
Pa 231.035882  U  238.028913
);
our @ATOM = qw (XX
            H                                                  He
            Li Be                               B  C  N  O  F  Ne
            Na Mg                               Al Si P  S  Cl Ar
            K  Ca Sc Ti V  Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
            Rb Sr Y  Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I  Xe
            Cs Ba La
                    Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb
                  Lu Hf Ta W  Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
            Fr Ra Ac
                    Th Pa U  Np Pu Am Cm Bk Cf Es Fm Md No
                  Lr Rf Db Sg Bh Hs Mt Ds Rg
           );
our $BA = 0.52917721;
  
my ($mol,$energy,$e_nuc,@hess) = read_HESS_file($ARGV[0]);
$energy ||= 0;
$e_nuc ||= 0;
my $N = @hess;

if ($orca) {
  print "\n\$orca_hessian_file\n\n";
  print "\$hessian\n";
  print "$N\n";
  my @num;
  # (0..8) => ([0..4],[5..8])
  foreach (0 .. $N-1) {
    push @num, [] if $_%5==0;
    push @{$num[-1]}, $_;
  }
  foreach my $r (@num) {
    print "  ";
    foreach (@$r) {
      printf "%19d", $_;
    }
    print "\n";
    for (my $i=0; $i<$N; $i++) {
      foreach my $j (@$r) {
        printf("%5d   ", $i) if $j == $r->[0];
        printf "%19.10E", $hess[$i][$j];
      }
      print "\n";
    }
  }
  print "\n";
  print "#\n# The atoms: label  mass x y z (in bohrs)\n#\n";
  print "\$atoms\n";
  for (my $i=1; $i<@$mol; $i++) {
    printf "%2s %12.5f %19.12f %18.12f %18.12f\n", 
    $mol->[$i][0],$ISOTOP{$mol->[$i][0]},
    $mol->[$i][1]/$BA, $mol->[$i][2]/$BA, $mol->[$i][3]/$BA;
  }
  print "\n";
}
else
{
  print " \$HESS\n";
  printf "ENERGY IS %19.10f E(NUC) IS %19.10f\n", $energy, $e_nuc;
  for (my $i=0; $i<$N; $i++) {
    my $k = 1;
    for (my $j=0; $j<$N; $j++) {
      printf("%2d%3d", ($i+1)%100, $k%100) if $j%5 == 0;
      printf "%15.8E", $hess[$i][$j];
      if ($j%5 == 4) {
        $k++;
        print "\n" if $j != $N-1;
      }
    }
    print "\n";
  }
  print " \$END\n";
}

sub read_HESS_file {
	my $file = shift;
  my $d = qr/-?\d+(?:\.\d+)?/;
  my $outtype = 'unknown';
  my $version = '';
  my @mol = ({});
  my $debug = 0;
  
  open PRI, $file or die "Can't open $file: $!\n";
  # determine the creator of the file
  while (<PRI>) {
    if (/Priroda(?:\s+version)?\s+(\d+)/) { # for Priroda v. 2-13
      $version = $1;
      last if eof(PRI); $_ = <PRI>;
      next unless (/^\s*copyright\s+\(c\)\s.*\sLaikov$/); # for Priroda v. 2-13
      $outtype = "priroda";
      last;
    }
    elsif (/^\s*\*\s+GAMESS\s+VERSION\s+=\s+(.+?)\s+\*/) {
      $version = $1;
      $outtype = 'gamessUS';
      last;
    }
    elsif (/^\s*\*\s+Firefly\s+version\s+([\d\.]+)/) {
      $version = $1;
      $outtype = 'Firefly';
      last;
    }
#    elsif (/^\s{8,}\*{30,}$/) {
#      last if eof(PRI); $_ = <PRI>;
#      next unless (/^\s*\*\s+GAMESS\s+VERSION\s+=\s+(.+?)\s+\*/);
#      $version = $1;
#      last if eof(PRI); $_ = <PRI>;
#      next unless (/^\s*\*\s+FROM\s+IOWA\s+STATE\s+UNIVERSITY\s+\*$/);
#      $outtype = 'gamessUS';
#      last;
#    }
    elsif (/\Q[Molden Format]/) {
      $outtype = 'molden';
      last;
    }
    elsif ($.==1 and m/^ \$DATA\s*$/) {
      $outtype = 'gamess_punch';
      last;
    }
    elsif (m/\$orca_hessian_file/) {
      $outtype = 'orca.hess';
      last;
    }
    
  }
  warn "$file: $outtype version $version\n" if $debug;
  
  my ($N,$energy,$e_nuc,@hess);
	
  if ($outtype eq 'priroda') {
    my @H;
    while (<PRI>) {
		  if (/Atomic Coordinates:/) {
        while (<PRI>) {
          last if m/#/;
          last if !/^\s+\d+\s+/;
          my ($atN,$x,$y,$z) = split;
          push @mol, [ $ATOM[$atN],$x,$y,$z ];
        }
      }
      !$N and m/Number of atoms\s+(\d+)/ and $N = $1;
      m/^eng> E=\s*(\S+)/ and $energy = $1;
      $mol->[0]{Energy} = $energy if $energy;
      if (m/^eng> H=/ .. m/^eng>\$end/) {
        last if m/^eng>\$end/;
        s/^eng> (?:H=)?//;
        push @H, split;
  		}
	  }
    #print "$N\n@H\n";
    die "Invalid hessian\n" if @H != 3*$N*(3*$N+1)/2;
    for (my $i=0; $i<3*$N; $i++) {
      for (my $j=0; $j<=$i; $j++) {
        $hess[$i][$j] = $hess[$j][$i] = shift @H;
      }
    }
	}
  
  elsif ($outtype eq 'gamessUS' or $outtype eq 'Firefly') {
	  #print "$outtype\n"; exit;
    LOOP:
    while (<PRI>) {
      m/TOTAL NUMBER OF ATOMS\s+=\s*($d)/ and $N = $1;
      m/NUCLEAR REPULSION ENERGY =\s+($d)/ and $e_nuc = $1;
      m/TOTAL ENERGY =\s+($d)/             and $energy = $1;
      if (/\$CONTRL OPTIONS/) {
        while (<PRI>) {
          if (/RUNTYP=(\w+)/) {
            die "Only RUNTYP=HESSIAN is supported\n" if $1 ne 'HESSIAN';
            last;
          }
        }
      }
      if (/COORDINATES \(BOHR\)/) {
        <PRI>;
        while (<PRI>) {
          last if m/^\s*$/;
          my ($at,$ch,$x,$y,$z) = split;
          push @mol, [ $ATOM[$ch],$x*$BA,$y*$BA,$z*$BA];
        }
      }
      $mol->[0]{Energy} = $energy if $energy;
      if (m/CARTESIAN FORCE CONSTANT MATRIX/) {
        <PRI>; <PRI>;
        while (<PRI>) {
          last LOOP if m/------/;
          #push @hess, $_;
          if (s/^\s*(\d+)\s+\w+\s+X//) {
            my $n = $1-1;
            push @{$hess[3*$n]}, unpack("A9" x (length($_)/9), $_);
            $_ = <PRI>; s/^\s+Y//;
            push @{$hess[3*$n+1]}, unpack("A9" x (length($_)/9), $_);
            $_ = <PRI>; s/^\s+Z//;
            push @{$hess[3*$n+2]}, unpack("A9" x (length($_)/9), $_);
          }
        }
      }
	  }
    for (my $i=0; $i<3*$N; $i++) {
      for (my $j=0; $j<=$i; $j++) {
        $hess[$j][$i] = $hess[$i][$j];
      }
    }
  }

  elsif ($outtype eq 'gamess_punch') {
    my @H;
    while (<PRI>) {
      if (m/^ \$HESS/) {
        @H = ();
        while (<PRI>) {
          last if m/\$END/;
          if (/ENERGY IS\s+($d) E\(NUC\) IS\s+($d)/) {
            $energy = $1;
            $e_nuc  = $2;
          }
          else {
            #next if /\$HESS/;
            #print 'x5'.("A15" x ((length($_)-5)/15)),"\n";
            push @H, unpack('x5'.("A15" x ((length($_)-5)/15)), $_);
          }
        }
  		}
	  }
    $N = sqrt(@H)/3;
    #print "$N\n@H\n";
    die "Invalid hessian\n" if $N !~ /^\d+$/;
    for (my $i=0; $i<3*$N; $i++) {
      $hess[$i] = [splice(@H,0,3*$N)];
    }
	}
  elsif ($outtype eq 'orca.hess') {
    my (@H,$NNN);
    while (<PRI>) {
      if (m/\$hessian/) {
        $_ = <PRI>;
        m/^\s*(\d+)\s*$/ ? $NNN = $1 : die "Can not find dimension of hessian\n";
        while (<PRI>) {
          last if m/^\s*$/;
          if (/^\s{15,}\d+/) {
            my @ind = split;
            for (my $i=0; $i<$NNN; $i++) {
              $_ = <PRI>;
              my @ar = split;
              my $ii = shift @ar;
              die "Invalid indexes if 1st column: $ii != $i\n" if $ii != $i;
              push @{$H[$i]}, @ar;
            }
          }
          else {
            die "Invalid indexes if string $_\n";
          }
        }
  		}
	  }
    @hess = @H;
	}
  
	close PRI;
  
  if ($debug) {
    warn "Hessian matrix:\n";
    for (my $i=0; $i<3*$N; $i++) {
      for (my $j=0; $j<=$i; $j++) {
        printf STDERR "%6.2f", $hess[$i][$j];
      }
      warn "\n";
    }
  }
  
  return (\@mol,$energy,$e_nuc,@hess);
}