/* Written by Jan Labanowski (jkl@ccl.net) in 1997 for the common good.
 * converts files in the Xmol (XYZ) format to the format required by xbs
 * just compile by giving the following command:
 *    cc -o xyz2bs xyz2bs.c -lm
 * and then use as
 *    xyz2bs myfile.xyz
 * which will create myfile.bs and eventually myfile.mv (if more than one frame
 */

#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <math.h>

#define BONDSCALE 0.1   /* bond radius is BONDSCALE * sum of rvdw for atoms */

#define ATOMSCALE 0.3   /* atom radius = ATOMSCALE * rvdw */

#define BONDMINLENGTH 0.1 /* bond not showed if shorter than BONDMINLENGTH
                             times sum of covalent radii */

#define BONDMAXLENGTH 1.5 /* bond not showed if longer than bondmaxlenghth
                             times sum of covalent radii */

typedef struct ATOM_DAT {  /* structure to hold data about atom types */
  char symb[4];           /* chemical symbol */
  float rvdw;             /* van der Waals radius */
  float rcov;             /* covalent radius */
  float r;                /* Red color intensity in RGB values (0 to 1) */
  float g;                /* Green color intensity in RGB values (0 to 1) */
  float b;                /* Blue color intensity in RGB values (0 to 1) */
  } ATOM_DAT;

static int n_at_types = 0;


/*symb rvdw rcov    R      G      B */
static ATOM_DAT atom_data[] = {
"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,
"LST",0.00, 0.00, 0.000, 0.000, 0.000}; /* do not delete this line ever !!! */


/* ====================================================================== */
/* given the chemical symbol, returns the element of atom_data table,
 *  or -1 if match was not found
 */
int find_type(symbol)
char *symbol;
{
 int i, type;

 for(i = 0; i < n_at_types; i++) {
   if(strcmp(symbol, atom_data[i].symb) == 0) {
     return(i);
     }
   }
 return(-1);
}


/* ====================================================================== */
/* skips blanks at the beginning of the string and moves chars to the left,
 *  Returns number of blanks found
 */
int skip_front_blanks(str)
char *str;
{
 int i, k;
 k = 0;
 while((isspace(str[k]) != 0) && (str[k] != '\0')) k++;
 
 if(k != 0) {
   i = 0;
   while(str[i] = str[i+k]) i++;
   }
 return(k);
}  

/* ====================================================================== */
/* skips blanks at the end of the string and returns number of blanks found
 */
int skip_end_blanks(str)
char *str;
{
 int i, k, l;
 l = strlen(str);
 k = l-1;
 for(k = l-1; k >= 0; k--) {
   if(isspace(str[k]) == 0) {
     break;
     }
   }
 str[++k] = '\0';
 return(l-k);
 }

/* ====================================================================== */
/* extract a word (i.e., sequence of nonblank characters) from the string,
 * starting at character pointer by start. start is then automatically
 * advanced to point at the next character after the word. Returns length
 * of extracted word, or zero if nothing extracted
 */
int extract_word(str, k, word)
char *str, *word;
int *k;
{
 int i, j;
 char ch;
 i = *k;
 word[0] = '\0';
 j = 0;
 while((ch = str[i++]) != '\0') {
   if(isspace(ch) == 0) {
     word[j++] = ch;
     word[j] = '\0';
     if(isspace(str[i]) != 0) {
       *k = i;
       return(j);
       }
     }
   }
  *k = i-1;
  return(j);
}

