/* this program produces XYZ animation file for molecular vibrations 
   represented by normal modes taken from output files from
   Aces2, Gamess, and Gaussian. This program autosenses what output file is
   given, so you do not have to worry to give some parameters. However
   if file is none of the listed, the results are unpredictable (i.e, 
   it will probably die with some last gasp messsage. You are welcome to
   add support for other packages... The program is not copyrighted as you
   see, so you are welcome to expand it and add your name to the author's list.
 */

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

/*
 <math.h> not needed here because of explicit declaration of
sin() below
*/

/* Some function delclarations, so compilers do not print these
   stupid warnings
 */
void wrcors();
void wrat();
void readln();

#define MAXAT 100

#define TOANG 0.529177249
#define TODEG 57.295779513082321
#define TWOPI 6.28318530717958647688

/* define codes for various quantum packages */
#define G94            1   /* when G94 file */
#define G88            0   /* when G88 file */
#define GAMESST        2   /* when GAMESS file */
#define ACES2T         3   /* ACES2 file */
#define UNKNOWN       -1   /* When could not decide what file is it */
#define ALLVIBS        1    /* When Program reports all 3*N modes */
#define NORTVIBS       0    /* When Program reports only 3*N-6 modes */

/* global variables */
int vibtotal = ALLVIBS;     /* initialize to all vibrations */
int debug = 0;              /* set to 1 if you want a debug output */

/**********************************************************

Program extracts normal modes from Gaussian, GAMESSS, and ACES2 output file.
It writes separate files, one simple cosine trajectory for
particular normal mode. Output is in Xmol XYZ format.
Click Animate button and enjoy!
Vectors representing the normal coordinate are attached to the
atoms, which is particularly useful for a static picture.
They can be toggled off/on in xmol.

This is version 3.1 (July 1997)

Compiling:

cc -o xvibs xvibs.c -lm

Running:

xvibs mymolecule.out {all|1-(3*N-6)}

Example:

g90 <benzene.inp > benzene.out (or equivalent)

xvibs benzene.out all

Output files are named benzene.{1-30}

or

xvibs benzene.out 15

for mode 15 only.

Otput will be benzene.15


-------------------------------------------
Written in 1990 by Milan Hodoscek
Institute of Chemistry, Ljubljana, Slovenia
Modified for XMOL 1992 @NIH, Bethesda, MD
Modified again Fred Brouwer Univ. of Amsterdam 1994
to include vectors (fred@org.chem.uva.nl (Fred Brouwer))
Modified to accept G94, Yos Ginting (Yos.Ginting@chem.utas.edu.au)
Modified by Jan Labanowski, jkl@ccl.net to accept input from GAMESS and ACESS

Bug report to: milan@kihp6.ki.si
-------------------------------------------

History:

Version 2 (June 1994)

-- Fixed bug to accept file names without point
-- Maximum number of atoms is now define with MAXAT
-- Added data for drawing vectors representing normal modes
-- Prints energy for G90 and G92 output

Version 3 (July 1996)
Yos Ginting <Yos.Ginting@chem.utas.edu.au>

-- Add capability to parse G94 output.
   The folowing two function were added:
     is_g94 , return 1 if input is in G94 format
              return 0 if it is in G92 format
              otherwise, display error messages
     vib_g94  is a modification of vib_g88.
              It now properly parses g94 output

Version 3.1 (July 1997)
Jan K. Labanowski <jkl@ccl.net>

-- replaced is_94 with which_program, and g94 with file_type
   The which_program returns the code (look up the #define statements) 
   which represents the type of format of OUTPUT file.
   Added global variables:
      vibtotal -- provides info on how many vibrations are listed
                  in output file (ALLVIBS -- 3N, i.e., rotations and 
                  translations included, NORTVIBS -- 3N-6, rotations
                  and translations removed
      debug -- for program debugging purposes, if set to 1, prints
               additional information to stderr
   Added routines for processing GAMESS and ACES2 output files


************************************************************/


