#include 
#include 

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

-- 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;i0) && (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