#!/usr/local/bin/perl

# Written by Jan Labanowski (jkl@ccl.net) in 1997 for the common good.

$BONDSCALE = 0.1;    #bond radius is $BONDSCALE * sum of rvdw for atoms
$ATOMSCALE = 0.3;     # atom radius = $ATOMSCALE * rvdw
$BONDMINLENGTH = 0.1;  # bond not showed if shorter than $BONDMINLENGTH times
                        # sum of covalent radii
$BONDMAXLENGTH = 1.5;   # bond not shoed if longer than $bondmaxlenghth times
                        # sum of covalent radii

#symb rwdw  rcov    R      G      B
@atom_data = (  # to mimick the Xmol assignements for Gaussian.
"C",  1.65, 0.77, 0.184, 0.310, 0.310,
"H",  1.25, 0.32, 0.973, 0.973, 1.000,
"Li", 1.86, 1.23, 0.698, 0.133, 0.133,
"B",  0.00, 0.82, 0.000, 0.392, 0.000,
"N",  1.50, 0.75, 0.000, 0.749, 1.000,
"O",  1.35, 0.73, 0.804, 0.361, 0.361,
"F",  1.40, 0.72, 0.933, 0.910, 0.667,
"Na", 1.12, 1.54, 0.941, 0.973, 1.000,
"Mg", 0.87, 1.36, 0.133, 0.545, 0.133,
"Al", 1.70, 1.18, 0.184, 0.310, 0.310,
"Si", 1.93, 1.16, 0.933, 0.910, 0.667,
"P",  1.75, 1.06, 1.000, 0.647, 0.000,
"S",  1.85, 1.02, 0.678, 1.000, 0.184,
"Cl", 1.80, 0.99, 0.000, 0.392, 0.000,
"K",  1.44, 2.03, 1.000, 0.078, 0.576,
"Ca", 1.18, 1.74, 0.184, 0.310, 0.310,
"Sc", 1.37, 1.44, 0.184, 0.310, 0.310,
"Ti", 1.37, 1.32, 0.184, 0.310, 0.310,
"V",  1.37, 1.22, 0.184, 0.310, 0.310,
"Cr", 1.37, 1.18, 0.184, 0.310, 0.310,
"Mn", 1.37, 1.17, 0.184, 0.310, 0.310,
"Fe", 0.90, 1.17, 1.000, 0.647, 0.000,
"Co", 0.88, 1.16, 0.737, 0.561, 0.561,
"Ni", 0.69, 1.15, 0.737, 0.561, 0.561,
"Cu", 0.72, 1.17, 0.737, 0.561, 0.561,
"Zn", 0.74, 1.25, 0.737, 0.561, 0.561,
"Ga", 1.37, 1.26, 0.737, 0.561, 0.561,
"Br", 1.95, 1.14, 0.737, 0.561, 0.561,
"Rb", 1.58, 2.16, 0.737, 0.561, 0.561,
"Ag", 1.27, 1.34, 0.184, 0.310, 0.310,
"I",  2.15, 1.33, 0.627, 0.125, 0.941,
"Cs", 1.84, 2.35, 0.737, 0.561, 0.561,
"Au", 1.37, 1.34, 0.933, 0.910, 0.667,
"Bq", 1.00, 0.00, 1.000, 0.078, 0.576,
"X",  1.00, 0.00, 1.000, 0.078, 0.576,
"Xx", 1.00, 0.00, 1.000, 0.078, 0.576);

for($i = 0; $i <= $#atom_data; $i=$i+6) {  #distribute atom_data to tables
  $j = int($i/6);                          # table row number
  $sym2ord{$atom_data[$i]} = int($i/6);    # needed to get pointer to data
  $symb[$j] = $atom_data[$i];              # atomic symbol
  $rvdw[$j] = $atom_data[$i+1];            # van der Waals radii
  $rcov[$j] = $atom_data[$i+2];            # covalent radii
  $color_r[$j] = $atom_data[$i+3];         # Red color component.
  $color_g[$j] = $atom_data[$i+4];         # Green color component.
  $color_b[$j] = $atom_data[$i+5];         # Blue color component.
  }

if($#ARGV != 0) {
  printf STDERR "Usage: xyz2bs filename.[xyz]\n";
  exit (2);
  }

$outroot = $inpname = $ARGV[0];

$outroot =~ s/\.xyz$//i;

