#!/usr/local/bin/perl

# Script to extract nomrmal modes from the DMOL 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 DMOL is given on the command line.
#
# search for line 'df              ATOMIC  COORDINATES' and save its position
# in the file with tell.
# In the same time look for line 'Eigenvalues, occupations, and MO coefficients'
# and if found, reset file pointer to last ATOMIC  COORDINATES with seek.
# Read in coordinates after skipping 2 lines. Read coordinates until line
# contains 'self consistent binding energy'
# Now search for line: 'Frequencies (CM-1) and normal modes'
# Then collect modes which are given as:
#        Frequencies (CM-1) and normal modes
#      0.0      0.0      0.0      0.0      0.0      0.0    193.9    268.7    363.2    709.9
#
#   -0.112    0.147    0.015    0.409   -0.435   -0.483    0.068   -0.047   -0.171    0.066




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

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

#now search for coordinates line
while (<DMOLOUT>) {
    if(/ATOMIC  COORDINATES                        DERIVATIVES/) {
      $coor_pointer = tell(DMOLOUT);
      }
  last if /Eigenvalues, occupations, and MO coefficients:/;
  }

#reset file pointer to last COORDINATES
seek(DMOLOUT,$coor_pointer,0);

#skip 2 lines;
$line = <DMOLOUT>;
$line = <DMOLOUT>;

#now collect cartesian coordinates in an array
$n_at = 0;
while ($line = <DMOLOUT>)  {
  last if ($line =~ /self consistent binding energy/);
  $n_at++;
  $line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;

  $at_symb[$n_at] = $2;
  $x_coor[$n_at] = $3 * 0.529177249;    # to Angstroms
  $y_coor[$n_at] = $4 * 0.529177249;
  $z_coor[$n_at] = $5 * 0.529177249;
  }

#search for line: "Frequencies (CM-1) and normal modes"

while ($line = <DMOLOUT>) {
  last if ($line =~ /Frequencies \(CM-1\)/);
  }

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

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

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


  $i = 0;

  #get frequences
  $line = <DMOLOUT>;
  chop($line);
  @vib_freq = split(/[ \t]+/, $line);
  if($first_pass == 1) {  #skip first 6 modes
    for ($j = 1; $j <= 6; $j++) {
      shift(@vib_freq);
      }
    }

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

  $line = <DMOLOUT>;   #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]);
      }

    if($n_curr_cols >= 10) {
      printf(NU10F "  %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 = <DMOLOUT>;
      chop($line);
      @modes = split(/[ \t]+/, $line);

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

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

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

  if($n_curr_cols >= 10) {
    close(NU10F);
    }

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

  #skip 2 lines to take next batch
  $line = <DMOLOUT>;
  $line = <DMOLOUT>;

  $first_pass = 0;
     
  }  # while

close(DMOLOUT);