Crowded

#! /usr/bin/perl 
# Crowded program can check for the geometry of input file, and determine how crowdy
# the geometry is
# It reads in both zmatrix and xyz input
# USAGE: Crowded $mol.com


use Term::ANSIColor;
my $maxBadness = 2;
my $tol = 0.3;

#atomic radii(empirical) (picometres)
#http://en.wikipedia.org/wiki/Atomic_radii_of_the_elements_%28data_page%29
my @symbo=("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");

my @radii= (25,  0,145,105, 85, 70, 65, 60, 50,  0,
        180,150,125,110,100,100,100, 71,220,180,
        160,140,135,140,140,140,135,135,135,135,
        130,125,115,115,115,  0,235,200,180,155,
        145,145,135,130,135,140,160,155,155,145,
        145,140,140,  0,260,215,195,185,185,185,
        185,185,185,180,175,175,175,175,175,175,
        175,155,145,135,135,130,135,135,135,150,
        190,180,160,190,  0,  0,  0,215,195,180,
        180,175,175,175,175,  0,  0,  0,  0,  0,
        0,  0,  0); 

my @coord= ( 1,  1,  1,  2,  3,  4,  3,  2,  1,  1,
        1,  2,  3,  4,  5,  6,  1,  1,  1,  2,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8,  8,  8,  8,  8,  8,  8,  8,
        8,  8,  8); 

if (isBadComFile($ARGV[0])) {
    print color 'blue';
    print "$ARGV[0] file too crowded\n";
    print color 'reset';
}

##############################################################

# Usage isBadComFile(nameOfMolecule)
sub isBadComFile {
    #my $molecule = shift;

    my @atom;
    my @X, @Y, @Z;

    my $badnessCount = 0;

    my $molecule = $ARGV[0];
    chomp($ARGV[0]);
    $cominput=file2string($ARGV[0]);
    stripcom($cominput);
    
    if ($iszmat){
        zmtCom2xyz($molecule, \@atom, \@X, \@Y, \@Z);
    }
    else{
        #print $cominput;
        @lines= split(/\n/,$cominput);
        for ($i =0 ; $i <= $#cominput ; $i++){
            ($Z[$i],$Y[$i],$X[$i],$atom[$i]) = reverse(split(/\s+/,$lines[$i]));
        }
    }

#    zmtCom2xyz($molecule, \@atom, \@X, \@Y, \@Z);

    for ($i=0;$i <= $#atom; $i++){
        for ($k=0; $k <= $#symbo; $k++){
            if (lc($atom[$i]) eq lc($symbo[$k])){
                $allow[$i] = $radii[$k]*0.01;
                $coordN[$i] = $coord[$k];
            }
        }
        for ($j=0;$j < $i; $j++){
            $dist[$i][$j] = sqrt(($X[$i] - $X[$j])**2 + ($Y[$i] - $Y[$j])**2 +($Z[$i] - $Z[$j])**2) ;
            $sbond[$i][$j] = $allow[$i] + $allow[$j] -$tol;
            $lbond[$i][$j] = $allow[$i] + $allow[$j] +$tol;

            # test if one atom is too close to another
            if ($sbond[$i][$j] > $dist[$i][$j]){
            printf "atom %s and %s has a bond length %f < allowance %f\n", $i+1, $j+1, $dist[$i][$j], $sbond[$i][$j] ;
                $badnessCount++;
            }
            # test if an atom has too many neighbours
            if ( $dist[$i][$j] < $lbond[$i][$j]){
                $bonding[$i] ++;      
                $bonding[$j] ++;      
                #printf "bonding site: %s, %s\n",$i+1,$j+1;
            }
        }
    }
    for ($i=0;$i <= $#atom; $i++){
        #print "$i+1, $bonding[$i]\n";
        if ($bonding[$i] > $coordN[$i]) {
            printf "atom %s has %s bonding > allowance %s\n",$i+1,$bonding[$i], $coordN[$i];
            $badnessCount++;
        }
    }

    print "$molecule file badness rating = $badnessCount \n";

    if ($badnessCount > $maxBadness) {
        #fakeGaussian($molecule);
        return 1;
    }
    return 0;
}

