#!/usr/local/bin/perl

# Script to extract nomrmal modes from the G90 output file.
# The OUT file may be the result of single calculations with
# link line (like example 6.15 in G90 manual) or a result of
# a separate restart with FREQ keyword from checkpoint file.
# The first task is to find out if OUT contains two runs
# (geometry optimization followed by frequences) or just one.
# Then if 2 runs, skip the portion of file up to and including
# first @ (it ends the result summary portion for geometry optimization).
# Then search for "Z-Matrix orientation:" string (if C1 symmetry
# and for "Standard orientation:" if symmetry used. Skip
# 4 lines and collect cartesian coordinates in the array
# until you reach line of "--------------". Count the atoms.
# Then search for the "Frequencies ---" string, which
# starts first portion of the table of frequences and
# normal modes. It looks like:
#       Frequencies ---   278.9971  314.8132  452.1797  864.0812  962.6147
#    Reduced masses ---     1.0702    1.1263    2.7343    1.0824    2.1573
#   Force constants ---     0.0491    0.0658    0.3294    0.4761    1.1778
#    IR Intensities ---    16.3380  137.0003   14.8384    2.4492    9.6849
#  Raman Activities ---     0.4115    2.3307    0.7819    0.2060   11.2770
#   Depolarizations ---     0.7439    0.7405    0.4521    0.6779    0.2744
# Coord Atom Element:
#   6     1     6          0.00129  -0.00502   0.07478  -0.03799  -0.05805
#   6     1     6         -0.00805   0.03163   0.11205   0.02176  -0.08881
#   6     1     6          0.00304  -0.01473  -0.15222   0.00036   0.19568
#  etc.
# The Coord column is probably the (frequency number - 1) for first column
# but beats me. There are 5 columns in each portion. Skip first column. Next is
# an atom number (X, Y, Z for each atom, correspond to succesive rows, I
# presume), then the atomic number), then normal modes. These need to be
# collected in arrays. The 2nd and 3rd column can be used for consistency
# checking. There should be 3N rows in each portion of the matrix.
# Then you search for another "Frequencies ---" and check if you have
# all 5 columns or if this is a last portion with less columns.
# There is a second version of this table in G90 output file, which
# lists 3 frequences and XYZ modes for each atom, but normal modes
# have only two digits after decimal point.
# Now how the data are printed? They are printed as required by xmol
# to display frequences. Each normal mode is printed as a separate file
# with an extension nu1 ... nu(3N-6). The root of file name is taken from
# as everything from the beginning of output file name from GAUSSIAN up to
# first period. The OUT file name from G90 is given on the command line.
#

die "You did not give G90 output file name as argument\n" if $#ARGV < 0;
die "You need only one parameter --- G90 output file name\n" if $#ARGV > 0;