main(argc, argv)
     int argc ;
     char *argv[] ;
{

    double xs[MAXAT], ys[MAXAT], zs[MAXAT] ;
    double x[MAXAT], y[MAXAT], z[MAXAT] ;
    int atnum[MAXAT];
    char vibnumb[10];

    double e_g88(), vib_g88(), vib_g94(), vib_gamess(), e_gamess(),
      vib_aces2(), e_aces2(),  e ;
    int atoi(), vibn, indx, all, file_type, which_program() ;
    char *pindx ;
    double sin() ;

    FILE *out, *dt1;
    int i, j, k, l, fact ;
    char ime[MAXAT] ;
    double scale = 0.7 ;
    double xx, yy, zz, ffact ;

    if (argc<3)
    {
        fprintf(stderr,"Too few parameters\n"); exit(1);
    }
    
    dt1 = fopen(argv[1],"r");
    if( dt1 == NULL )
    {
        fprintf(stderr,"\t<--- file error on %s\n",argv[1]);
        exit(2) ;
    }

    /* Get information about the input file format. */
    file_type=which_program(dt1);
    if(debug == 1) fprintf(stderr, "File type is = %d\n", file_type);
    /* Back to the start of the file */
    rewind(dt1);
    all = 0 ;
    if (argv[2][0] == 'a' )
    {
     all = 1 ; vibn=0 ;  /* all vibrations will be processed */
    }
    else vibn = atoi(argv[2]) - 1 ;  /* one vibration given on command line
                                      * subtract 1 before advancing it */
    
    if ( ( pindx = strchr(argv[1],'.') ) == NULL )
        indx = strlen(argv[1]) ;   /* length of the file name root */
    else
        indx = strlen(argv[1]) - strlen(pindx) ;

    ime[0] = 0 ;  
    strncpy(ime,argv[1],indx) ;  /* copy the file name root to ime */

    allvib:     /* LOOP OVER VIBRATIONS */
    vibn++ ;    /* next vibration */

    sprintf(vibnumb,"%s.%03d",ime,vibn);   /* create file name */
    if(debug == 1) fprintf(stderr, "Creating file %s\n", vibnumb);
    out = fopen(vibnumb,"w");            /* open file name */
    if( out == NULL )
    {
        fprintf(stderr, "\t<--- file error on %s",vibnumb);
        exit(3) ;
    }
    rewind(dt1) ;                        /* revind file name */
    if((file_type == G94) || (file_type == G88)) {
      l = coor_g88(dt1,atnum,x,y,z);       /* read in coordinates */
      }
    else if(file_type == GAMESST) {
      l = coor_gamess(dt1, atnum, x, y, z);
      }
    else if(file_type == ACES2T) {
      l = coor_aces2(dt1, atnum, x, y, z);
      }
    else {
      fprintf(stderr, "Unknown input file type at coor\n");
      exit(2);
      }
    
    if(debug == 1) fprintf(stderr, "Natoms=%d\n",l); 
    rewind(dt1);                         /* revind file */
    if((file_type == G94) || (file_type == G88)) {
      e = e_g88(dt1);       /* read energy */
      }
    else if(file_type == GAMESST) {
      e = e_gamess(dt1);
      }
    else if(file_type == ACES2T) {
      e = e_aces2(dt1);
      }
    else {
      fprintf(stderr, "Unknown input file type\n");
      exit(2);
      }
    
    if(debug == 1) fprintf(stderr,"Energy=%f\n",e);
    /*
        Write the header for the equilibrium structure:
         (this is done now because e gets a new value with the call
        to vib_g88)
    */
    fprintf(out,"    %d\n",l);          /* l is number of real atoms */
    fprintf(out,"* %s %15.8f\n",argv[1],e);  /* print title */

    /*
        Get normal coordinates:
    */
    rewind(dt1) ;
    if (file_type == G94)
    {
        e=vib_g94(dt1, vibn, l, xs, ys, zs);
    }
    else if(file_type == G88)
    {
        e=vib_g88(dt1, vibn, l, xs, ys, zs);
    }
    else if(file_type == GAMESST)
    {
        e=vib_gamess(dt1, vibn, l, xs, ys, zs);
    }
    else if(file_type == ACES2T)
    {
        e=vib_aces2(dt1, vibn, l, xs, ys, zs);
    }
    if(debug == 1) fprintf(stderr, "Frequency = %f\n",e);

    /*
        Write the coordinates of the equilibrium structure:
    */
    for(i=0;i<l;i++)
    {
        /* try vector length = 2*normal coordinate */
        wrcors(out,atnum[i],x[i],y[i],z[i],2*xs[i],2*ys[i],2*zs[i]);
    }

    /* And write all other vibrational structures:
      (skip j=0 and j=fact because they are equal to the equilibrium geometry=)
    */
    fact = 20 ;  /* 20 frames will be generated */
    ffact = (double)fact ;  
    for ( j = 1; j < fact; j++ )  /* write each frame */
    {
        fprintf(out,"    %d\n",l);  /* no of real atoms */
        fprintf(out,"* %s %15.8f\n",vibnumb,e);   /* title */

        for(i=0;i<l;i++)  /* modify coordinates along normal modes */
        {
         xx = x[i]+scale*sin(TWOPI*j/ffact)*xs[i] ;
         yy = y[i]+scale*sin(TWOPI*j/ffact)*ys[i] ;
         zz = z[i]+scale*sin(TWOPI*j/ffact)*zs[i] ;
         wrcors(out,atnum[i],xx,yy,zz,2*xs[i],2*ys[i],2*zs[i]);
        }
    }
        fclose(out);
    if(vibtotal == ALLVIBS) {k = 3*l;}
    else {k = 3*l-6;}
    if ( (all == 1) && ( vibn < k )) goto allvib ;
}