$outname = $outroot . ".bs";    # filename for definitions and 1st frame
$outmove = $outroot . ".mv";    # filename for frames (if exist)

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

$n = -1;
$n_spec = -1;
$n_atoms = 1000;
$n_frames = -1;
while ($line = ) {
  chop($line);                           #new line chopped
  $n++;
  if(($n == 0) || ($n == ($n_frames+1)*($n_atoms + 2))) {
    if($line !~ /\S/) {  # terminate if empty line found where count should be
      $n--;
      last;
      }
    $k = $line + 0;
    if(($n_frames > 0) && ($k != $n_atoms)) { #get  no. of atoms
      fprintf STDERR,
      "Number of atoms in frame %d does not match previous frame\n",
      $n_frames + 1;
      exit(2);
      }
    $n_atoms = $k;
    $n_frames++;
    }
  elsif(($n == 1) || ($n == $n_frames*($n_atoms + 2) + 1)) { #get titles
    $title[$n_frames] = $line;     # get title
    }
  else {
    $line =~ s/^\s+//;                     #front blanks chopped
    $line =~ s/\s+$//;                     #trailing blanks chopped
    if($line =~ /\S/) {                   #if nonblank line
      $k = $n % ($n_atoms + 2) - 1;        #atom number
      @fields = split(/\s+/,$line);        #split into individual columns
      if($#fields < 3) {                   #if not all coordinates given
        printf STDERR
        "Wrong line in XYZ file for atom %d in frame %d\n", $k, $n_frames + 1;
        exit(2);
        }
      $sym = $fields[0];                   #get atom symbol;
      $sym =~ s/[0-9]+$//;                 #skip digits if appended to symbol
      $sym =~ tr/A-Z/a-z/;                 #convert symbol to uppercase
      substr($sym, 0, 1) =~ tr/a-z/A-Z/;   #capitalize the first letter
      for($i = 1; $i <= 3; $i++) {
        if($fields[$i] !~ /^[-+]?[0-9]*\.[0-9]+$/) {
          printf STDERR
          "Line for atom %d in frame %d has wrong number\n",
          $k, $n_frame + 1;
          exit(2);
          }
        }
      # check if atom specs were already saved
      if(defined $sym2ord{$sym}) {
        $ord = $sym2ord{$sym};
        }
      else {

        printf STDERR < 0.001) &&
       ($rcov[$ord2] > 0.001)) {
      $dcov = $rcov[$ord1] + $rcov[$ord2];
      }
    else {
      $dcov = 0.0;
      }

    if($d > 1.5*$dcov) { #skip nonbonded atom pairs
      next;
      }

    # skip bonds already saved
    $saved = 0;
    for($l = 0; $l <= $n_bonds; $l++) {
      if((($b_at1[$l] == $ord1) && ($b_at2[$l] == $ord2)) ||
         (($b_at1[$l] == $ord2) && ($b_at2[$l] == $ord1))) {
        $saved = 1;
        next;
        }
      }

    if($saved == 1) {  # do not save is alrderady there
      next;
      }

    #save unique bond
    $n_bonds++;
    $b_at1[$n_bonds] = $ord1;
    $b_at2[$n_bonds] = $ord2;
    $b_min[$n_bonds] = $BONDMINLENGTH * $dcov;
    $b_max[$n_bonds] = $BONDMAXLENGTH * $dcov;
    $b_rad[$n_bonds] = $BONDSCALE * $dcov;
    $b_col[$n_bonds] = 0.5;
    }
  }
  
# Now print the stuff

open(OUTBS, ">$outname") || die "Cannot create $outname\n";

print OUTBS "* $title[0]\n";

# print coordinates for the 1st frame
for($i = 1; $i <= $n_atoms; $i++) {
  printf OUTBS "atom   %s   %11.6f %11.6f %11.6f\n", $symb[$at_ord[$i+1]], 
        $x[$i+1], $y[$i+1], $z[$i+1];
  }
print OUTBS "\n";

# print specs for atoms
for ($i = 0; $i <= $n_spec; $i++) {
 $ord = $spec_ord[$i];
 printf OUTBS "spec   %s  %6.3f  %5.2f %5.2f %5.2f\n", $symb[$ord],
        $ATOMSCALE*$rvdw[$ord],  $color_r[$ord], $color_g[$ord],
        $color_b[$ord];
 }
print OUTBS "\n";