/* ====================================================================== */
main(argc, argv)
int argc;
char **argv;
{
 char inpname[200], outname[200], outmove[200];
 FILE *inpf, *outbs, *outmv;
 int i, j, k, l, n, n_spec, n_frames, n_atoms, n_bonds, ord, ord1, ord2, saved;
 int at_ord[10000];
 int spec_ord[100];
 float d, dx, dy, dz, dcov;
 float b_min[100], b_max[100], b_rad[100], b_col[100];
 float x[10000], y[10000], z[10000];
 int b_at1[10000], b_at2[10000];
 unsigned ip;
 char *pch;
 char line[500];
 char sym[3];
 char *title[1000];
 char fields[20][20];


 /* find how many atom types are available */
 n_at_types = 0;
 while(strcmp(atom_data[n_at_types++].symb, "LST") != 0);

 if(argc != 2) {
   fprintf(stderr, "Usage: xyz2bs filename.[xyz]\n");
   exit (2);
   }

 strcpy(inpname, argv[1]);  /* copy file name to inpfilename */
 skip_front_blanks(inpname);
 skip_end_blanks(inpname);
 strcpy(outname, inpname);  /* copy file name to inpfilename */
 strcpy(outmove, inpname);  /* copy file name to inpfilename */
 

 pch = strchr(inpname, '.');   /* find position of a period */
 if(pch == NULL) {
   pch = strchr(inpname, '\0');
   }

 ip = pch - inpname;  /* point at period */
 strcpy(inpname + ip, ".xyz");  /* input file */
 strcpy(outname + ip, ".bs");   /* output file for 1st frame */
 strcpy(outmove + ip, ".mv");   /* output file for frames 2 to last */

 if(((inpf = fopen(inpname, "r")) == NULL) &&
    ((inpf = fopen(argv[1], "r")) == NULL)) {
   fprintf(stderr,"Could open neither %s nor %s\n", inpname, argv[1]);
   exit(2);
   }

 n = -1;
 n_spec = -1;
 n_atoms = 1000;
 n_frames = -1;
 while (fgets(line, 490, inpf) != NULL) {
   n++;
   skip_front_blanks(line);
   skip_end_blanks(line);
fprintf(stderr,"%d: |%s|\n", n, line);

   if((n == 0) || (n == (n_frames+1)*(n_atoms + 2))) {
     if(line[0] == '\0') {  /* terminate, empty line where atom count */
       n--;
       break;
       }
     sscanf(line, "%d", &k); /* get number of atoms */
fprintf(stderr, "Nat=%d\n", k);
     if((n_frames > 0) && (k != n_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++;
     }
   else if((n == 1) || (n == n_frames*(n_atoms + 2) + 1)) { /* get titles */
     l = strlen(line) + 1;
     title[n_frames] = malloc(sizeof(char)*l);
     strcpy(title[n_frames], line);    /* save title */
fprintf(stderr,"title[%d]=|%s|\n", n_frames, line);
     }
   else {
     k = n % (n_atoms + 2) - 1;        /* atom number */
     l = 0;
     j = 0;
     while(extract_word(line, &j, fields[l]) != 0) {
fprintf(stderr,"Word[%d]=|%s|\n", l, fields[l]);
       l++; /* extract columns */
       }
     if(l < 4) {   /* if line has less columns than symbol and coordinates */
       fprintf(stderr, 
       "Wrong line in XYZ file for atom %d in frame %d\n", k, n_frames + 1);
       exit(2);
       }
     sym[0] = fields[0][0];             /* get atom symbol */;
     sym[1] = fields[0][1];
     sym[2] = '\0';
     if(isalpha(sym[0]) == 0) {
       fprintf(stderr, "Wrong atom symbol at line %d\n", n);
       exit(2);
       }
     if(isalpha(sym[1]) == 0) {
       sym[1] = '\0';
       }
     if(islower(sym[0]) != 0) {
       sym[0] = toupper(sym[0]);
       }
     if(isupper(sym[1]) != 0) {
       sym[1] = tolower(sym[1]);
       }
     
     x[n] = atof(fields[1]);
     y[n] = atof(fields[2]);
     z[n] = atof(fields[3]);
     if((ord = find_type(sym)) < 0) {
       fprintf (stderr,
"Atom symbol %s is not defined in this perl script. Add this symbol to\n",
       sym);
       fprintf(stderr,
       "the table at the top of the program, recompile,  and rerun.\n");
       exit(2);
       }
     at_ord[n] = ord;
     if(n_frames == 0) {
       saved = 0;
       for(i = 0; i <= n_spec; i++) {
         if(spec_ord[i] == ord) {
           saved = 1;
           break;
           }
         }
       if(saved == 0) {
         n_spec++;
         spec_ord[n_spec] = ord;
         }
       }
     }
   }

 if(n != (n_frames+1)*(n_atoms + 2) - 1) {
   fprintf(stderr, 
   "Number of atoms and number of coordinates do not much in last frame\n");
   exit(2);
   }
    
/* find bonds which are present in the 1st frame. */

 k = 0;
 n_bonds = -1;
 for (i = 2; i < n_atoms + 1; i++) {      /* scan all atom pairs */
   for(j = i+1; j < n_atoms + 2; j++) {
 
     /*  calculate distance between atoms; */
     dx = x[i] - x[j];
     dy = y[i] - y[j];
     dz = z[i] - z[j];
     d = sqrt(dx*dx + dy*dy + dz*dz);

     /* get atom pointers in the table */
     ord1 = at_ord[i];
     ord2 = at_ord[j];

     /* calculate sum of covalent radii for the atoms */
     if((atom_data[ord1].rcov > 0.001) &&
        (atom_data[ord2].rcov > 0.001)) {
        dcov = atom_data[ord1].rcov + atom_data[ord2].rcov;
        }
     else {
        dcov = 0.0;
        }

     if(d > 1.5*dcov) { /* skip nonbonded atom pairs */
       continue;
       }

     /* 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;
         continue;
         }
       }

     if(saved == 1) {  /* do not save is alrderady there */
       continue;
       }

     /* 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 */

 if((outbs = fopen(outname, "w")) == NULL) {
   fprintf(stderr, "Cannot create %s\n", outname);
   exit(2);
   }

 fprintf(outbs, "* %s\n",title[0]);

 /* print coordinates for the 1st frame */
 for(i = 1; i <= n_atoms; i++) {
   fprintf(outbs, "atom   %s   %11.6f %11.6f %11.6f\n", atom_data[at_ord[i+1]].symb,
        x[i+1], y[i+1], z[i+1]);
   }
 fprintf(outbs, "\n");

 /*  print specs for atoms */
 for (i = 0; i <= n_spec; i++) {
   ord = spec_ord[i];
   fprintf(outbs, "spec   %s  %6.3f  %5.2f %5.2f %5.2f\n", atom_data[ord].symb,
        ATOMSCALE*atom_data[ord].rvdw,  atom_data[ord].r, atom_data[ord].g,
        atom_data[ord].b);
   }
 fprintf(outbs, "\n");

 /* print bond specs */
 for(i = 0; i <= n_bonds; i++) {
   fprintf(outbs, "bonds  %s   %s  %6.3f  %6.3f  %5.2f  %5.2f\n",
         atom_data[b_at1[i]].symb, atom_data[b_at2[i]].symb, b_min[i],
         b_max[i], b_rad[i], b_col[i]);
    }

 fclose(outbs);

 /* check if additional frames are present in XYZ file, and print them */
 if(n_frames > 0) {
   if((outmv = fopen(outmove, "w")) == NULL) {
     fprintf(stderr, "Cannot create %s\n", outmove);
     exit(2);
     }
   for(i = 1; i <= n_frames; i++) {
     fprintf(outmv, "frame    %s\n", title[i]);
     for(j = i * (n_atoms + 2) + 2;  j <= (i+1) * (n_atoms + 2)-1; j++) {
       fprintf(outmv, "%9.4f  %9.4f  %9.4f\n", x[j], y[j], z[j]);
       }
     }
   fclose(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
#
*/