/*************************************************/
double vib_g88(fil, vibn, nat, x, y, z)
     FILE *fil ;
     int vibn, nat ;
     double x[], y[], z[] ;
{
  char line[160] ;
  char num[50] ;
  int i, ii, j, k, l, m, n ;
  double omega, atof() ;

  while( fscanf(fil,"%s",line) != EOF )
  {
    if ( strcmp(line,"(mDyne/A)") == 0 ) 
    {
    /* which column to read. vibn=2, read col 2, vibn=6 read col 1 etc*/
       k = vibn % 5 ; 
       n = (int)(vibn/5) ; /* 
       if ( k == 0 )
         {
          k = 5 ; n-- ;
         }
        
       readln(fil) ; /* advance to the end of line contaning the (mDyne/A) */
       readln(fil) ; readln(fil) ;
        
       /* How many lines are we going to skip */
       i = 0 ;
       j = 3*nat ;
       while ( i < n )
         {
         for (ii=0;ii<7;ii++) readln(fil) ;
         for(l=0;l<j;l++) readln(fil) ;
         for (ii=0;ii<2;ii++) readln(fil) ;
           i++;
         }
       /* We got the right frequency! */
       fscanf(fil,"%s",num);
       fscanf(fil,"%s",num);
       m = 1 ;
       while (m<k)
         {
         fscanf(fil,"%s",num);
         m++;
         }
       /* read the frequency */
       fscanf(fil,"%s",num);
       omega = atof(num) ; /* convert it as a floating point number */
        
       /* Go to normal mode */
       for (ii=0;ii<7;ii++) readln(fil) ;
       
       for (i=0; i<nat; i++)
          {
          fscanf(fil,"%d %d %d",&l,&l,&l);
          m = 1 ;
          while (m<k)
            {
            fscanf(fil,"%s",num);
            m++;
            }
          fscanf(fil,"%s",num); x[i] = atof(num) ;
          readln(fil);
          fscanf(fil,"%d %d %d",&l,&l,&l);
          m = 1 ;
          while (m<k) 
             {
             fscanf(fil,"%s",num);
             m++;
             }
          fscanf(fil,"%s",num); y[i] = atof(num) ;
          readln(fil);
          fscanf(fil,"%d %d %d",&l,&l,&l);
          m = 1 ;
          while (m<k)
            {
            fscanf(fil,"%s",num);
            m++;
            }
          fscanf(fil,"%s",num); z[i] = atof(num) ;
          readln(fil);
          if(debug == 1) 
            fprintf(stderr,"%d = %10.5f %10.5f %10.5f\n",i, x[i], y[i], z[i]);
          }/* End For */
       return(omega) ;
    }/* End If */
  }/* End While */
  return(omega);     /* to keep some compilers happy */
}

