#!/usr/local/bin/perl

@atom = ("",
         "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", ""
        );

#$infile = "h2o_frq.out";
if ($#ARGV < 0){
    print "USAGE: freq.pl Gaussian-output [-d]\n";
    print "       -d   DOS mode (with CRLF)\n";
    exit;
}
$infile = "$ARGV[0]";
$outfn  = "$ARGV[0]";
$outfn  =~ s/[.].*$//;
$le = "\n";
if ($#ARGV == 1){
  if ($ARGV[1] =~ /-d/i){
    $le = "\r\n";
  }
}

$ct  = 0;
$mct = 0;
$vct = 0;
$fct = 0;
$gvr = 0;
open(LST, "$infile");
while(<LST>){
  chop;
  s/^ //;
  if (/---/){
    $ct++;
    next;
  }
  if ($ct == 3){
    $molname = $_;
    $ct++;
    next;
  } 
  if (/ Gaussian 86\(TM\) / && $gvr < 1986){ $gvr = 1986; next; }
  if (/ Gaussian 88\(TM\) / && $gvr < 1988){ $gvr = 1988; next; }
  if (/ Gaussian 90\(TM\) / && $gvr < 1990){ $gvr = 1990; next; }
  if (/ Gaussian 92\(TM\) / && $gvr < 1992){ $gvr = 1992; next; }
  if (/ Gaussian\(R\) 94 /  && $gvr < 1994){ $gvr = 1994; next; }
  if (/ Gaussian\(R\) 98 /  && $gvr < 1998){ $gvr = 1998; next; }
  if (/ Gaussian\(R\) 03 /  && $gvr < 2003){ $gvr = 2003; next; }
  if (/standard orientation/i){
    if ($mct != 0){
      $mct = 0;
    }
    $_ = <LST>;
    $_ = <LST>;
    $_ = <LST>;
    $_ = <LST>;
    while(<LST>){
      chop;
      s/^ //;
      last if (/^----/);
      @a = split(" ");
      $mol[$mct]  = $atom[$a[1]];
      if ($gvr < 1998){
        $molx[$mct] = $a[2];
        $moly[$mct] = $a[3];
        $molz[$mct] = $a[4];
      }
      else {
        $molx[$mct] = $a[3];
        $moly[$mct] = $a[4];
        $molz[$mct] = $a[5];
      }
      $mct++;
    }
  }
  if (/Harmonic frequencies/){
    $_ = <LST>;
    $_ = <LST>;
    if ($gvr > 1998){ $_ = <LST>; }
    while(<LST>){
      chop;
      s/^ //;
      last if (length($_) == 0);
      @c = ();
      @c = split(" ");
      $_ = <LST>; chop; s/^ //; @d1 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d2 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d3 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d4 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d5 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d6 = split(" ");
      $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d7 = split(" ");
      if ($gvr > 1998){ $_ = <LST>; chop; s/^ //; s/^[^-]+--//; @d8 = split(" "); }
      for($i = 0; $i < @c; $i++){
        $sym[$fct + $i] = $d1[$i];
        $frq[$fct + $i] = $d2[$i];
        $red[$fct + $i] = $d3[$i];
        $frc[$fct + $i] = $d4[$i];
        $ir[$fct + $i] = $d5[$i];
        $rom[$fct + $i] = $d6[$i];
        $dep[$fct + $i] = $d7[$i];
        if ($gvr > 1998){ $depu[$fct + $i] = $d8[$i]; }
      }
      $_ = <LST>;
      for ($j = 0; $j < $mct; $j++){
        $_ = <LST>; chop; s/^ //; @e = split(" ");
        for ($i = 0; $i < @c; $i++){
          $num = sprintf("%d.%d", $j, $fct + $i);
          $frqx{$num} = $e[2 + $i * 3];
          $frqy{$num} = $e[2 + $i * 3 + 1];
          $frqz{$num} = $e[2 + $i * 3 + 2];
        }
      }
      $fct += 3;
    }
  }
}
close(LST);
if ($mct == 0 || $fct == 0){
  print "ERROR: strange format.\n";
  exit;
}
$gvrs = sprintf("%04d", $gvr);
$gvrs =~ s/^[12][0-9]//;
print "Gaussian $gvrs output file.\n";
print "Molecular count = $mct\n";
for ($i = 0; $i < $mct; $i++){
  print "$mol[$i] : $molx[$i], $moly[$i], $molz[$i]\n";
}
print "Frequency mode = $fct\n";
$outfile = "$outfn";
$outfile .= "_table.txt";
open(LST, "> $outfile");
print LST "Num. ";
print LST "Sym. ";
print LST "Frequencies ";
print LST " Red.masses ";
print LST " Frc_consts ";
print LST "   IR_Inten ";
print LST "Raman_Activ ";
if ($gvr < 1998){
  print LST "    Depolar";
}
else {
  print LST " Depolar(P) ";
  print LST " Depolar(U)";
}
print LST "$le";
for ($i = 0; $i < $fct; $i++){
  printf LST "%04d ",$i+1;
  printf LST "%4s ",$sym[$i];
  printf LST "%11.4f ",$frq[$i]; 
  printf LST "%11.4f ",$red[$i]; 
  printf LST "%11.4f ",$frc[$i]; 
  printf LST "%11.4f ",$ir[$i]; 
  printf LST "%11.4f ",$rom[$i]; 
  if ($gvr < 1998){
    printf LST "%11.4f",$dep[$i]; 
  }
  else {
    printf LST "%11.4f ",$dep[$i]; 
    printf LST "%11.4f", $depu[$i]; 
  }
  printf LST "$le"; 
}
close(LST);
for ($i = 0; $i < $fct; $i++){
  $ii = $i + 1;
  $outfile = "$outfn";
  $outfile .= sprintf("_f%04d.xyz", $ii);
  open(LST, "> $outfile");
  print LST "$mct$le";
  print LST "$molname No.$ii $sym[$i] Freq=$frq[$i] IR=$ir[$i]$le$le";
  for ($j = 0; $j < $mct; $j++){
    printf LST "%-2s ",$mol[$j];
    print LST "$molx[$j] ";
    print LST "$moly[$j] ";
    print LST "$molz[$j] ";
    $num = sprintf("%d.%d", $j, $i);
    print LST "$frqx{$num} ";
    print LST "$frqy{$num} ";
    print LST "$frqz{$num}$le";
  } 
  close(LST);
}