$G90OUT = $ARGV[0];
$G90OUT =~ /([\w\/+#-]+)/;
$OUT_root = $1;

@at_symbols = ('H ', 'He',
               'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne',
               'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar',
               'K ', 'Ca',
               'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn',
                           'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'X ');



open(G90OUT,"<$G90OUT") || die "Could not open $G90OUT\n";

#check if 1 or 2 @ signes in the file. Check is symmetry
$n_runs = 0;
$symmetry = 1;   #default is symmetry
while(<G90OUT>) {
  $n_runs++ if /@/;
  $symmetry = 0 if /FULL\s+POINT\s+GROUP\s+C1/;
  }

die "No @ sign found in the G90 output file\n" if $n_runs == 0;

seek(G90OUT, 0, 1);  #clean EOF status
seek(G90OUT, 0, 0);  #rewind G90OUT file

if($n_runs == 2) {  #skipp after 1st @ if 2 runs
  while ($_ = <G90OUT>) {
    last if ($_ =~ /@/);
    }
  }

#now search for "Z-Matrix orientation:" line if no symmetry
#or Standard orientation: is symmetry
while (<G90OUT>) {
  if($symmetry == 0) {
    last if /Z-Matrix orientation:/;
    }
  else {
    last if /Standard orientation:/;
    }
  }

#skip 4 lines
for ($i = 0; $i < 4; $i++ ) {
  $line = <G90OUT>;
  }

#now collect cartesian coordinates in an array
$n_at = 0;
while ($line = <G90OUT>) {
  last if ($line =~ /--------------------/);  # when last line
  $n_at++;
  $line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;
  die "Unstandard format of G90 output file\n" if $1 != $n_at;
  if($2 < $#at_symbols) {
    $at_symb[$n_at] = $at_symbols[$2-1];
    }
  else {
    $at_symb[$n_at] = $at_symbols[$#at_symbols];
    }
  $at_num[$n_at] = $2;
  $x_coor[$n_at] = $3;
  $y_coor[$n_at] = $4;
  $z_coor[$n_at] = $5;
  }

#search for line: "Low frequencies ---" which starts frequences block
while ($line = <G90OUT>) {
  last if ($line =~ /Low frequencies ---/);
  }

$n_cols = 3*$n_at-6;  #will be decremented by number of processed modes
$n_col_num = 0;
while($n_cols > 0) {
  if($n_cols >= 5) {
    $n_curr_cols = 5;
    }
  else {
    $n_curr_cols = $n_cols;
    }
  #open files for modes in the xmol format
  if($n_curr_cols >= 1) {
    $k = $n_col_num + 1;    # Which frequency we are"
    $temp = sprintf("%d",$k);
    $NU1F = $OUT_root . '.nu' . $temp;
    open(NU1F,">$NU1F") || die "Could not open $NU1F\n";
    }

  #open files for modes in the xmol format
  if($n_curr_cols >= 2) {
    $k = $n_col_num + 2;    # Which frequency we are"
    $temp = sprintf("%d",$k);
    $NU2F = $OUT_root . '.nu' . $temp;
    open(NU2F,">$NU2F") || die "Could not open $NU2F\n";
    }

  #open files for modes in the xmol format
  if($n_curr_cols >= 3) {
    $k = $n_col_num + 3;    # Which frequency we are"
    $temp = sprintf("%d",$k);
    $NU3F = $OUT_root . '.nu' . $temp;
    open(NU3F,">$NU3F") || die "Could not open $NU3F\n";
    }

  #open files for modes in the xmol format
  if($n_curr_cols >= 4) {
    $k = $n_col_num + 4;    # Which frequency we are"
    $temp = sprintf("%d",$k);
    $NU4F = $OUT_root . '.nu' . $temp;
    open(NU4F,">$NU4F") || die "Could not open $NU4F\n";
    }

  #open files for modes in the xmol format
  if($n_curr_cols >= 5) {
    $k = $n_col_num + 5;    # Which frequency we are"
    $temp = sprintf("%d",$k);
    $NU5F = $OUT_root . '.nu' . $temp;
    open(NU5F,">$NU5F") || die "Could not open $NU5F\n";
    }

  $search_expr = 'Frequencies ---' . ('\s+(\S+)' x $n_curr_cols);
  $i = 0;
  while ($line = <G90OUT>) {
    $i++;
    if($i > 10) {
      die "Could not find Frequencies --- in OUT file\n";
      }
    if($line =~ /$search_expr/) {
      if($n_curr_cols >= 1) {
        $vib_freq[1] = $1; 
        printf NU1F "%d\nNU=%f\n",$n_at, $vib_freq[1];
        }
      if($n_curr_cols >= 2) {
        $vib_freq[2] = $2;
        printf NU2F "%d\nNU=%f\n",$n_at, $vib_freq[2];
        }
      if($n_curr_cols >= 3) {
        $vib_freq[3] = $3;
        printf NU3F "%d\nNU=%f\n",$n_at, $vib_freq[3];
        }
      if($n_curr_cols >= 4) {
        $vib_freq[4] = $4;
        printf NU4F "%d\nNU=%f\n",$n_at, $vib_freq[4];
        }
      if($n_curr_cols >= 5) {
        $vib_freq[5] = $5;
        printf NU5F "%d\nNU=%f\n",$n_at, $vib_freq[5];
        }
      last;
      }
    }

  $i = 0;
  while ($line = <G90OUT>) {
    $i++;
    if($i > 10) {
      die "Could not find Coord Atom Element: in OUT file\n";
      }
    if($line =~ /Coord Atom Element:/) {
      last;
      }
    }

  $search_expr = '\S+\s+(\S+)\s+(\S+)' . ('\s+(\S+)' x $n_curr_cols);

  for ($n = 1; $n <= $n_at; $n++)  {
    if($n_curr_cols >= 1) {
      printf(NU1F "  %2s  %11.5f %11.5f %11.5f", $at_symb[$n],
             $x_coor[$n], $y_coor[$n], $z_coor[$n]);
      }

    if($n_curr_cols >= 2) {
      printf(NU2F "  %2s  %11.5f %11.5f %11.5f", $at_symb[$n],
             $x_coor[$n], $y_coor[$n], $z_coor[$n]);
      }

    if($n_curr_cols >= 3) {
      printf(NU3F "  %2s  %11.5f %11.5f %11.5f", $at_symb[$n],
             $x_coor[$n], $y_coor[$n], $z_coor[$n]);
      }

    if($n_curr_cols >= 4) {
      printf(NU4F "  %2s  %11.5f %11.5f %11.5f", $at_symb[$n],
             $x_coor[$n], $y_coor[$n], $z_coor[$n]);
      }

    if($n_curr_cols >= 5) {
      printf(NU5F "  %2s  %11.5f %11.5f %11.5f", $at_symb[$n],
             $x_coor[$n], $y_coor[$n], $z_coor[$n]);
      }

    for ($i = 1; $i <= 3; $i++) {
      $line = <G90OUT>;
      if( !($line =~ $search_expr)) {
	  die "Unstandard format of G90 OUTPUT file\n";
        }

      if($n != $1) {
        die "Unstandard format of G90 OUTPUT file --- mode line\n";
        }

      if($at_num[$n] != $2) {
        die "Unstandard format of G90 OUTPUT file --- at. num. wrong\n";
        }

      if($n_curr_cols >= 1) {
        $f =  $3;
        printf NU1F "  %11.5f", $f;
        }
      
    
      if($n_curr_cols >= 2) {
        $f =  $4;
        printf NU2F "  %11.5f", $f;
        }
      
    
      if($n_curr_cols >= 3) {
        $f =  $5;
        printf NU3F "  %11.5f", $f;
        }
      
    
      if($n_curr_cols >= 4) {
        $f =  $6;
        printf NU4F "  %11.5f", $f;
        }

      if($n_curr_cols >= 5) {
        $f = $7;
        printf NU5F "  %11.5f", $f;
        }
      }  # for $i

    if($n_curr_cols >= 1) {
      printf NU1F "\n";
      }

    if($n_curr_cols >= 2) {
      printf NU2F "\n";
      }

    if($n_curr_cols >= 3) {
      printf NU3F "\n";
      }

    if($n_curr_cols >= 4) {
      printf NU4F "\n";
      }

    if($n_curr_cols >= 5) {
      printf NU5F "\n";
      }

    }  # for $n

  if($n_curr_cols >= 1) {
    close(NU1F);
    }

  if($n_curr_cols >= 2) {
    close(NU2F);
    }

  if($n_curr_cols >= 3) {
    close(NU3F);
    }

  if($n_curr_cols >= 4) {
    close(NU4F);
    }

  if($n_curr_cols >= 5) {
    close(NU5F);
    }

  $n_col_num += $n_curr_cols;
  $n_cols -= $n_curr_cols;
  }  # while

close(G90OUT);