/*============================================
 * Find total energy from GAMESS output -- jkl's hack
 *
 */

double e_gamess(fil)
FILE *fil;
{
 char line[300];
 double ene;

 while(fscanf(fil, "%s", line) != EOF) {
   if(strcmp(line, "TOTAL") == 0) {
     fscanf(fil, "%s", line);
     if(strcmp(line, "ENERGY") == 0) {
       fscanf(fil, "%s", line);
       if(strcmp(line, "=") == 0) {
         fscanf(fil, "%s", line);
         ene = atof(line);
         }
       }
     }
   }
 return(ene);
}

/*============================================
 * Find total energy from ACES2 output 
 *
 */

double e_aces2(fil)
FILE *fil;
{
 char line[300];
 double ene;

 while(fscanf(fil, "%s", line) != EOF) {
   if(strcmp(line, "E(SCF)=") == 0) {
     fscanf(fil, "%s", line);
     ene = atof(line);
     }
   }
 return(ene);
}


/* ==================================================
 *  Read GAMESS modes  -- jkl's hack
 *
 */

double vib_gamess(fil, vibn, nat, x, y, z)
FILE *fil;
int vibn, nat;
double x[], y[], z[];
{
 char line[300];
 double freq;
 int i, j, k, n;

 rewind(fil);
 while( fscanf(fil,"%s",line) != EOF ) {
   if(strcmp(line, "DEBYE**2/AMU-ANGSTROM**2") == 0) { /* vibrations start */
     k = vibn%9;  /* get column number */
     n = (int)(vibn/9);   /* get block number */
     if(k == 0) {k = 9; n--;}
     for(i = 0; i < n; i++) {  /* skip n blocks of vibrations */
       while(fscanf(fil,"%s",line) != EOF) {
         if(strcmp(line, "FREQUENCY:") == 0) {
           break;
           }
         }
       for(j = 0; j < 4; j++) readln(fil);  /* skip numbers, freqs, ints */
       for(j = 0; j < 3*nat; j++) readln(fil); /* skip modes */
       for(j = 0; j < 13; j++) readln(fil); /* skip TRANS/ROT. SAYVETZ */
       }
     while(fscanf(fil,"%s",line) != EOF) {
       if(strcmp(line, "FREQUENCY:") == 0) {
         break;
         }
       }
     for(i = 0; i < k; i++) fscanf(fil, "%s", line);
     freq = atof(line); /* save frequency value */
     readln(fil); /* skip rest of the line */
     readln(fil); /* skip next line line */
     for(j = 0; j < nat; j++) {
       /*read first atom line */
       fscanf(fil, "%*s %*s %*s");
       for(i = 0; i < k; i++) {fscanf(fil, "%s", line);}  /* get to the k-th */
       x[j] = atof(line);
       readln(fil);  /* skip rest of the line */
       fscanf(fil, "%*s");
       for(i = 0; i < k; i++) {fscanf(fil, "%s", line);}  /* get to the k-th */
       y[j] = atof(line);
       readln(fil);  /* skip rest of the line */
       fscanf(fil, "%*s");
       for(i = 0; i < k; i++) {fscanf(fil, "%s", line);}  /* get to the k-th */
       z[j] = atof(line);
       readln(fil);  /* skip rest of the line */
       }
     }
   }
 return(freq);
 }
 
