|
|
/* 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
#include
#include
/*
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.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
-- 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
-- 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;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*/
}
/***************************************************/
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 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);
}
|