# usage: zmtCom2xyz($molecule, \@atom, \@X, \@Y, \@Z);
sub zmtCom2xyz {
    my $mol = shift;
    $atom = shift;
    my ($x,$y,$z) = @_;
    
    open COMFILE, "< $mol" or die "Couldn't open file $mol\n$!";
    
    # this takes Gaussian input file with geometry given as
    # Zmat, and converts it to the XYZ format (cartesians)
    # This is an educational script, hence many comments.
    # The script is a filter, i.e., should be used as
    #   cat name.inp | g9xinp2xyz.pl > name.xyz
    # where name.inp is the Gaussian input file and
    # name.xyz is a file with cartesian coordinates a la XMol
    # It is assumed that G9x input file has a format:
    #   % part (e.g., checkpoint specs)
    #   # keywords/links part
    #     empty line
    #   title lines
    #     empty line
    #   charge multiplicity line
    #   zmatrix lines
    #     empty line or "Variables"
    #   symbolic variables for Z-matrix (if ZMAT has symbolic parameters)
    #     empty line or "Constants"
    #   symbolic constants for Z-matriz (if all references were not resolved yet)
    #     empty line or line starting with "-"
    #
    #  Written by Jan Labanowski in July 1997 <jkl@ccl.net>

    #

    $keywords = "";
    while($line = <COMFILE>) { # read next lines of keywords (route section)
      last if ($line =~ /\s*\d\s+\d\s*/);
    }


    $natoms = 0;    #number of atoms found initialized
    while ($line = <COMFILE>) {  # read molecule description lines
      chop($line);
      $line =~ tr/a-z/A-Z/;
      $line =~ s/^\s+//;   # get rid of front and back spaces
      $line =~ s/$\s+//;
      if(($line !~ /\S/) ||
         ($line =~ /variab/i) ||
         ($line =~ /consta/i)) {  # if blank line or end of Zmat
        last;              # exit the loop
        }
      if($natoms == 0) {   # first atom is only an atomic symbol
        $symbol[$natoms] = $line;
        # no parameters necessary for atom 1
        ($bf[$natoms], $af[$natoms], $tf[$natoms]) =  (1, 1, 1);
        }
      else {
        # split the line into individual fields:
        (@fields) =  split(/\s+/, $line);   # split line into fields

        # assumed most standard z-matrix form, so check number of entries
        if((($#fields != 2) && ($natoms == 1)) &&
           (($#fields != 4) && ($natoms == 2)) &&
           (($#fields != 6) && ($natoms >= 3))) {
          printf STDERR "Error in Z-Matrix for atom %d\n", $natoms+1;
          exit(2);   # abort 
          }

        for($i = $#fields+1; $i <8; $i++) {  # fill up with empty strings
          $fields[$i] = "";
          }

        # parse fields  symbol_i  j  d_ij   k   theta_ijk  l  phi_ijkl 
        ($symbol[$natoms],  # atom symbol
         $bondat[$natoms],  # bonded atom
         $blen[$natoms],    # bond length
         $angat[$natoms],   # atom in an angle
         $angle[$natoms],   # angle
         $torat[$natoms],   # last atom in torsion angle
         $tors[$natoms])    # torsion angle
                              = (@fields);
        # initialize flags for required parameters to "NotFound" status (i.e., 0)
        if($natoms == 1) {
          ($bf[$natoms], $af[$natoms], $tf[$natoms]) =  (0, 1, 1);
          }
        elsif($natoms == 2) {
          ($bf[$natoms], $af[$natoms], $tf[$natoms]) =  (0, 0, 1);
          }
        else {
          ($bf[$natoms], $af[$natoms], $tf[$natoms]) =  (0, 0, 0);
          }
        
        # check if parameters are given as numbers (i.e., sign, possibly followed
        # by digits, possibly followed by period, possibly followed by digits
        # and containing at least one digit total. Set flag if number.
        if(($blen[$natoms] =~ /^[-+]?[0-9]*\.?[0-9]*$/) &&
           ($blen[$natoms] =~ /[0-9]+/)) {
          $bf[$natoms] = 1;
          }
        if(($angle[$natoms] =~ /^[-+]?[0-9]*\.?[0-9]*$/) &&
           ($angle[$natoms] =~ /[0-9]+/)) {
          $af[$natoms] = 1;
          }
        if(($tors[$natoms] =~ /^[-+]?[0-9]*\.?[0-9]*$/) &&
           ($tors[$natoms] =~ /[0-9]+/)) {
          $tf[$natoms] = 1;
          }

        # Atom numbers must be all digits
        if(($bondat[$natoms] !~ /^[0-9]+$/) ||
           (($angat[$natoms] !~ /^[0-9]+$/) && ($natoms > 1)) ||
           (($torat[$natoms] !~ /^[0-9]+$/) && ($natoms > 2))) {
          printf STDERR "Error in atom definition for %d\n", $i+1;
          exit(2); 
          }
        }
      $natoms++;
      }

    # check any symbolic parameters are given in the Z-mat
    $parms_given = 0;
    for($i = 0; $i < $natoms; $i++) {
      if(($bf[$i] == 0) || 
         ($af[$i] == 0) || 
         ($tf[$i] == 0)) {
        $parms_given = 1;     # symbolic parameters present
        last;
        }
      }

    if($parms_given == 1) { # read in symbolic parameters (variables) when present
      while($line = <COMFILE>) {
        $line =~ s/^\s//;        # trailing blanks out!
        $line =~ s/\s$//;
        $line =~ tr/a-z/A-Z/;    # convert to uppercase
        if(($line !~ /\S/) ||
           ($line =~ /const/i)) {  # if blank line or end of Variables
          last;              # exit the loop
          }

        ($name, $value) = @fields = split(/[ =,]+/, $line);  # retrieve fields
        if($#fields != 1) {
          die "Wrong format on Variables assignment line\n";
          }
        if(($value !~ /^[-+]?[0-9]*\.?[0-9]*$/) ||
           ($value !~ /[0-9]+/)) {
          die "Wrong number format of Variables assignment line\n";
          }


        # substitute parameter in Z-Matrix. The symbolic parameter name is
        # allowed to be preceded with a sign, so take this into account
        $parameter_found = 0;
        for($i = 0; $i < $natoms; $i++) {
          if($bf[$i] == 0) {
            if($blen[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $blen[$i] = $sign * $value;
              $parameter_found = 1;
              $bf[$i] = 1;
              }
            }

          if($af[$i] == 0) {
            if($angle[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $angle[$i] = $sign * $value;
              $parameter_found = 1;
              $af[$i] = 1;
              }
            }

          if($tf[$i] == 0) {
            if($tors[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $tors[$i] = $sign * $value;
              $parameter_found = 1;
              $tf[$i] = 1;
              }
            }
          }
        # parameter, if given, must appear at least once in Z-Matrix
        if($parameter_found == 0) {
          die "Parameter $name = $value does not appear in Z-Matrix\n";
          }

        }
      }

    # check if all symbolic paramemters were reslolved
    $parms_given = 0;
    for($i = 0; $i < $natoms; $i++) {
      if(($bf[$i] == 0) || 
         ($af[$i] == 0) || 
         ($tf[$i] == 0)) {
        $parms_given = 1;    # symbolic parameters present
        last;
        }
      }
    # if all parameters are not resloved, read in the Constants section
    if($parms_given == 1) { # read in symbolic parameters (variables) when present
      while($line = <COMFILE>) {
        $line =~ s/^\s//;        # trailing blanks out!
        $line =~ s/\s$//;
        $line =~ tr/a-z/A-Z/;    # convert to uppercase
        if($line =~ /const/i) {
          next;
          }
        if(($line !~ /\S/) ||
           ($line =~ /^\s*\-/)) {  # if blank line or - at the front
          last;              # exit the loop
          }
        ($name, $value) = @fields = split(/[ =,]+/, $line);  # retrieve fields
        if($#fields != 1) {
          die "Wrong format of Constants assignment line\n";
          }
        if(($value !~ /^[-+]?[0-9]*\.?[0-9]*$/) ||
           ($value !~ /[0-9]+/)) {
          die "Wrong number format on Constants assignment line\n";
          }

        # substitute parameter in Z-Matrix. The symbolic parameter name is
        # allowed to be preceded with a sign, so take this into account
        $parameter_found = 0;
        for($i = 0; $i < $natoms; $i++) {
          if($bf[$i] == 0) {
            if($blen[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $blen[$i] = $sign * $value;
              $parameter_found = 1;
              $bf[$i] = 1;
              }
            }

          if($af[$i] == 0) {
            if($angle[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $angle[$i] = $sign * $value;
              $parameter_found = 1;
              $af[$i] = 1;
              }
            }

          if($tf[$i] == 0) {
            if($tors[$i] =~ /^([-+]?)$name$/) {
              $sign = $1 . "1.0";
              $tors[$i] = $sign * $value;
              $parameter_found = 1;
              $tf[$i] = 1;
              }
            }
          }

        # parameter, if given, must appear at least once in Z-Matrix
        if($parameter_found == 0) {
          die "Parameter $name = $value does not appear in Z-Matrix\n";
          }
        }
      }

    # check if all symbolic paramemters were reslolved
    $parms_given = 0;
    for($i = 0; $i < $natoms; $i++) {
      if(($bf[$i] == 0) || 
         ($af[$i] == 0) || 
         ($tf[$i] == 0)) {
        $parms_given = 1;    # symbolic parameters present
        last;
        }
      }
    if($parms_given == 1) {
      print STDERR
      "The follwing symbolic parameters are not resolved in Z-matrix:\n";
      for($i = 0; $i < $natoms; $i++) {
        if($bf[$i] == 0) {
          printf STDERR "   %s for bond length for atom %d\n", $blen[$i], $i;
          }
        if($af[$i] == 0) {
          printf STDERR "   %s for angle value for atom %d\n", $angle[$i], $i;
          }
        if($tf[$i] == 0) {
          printf STDERR "   %s for torsion for atom %d\n", $tors[$i], $i;
          }
        }
      exit(2);
      }


    # output cartesian coordinates in XYZ XMol format
    for ($i = 0; $i < $natoms; $i++) {
        oneatom($i+1, $bondat[$i], $blen[$i], $angat[$i], $angle[$i],
              $torat[$i], $tors[$i]);

        $$atom[$i] = $symbol[$i];
        $$z[$i] = $Zcoor[$i+1];
        $$x[$i] = $Xcoor[$i+1];
        $$y[$i] = $Ycoor[$i+1];

#        printf STDOUT "%-3s  %12.6f  %12.6f  %12.6f\n",$$atom[$i],$$z[$i],$$x[$i],$$y[$i];

    }

    close COMFILE;

}

# -- Takes Z-mat line, i.e., 
#        atomnumber, ibond, bondlen, iangle, bondangle, itors, torsangle
#    and updates NONLOCAL arrays Xcoor, Ycoor, Zcoor which hold cartesian
#    coordinates. It is assumed that atom numbering starts from 1
#    (not from 0 as usually in C and perl)
#    USES nonlocal array Xcoor, Ycoor, Zcoor where it stores coordinates
#    It is assumed that angles are given in DEGREEs and bondlen in the same
#    units in which cartesian coordinates will be given.
#                  
#             ibond -- iangle
#              /             \ 
#          atnum             itors
#
#  CAUTIONS -- contrary to the UNIX/PERL/C tradition, the
#  arrays XYZ start from subscript 1 (not 0), and XYZ[0] are left unused,
#  in other words, array element numbers correspond to atom numbers.

sub oneatom {
  # variables are declared local, so they do not interfere with
  # variables outside this subroutine
  local ($atnum,
         $ibond, $bondlen,
         $iangle, $bondangle,
         $itors, $torsangle);
  local ($deg2r);
  local ($rinv, $aux);
  local ($xa, $ya, $za);
  local ($xb, $yb, $zb);
  local ($xc, $yc, $zc);
  local ($xx, $yy, $zz);

  $atnum = $_[0];

  if($atnum == 2) {
    ($atnum, $ibond, $bondlen) = @_;
    }
  elsif ($atnum == 3) {
    ($atnum, $ibond, $bondlen, $iangle, $bondangle) = @_;
    }
  elsif ($atnum > 3) {
    ($atnum, $ibond, $bondlen, $iangle, $bondangle, $itors, $torsangle) = @_;
    }

  $deg2r = atan2(1.0, 1.0)*4.0/180.0;

  # first atom in the origin of coordinate system      
  if($atnum == 1) {
    $Xcoor[1] = 0.0;
    $Ycoor[1] = 0.0;
    $Zcoor[1] = 0.0;
    return;
    }

  # second atom on Z-axis, bondlen from the first atom
  if($atnum == 2) {
    $Xcoor[2] = 0.0;
    $Ycoor[2] = 0.0;
    $Zcoor[2] = $bondlen;
    return;
    }

  # bonded atom coordinates
  $x0 = $Xcoor[$ibond];
  $y0 = $Ycoor[$ibond];
  $z0 = $Zcoor[$ibond];

  if($atnum == 3) {    # third atom placed on the Oxz plane
    $Xcoor[3] = $x0+$bondlen*sin($bondangle*$deg2r);
    $Ycoor[3] = 0.0;
    $aux = $bondlen*cos($bondangle*$deg2r);
    if($ibond == 1) {
      $Zcoor[3] = $z0+$aux;
      }
    elsif ($ibond == 2) {
      $Zcoor[3] = $z0-$aux;
      }
    return;
    }


  # vector from j-->k
  $xx = $Xcoor[$iangle] - $Xcoor[$ibond];
  $yy = $Ycoor[$iangle] - $Ycoor[$ibond];
  $zz = $Zcoor[$iangle] - $Zcoor[$ibond];

  # inverted vector length
  $rinv = 1.0/sqrt($xx*$xx + $yy*$yy + $zz*$zz);

  # unit vector ibond--->iangle
  $xa = $xx*$rinv;
  $ya = $yy*$rinv;
  $za = $zz*$rinv;

  # vector from ibond--->itors  
  $xb = $Xcoor[$itors] - $Xcoor[$ibond];
  $yb = $Ycoor[$itors] - $Ycoor[$ibond];
  $zb = $Zcoor[$itors] - $Zcoor[$ibond];

# unit vector ibond--->itors
  $xc = $ya*$zb - $za*$yb;
  $yc = $za*$xb - $xa*$zb;
  $zc = $xa*$yb - $ya*$xb;

  $rinv = 1.0/sqrt($xc*$xc + $yc*$yc + $zc*$zc);

  $xc = $xc*$rinv;
  $yc = $yc*$rinv;
  $zc = $zc*$rinv;

  $xb = $yc*$za - $zc*$ya;
  $yb = $zc*$xa - $xc*$za;
  $zb = $xc*$ya - $yc*$xa;

  $zz = $bondlen*cos($bondangle*$deg2r);
  $xx = $bondlen*sin($bondangle*$deg2r)*cos($torsangle*$deg2r);
  $yy = $bondlen*sin($bondangle*$deg2r)*sin($torsangle*$deg2r);

  $Xcoor[$atnum] = $x0 + $xa*$zz + $xb*$xx + $xc*$yy;
  $Ycoor[$atnum] = $y0 + $ya*$zz + $yb*$xx + $yc*$yy;
  $Zcoor[$atnum] = $z0 + $za*$zz + $zb*$xx + $zc*$yy;

  return;
  }


################################################################################
sub file2string {
    # file2string("filename") : gets a txt file, turns it into a string.
    my ($i,$file,$temp);
    local $/;

    open(INPUTFILE, $_[0]) ||die "Can't open $_[0] for reading.\n";
    undef $/;
    $file = <INPUTFILE>;
    close(INPUTFILE)|| die "Can't close $_[0].\n";
    return $file;
}

sub stripcom{
    #stripcom("geometry_string"), returns stripped geometry string.
    my($geometry, @file, $cm_flag, $bk, $cm_index, $i);
    my($temp2, @temp);

    $geometry = $_[0];
    @file = split(/\n/,$geometry);
    $cm_flag = 0;   # whether we've hit the (first) charge/multiplicity section or not.
    $bk = 0;        #blank line counter for zmat section
    @cominput = ();
    $cm_index = 9999999;

    for($i = 0; $i <= $#file; $i++){
        #match the charge/multiplicity section.
        if ($file[$i] =~/#/){
            $theory =${theory}.${file[$i]};
        }

        if ($file[$i] =~/^[ \t]*-*\+*[ \t]*[0-9]{1,2}[ \t]*\+*[ \t]*[0-9]{1,2}[ \t]*$/){
            $cm_flag +=1;         #means you only get the top geometry.
            if ($cm_flag == 1){
                #   push(@cominput, $file[$i]);       #put cm into array
                $cm_index = $i;                #mark where cm is
                $temp2 = $file[$i+1];              #get next line
                $temp2 =~s/^[ \t]*//;             #strip leading whitespace
                $temp2 =~s/[ \t]*$//;              #strip trailing whitespace
                @temp = split(/[ \t]+/,$temp2);        #split
                if ($temp[0] =~ /\b[a-z]{1,2}\b/i){        #see if zmat
                    if (@temp == 1){
                        $iszmat = 1;
                    }else{
                        if (@temp == 4){
                            $iszmat = 0;
                        }else{
                            die "$temp2 : Not xyz or xmat!\n";
                        }
                    }
                }else{
                    die "$temp[0] : not an atom!!\n";
                }
            }
        }
        if ($i > $cm_index && $cm_flag == 1){    #discard useless header and footer
            if($iszmat == 0 && $file[$i] =~/^[ \t]*[a-zA-Z]{1,2}[ \t]{1,}-?\+? *\.?[0-9]{1,}\.?[0-9]*[ \t]{1,}-?\+? *\.?[0-9]{1,}\.?[0-9]*[ \t]{1,}-?\+? *\.?[0-9]{1,}\.?[0-9]*[ \t]*$/){     #if zmat and a coordinate
                push(@cominput, $file[$i]);
            }
            if($iszmat == 1){
                if($file[$i] =~ /^[ \t]*$/){
                    $bk+=1;             #count blank lines in zmat
                }
                if ($bk < 2){      # if zmat and before the end of it, then
                    push(@cominput, $file[$i]);
                }
            }
        }
    }
    $cominput = join("\n",@cominput);
    return $cominput;
}