/***************************************************************/

double vib_g94(fil, vibn, nat, x, y, z)
FILE *fil ;
int vibn, nat ;
double x[], y[], z[] ;
{
 char line[160] ;
 char num[50],numx[50],numy[50],numz[50] ;
 int i, ii, j, k, l, m, n, linecount ;
 double omega, atof() ;

 while( fscanf(fil,"%s",line) != EOF )
   {
/* G94 output is in 3 column set of XYZ as showm in the example below
 Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
   1   6    -0.03  -0.02   0.00    -0.02   0.05  -0.08     0.00  -0.12   0.00
   2   6    -0.01  -0.01  -0.05    -0.04  -0.02   0.12    -0.11  -0.07   0.11
   3   6    -0.03  -0.03   0.06    -0.05   0.05  -0.17     0.14  -0.06  -0.12
   4   6    -0.01   0.00  -0.05    -0.07  -0.09   0.22    -0.10   0.04   0.12
*/                
   if ( strcmp(line,"(mDyne/A)") == 0 ) 
     {
     /* which column to read.vibn=2 read col 2, vibn=4 read col 1 etc */
     k = vibn % 3 ;  /* column number */
     n = (int)(vibn/3) ;   /* block number */
     if ( k == 0 )
       {
       k = 3 ; n-- ;
       }
        
     /* advance to the end of the line containing "(mDyne/A)" */
     readln(fil) ;
     readln(fil) ; readln(fil) ;
        
     /* How many lines are we going to skip */
     i = 0 ;
     while ( i < n )
       {
      /* Skip the following 7 lines (shown as an example only):
                     Frequencies --  2197.1392              3472.2818              3726.1166
                     Red. masses --     1.0779                 1.0497                 1.0794
                     Frc consts  --     3.0657                 7.4567                 8.8293
                     IR Inten    --     5.1595                42.9199                34.5097
                     Raman Activ --    11.2044                64.3053                26.5601
                     Depolar     --     0.7474                 0.1885                 0.7500
                     Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
                   */
       for (ii=0;ii<7;ii++) readln(fil) ;
         /* Skips num of Atom (nat) of coordinates set. See typical example below:
                       1   6    -0.03  -0.02   0.00    -0.02   0.05  -0.08     0.00  -0.12   0.00
                       2   6    -0.01  -0.01  -0.05    -0.04  -0.02   0.12    -0.11  -0.07   0.11
                       3   6    -0.03  -0.03   0.06    -0.05   0.05  -0.17     0.14  -0.06  -0.12
                       4   6    -0.01   0.00  -0.05    -0.07  -0.09   0.22    -0.10   0.04   0.12
                       .   .      .      .      .        .      .      .
                       .   .      .      .      .        .      .      .
                       nat .      .      .      .        .      .      .
                    */                
       for(l=0;l<nat;l++) readln(fil) ;
       /* Skips the following two lines
                                        vibn                 vibn+1                 vibn+2
                                        ?A                     ?A                     ?A
                    */
       for (ii=0;ii<2;ii++) readln(fil) ;
          i++;
       }
     /* We got the right frequency! */
     /* this one reads the string " Frequencies " */        
     fscanf(fil,"%s",num);
     /* this one reads the string " --" */
     fscanf(fil,"%s",num);
     m = 1 ;
     while (m<k)
       {
       fscanf(fil,"%s",num);
       m++;
       }
     /* read the frequency */
     fscanf(fil,"%s",num);
     omega = atof(num) ; /* convert it to a floating point number */
        
     /* Go to normal mode */
     /* Skip the following 7 lines (shown as an example only):
                 Frequencies --  2197.1392              3472.2818              3726.1166
                 Red. masses --     1.0779                 1.0497                 1.0794
                 Frc consts  --     3.0657                 7.4567                 8.8293
                 IR Inten    --     5.1595                42.9199                34.5097
                 Raman Activ --    11.2044                64.3053                26.5601
                 Depolar     --     0.7474                 0.1885                 0.7500
                 Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
                */
     for (ii=0;ii<7;ii++) readln(fil) ;
     /* Read the first two columns under the following heading */
     /* Atom AN */

     for (i=0; i<nat; i++)
       {
       fscanf(fil,"%d %d",&l,&l);
       m = 1 ;
       while (m<k)  /* loop for each column m < k */
         {
          fscanf(fil,"%s %s %s",numx,numy,numz);
          m++;
          }
       fscanf(fil,"%s %s %s",numx,numy,numz); /* read k-th column */
       x[i] = atof(numx) ;
       y[i] = atof(numy) ;
       z[i] = atof(numz) ;
       if(debug == 1) 
          fprintf(stderr, "%d = %10.5f %10.5f %10.5f\n",i, x[i], y[i], z[i]);
       readln(fil);                    
       }/* End For */
     return(omega) ;
     }/* End If */
  }/* End While */
  return(omega) ;  /* to keep some compilers happy */
}


