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

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

#define MAXAT 100

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

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

Program extracts normal modes from Gaussian 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 (July 20, 1996)

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)

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

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

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(), e, atof() ;
	int atoi(), vibn, indx, all, g94, is_g94() ;
	char *pindx ;
	double sin() ;

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

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

	/* Get information about the input file format. Is it g92 or g94 */
	g94=is_g94(dt1);
	/* Back to the start of the file */
	rewind(dt1);
	all = 0 ;
	if (argv[2][0] == 'a' )
	{
		all = 1 ; vibn=0 ;
	}
	else vibn = atoi(argv[2]) - 1 ;
	
	if ( ( pindx = strchr(argv[1],'.') ) == NULL )
		indx = strlen(argv[1]) ;
	else
		indx = strlen(argv[1]) - strlen(pindx) ;

	ime[0] = 0 ;
	strncpy(ime,argv[1],indx) ;

	allvib:
	vibn++ ;

	sprintf(vibnumb,"%s.%d",ime,vibn);
	printf("Creating file %s\n", vibnumb);
	out = fopen(vibnumb,"w");
	if( out == NULL )
	{
		printf("\t<--- file error on %s",vibnumb);
		exit(3) ;
	}
	rewind(dt1) ;
	l = coor_g88(dt1,atnum,x,y,z);
	printf("Natoms=%d\n",l);
	rewind(dt1);
	e=e_g88(dt1);
	printf("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);
	fprintf(out,"* %s %15.8f\n",argv[1],e);

	/*
		Get normal coordinate:
	*/
	rewind(dt1) ;
	if (g94)
	{
		e=vib_g94(dt1, vibn, l, xs, ys, zs);
	}
	else
	{
		e=vib_g88(dt1, vibn, l, xs, ys, zs);
	}
	printf("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 ;
	ffact = (double)fact ;
	for ( j = 1; j < fact; j++ )
	{
		fprintf(out,"	%d\n",l);
		fprintf(out,"* %s %15.8f\n",vibnumb,e);

		for(i=0;i<l;i++)
		{
			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 ( (all == 1) && ( vibn < 3*l-6 )) 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 ) 
		{
				k = vibn % 5 ; /* which column to read. vibn=2, read col 2, vibn=6 read col 1 etc*/
		      	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);
			
					printf("%d = %10.5f %10.5f %10.5f\n",i, x[i], y[i], z[i]);
		
				}/* End For */
		      	return(omega) ;
		}/* End If */
  }/* End While */
}

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 )
  {
	    if ( strcmp(line,"(mDyne/A)") == 0 ) 
		{
/* 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
*/				
				k = vibn % 3 ; /* which column to read.vibn=2 read col 2, vibn=4 read col 1 etc */
		      	n = (int)(vibn/3) ;
		      	if ( k == 0 )
		      	{
					k = 3 ; n-- ;
		      	}
		
				readln(fil) ; /* advance to the end of the line containing "(mDyne/A)" */
		
		      	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)
					{
				  			fscanf(fil,"%s %s %s",numx,numy,numz);
				  		m++;
					}
					fscanf(fil,"%s %s %s",numx,numy,numz);
					x[i] = atof(numx) ;
					y[i] = atof(numy) ;
					z[i] = atof(numz) ;
					printf("%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 */
}

wrcors(out,j,x,y,z,xs,ys,zs)
     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*/
}

wrat(out, c, x,y,z,xs,ys,zs)
     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);
}

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 )
				{
	  				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]);*/
	  				printf("l=%d %d %f %f %f\n",l,atnum[l],x[l],y[l],z[l]);
	  				l++ ;
				}
				readln(fil); getc(fil) ; /* unix */
	   		}
		}
	}
	return(l);
}

readln (fil)
     FILE *fil ;
{
  int i=0;
  while(getc(fil) != '\n') i++ ;
}

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);
}

/* 
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 is_g94 ( fil )
     FILE *fil ;
{
  	char line[160] ;
	int count, flag ;
	count=0;
	flag=3;
	while( fscanf(fil,"%s",line) != EOF )
	{
		if ( strcmp(line,"(mDyne/A)") == 0 )
		{
			count++ ;
		}
	}
	if (count==1)
	{
		flag=1;
	}
	else if (count==2)
	{
		flag=0;
	}
	else
	{
		printf("Error.The file does not look like a Gaussian output\n"); exit(1);
	}
	return(flag);
}