#!/usr/local/bin/perl

# Script to extract nomrmal modes from the DGAUSS output file.
# 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 DGAUSS is given on the command line.
#
# search for line '---------- Optimization has converged  ----------'
# search for line 'Final Coordinates File: DOptx'
# next line is no of atoms
# next n_at lines are:
#  At. numb    X  Y  Z (in angastroms)
# search for line 'Mass weighted Cartesian displacement vectors'
# what follows has a following structure:
#
#          Mass weighted Cartesian displacement vectors
#          Vibrational frequencies  cm-1
#
#
#           1       2       3       4       5       6       7       8       9
#
#          0.00    0.00    0.00    0.00    0.00    0.00  274.20  277.21  410.79
#
#   1    0.0000  0.0000  0.0000  0.0000  0.0000  0.2262 -0.0174 -0.0101  0.0015
#   2   -0.0061  0.0375 -0.1292  0.0503 -0.0626  0.0314 -0.0316 -0.0150 -0.0987
#   3   -0.0820 -0.0799  0.0230  0.0333 -0.1856  0.0063  0.0115  0.0028  0.0619
#   4    0.0299 -0.1387 -0.0505 -0.0420  0.0616  0.0741  0.0628 -0.0249  0.0649
#   5    0.0100 -0.0272 -0.1520  0.0302 -0.0423 -0.0398  0.0124 -0.0182 -0.0391
#   6   -0.1327 -0.0971 -0.0034  0.0411  0.0279 -0.0052  0.0169 -0.0347 -0.0831
#   7   -0.0595 -0.2769 -0.1644 -0.1831  0.1866  0.0080  0.1378 -0.0920  0.0467
#      -------- etc----
#
#  i.e., there is a row of frequences number, then new line, then frequences,
#  then columns of normal modes (3*N rows, 9 columns)

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

$DGOUT = $ARGV[0];
$DGOUT =~ /([\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(DGOUT,"<$DGOUT") || die "Could not open $DGOUT\n";

while ($_ = <DGOUT>) {
  last if ($_ =~ /---------- Optimization has converged  ----------/);
  }

#now search for "Final Coordinates File: DOptx" line
while (<DGOUT>) {
  last if /Final Coordinates File: DOptx/;
  }


#now collect cartesian coordinates in an array
$n_at = <DGOUT>;
for ($i = 1; $i <= $n_at; $i++) {
  $line = <DGOUT>;
  $line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;
  if($1 < $#at_symbols) {
    $at_symb[$i] = $at_symbols[$1-1];
    }
  else {
    $at_symb[$i] = $at_symbols[$#at_symbols];
    }
  $at_num[$i] = $1;
  $x_coor[$i] = $2;
  $y_coor[$i] = $3;
  $z_coor[$i] = $4;
  }

#search for line: "Low frequencies ---" which starts frequences block
while ($line = <DGOUT>) {
  last if ($line =~ /Mass weighted Cartesian displacement vectors/);
  }

#skip 5 lines
for($i = 1; $i <= 5; $i++) {
  $line = <DGOUT>;
  }

$n_cols = 3*$n_at-6;  #will be decremented by number of processed modes
$n_col_num = 0;
$first_pass = 1;
while($n_cols > 0) {
  if($n_cols >= 9) {
    $n_curr_cols = 9;
    }
  else {
    $n_curr_cols = $n_cols;
    }

  if($first_pass == 1) {   #skip first 6 modes
    $n_curr_cols = $n_curr_cols - 6;
    }

  #get frequences
  $line = <DGOUT>;
  chop($line);
  @vib_freq = split(/[ \t]+/, $line);

  if($first_pass == 1) {   #skip first 6 modes
    for($i = 1; $i <= 6; $i++) {
      shift(@vib_freq);
      }
    }

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

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

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

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

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

  if($n_curr_cols >= 1) {
    printf NU1F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[1];
    }
  if($n_curr_cols >= 2) {
    printf NU2F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[2];
    }
  if($n_curr_cols >= 3) {
    printf NU3F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[3];
    }
  if($n_curr_cols >= 4) {
    printf NU4F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[4];
    }
  if($n_curr_cols >= 5) {
    printf NU5F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[5];
    }
  if($n_curr_cols >= 6) {
    printf NU6F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[6];
    }
  if($n_curr_cols >= 7) {
    printf NU7F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[7];
    }
  if($n_curr_cols >= 8) {
    printf NU8F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[8];
    }
  if($n_curr_cols >= 9) {
    printf NU9F "%d\nNUDGAUSS=%f\n",$n_at, $vib_freq[9];
    }

  $line = <DGOUT>;   #skip one line;

  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]);
      }

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

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

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

    if($n_curr_cols >= 9) {
      printf(NU9F "  %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 = <DGOUT>;
      chop($line);
      @modes = split(/[ \t]+/, $line);
      shift(@modes);   # skip numbering
      if($first_pass == 1) {   #skip first 6 modes
        for ($j = 1; $j <= 6; $j++) {
          shift(@modes);
          }
        }

      if($n_curr_cols >= 1) {
        printf NU1F "  %10.4f", $modes[1];
        }

      if($n_curr_cols >= 2) {
        printf NU2F "  %10.4f", $modes[2];
        }

      if($n_curr_cols >= 3) {
        printf NU3F "  %10.4f", $modes[3];
        }

      if($n_curr_cols >= 4) {
        printf NU4F "  %10.4f", $modes[4];
        }

      if($n_curr_cols >= 5) {
        printf NU5F "  %10.4f", $modes[5];
        }

      if($n_curr_cols >= 6) {
        printf NU6F "  %10.4f", $modes[6];
        }

      if($n_curr_cols >= 7) {
        printf NU7F "  %10.4f", $modes[7];
        }

      if($n_curr_cols >= 8) {
        printf NU8F "  %10.4f", $modes[8];
        }

      if($n_curr_cols >= 9) {
        printf NU9F "  %10.4f", $modes[9];
        }
      }  # 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";
      }

    if($n_curr_cols >= 6) {
      printf NU6F "\n";
      }

    if($n_curr_cols >= 7) {
      printf NU7F "\n";
      }

    if($n_curr_cols >= 8) {
      printf NU8F "\n";
      }

    if($n_curr_cols >= 9) {
      printf NU9F "\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);
    }

  if($n_curr_cols >= 6) {
    close(NU6F);
    }

  if($n_curr_cols >= 7) {
    close(NU7F);
    }

  if($n_curr_cols >= 8) {
    close(NU8F);
    }

  if($n_curr_cols >= 9) {
    close(NU9F);
    }

  $n_col_num += $n_curr_cols;
  $n_cols -= $n_curr_cols;

  #skip 3 lines to take next batch
  $line = <DGOUT>;
  $line = <DGOUT>;
  $line = <DGOUT>;

  $first_pass = 0;
     
  }  # while

close(DGOUT);