/***********************************************************/
/* Read ACES2 modes -- jkl's hack
 */
double vib_aces2(fil, vibn, nat, x, y, z)
FILE *fil ;
int vibn, nat ;
double x[], y[], z[] ;
{
 char line[160] ;
 char num[50],numx[50],numy[50],numz[50] ;
 int i, ii, j, k, l, m, n, linecount ;
 double omega, atof() ;

 while( fscanf(fil,"%s",line) != EOF )
   {
   if ( strcmp(line,"Mass**-1/2") == 0 ) 
     {
      if(debug == 1) fprintf(stderr,"Line: |%s| read\n", line);
     /* which column to read.vibn=2 read col 2, vibn=4 read col 1 etc */
     k = vibn % 3 ;  /* column number */
     n = (int)(vibn/3) ;   /* block number */
     if ( k == 0 )
       {
       k = 3 ; n-- ;
       }
        
     /* advance to the representation line */
        
     /* Skip bloks of frequencies */
     i = 0 ;
     while ( i < n )
       {
       for(j = 0; j < nat + 5; j++) {  /* 5 + nat lines */
         readln(fil);
         }
       i++;
       }
     /* We got the right frequency! */
     readln(fil) ;
     readln(fil);
     readln(fil);
fprintf(stderr, "Where it wants frequency\n");
     for(i = 0; i < k; i++) {
       fscanf(fil, "%s", num);
       }
     omega = atof(num) ; /* convert it to a floating point number */

     /* skip next 3 lines */
     readln(fil); readln(fil); readln(fil);

     for (i=0; i<nat; i++)
       {
       fscanf(fil,"%*s");  /* skip element symbol */
       for(m = 0; m < k; m++) {
         fscanf(fil,"%s %s %s",numx,numy,numz);
         }
       x[i] = atof(numx) ;
       y[i] = atof(numy) ;
       z[i] = atof(numz) ;
       if(debug == 1) 
         fprintf(stderr, "%d = %10.5f %10.5f %10.5f\n",i, x[i], y[i], z[i]);
       readln(fil);                    
       }/* End For */
     return(omega) ;
     }/* End If */
  }/* End While */
  return(omega);  /* to keep some compilers happe -- this is never reached */
}
/*********************************************************/
void wrcors(out,j,x,y,z,xs,ys,zs) /* find atom symbol from atomic number and write*/
     FILE *out;
     int j;
     double x,y,z,xs,ys,zs ;
{
    /*
    XMol doesn't understand all the periodic table of elements, so 
    everything above Kr is set to X
    */
    static char atom[37][3] =
      {"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", "X "} ;

      if ( (j>0) && (j<37) ) wrat(out,atom[j-1],x,y,z,xs,ys,zs);
      else wrat(out,atom[37],x,y,z,xs,ys,zs) ;
    /* perhaps xs,ys and zs should be set to zero for too small amplitudes*/
}