# print bond specs
for($i = 0; $i <= $n_bonds; $i++) {
  printf OUTBS "bonds  %s   %s  %6.3f  %6.3f  %5.2f  %5.2f\n",
         $symb[$b_at1[$i]], $symb[$b_at2[$i]], $b_min[$i],
         $b_max[$i], $b_rad[$i], $b_col[$i];
  }

close(OUTBS);

# check if additional frames are present in XYZ file, and print them
if($n_frames > 0) {
  open(OUTMV, ">$outmove") || die "Cannot create $outmove\n";
  for($i = 1; $i <= $n_frames; $i++) {
    print OUTMV "frame    $title[$i]\n";
    for($j = $i * ($n_atoms + 2) + 2;  $j <= ($i+1) * ($n_atoms + 2)-1; $j++) {
      printf OUTMV "%9.4f  %9.4f  %9.4f\n", $x[$j], $y[$j], $z[$j];
      }
   }
  close(OUTMV);
  }

# COMMENTS...
# The Xmol is a program originally conceived in the Minnesota Supercomputer
# center and available from anon.ftp ftp.msc.edu (pub/xmol).
# The xbs is the program written by Michael Methfessel, and retrieved from
# http://ihp02.ihp-ffo.de/~msm/.
# the xbs is an X-windows program written in C language which displays
# molecules. The xbs uses its own format for representing molecules.
# For this reason this script is used to convert from the popular XYZ format
# (which is more populer) to the format required by xbs.
# 
# Short example of XYZ format:
#  6
#  This is methanol with normal modes
#  C      0.050065    0.665047    0.000000    0.200000   -0.020000    0.000000
#  O      0.050065   -0.767876    0.000000   -0.140000    0.100000    0.000000
#  H     -0.912101   -1.005927    0.000000    0.340000   -1.740000    0.000000
#  H      1.090132    0.995415    0.000000   -0.040000    0.800000    0.000000
#  H     -0.439468    1.081620    0.886605   -0.180000   -0.080000   -0.180000
#  H     -0.439468    1.081620   -0.886605   -0.180000   -0.080000    0.180000
#
#   1st line -- number of atoms
#   2nd line -- some title (or empty line)
#   3rd line to the end -- atom symbol and XYZ cartesian coordinates in the
#            first columns.  If there are more than 4 columns, the columns
#            additional columns represent:
#                if 5 columns, last column is atomic charges
#                if 7 columns, last column is normal modes vectors
# The XYZ file may contain more then one frame for animations. In this case
# the subsequent frames follow exactly the same format as above.
# 
# 
# The xbs can use one or 2 files which have extensions .bs and .mv.
# The .bs represents the information about the 1st frame and this file
# has to be present. It has the following format
# Short example of xbs format: for a first frame
#  * methanol
# atom   C      0.050065    0.665047    0.000000
# atom   O      0.050065   -0.767876    0.000000
# atom   H     -0.912101   -1.005927    0.000000
# atom   H      1.090132    0.995415    0.000000
# atom   H     -0.439468    1.081620    0.886605
# atom   H     -0.439468    1.081620   -0.886605
# 
# spec   C   0.495   0.70 
# spec   O   0.405   0.40
# spec   H   0.375   1.00
# 
# bonds  C   O   0.150   2.250   0.15   0.70
# bonds  C   H   0.109   1.635   0.11   0.70
# bonds  O   H   0.105   1.575   0.11   0.70
#
# Only lines which start with the words "atom", "spec", and "bonds".
# 
# It is assumed that length is given in Angstroms.
# are used by the program. The other lines are skipped over and not used.
# What the tagged lines mean
#   atom symbol  xcoor ycoor zcoor
#   spec symbol  radius color
#   bonds symbol1 symbol2 min-length max-length radius color.
# The color can be either one number or 3 numbers. One number is the
# degree of gray between 0.0 -- 1.0 (0.0 black, 1.0 white).
# Gray is nice to print to the b/w laser printer. You can also use the
# RGB values, i.e., 3 values between 0.0 and 1.0 which will give the
# saturation with Red, Green, and Blue. E.g. using 1.0 0.0 0.0 will give
# you red color, 1.0 1.0 0.0 bright yellow, 1.0 0.65 0.0 orange,
# and 0.65 0.17 0.17 brown. You can also use a gualified (i.e., official)
# Xwindows color name which is defined in the file /usr/lib/X11/rgb.txt, 
# e.g., spec 0.5 blue
#