#!/usr/local/bin/perl
# 04.11.92
# Martin Schuetz, Institute of physical chemistry,
# University of Berne, Switzerland
# script to extract normal modes from Gaussian90 or CADPAC 4.0 frequency
# job outputfiles. A set of 3n-6 input files (G90) or 3n input files (CADPAC)
# are created for the SCHAKAL88B program of Egbert Keller, corresponding
# to 3n-6 normal modes (or 3n, if translations anr rotationas are included as
# in CADPAC). The displacement vectors of the normal modes are visualized
# in SCHAKAL88 as sticks, pointing to dummy atoms using a special bond type.

die "You did not give G90 or CADPAC output file name as argument\n" if $#ARGV < 0;
die "You need at most two parameters: filename scalingfactor\n" if $#ARGV > 1;

$GCOUT = $ARGV[0];
$GCOUT =~ /([\w\/+#-]+)/;
$scal_fac = 1.0;
if ($#ARGV == 1) {
  $scal_fac = $ARGV[1];
  }
$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 ');

@at_radii   = ( 1.20, 1.15,
                1.60, 1.60, 1.65, 1.60, 1.50, 1.40, 1.35, 1.30,
                2.00, 2.00, 2.05, 2.00, 1.90, 1.85, 1.80, 1.75,
                2.00, 2.00,
                2.05, 2.10, 2.20, 2.20, 2.15, 2.10,2.10, 2.10, 2.10,
                            2.15, 2.10, 2.00, 2.00, 1.95, 1.90, 0.01 );
@at_colors  = ( 0122, 5122, 
                2123, 3124, 1001, 1026, 3005, 2005, 3064, 1063,
                7064, 1001, 1001, 1064, 3026, 5004, 3025, 2003,
                4001, 1001,
                1001, 4001, 3001, 7001, 4001, 4001, 3001, 2001, 1001,
                            1001, 1001, 2026, 2026, 2026, 1001, 8001 );
                
@at_weights = ( 1, 4,
                7, 9, 11, 12, 14, 16, 19, 20,
                23, 24, 27, 28, 31, 32, 35, 40 );

$at_dark = 1003;
$at_scalef = 0.2; # internal scaling factor of Schakal

@ab_initio = ('unknown', 'Gaussian', 'CADPAC');

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

# check if Gaussian or CADPAC
$outf_type = do CHECK();
die "$GCOUT is neither Gaussian90 nor CADPAC outputfile" if !$outf_type;
printf "Hmm, this is a @ab_initio[$outf_type] output file\n";

# process output file...
if ($outf_type == 1) { do process_G90(); }
if ($outf_type == 2) { do process_CADPAC(); }
close(GCOUT);
printf "all done\n";

#
#
# -- subroutines
#
#

sub CHECK {
  local($chk) = 0;
  while (<GCOUT>) {
    if (/Gaussian/) { $chk = 1; }
    if (/CAMBRIDGE ANALYTIC DERIVATIVES PACKAGE/) { $chk = 2; }
    last if ($chk != 0);
    }
  $chk;
  }


sub SphereOffs {
  local($x1,$y1,$z1,$dx,$dy,$dz,$rad) = @_;
  local($norm) = sqrt($dx*$dx + $dy*$dy + $dz*$dz);
  if ($norm <= 0.0) {
    $norm = 1.0
  }
  $rad *= $at_scalef;
  local($ndx) = $dx/$norm;
  local($ndy) = $dy/$norm;
  local($ndz) = $dz/$norm;
  local($newpos[0]) = $x1 + $rad*$ndx + $dx;
  local($newpos[1]) = $y1 + $rad*$ndy + $dy;
  local($newpos[2]) = $z1 + $rad*$ndz + $dz;
  @newpos;
  }


sub process_G90 {
  #check if 1 or 2 @ signes in the file. Check is symmetry
  $n_runs = 0;
  $symmetry = 1;   #default is symmetry
  while(<GCOUT>) {
    $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(GCOUT, 0, 1);  #clean EOF status
  seek(GCOUT, 0, 0);  #rewind G90OUT file

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

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

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

  #now collect cartesian coordinates in an array
  $n_at = 0;
  while ($line = <GCOUT>) {
    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];
      $at_rad[$n_at] = $at_radii[$2-1];
      $at_col[$n_at] = $at_colors[$2-1];
      }
    else {
      $at_symb[$n_at] = $at_symbols[$#at_symbols];
      $at_rad[$n_at] = $at_radii[$#at_symbols];
      $at_col[$n_at] = $at_colors[$#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 = <GCOUT>) {
    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 schakal 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 schakal 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 schakal 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 schakal 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 schakal 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 = <GCOUT>) {
      $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 "TITL NU=%fcm-1\n", $vib_freq[1];
          }
        if($n_curr_cols >= 2) {
          $vib_freq[2] = $2;
          printf NU2F "TITL NU=%fcm-1\n", $vib_freq[2];
          }
        if($n_curr_cols >= 3) {
          $vib_freq[3] = $3;
          printf NU3F "TITL NU=%fcm-1\n", $vib_freq[3];
          }
        if($n_curr_cols >= 4) {
          $vib_freq[4] = $4;
          printf NU4F "TITL NU=%fcm-1\n", $vib_freq[4];
          }
        if($n_curr_cols >= 5) {
          $vib_freq[5] = $5;
          printf NU5F "TITL NU=%fcm-1\n", $vib_freq[5];
          }
        last;
        }
      }

    $i = 0;
    while ($line = <GCOUT>) {
      $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 "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n
",
               $at_symb[$n], $x_coor[$n], $y_coor[$n], $z_coor[$n],
               $at_rad[$n], $at_col[$n], $at_dark);
        }

      if($n_curr_cols >= 2) {
        printf(NU2F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n
",
               $at_symb[$n], $x_coor[$n], $y_coor[$n], $z_coor[$n],
               $at_rad[$n], $at_col[$n], $at_dark);
        }

      if($n_curr_cols >= 3) {
        printf(NU3F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n
",
               $at_symb[$n], $x_coor[$n], $y_coor[$n], $z_coor[$n],
               $at_rad[$n], $at_col[$n], $at_dark);
        }

      if($n_curr_cols >= 4) {
        printf(NU4F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
               $at_symb[$n], $x_coor[$n], $y_coor[$n], $z_coor[$n],
               $at_rad[$n], $at_col[$n], $at_dark);
        }

      if($n_curr_cols >= 5) {
        printf(NU5F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n
",
               $at_symb[$n], $x_coor[$n], $y_coor[$n], $z_coor[$n],
               $at_rad[$n], $at_col[$n], $at_dark);
        }

      for ($k = $n + 1; $k <= $n_at; $k++)  {  # loop for BONDS

        $dx = $x_coor[$n] - $x_coor[$k];
        $dy = $y_coor[$n] - $y_coor[$k];
        $dz = $z_coor[$n] - $z_coor[$k];
        $cut_off = 2.5;
        if (($at_symb[$n] eq "H ") || ($at_symb[$k] eq "H ")) {
          $cut_off = 1.2;
          }

        if ($dx*$dx + $dy* $dy + $dz * $dz < $cut_off) {

          if($n_curr_cols >= 1) {
            printf(NU1F "CONN     %03d%03d    6004    .068\n", $n, $k);
            }

          if($n_curr_cols >= 2) {
            printf(NU2F "CONN     %03d%03d    6004    .068\n", $n, $k);
            }

          if($n_curr_cols >= 3) {
            printf(NU3F "CONN     %03d%03d    6004    .068\n", $n, $k);
            }

          if($n_curr_cols >= 4) {
            printf(NU4F "CONN     %03d%03d    6004    .068\n", $n, $k);
            }

          if($n_curr_cols >= 5) {
            printf(NU5F "CONN     %03d%03d    6004    .068\n", $n, $k);
            }

          }

        }  # for $k

      }  # for $n

    for ($n = 1; $n <= $n_at; $n++)  {
      
      @at_pos = ($x_coor[$n], $y_coor[$n], $z_coor[$n]);

      for ($i = 0; $i <= 2; $i++) {
        $line = <GCOUT>;
        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) { $vec1[$i] = $scal_fac*$3; }
        if ($n_curr_cols >= 2) { $vec2[$i] = $scal_fac*$4; }
        if ($n_curr_cols >= 3) { $vec3[$i] = $scal_fac*$5; }
        if ($n_curr_cols >= 4) { $vec4[$i] = $scal_fac*$6; }
        if ($n_curr_cols >= 5) { $vec5[$i] = $scal_fac*$7; }
       
        }  # for $i

      if($n_curr_cols >= 1) {
        @dap = &SphereOffs(@at_pos, @vec1, $at_rad[$n]);
        printf NU1F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2]; 
        printf NU1F "    ##  0.020   8001  1000         0\n";
        printf NU1F "CONN     %03d%03d    8011    .034\n", $n, $n_at + $n;
        }

      if($n_curr_cols >= 2) {
        @dap = &SphereOffs(@at_pos, @vec2, $at_rad[$n]);
        printf NU2F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
        printf NU2F "    ##  0.020   8001  1000         0\n";
        printf NU2F "CONN     %03d%03d    8011    .034\n", $n, $n_at + $n;
        }

      if($n_curr_cols >= 3) {
        @dap = &SphereOffs(@at_pos, @vec3, $at_rad[$n]);
        printf NU3F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
        printf NU3F "    ##  0.020   8001  1000         0\n";
        printf NU3F "CONN     %03d%03d    8011    .034\n", $n, $n_at + $n;
        }

      if($n_curr_cols >= 4) {
        @dap = &SphereOffs(@at_pos, @vec4, $at_rad[$n]);
        printf NU4F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
        printf NU4F "    ##  0.020   8001  1000         0\n";
        printf NU4F "CONN     %03d%03d    8011    .034\n", $n, $n_at + $n;
        }

      if($n_curr_cols >= 5) {
        @dap = &SphereOffs(@at_pos, @vec5, $at_rad[$n]);
        printf NU5F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
        printf NU5F "    ##  0.020   8001  1000         0\n";
        printf NU5F "CONN     %03d%03d    8011    .034\n", $n, $n_at + $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
  }  # end of subroutine process_G90


sub process_CADPAC {
  
  # search for "MOLECULAR GEOMETRY (ANGSTROM)"
  while (<GCOUT>) {
    last if /MOLECULAR GEOMETRY \(ANGSTROM\)/;
    }

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

  # collect cartesian coordinates in an array
  $n_at = 0;
  while ($line = <GCOUT>) {
    last if ($line =~ /^$/);  # when last line
    $n_at++;
    $line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;
    $at_weight = int ($5 + 0.5); 
    for ($i = 0; $i < $#at_weights; $i++) {
      last if ($at_weight == $at_weights[$i]);
      } 
    $at_symb[$n_at] = $at_symbols[$i]; 
    $at_rad[$n_at] = $at_radii[$i];
    $at_col[$n_at] = $at_colors[$i];
    $at_num[$n_at] = $1;
    $x_coor[$n_at] = $2;
    $y_coor[$n_at] = $3;
    $z_coor[$n_at] = $4;
    }

  # search for line: 
  # "Transformation matrix from cartesians to normal mode coordinates"
  while ($line = <GCOUT>) {
    last if ($line =~ /Transformation matrix from cartesians to normal mode coordinates/);
    }

  $n_cols = 3*$n_at;  # translations/rotations included in CADPAC
  $n_col_num = 0;
  while($n_cols > 0) {
    if($n_cols >= 8) {
      $n_curr_cols = 8;
      }
    else {
      $n_curr_cols = $n_cols;
      }
    # open files for modes in the schakal 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";
      }
    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";
      }
    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";
      }
    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";
      }
    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";
      }
    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";
      }
    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";
      }
    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";
      }
    # skip 5 lines
    for ($i = 0; $i < 5; $i++ ) {
      $line = <GCOUT>;
      }
    # read frequencies
    $search_expr = ('\s+(\S+)' x $n_curr_cols);
    $line = <GCOUT>;
    # no loop, arrays of filehandles seem to be impossible in perl (ARGHH)
    if($line =~ /$search_expr/) {
      if($n_curr_cols >= 1) { 
        $vib_freq[1] = $1;
        printf(NU1F "TITL NU=%fcm-1\n", $vib_freq[1]);
        }
      if($n_curr_cols >= 2) { 
        $vib_freq[2] = $2; 
        printf(NU2F "TITL NU=%fcm-1\n", $vib_freq[2]);
        }
      if($n_curr_cols >= 3) { 
        $vib_freq[3] = $3; 
        printf(NU3F "TITL NU=%fcm-1\n", $vib_freq[3]);
        }
      if($n_curr_cols >= 4) { 
        $vib_freq[4] = $4; 
        printf(NU4F "TITL NU=%fcm-1\n", $vib_freq[4]);
        }
      if($n_curr_cols >= 5) { 
        $vib_freq[5] = $5; 
        printf(NU5F "TITL NU=%fcm-1\n", $vib_freq[5]);
        }
      if($n_curr_cols >= 6) { 
        $vib_freq[6] = $6; 
        printf(NU6F "TITL NU=%fcm-1\n", $vib_freq[6]);
        }
      if($n_curr_cols >= 7) { 
        $vib_freq[7] = $7; 
        printf(NU7F "TITL NU=%fcm-1\n", $vib_freq[7]);
        }
      if($n_curr_cols >= 8) { 
        $vib_freq[8] = $8; 
        printf(NU8F "TITL NU=%fcm-1\n", $vib_freq[8]);
        }

      # write previously read coordinates and appropriate bonds
      for  ($i = 1; $i <= $n_at; $i++) {
        if($n_curr_cols >= 1) {
          printf(NU1F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 2) {
          printf(NU2F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 3) {
          printf(NU3F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 4) {
          printf(NU4F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 5) {
          printf(NU5F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 6) {
          printf(NU6F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 7) {
          printf(NU7F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }
        if($n_curr_cols >= 8) {
          printf(NU8F "ATOM  %2s  %9.4f %9.4f %9.4f   ## %6.3f %6d %5d         0\n",
                 $at_symb[$i], $x_coor[$i], $y_coor[$i], $z_coor[$i],
                 $at_rad[$i], $at_col[$i], $at_dark);
          }

        for  ($j = $i + 1; $j <= $n_at; $j++) {  # loop for BONDS
          $dx = $x_coor[$i] - $x_coor[$j];
          $dy = $y_coor[$i] - $y_coor[$j];
          $dz = $z_coor[$i] - $z_coor[$j];
          $cut_off = 2.5;
          if (($at_symb[$i] eq "H ") || ($at_symb[$j] eq "H ")) {
            $cut_off = 1.2;
            }
          if ($dx*$dx + $dy* $dy + $dz * $dz < $cut_off) {
            if($n_curr_cols >= 1) {
              printf(NU1F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 2) {
              printf(NU2F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 3) {
              printf(NU3F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 4) {
              printf(NU4F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 5) {
              printf(NU5F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 6) {
              printf(NU6F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 7) {
              printf(NU7F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            if($n_curr_cols >= 8) {
              printf(NU8F "CONN     %03d%03d    6004    .068\n", $i, $j);
              }
            }
          }
        }
      # skip 2 lines
      for ($i = 0; $i < 2; $i++ ) {
        $line = <GCOUT>;
      }
      # now read displacement vectors
      for  ($i = 1; $i <= $n_at; $i++) {

        @at_pos = ($x_coor[$i], $y_coor[$i], $z_coor[$i]);

        for ($l = 0; $l <= 2; $l++) {
          if ($l == 0) {
            $search_expr =  '\S+\s+\S+\s+\S+' . ('\s+(\S+)' x $n_curr_cols);
            }
          else {
            $search_expr =  '\S+' . ('\s+(\S+)' x $n_curr_cols);
            }
          $line = <GCOUT>;
          if( !($line =~ $search_expr)) {
            die "Unstandard format of CADPAC OUTPUT file\n";
            }
          if ($n_curr_cols >= 1) { $vec1[$l] = $scal_fac*$1; }
          if ($n_curr_cols >= 2) { $vec2[$l] = $scal_fac*$2; }
          if ($n_curr_cols >= 3) { $vec3[$l] = $scal_fac*$3; }
          if ($n_curr_cols >= 4) { $vec4[$l] = $scal_fac*$4; }
          if ($n_curr_cols >= 5) { $vec5[$l] = $scal_fac*$5; }
          if ($n_curr_cols >= 6) { $vec6[$l] = $scal_fac*$6; }
          if ($n_curr_cols >= 7) { $vec7[$l] = $scal_fac*$7; }
          if ($n_curr_cols >= 8) { $vec8[$l] = $scal_fac*$8; }

          }  # for $l

        if($n_curr_cols >= 1) { 
          @dap = &SphereOffs(@at_pos, @vec1, $at_rad[$i]);
          printf NU1F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU1F "    ##  0.020   8001  1000         0\n";
          printf NU1F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 2) {
          @dap = &SphereOffs(@at_pos, @vec2, $at_rad[$i]);
          printf NU2F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU2F "    ##  0.020   8001  1000         0\n";
          printf NU2F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 3) {
          @dap = &SphereOffs(@at_pos, @vec3, $at_rad[$i]);
          printf NU3F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU3F "    ##  0.020   8001  1000         0\n";
          printf NU3F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 4) {
          @dap = &SphereOffs(@at_pos, @vec4, $at_rad[$i]);
          printf NU4F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU4F "    ##  0.020   8001  1000         0\n";
          printf NU4F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 5) {
          @dap = &SphereOffs(@at_pos, @vec5, $at_rad[$i]);
          printf NU5F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU5F "    ##  0.020   8001  1000         0\n";
          printf NU5F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 6) {
          @dap = &SphereOffs(@at_pos, @vec6, $at_rad[$i]);
          printf NU6F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU6F "    ##  0.020   8001  1000         0\n";
          printf NU6F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 7) {
          @dap = &SphereOffs(@at_pos, @vec7, $at_rad[$i]);
          printf NU7F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU7F "    ##  0.020   8001  1000         0\n";
          printf NU7F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        if($n_curr_cols >= 8) {
          @dap = &SphereOffs(@at_pos, @vec8, $at_rad[$i]);
          printf NU8F "ATOM  Ar %9.4f %9.4f %9.4f", $dap[0], $dap[1], $dap[2];
          printf NU8F "    ##  0.020   8001  1000         0\n";
          printf NU8F "CONN     %03d%03d    8011    .034\n", $i, $n_at + $i;
          }
        }  # for $i
      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); }
      $n_col_num += $n_curr_cols;
      $n_cols -= $n_curr_cols;
      }
    }
  }