/***************************************************/
void wrat(out, c, x,y,z,xs,ys,zs)  /* write single atom line */
     FILE *out;
     char *c;
     double x,y,z,xs,ys,zs ;
{
  fprintf(out,"%s   %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n",c,x,y,z,xs,ys,zs);
}

/******************************************************/
/* READ GAMESS coordinates  -- jkl's hack
 */ 
int coor_gamess(fil, atnum, x, y, z)
FILE *fil;
int atnum[];
double x[], y[], z[];
{
 char line[300];
 char s1[30], s2[30], s3[30], s4[30];
 int nat;

 nat = 0; 
 while(fscanf(fil, "%s", line) != EOF) {
   if(strcmp(line, "COORDINATES") == 0) {
     fscanf(fil, "%s", line);
     if(strcmp(line, "(BOHR)") == 0) { /* we found the coordinate table */
       readln(fil);
       readln(fil);
       while(fgets(line, 299, fil) != NULL) {
         if(debug == 1) fprintf(stderr, "Line=|%s|\n", line);
         if(strlen(line) < 30) { /* end of coordinates */
           break;
           }
         sscanf(line, "%*s %s %s %s %s", s1, s2, s3, s4);
         atnum[nat] = (int)atof(s1);
         x[nat] = atof(s2)*0.529177249;
         y[nat] = atof(s3)*0.529177249;
         z[nat] = atof(s4)*0.529177249;
         if(debug == 1) 
           fprintf(stderr, "Coors[%d] = %.4f %.4f %.4f\n", nat, x[nat],
                  y[nat], z[nat]);
         nat++;
         }
       }
     }
   }
 return(nat);
}

/******************************************************/
/* read ACES2 coordinates  -- jkl's hack
 */
int coor_aces2(fil, atnum, x, y, z)
FILE *fil;
int atnum[];
double x[], y[], z[];
{
 char line[300];
 char s1[30], s2[30], s3[30], s4[30];
 int i, j, k, nat;
 double coor;

 nat = 0; 
 while(fscanf(fil, "%s", line) != EOF) {
   if(strcmp(line, "coordinate") == 0) {
     fscanf(fil, "%s", line);
     if(strcmp(line, "input") == 0) {
       fscanf(fil, "%s", line);
       if(strcmp(line, "(Bohr)") == 0) { /* we found the coordinate table */
         readln(fil);
         readln(fil);
         readln(fil);
         readln(fil);
         readln(fil);
         /* This are initial coordinates, we use them to find nat and 
            atomic numbers for atoms */
         while(fgets(line, 299, fil) != NULL) {
           if(debug == 1) fprintf(stderr, "Line=|%s|\n", line);
           sscanf(line, "%s %s %s %s %s", line, s1, s2, s3, s4);
           if((line[1] == '-') && (line[2] == '-')) { /* we reached end line */
             return(nat);
             }
           if(atoi(s1) != 0) {  /* if not dummy atom */
             atnum[nat] = atoi(s1);
             x[nat] = atof(s2)*0.529177249;
             y[nat] = atof(s3)*0.529177249;
             z[nat] = atof(s4)*0.529177249;
             if(debug == 1) 
               fprintf(stderr, "Coors[%d] = %.4f %.4f %.4f\n",
                       nat, x[nat], y[nat], z[nat]);
             nat++;
             }
           }
         }
       }
     }
   }
 return(nat);
}


/******************************************************/
int coor_g88 ( fil, atnum, x, y, z )
     FILE *fil ;
     int atnum[] ;
     double x[],y[],z[] ;
{
      char line[160] ;
    char num1[30], num2[30],num3[30] ;
    int i, j, l, flag ;
    double atof() ;
    flag=0;
    while( fscanf(fil,"%s",line) != EOF )
    {
     if ( (strcmp(line,"Standard") == 0) || (strcmp(line,"Z-Matrix") == 0) )
       flag=1 ;
       if ((strcmp(line,"Coordinates") == 0) && (flag==1))
        {
            flag=0;
            readln(fil); readln(fil); readln(fil);
            l=0;
            while ( (getc(fil) != '-') && (l<MAXAT))
            {
                fscanf(fil,"%d %d",&i,&j) ;
                if( j != -1 ) /* if real atom */
                {
                      fscanf(fil,"%s %s %s",num1, num2, num3);
                      atnum[l]=j ;
                      x[l]=atof(num1); y[l]=atof(num2); z[l]=atof(num3);
                      /*fscanf(fil," %f %f %f",&x[l],&y[l],&z[l]);*/
                      if(debug == 1) 
                        fprintf(stderr, "l=%d %d %f %f %f\n",
                        l,atnum[l],x[l],y[l],z[l]);
                      l++ ;
                }
                readln(fil); getc(fil) ; /* unix */
               }
        }
    }
    return(l);  /* l = number of real atoms */
}

/****************************************************/

void readln (fil)  /* read/skip a line */
     FILE *fil ;
{
  int i=0;
  char line[300];
  i = 0;
  while((line[i] = getc(fil)) != '\n') {
    i++ ;
    if(i > 299) i--;
    }
  line[i] = '\0';
  if(debug == 1) {
    fprintf(stderr, "Skipping:|%s|\n", line);
    }
}

/****************************************************/

double e_g88 ( fil )
     FILE *fil ;
{
  double energy, atof() ;
  char line[132], num[40] ;
  energy = 0 ;
  while( fscanf(fil,"%s",line) != EOF ) {
    if ( (strcmp(line,"DONE:") == 0) || (strcmp(line,"Done:") == 0 ) ) {
      fscanf(fil,"%s",line) ; fscanf(fil,"%s",line) ;
      fscanf(fil,"%s",num); energy=atof(num);
      /* printf("%20.10f\n",energy); */
    }
  }
  return(energy);
}
/****************************************************/

/*  redone by jkl
Checks file format
If the file has two lines containing "mDyne/A" then it is in G88 or G92 format.
If the file has only one line containing "mDyne/A" then it is in G94 format.
NOTE: This is the most obvious thing I observed.
      There could be a better way to check for different file formats.
*/
int which_program ( fil )
FILE *fil ;
{
 char line[300], line1[300];
 char num1[30], num2[30], num3[30];
 int count, ftype, natoms, i, j;
 float x, y, z;

 ftype = UNKNOWN;
 rewind(fil);
 while(fscanf(fil, "%s", line) != EOF) {
   if(strcmp(line, "Gaussian,") == 0) {
     vibtotal = NORTVIBS;
     count=0;
     while( fscanf(fil,"%s",line) != EOF ) {
       if (strcmp(line,"(mDyne/A)") == 0 ) {
         count++ ;
         }
       }
     if (count==1) {
       ftype=G94;
       }
     else if (count==2) {
       ftype=G88;
       }
     else {
       fprintf(stderr, 
          "Error.The file does not look like a Gaussian output\n");
       exit(1);
       }
     return(ftype);
     }
   else if(strcmp(line, "GAMESS") == 0) {
     ftype=GAMESST;
     if(debug == 1) fprintf(stderr, "The GAMESS was found\n");
     vibtotal = ALLVIBS;

     return(ftype);
     }
   else if(strcmp(line, "ACES2:") == 0) {
     vibtotal = NORTVIBS;
     ftype=ACES2T;
     return(ftype);
     }

   }
 return(ftype);
}