gamessplugin.c

/* MACHINE GENERATED FILE, DO NOT EDIT! */

#define VMDPLUGIN molfile_gamessplugin
#define STATIC_PLUGIN 1

/***************************************************************************
 *cr
 *cr            (C) Copyright 1995-2009 The Board of Trustees of the
 *cr                        University of Illinois
 *cr                         All Rights Reserved
 *cr
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: gamessplugin.c,v $
 *      $Author: saam $       $Locker:  $             $State: Exp $
 *      $Revision: 1.144 $       $Date: 2009/06/27 00:47:00 $
 *
 ***************************************************************************/

/* *******************************************************
 *
 *          G A M E S S     P L U G I N 
 *
 * This plugin allows VMD to read GAMESS log files
 * currently only single point geometries and trajectories
 * for optimizations, saddle point runs are supported 
 *
 * ********************************************************/

#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <errno.h>
#include <time.h>
#include <math.h>

#include "gamessplugin.h"
#include "unit_conversion.h"
 
#define ANGSTROM 0
#define BOHR     1
#define SPIN_ALPHA  0
#define SPIN_BETA   1

/* #define DEBUGGING 1 */

/*
 * Error reporting macro for use in DEBUG mode
 */

#define GAMESS_DEBUG
#ifdef GAMESS_DEBUG
#define PRINTERR fprintf(stderr, "\n In file %s, line %d: \n %s \n \n", \
                            __FILE__, __LINE__, strerror(errno))
#else
#define PRINTERR (void)(0)
#endif

/*
 * Error reporting macro for the multiple fgets calls in
 * the code
 */
#define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE

#define UNK_SHELL -666
#define SPD_D_SHELL -5
#define SPD_P_SHELL -4
#define SPD_S_SHELL -3
#define SP_S_SHELL -2
#define SP_P_SHELL -1
#define S_SHELL 0
#define P_SHELL 1
#define D_SHELL 2
#define F_SHELL 3
#define G_SHELL 4
#define H_SHELL 5

#define FOUND   1
#define STOPPED 2


/* ######################################################## */
/* declaration/documentation of internal (static) functions */
/* ######################################################## */

/* this routine is the main gamess log file
 * parser responsible for static, i.e. 
 * non-trajectory information */
static int parse_static_data(gamessdata *, int *);

static void print_input_data(gamessdata *);

/* this routine checks if the current run is an
 * actual GAMESS run; returns true/false */
static int have_gamess(gamessdata *);


/* this function reads the number of processors requested */
static int get_proc_mem(gamessdata *);


/* Parse the $BASIS options*/
static int get_basis_options(gamessdata *);


/* Determine the run title line */
static int get_runtitle(gamessdata *);

/* Read the input atom definitions and geometry */
static int get_input_structure(gamessdata *data);

/* Read basis set and orbital statistics such as
 * # of shells, # of A/B orbitals, # of electrons, 
 * multiplicity and total charge */
static int get_basis_stats(gamessdata *);

/* Read the contrl group and check for
 * supported RUNTYPes. Terminate the plugin
 * if an unsupported one is encountered. */
static int get_contrl(gamessdata *);

/* Read input parameters regarding calculation of 
 * certain molecular properties such as electrostatic
 * moments and the MEP. */
static int get_properties_input(gamessdata *);

/* Read symmetry point group and highest axis */
static int get_symmetry(gamessdata *);

/* read in the $GUESS options */
static int get_guess_options(gamessdata *);

/* the function get_initial_info provides the atom number,
 * coordinates, and atom types and stores them
 * temporarily. */ 
static int get_final_properties (gamessdata *);

static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
                           int *numatoms);


/* the function get_basis we also parse the basis function section to
 * determine the number of basis functions, contraction
 * coefficients. For Pople/Huzinga style basis sets
 * this numbers are in principle fixed, and could hence
 * be provided by the the plugin itself; however, the user might
 * define his own basis/contraction coeffients and hence reading
 * them from the input file seem to be somewhat more general. */
static int get_basis (gamessdata *);


/* read all primitives for the current shell */
static int read_shell_primitives(gamessdata *, prim_t **prim,
                                 char *shellsymm, int icoeff);

/* convert shell symmetry type from char to int */
static int shellsymm_int(char symm);

/* Populate the flat arrays containing the basis set data */
static int fill_basis_arrays(gamessdata *);

static int read_first_frame(gamessdata *);

/* this subroutine scans the output file for
 * the trajectory information */
static int get_traj_frame(gamessdata *, qm_atom_t *, int);


/* returns 1 if the optimization has converged */
static int analyze_traj(gamessdata *);

/* read the number of scf iterations and the scf energies
 * for the current timestep. */
static int get_scfdata(gamessdata *, qm_timestep_t *);


/* this function parses the input file for the final
 * wavefunction and stores it in the appropriate arrays; */
static int get_wavefunction(gamessdata *, qm_timestep_t *, qm_wavefunction_t *);
static int read_coeff_block(FILE *file, int wavef_size,
                              float *wave_coeff, int *angular_momentum);

static int read_localized_orbitals(gamessdata *data);

static int get_population(gamessdata *, qm_timestep_t *);

static int get_gradient(gamessdata *, qm_timestep_t *ts);

/* Read ESP charges. */
static int get_esp_charges(gamessdata *data);

/* For runtyp=HESSIAN, this subroutine scans the file for 
 * the hessian matrix in internal coordinates 
 * as well as the internal coordinate information */
static int get_int_coords(gamessdata *);


/* For runtyp=HESSIAN, this subroutine scans the file for 
 * the cartesian hessian matrix */ 
static int get_cart_hessian(gamessdata *);


/* For runtyp=HESSIAN, this subroutine reads the frequencies
 * and intensities of the normal modes */
static int get_normal_modes(gamessdata *);


/* helper routine to chop spaces/newlines off
 * a C character string 
 *
 * TODO: This function is horrible and should
 *       be replaced by a cleaner solution */
static char* chop_string_all(char *);

static char* trimleft(char *);
static char* trimright(char *);
static void eatwhitelines(FILE *fd);

static int goto_keyline(FILE *file, const char *keystring,
                        const char *stopstring);
static int pass_keyline(FILE *file, const char *keystring,
                          const char *stopstring);
static int goto_keystring2(FILE *file, const char *keystring,
        const char *stopstring1, const char *stopstring2);
static void whereami(FILE *file);

static void thisline(FILE *file) {
  char buffer[BUFSIZ];
  long filepos;
  filepos = ftell(file);
  if (!fgets(buffer, sizeof(buffer), file)) {
    if (feof(file)) printf("HERE) EOF\n");
    else printf("HERE) ????\n");
    return;
  }
  printf("HERE) %s\n", buffer);
  fseek(file, filepos, SEEK_SET);
}

/* helper routine to chop newlines off
 * a C character string 
 * 
 * TODO: This function is horrible and should
 *       be replaced by a cleaner solution */
static char* chop_string_nl(char *);


/* skip n line at a time */
static void eatline(FILE * fd, int n)
{
  int i;
  for (i=0; i<n; i++) {
    char readbuf[1025];
    fgets(readbuf, 1024, fd);
  }
}


/* Increase wavefunction array in ts by one. */
static qm_wavefunction_t* add_wavefunction(qm_timestep_t *ts) {
  if (ts->numwave) {
    /* Add a new wavefunction */
    ts->wave = (qm_wavefunction_t *)realloc(ts->wave,
                        (ts->numwave+1)*sizeof(qm_wavefunction_t));
    ts->numwave++;
  } else {
    /* We have no wavefunction for this timestep so create one */
    ts->wave = (qm_wavefunction_t *)calloc(1, sizeof(qm_wavefunction_t));
    ts->numwave = 1;
  }
  memset(&ts->wave[ts->numwave-1], 0, sizeof(qm_wavefunction_t));

  return &ts->wave[ts->numwave-1];
}

/* Replace the n-th wavefunction in ts with the last
 * one and decrease the array length by one. */
static void replace_wavefunction(qm_timestep_t *ts, int n) {
  if (ts->numwave>=2 && n>=0 && n<ts->numwave-1) {
    qm_wavefunction_t *w1, *w2;
    w2 = &ts->wave[n];
    w1 = &ts->wave[ts->numwave-1];
    free(w2->wave_coeffs);
    free(w2->orb_energies);
    free(w2->occupancies);
    memcpy(w2, w1, sizeof(qm_wavefunction_t));
    ts->wave = (qm_wavefunction_t *) realloc(ts->wave,
                   (ts->numwave-1)*sizeof(qm_wavefunction_t));
    ts->numwave--;
  }
}


/* Delete the last wavefunction in ts */
static void del_wavefunction(qm_timestep_t *ts) {
  if (ts->numwave) {
    qm_wavefunction_t *w;
    w = &ts->wave[ts->numwave-1];
    free(w->wave_coeffs);
    free(w->orb_energies);
    free(w->occupancies);
    ts->numwave--;
    ts->wave = (qm_wavefunction_t *)realloc(ts->wave,
                        ts->numwave*sizeof(qm_wavefunction_t));    
  }
}


/* ######################################################## */
/* Functions that are needed by the molfile_plugin          */
/* interface to provide VMD with the parsed data            */
/* ######################################################## */


/***************************************************************
 *
 * Called by VMD to open the GAMESS logfile and get the number
 * of atoms.
 * We are also reading all the static (i.e. non-trajectory)
 * data here since we have to parse a bit to get the atom count
 * anyway. These data will then be provided to VMD by
 * read_gamess_metadata() and read_gamess_rundata().
 *
 * *************************************************************/
static void *open_gamess_read(const char *filename, 
                  const char *filetype, int *natoms) {

  FILE *fd;
  gamessdata *data;
  
  /* open the input file */
  fd = fopen(filename, "rb");
 
  if (!fd) {
    PRINTERR;
    return NULL;
  }

  /* allocate memory for main data structure */
  data = (gamessdata *)calloc(1,sizeof(gamessdata));

  /* make sure memory was allocated properly */
  if (data == NULL) {
    PRINTERR;
    return NULL;
  }

  data->runtype = NONE;
  data->scftype = NONE;
  data->dfttype = NONE;
  data->citype  = NONE;
  data->num_shells = 0;
  data->num_basis_funcs = 0;
  data->num_basis_atoms = 0;
  data->num_frames = 0;
  data->num_frames_sent = 0;
  data->num_frames_read = 0;
  data->trajectory_done = FALSE;
  data->opt_status = MOLFILE_QM_STATUS_UNKNOWN;
  data->have_internals = FALSE;
  data->have_cart_hessian = FALSE;
  data->have_normal_modes = FALSE;
  data->nimag = 0;
  data->num_electrons = 0;

  data->nintcoords = 0;
  data->nbonds = 0;
  data->nangles = 0;
  data->ndiheds = 0;
  data->nimprops = 0;

  data->version = 0;
  data->have_pcgamess = 0;

  /* initialize some of the character arrays */
  memset(data->basis_string,0,sizeof(data->basis_string));
  memset(data->version_string,0,sizeof(data->version_string));
  memset(data->memory,0,sizeof(data->memory));

  /* store file pointer and filename in gamess struct */
  data->file = fd;


  /* check if the file is GAMESS format; if yes
   * parse it, if no exit */
  if (have_gamess(data)==TRUE) {
    /* if we're dealing with an unsupported GAMESS
     * version, we better quit */
    if (data->version==0) {
      printf("gamessplugin) GAMESS version %s not supported. \n",
             data->version_string);
      printf("gamessplugin) .... bombing out! Sorry :( \n");
      return NULL;
    }

    /* get the non-trajectory information from the log file */    
    if (parse_static_data(data, natoms) == FALSE) 
      return NULL;
  }
  else {
    return NULL;
  }

  return data;
}


/************************************************************
 * 
 * Provide VMD with the structure of the molecule, i.e the
 * atoms coordinates names, etc.
 *
 *************************************************************/
static int read_gamess_structure(void *mydata, int *optflags, 
                      molfile_atom_t *atoms) 
{
  gamessdata *data = (gamessdata *)mydata;
  qm_atom_t *cur_atom;
  molfile_atom_t *atom;
  int i = 0;
 
  *optflags = MOLFILE_ATOMICNUMBER; /* no optional data */
  if (data->have_mulliken) 
    *optflags |= MOLFILE_CHARGE;

  /* all the information I need has already been read in
   * via the initial scan and I simply need to copy 
   * everything from the temporary arrays into the 
   * proper VMD arrays.
   * Since there are no atom names in the GAMESS output
   * I use the atom type here --- maybe there is a better
   * way to do this !!?? */

  /* get initial pointer for atom array */
  cur_atom = data->initatoms;

  for(i=0; i<data->numatoms; i++) {
    atom = atoms+i;
    strncpy(atom->name, cur_atom->type, sizeof(atom->name)); 
    strncpy(atom->type, cur_atom->type, sizeof(atom->type));
    strncpy(atom->resname,"", sizeof(atom->resname)); 
    atom->resid = 1;
    atom->chain[0] = '\0';
    atom->segid[0] = '\0';
    atom->atomicnumber = cur_atom->atomicnum;
#ifdef DEBUGGING
    printf("gamessplugin) atomicnum[%d] = %d\n", i, atom->atomicnumber);
#endif

    if (data->have_mulliken)
      atom->charge = data->qm_timestep->mulliken_charges[i];

    cur_atom++;
  }
 
  return MOLFILE_SUCCESS; 
}


/*****************************************************
 *
 * provide VMD with the sizes of the QM related
 * data structure arrays that need to be made
 * available
 *
 *****************************************************/
static int read_gamess_metadata(void *mydata, 
    molfile_qm_metadata_t *metadata) {

  gamessdata *data = (gamessdata *)mydata;

  if (data->runtype == HESSIAN) {
    metadata->ncart = (3*data->numatoms);
    metadata->nimag = data->nimag;             
   
    if (data->have_internals) {
      metadata->nintcoords = data->nintcoords; 
    } else {
      metadata->nintcoords = 0;
    }
  }
  else {
    metadata->ncart = 0;
    metadata->nimag = 0;
    metadata->nintcoords = 0;
  }

  /* orbital data */
  metadata->num_basis_funcs = data->num_basis_funcs;
  metadata->num_basis_atoms = data->num_basis_atoms;
  metadata->num_shells      = data->num_shells;
  metadata->wavef_size      = data->wavef_size;  

#if vmdplugin_ABIVERSION > 11
  /* system and run info */
  metadata->have_sysinfo = 1;

  /* charges */
  metadata->have_esp = data->have_esp;

  /* hessian info */
  metadata->have_carthessian = data->have_cart_hessian;
  metadata->have_internals   = data->have_internals;

  /* normal mode info */
  metadata->have_normalmodes = data->have_normal_modes;
#endif

  return MOLFILE_SUCCESS;
}


/******************************************************
 * 
 * Provide VMD with the static (i.e. non-trajectory)
 * data. That means we are filling the molfile_plugin
 * data structures.
 *
 ******************************************************/
static int read_gamess_rundata(void *mydata, 
                               molfile_qm_t *qm_data) {

  gamessdata *data = (gamessdata *)mydata;
  int i, j;
  int ncart;
  molfile_qm_hessian_t *hessian_data = &qm_data->hess;
  molfile_qm_basis_t   *basis_data   = &qm_data->basis;
  molfile_qm_sysinfo_t *sys_data     = &qm_data->run;

  /* fill in molfile_qm_hessian_t */
  if (data->runtype == HESSIAN) {
    ncart = 3*data->numatoms;

    /* Hessian matrix in cartesian coordinates */
    if (data->have_cart_hessian) {
      for (i=0; i<ncart; i++) {
        for (j=0; j<=i; j++) {
          hessian_data->carthessian[ncart*i+j] = data->carthessian[ncart*i+j];
          hessian_data->carthessian[ncart*j+i] = data->carthessian[ncart*i+j];
        }
      }
    }

    /* Hessian matrix in internal coordinates */
    if (data->have_internals) {
      for (i=0; i<(data->nintcoords)*(data->nintcoords); i++) {
        hessian_data->inthessian[i] = data->inthessian[i];
      }
    }

    /* wavenumbers, intensities, normal modes */
    if (data->have_normal_modes) {
      for (i=0; i<ncart*ncart; i++) {
        hessian_data->normalmodes[i] = data->normal_modes[i];
      }
      for (i=0; i<ncart; i++) {
        hessian_data->wavenumbers[i] = data->wavenumbers[i];
        hessian_data->intensities[i] = data->intensities[i];
      }
    }

    /* imaginary modes */
    for (i=0; i<data->nimag; i++) {
      hessian_data->imag_modes[i] = data->nimag_modes[i];
    }
  }

  /* fill in molfile_qm_sysinfo_t */
  sys_data->runtype = data->runtype;
  sys_data->scftype = data->scftype;
  sys_data->nproc = data->nproc;
  sys_data->num_electrons = data->num_electrons;
  sys_data->totalcharge = data->totalcharge;
  sys_data->multiplicity = data->multiplicity;
  sys_data->num_occupied_A = data->num_occupied_A;
  sys_data->num_occupied_B = data->num_occupied_B;

#if vmdplugin_ABIVERSION > 11
  sys_data->status = data->opt_status;

  if (data->have_esp) {
    for (i=0; i<data->numatoms; i++) {
      sys_data->esp_charges[i] = data->esp_charges[i];
    }
  }
#endif

  strncpy(sys_data->basis_string, data->basis_string,
          sizeof(sys_data->basis_string));

  sys_data->memory = 0; /* XXX fixme */

  strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
  strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
  strncpy(sys_data->version_string, data->version_string,
          sizeof(sys_data->version_string));

#if vmdplugin_ABIVERSION > 11
  /* fill in molfile_qm_basis_t */
  if (data->num_basis_funcs) {
    for (i=0; i<data->num_basis_atoms; i++) {
      basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
      basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
    }
    
    for (i=0; i<data->num_shells; i++) {
      basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
      basis_data->shell_symmetry[i] = data->shell_symmetry[i];
    }
    
    for (i=0; i<2*data->num_basis_funcs; i++) {
      basis_data->basis[i] = data->basis[i];
    }

    for (i=0; i<3*data->wavef_size; i++) {
      basis_data->angular_momentum[i] = data->angular_momentum[i];
    }
  }
#endif
 
  return MOLFILE_SUCCESS;
}




#if vmdplugin_ABIVERSION > 11

/***********************************************************
 * Provide non-QM metadata for next timestep. 
 * Required by the plugin interface.
 ***********************************************************/
static int read_timestep_metadata(void *mydata,
                                  molfile_timestep_metadata_t *meta) {
  meta->count = -1;
  meta->has_velocities = 0;

  return MOLFILE_SUCCESS;
}

/***********************************************************
 * Provide QM metadata for next timestep. 
 * This actually triggers reading the entire next timestep
 * since we have to parse the whole timestep anyway in order
 * to get the metadata. So we store the read data locally
 * and hand them to VMD when requested by read_timestep().
 *
 ***********************************************************/
static int read_qm_timestep_metadata(void *mydata,
                                    molfile_qm_timestep_metadata_t *meta) {
  int have = 0;

  gamessdata *data = (gamessdata *)mydata;
  printf("gamessplugin) read_qm_timestep_metadata(): %i/%i/%i\n",
         data->num_frames, 
         data->num_frames_read,
         data->num_frames_sent);

  meta->count = -1; /* Don't know the number of frames yet */
  meta->has_gradient = 0;

  if (data->num_frames_read > data->num_frames_sent) {
    have = 1;
  }
  else if (data->num_frames_read < data->num_frames) {
    printf("gamessplugin) Probing timestep %i\n", data->num_frames_read);

    have = get_traj_frame(data, data->initatoms, data->numatoms);
  }

  if (have) {
    int i;
    qm_timestep_t *cur_ts;

    /* get a pointer to the current qm timestep */
    cur_ts = data->qm_timestep+data->num_frames_sent;
    printf("gamessplugin) Approved timestep %i\n", data->num_frames_sent);
    
    meta->num_scfiter  = cur_ts->num_scfiter;

    for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
      char typestr[MOLFILE_BUFSIZ];
      strcpy(typestr, cur_ts->wave[i].info);
      
      printf("gamessplugin) num_orbitals_per_wavef[%d/%d]=%3d type=%s\n",
             i+1, cur_ts->numwave, cur_ts->wave[i].num_orbitals, typestr);
      meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
    }
    meta->num_wavef  = cur_ts->numwave;
    meta->wavef_size = data->wavef_size;

  } else {
    meta->num_scfiter  = 0;
    meta->num_orbitals_per_wavef[0] = 0;
    meta->num_wavef = 0;
    meta->wavef_size = 0;

    data->trajectory_done = TRUE;
  }

  return MOLFILE_SUCCESS;
}


/***********************************************************
 *
 * This function provides the data of the next timestep.
 * Here we actually don't read the data from file, that had
 * to be done already upon calling read_timestep_metadata().
 * Instead we copy the stuff from the local data structure
 * into the one's provided by VMD.
 *
 ***********************************************************/
static int read_timestep(void *mydata, int natoms, 
       molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
                   molfile_qm_timestep_t *qm_ts) 
{
  gamessdata *data = (gamessdata *)mydata;
  qm_atom_t *cur_atom;
  int i = 0;
  qm_timestep_t *cur_ts;

  if (data->trajectory_done == TRUE) return MOLFILE_ERROR;

  printf("gamessplugin) Sending timestep %i\n", data->num_frames_sent);

  /* initialize pointer for temporary arrays */
  cur_atom = data->initatoms; 
  
  /* copy the coordinates */
  for(i=0; i<natoms; i++) {
    ts->coords[3*i  ] = cur_atom->x;
    ts->coords[3*i+1] = cur_atom->y;
    ts->coords[3*i+2] = cur_atom->z; 
    cur_atom++;
  }    
  
  /* get a convenient pointer to the current qm timestep */
  cur_ts = data->qm_timestep+data->num_frames_sent;

  /* store the SCF energies */
  for (i=0; i<cur_ts->num_scfiter; i++) {
    qm_ts->scfenergies[i] = cur_ts->scfenergies[i];
  }

  /* store the wave function and orbital energies */
  if (cur_ts->wave) {
    for (i=0; i<cur_ts->numwave; i++) {
      qm_wavefunction_t *wave = &cur_ts->wave[i];
      qm_ts->wave[i].type        = wave->type;
      qm_ts->wave[i].spin        = wave->spin;
      qm_ts->wave[i].excitation  = wave->exci;
      strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);

      if (wave->wave_coeffs && wave->orb_energies) {
        memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
               wave->num_orbitals*data->wavef_size*sizeof(float));
        memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
               wave->num_orbitals*sizeof(float));
      }
    }
  }

  if (data->runtype == ENERGY || data->runtype == HESSIAN) {
    /* We have only a single point */
    data->trajectory_done = TRUE;
  }

  data->num_frames_sent++;

  return MOLFILE_SUCCESS;
}
#endif



/**********************************************************
 *
 * clean up when done and free all the memory do avoid
 * memory leaks
 *
 **********************************************************/
static void close_gamess_read(void *mydata) {

  gamessdata *data = (gamessdata *)mydata;
  int i, j;
  fclose(data->file);

  free(data->initatoms);
  free(data->basis);
  free(data->shell_symmetry);
  free(data->atomicnum_per_basisatom);
  free(data->num_shells_per_atom);
  free(data->num_prim_per_shell);
  free(data->mulliken_charges);
  free(data->esp_charges);
  free(data->bonds);
  free(data->angles);
  free(data->dihedrals);
  free(data->impropers);
  free(data->internal_coordinates);
  free(data->bond_force_const);
  free(data->angle_force_const);
  free(data->dihedral_force_const);
  free(data->improper_force_const);
  free(data->inthessian);
  free(data->carthessian);
  free(data->wavenumbers);
  free(data->intensities);
  free(data->normal_modes);
  free(data->angular_momentum);
  free(data->filepos_array);

  if (data->basis_set) {
    for(i=0; i<data->num_basis_atoms; i++) {
      for (j=0; j<data->basis_set[i].numshells; j++) {
        free(data->basis_set[i].shell[j].prim);
      }
      free(data->basis_set[i].shell);
    } 
    free(data->basis_set);
  }

  for (i=0; i<data->num_frames; i++) {
    free(data->qm_timestep[i].scfenergies);
    free(data->qm_timestep[i].gradient);
    free(data->qm_timestep[i].mulliken_charges);
    for (j=0; j<data->qm_timestep[i].numwave; j++) {
      free(data->qm_timestep[i].wave[j].wave_coeffs);
      free(data->qm_timestep[i].wave[j].orb_energies);
      free(data->qm_timestep[i].wave[j].occupancies);
    }
    free(data->qm_timestep[i].wave);
  }
  free(data->qm_timestep);
  
  free(data);
}

/* ####################################################### */
/*             End of API functions                        */
/* The following functions actually do the file parsing.   */
/* ####################################################### */



/********************************************************
 *
 * Main gamess log file parser responsible for static,  
 * i.e. non-trajectory information.
 *
 ********************************************************/
static int parse_static_data(gamessdata *data, int *natoms) 
{
  /* Read # of procs and amount of requested memory */
  if (!get_proc_mem(data))        return FALSE;

  /* Read the $BASIS options */
  if (!get_basis_options(data))   return FALSE;

  /* Read the run title */
  if (!get_runtitle(data))        return FALSE;

  /* Read the input atom definitions and geometry */
  if (!get_input_structure(data)) return FALSE;

  /* Read the basis set */
  if (!get_basis(data))           return FALSE; 

  /* Read the number of orbitals, electrons, 
   * charge, multiplicity, ... */
  if (!get_basis_stats(data))     return FALSE;

  /* Read the $CONTRL group; here we determine the
   * RUNTYP and bomb out if we don't support
   * it; also read in some control flow related
   * info as well as the COORD type */
  if (!get_contrl(data))          return FALSE;

  /* Read input parameters regarding calculation of 
   * certain molecular properties such as electrostatic
   * moments and the MEP. */
  if (!get_properties_input(data)) return FALSE;

  /* Read symmetry point group and highest axis */
  if (!get_symmetry(data))         return FALSE;

  /* Read the $GUESS options */
  if (!get_guess_options(data))    return FALSE;

  /* Find the end of the trajectory and count the
   * frames on the way.
   * If no regular end is found we won't find any
   * properties to read and return. */
  if (!analyze_traj(data)) {
    printf("gamessplugin) WARNING: Truncated or abnormally terminated file!\n\n");
  }



  /* provide VMD with the proper number of atoms */
  *natoms = data->numatoms;

  /* Read the first frame*/
  read_first_frame(data);


  /* Read the properties at the end of a calculation */
  get_final_properties(data);

  printf("gamessplugin) num_frames_read = %i\n", data->num_frames_read);
  printf("gamessplugin) num_frames_sent = %i\n", data->num_frames_sent);

#ifdef DEBUGGING 
  /* Test print the parsed data in same format as logfile */
  print_input_data(data);
#endif

  return TRUE;
}


#define TORF(x) (x ? "T" : "F")

static void print_input_data(gamessdata *data) {
  int i, j, k;
  int primcount=0;
  int shellcount=0;

  printf("\nDATA READ FROM FILE:\n\n");
  printf(" %10s WORDS OF MEMORY AVAILABLE\n", data->memory);
  printf("\n");
  printf("     BASIS OPTIONS\n");
  printf("     -------------\n");
  printf("%s\n", data->basis_string);
  printf("\n\n\n");
  printf("     RUN TITLE\n");
  printf("     ---------\n");
  printf(" %s\n", data->runtitle);
  printf("\n");
  printf(" THE POINT GROUP OF THE MOLECULE IS %s\n", "XXX");
  printf(" THE ORDER OF THE PRINCIPAL AXIS IS %5i\n", 0);
  printf("\n");
  printf(" YOUR FULLY SUBSTITUTED Z-MATRIX IS\n");
  printf("\n");
  printf(" THE MOMENTS OF INERTIA ARE (AMU-ANGSTROM**2)\n");
  printf(" IXX=%10.3f   IYY=%10.3f   IZZ=%10.3f\n", 0.0, 0.0, 0.0);
  printf("\n");
  printf(" ATOM      ATOMIC                      COORDINATES (BOHR)\n");
  printf("           CHARGE         X                   Y                   Z\n");
  for (i=0; i<data->numatoms; i++) {
    printf(" %-8s %6d", data->initatoms[i].type, data->initatoms[i].atomicnum);
    
    printf("%17.10f",   ANGS_TO_BOHR*data->initatoms[i].x);
    printf("%20.10f",   ANGS_TO_BOHR*data->initatoms[i].y);
    printf("%20.10f\n", ANGS_TO_BOHR*data->initatoms[i].z);
  }
  printf("\n");
  printf("     ATOMIC BASIS SET\n");
  printf("     ----------------\n");
  printf(" THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED\n");
  printf(" THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY\n");
  printf("\n");
  printf("  SHELL TYPE  PRIMITIVE        EXPONENT          CONTRACTION COEFFICIENT(S)\n");
  printf("\n");

#if 0
  for (i=0; i<data->numatoms; i++) {
    printf("%-8s\n\n", data->initatoms[i].type);
    printf("\n");
    printf("nshells=%d\n", data->num_shells_per_atom[i]);

    for (j=0; j<data->num_shells_per_atom[i]; j++) {
      printf("nprim=%d\n", data->num_prim_per_shell[shellcount]);

      for (k=0; k<data->num_prim_per_shell[shellcount]; k++) {
        printf("%6d   %d %7d %22f%22f\n", j, data->shell_symmetry[shellcount],
               primcount+1, data->basis[2*primcount], data->basis[2*primcount+1]);
        primcount++;
      }

      printf("\n");
      shellcount++;
    }
  }
#endif
  printf("gamessplugin) =================================================================\n");
  for (i=0; i<data->num_basis_atoms; i++) {
    printf("%-8s (%10s)\n\n", data->initatoms[i].type, data->basis_set[i].name);
    printf("\n");

    for (j=0; j<data->basis_set[i].numshells; j++) {

      for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
        printf("%6d   %d %7d %22f%22f\n", j,
               data->basis_set[i].shell[j].symmetry,
               primcount+1,
               data->basis_set[i].shell[j].prim[k].exponent,
               data->basis_set[i].shell[j].prim[k].contraction_coeff);
        primcount++;
      }

      printf("\n");
      shellcount++;
    }
  }
  printf("\n");
  printf(" TOTAL NUMBER OF BASIS SET SHELLS             =%5d\n", data->num_shells);
  printf(" NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =%5d\n", data->wavef_size);
  printf(" NUMBER OF ELECTRONS                          =%5d\n", data->num_electrons);
  printf(" CHARGE OF MOLECULE                           =%5d\n", data->totalcharge);
  printf(" SPIN MULTIPLICITY                            =%5d\n", data->multiplicity);
  printf(" NUMBER OF OCCUPIED ORBITALS (ALPHA)          =%5d\n", data->num_occupied_A);
  printf(" NUMBER OF OCCUPIED ORBITALS (BETA )          =%5d\n", data->num_occupied_B);
  printf(" TOTAL NUMBER OF ATOMS                        =%5i\n", data->numatoms);
  printf("\n");
}



/**********************************************************
 *
 * this subroutine checks if the provided files is
 * actually a GAMESS file;
 *
 **********************************************************/
static int have_gamess(gamessdata *data) 
{
  char word[3][BUFSIZ];
  char buffer[BUFSIZ];
  int day, year;
  char month[BUFSIZ], rev[BUFSIZ];
  int i = 0;
 
  buffer[0] = '\0';
  for (i=0; i<3; i++) word[i][0] = '\0';


  /* check if the file is GAMESS format 
   * for now I just read line by line until 
   * the word GAMESS appears                */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],&word[2][0]);

  } while( strcmp(&word[1][0],"GAMESS") || 
           strcmp(&word[2][0],"VERSION") );

  /* extract the version number if possible; otherwise
   * return empty string */
  if (strstr(buffer,"=") != NULL) {
    strncpy(data->version_string,strstr(buffer,"=")+2,16);
    data->version_string[16] = '\0';
  }
  
  /* determine if we're dealing with pre-"27 JUN 2005"
   * version */
  sscanf(data->version_string,"%d %s %d %s",&day, month, &year, rev);
  
  if ( ( year >= 2006 ) ||
       ( year == 2005 && !strcmp(month,"JUN") ) ||
       ( year == 2005 && !strcmp(month,"NOV") ) ||
       ( year == 2005 && !strcmp(month,"DEC") ) )
  {
    data->version = 2;
  }
  else { 
    data->version = 1;
  }

  /* scan the next 20 lines and see if we're dealing with
   * PC Gamess output */
  for (i=0; i<20; i++) {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],&word[2][0]);

    if (!strcmp(&word[1][0],"PC") && !strcmp(&word[2][0],"GAMESS"))
    {
      data->have_pcgamess = 1;
      break;
    }
  }

  /* short messsage to stdout */
  if ( data->have_pcgamess == 1 ) {
    printf("gamessplugin) Detected PC GAMESS output :)\n");
  }
  else {
    printf("gamessplugin) Detected GAMESS format :)\n");
  }
  printf("gamessplugin) GAMESS version = %s \n", 
         data->version_string);

  return TRUE;
}


/**********************************************************
 *
 * this subroutine reads the number of procs and the amount
 * of memory requested
 *
 **********************************************************/
static int get_proc_mem(gamessdata *data) {

  char word[4][BUFSIZ];
  char buffer[BUFSIZ];
  char *temp;
  int nproc;
  int i;

  buffer[0] = '\0';
  for (i=0; i<3; i++) word[i][0] = '\0';

  rewind(data->file);

  
  /* scan for the number of processors; here we need
   * distinguish between vanilla Gamess and PC Gamess */
  if (data->have_pcgamess == 1) {
    /* for now we fake ncpu = 1 until we know exactly
     * how the output format looks like */
    nproc = 1;
  }
  else {
    do {
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %d %s",&word[0][0],&nproc,&word[1][0]);

      if (!strcmp(&word[0][0],"Initiating") &&
          !strcmp(&word[1][0],"compute")) {
        break;
      }

      /* Some versions */
      if (!strcmp(&word[0][0],"PARALLEL") &&
          !strcmp(&word[0][1],"RUNNING")) {
        sscanf(buffer,"%*s %*s %*s %*s %d %*s",&nproc);
        break;
      }

    } while (strcmp(&word[0][0],"ECHO") || 
             strcmp(&word[1][0],"THE") );
  }
  
  /* store the number of processors */
  data->nproc = nproc;

  
  /* scan for the amount of memory requested */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);

  } while( strcmp(&word[0][0],"$SYSTEM") || 
           strcmp(&word[1][0],"OPTIONS") );

  eatline(data->file, 1);


  /* next line contains the amount of memory requested,
   * vanilla Gamess and PC Gamess need separate treatment */
  if (data->have_pcgamess == 1) {
    GET_LINE(buffer, data->file);

    /* store it */
    if ((temp = strstr(buffer,"MEMORY=")+8)==NULL) return FALSE;
    strncpy(data->memory,chop_string_all(temp),sizeof(data->memory));
  }
  else {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],&word[2][0]);

    /* store it */
    strncpy(data->memory,&word[2][0],sizeof(data->memory));
  }

  printf("gamessplugin) GAMESS used %d compute processes \n", nproc);
  printf("gamessplugin) GAMESS used %s words of memory \n", data->memory);

  return TRUE;
}


/**********************************************************
 *
 * Extract the $BASIS options
 *
 **********************************************************/
static int get_basis_options(gamessdata *data) {

  char word[3][BUFSIZ];
  char buffer[BUFSIZ];
  char diffuse[BUFSIZ];
  char polarization[BUFSIZ];
  int i = 0;

  buffer[0] = '\0';
  diffuse[0] = '\0';
  polarization[0] = '\0';
  for (i=0; i<3; i++) word[i][0] = '\0';

  /* to be safe let's rewind the file */
  rewind(data->file);

  /* start scanning */
  if (!goto_keyline(data->file, "BASIS OPTIONS", "RUN TITLE")) {
    /* No Basis options section found
     * (basis was entered explicitly) */
    return TRUE;
  }

  eatline(data->file, 2);


  /* the first string in the current line contains the
   * GBASIS used; copy it over into the gbasis variable
   * of gamessdata */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%s %s %s", &word[0][0], &word[1][0], &word[2][0]);
 
  strncpy(data->gbasis, (&word[0][0])+7, sizeof(data->gbasis));


  /* in case we're using a pople style basis set, i.e. 
   * GBASIS=N311,N31,N21 or STO we also scan for the number 
   * of gaussians, as well as p,d,f and diffuse functions
   * and use this info to assemble a "basis set string" */
  if ( !strncmp(data->gbasis,"N311",sizeof(data->gbasis)) ||
       !strncmp(data->gbasis,"N31",sizeof(data->gbasis)) ||
       !strncmp(data->gbasis,"N21",sizeof(data->gbasis)) ||
       !strncmp(data->gbasis,"STO",sizeof(data->gbasis)) ) 
  {
    int ngauss, npfunc, ndfunc, nffunc;
    int diffs=FALSE, diffsp=FALSE;
    char torf;

    /* word[2] read previously should still contain the
     * number of gaussians used */
    ngauss = atoi(&word[2][0]);


    /* the next line gives us the d,f and diffuse sp
     * functions */
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%*s %d %*s %d %*s %c", &ndfunc, &nffunc, &torf);

    /* convert GAMESS' .TRUE./.FALSE. for DIFFSP into 1/0 */
    if (torf=='T') diffsp = TRUE;


    /* the next line gives us the p and diffuse s
     * functions */
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%*s %d %*s %c", &npfunc, &torf);

    /* convert GAMESS' .TRUE./.FALSE. for DIFFSP into 1/0 */
    if (torf=='T') diffs = TRUE;


    /* now we need some logic to assemble this info into
     * some nice looking string a la "6-31G*" */

    /* create the diffuse function string */
    if (diffs && diffsp) {
      strncpy(diffuse,"++",sizeof(diffuse));
    }
    else if (diffsp) {
      strncpy(diffuse,"+",sizeof(diffuse));
    }
    else {
      strncpy(diffuse,"",sizeof(diffuse));
    }


    /* create the polarization function string */
    if (npfunc>0 && ndfunc>0 && nffunc>0) {
      sprintf(polarization, "(%dp,%dd,%df)", npfunc, ndfunc, nffunc);
    } else if (npfunc>0 && ndfunc>0) {
      sprintf(polarization, "(%dp,%dd)", npfunc, ndfunc);
    }
    else if (npfunc>0) {
      sprintf(polarization, "(%dp)", npfunc);
    }
    else if (ndfunc>0) {
      sprintf(polarization, "(%dd)", ndfunc);
    } 
    else {
      strncpy(polarization, "", sizeof(polarization));
    } 

    /* assemble the bits */ 
    if (!strcmp(data->gbasis, "STO")) {
      sprintf(data->basis_string, "STO-%dG%s%s",
              ngauss, diffuse, polarization);
    }
    else {
      sprintf(data->basis_string, "%d-%s%sG%s",
              ngauss, (data->gbasis+1), diffuse, 
              polarization);
    }      
  }

  /* cc-pVnZ and cc-pCVnZ */
  else if (!strncmp(data->gbasis, "CC",  2)) {
    strcpy(data->basis_string, "cc-p");
    if (strlen(data->gbasis)==4 && data->gbasis[3]=='C') {
      strcat(data->basis_string, "C");
    }
    strcat(data->basis_string, "V");
    strncat(data->basis_string, &data->gbasis[2], 1);
    strcat(data->basis_string, "Z");
  }

  /* aug-cc-pVnZ and aug-cc-pCVnZ */
  else if (!strncmp(data->gbasis, "ACC", 3)) {
    strcpy(data->basis_string, "aug-cc-p");
    if (strlen(data->gbasis)==5 && data->gbasis[4]=='C') {
      strcat(data->basis_string, "C");
    }
    strcat(data->basis_string, "V");
    strncat(data->basis_string, &data->gbasis[3], 1);
    strcat(data->basis_string, "Z");
  }

  /* for non Pople style basis sets we just use the GBASIS
   * for the basis string;
   * TODO: make the basis_string more comprehensive for non
   *       pople-style basis sets */
  else {
    strncpy(data->basis_string,data->gbasis,
            sizeof(data->basis_string));
  }

  return TRUE;
}



/**********************************************************
 *
 * Extract the run title line
 *
 **********************************************************/
static int get_runtitle(gamessdata *data) {

  char word[2][BUFSIZ];
  char buffer[BUFSIZ];

  /* initialize arrays */
  word[0][0] = '\0';
  word[1][0] = '\0';
  buffer[0] = '\0';

  /* look for RUN TITLE section */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);

  } while (strcmp(&word[0][0],"RUN") || 
           strcmp(&word[1][0],"TITLE"));

  eatline(data->file, 1);

  GET_LINE(buffer, data->file);
  strncpy(data->runtitle,chop_string_nl(buffer),sizeof(data->runtitle));

  return TRUE;
} 


/* Read the input atom definitions and geometry */
static int get_input_structure(gamessdata *data) {
  char buffer[BUFSIZ];
  char word[4][BUFSIZ];
  int i, numatoms;
  int bohr;

  buffer[0] = '\0';
  for (i=0; i<4; i++) word[i][0] = '\0';

  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s %s",&word[0][0],&word[1][0],
           &word[2][0],&word[3][0]);

  } while(strcmp(&word[0][0],"ATOM") || 
          strcmp(&word[1][0],"ATOMIC"));

  
  /* test if coordinate units are Bohr */
  bohr = (!strcmp(&word[3][0],"(BOHR)"));

  eatline(data->file, 1);

  /* we don't know the number of atoms yet */
  numatoms=-1;  

  /* Read the coordinate block */
  if (get_coordinates(data->file, &data->initatoms, bohr, &numatoms))
    data->num_frames_read = 0;
  else return FALSE;


  /* store number of atoms in data structure */
  data->numatoms = numatoms;

  return TRUE; 
}


/**********************************************************
 *
 * Read data from the $CONTRL card
 *
 **********************************************************/
static int get_contrl(gamessdata *data) {

  char word[3][BUFSIZ];
  char buffer[BUFSIZ];
  char *temp;

  word[0][0] = '\0';
  word[1][0] = '\0';
  word[2][0] = '\0';
  buffer[0] = '\0';


  /* start scanning; currently we support
   * RUNTYP = ENERGY, OPTIMIZE, SADPOINT, HESSIAN, SURFACE */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);

  } while( strcmp(&word[0][0],"\044CONTRL") || 
           strcmp(&word[1][0],"OPTIONS") );

  eatline(data->file, 1);

  /* current line contains SCFTYP, RUNTYP, EXETYP info; scan it */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);

  /* check for supported RUNTYPs */
  if      (!strcmp(&word[1][0],"RUNTYP=ENERGY")) {
    data->runtype = ENERGY;
  }
  else if (!strcmp(&word[1][0],"RUNTYP=OPTIMIZE")) {
    data->runtype = OPTIMIZE;
  }
  else if (!strcmp(&word[1][0],"RUNTYP=SADPOINT")) {
    data->runtype = SADPOINT;
  }
  else if (!strcmp(&word[1][0],"RUNTYP=HESSIAN")) {
    data->runtype = HESSIAN;
  }
  else if (!strcmp(&word[1][0],"RUNTYP=SURFACE")) {
    data->runtype = SURFACE;
  }
  else if (!strcmp(&word[1][0],"RUNTYP=GRADIENT")) {
    data->runtype = GRADIENT;
  }
  else {
    printf("gamessplugin) The %s is currently not supported \n",
           &word[1][0]);
  /*   return FALSE; */
  }
  printf("gamessplugin) File generated via %s \n",&word[1][0]);


  /* determine SCFTYP */
  if (!strcmp(&word[0][0],"SCFTYP=RHF")) {
    data->scftype = RHF;
  }
  else if (!strcmp(&word[0][0],"SCFTYP=UHF")) {
    data->scftype = UHF;
  }
  else if (!strcmp(&word[0][0],"SCFTYP=ROHF")) {
    data->scftype = ROHF;
  }
  else if (!strcmp(&word[0][0],"SCFTYP=GVB")) {
    data->scftype = GVB;
  }
  else if (!strcmp(&word[0][0],"SCFTYP=MCSCF")) {
    data->scftype = MCSCF;
  }
  else if (!strcmp(&word[0][0],"SCFTYP=NONE")) {
    data->scftype = NONE;
  }
  else {
    /* if we don't find a supported SCFTYP we bomb out; this
     * might be a little drastic */
    printf("gamessplugin) %s is currently not supported \n",
           &word[0][0]);
    return FALSE;
  }
  printf("gamessplugin) Type of wavefunction used %s \n",
         &word[0][0]);

  /* scan for MPLEVL, CITYP, CCTYP, VBTYP info; */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%s %s %*s %s",&word[0][0],&word[1][0],&word[2][0]);

  if (!strcmp(&word[0][0],"MPLEVL=")) {
    /* Moller-Plesset perturbation level */
    printf("gamessplugin) MP perturbation level %s \n",&word[1][0]);
    data->mplevel = atoi(&word[1][0]);

    /* determine CITYP */
    if      (!strcmp(&word[2][0],"=NONE"))  data->citype = NONE;
    else if (!strcmp(&word[2][0],"=CIS"))   data->citype = CIS;
    else if (!strcmp(&word[2][0],"=ALDET")) data->citype = ALDET;
    else if (!strcmp(&word[2][0],"=ORMAS")) data->citype = ORMAS;
    else if (!strcmp(&word[2][0],"=GUGA"))  data->citype = GUGA;
    else if (!strcmp(&word[2][0],"=FSOCI")) data->citype = FSOCI;
    else if (!strcmp(&word[2][0],"=GENCI")) data->citype = GENCI;
    else                                    data->citype = UNKNOWN;
    printf("gamessplugin) CI method %s \n",&word[2][1]);

    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);
  }

  /* scan for DFTTYP, TDDFT info; */
  if (!strncmp(&word[0][0],"DFTTYP=", 7)) {
    printf("gamessplugin) Density functional used is %s \n",&word[0][7]);
    GET_LINE(buffer, data->file);
  }


  /* find the coordinate type in next line */
  while ( (temp=strstr(buffer,"COORD =")) == NULL ) {
    GET_LINE(buffer, data->file);;
  }
  strncpy(data->geometry,chop_string_all(temp+7),
          sizeof(data->geometry)); 
  printf("gamessplugin) Coordinate type used is %s \n", data->geometry);

  return TRUE;
}

/* Read input parameters regarding calculation of 
 * certain molecular properties such as electrostatic
 * moments and the MEP. */
static int get_properties_input(gamessdata *data) {
  /* TODO!! */
  return TRUE;
}

/* Read symmetry point group and highest axis */
static int get_symmetry(gamessdata *data) {
  /* TODO!! */
  /* This could be lumped together with 
   * get_guess_options() where we deal with orbital
   * symmetry and which comes right after the pointgroup
   * stuff */
  return TRUE;
}


static int read_first_frame(gamessdata *data) {
  /* the angular momentum is populated in get_wavefunction 
   * which is called by get_traj_frame(). We have obtained
   * the array size wavef_size already from the basis set
   * statistics */
  data->angular_momentum = (int*)calloc(3*data->wavef_size, sizeof(int));

  /* Try reading the first frame. 
   * If there is only one frame then also read the
   * final wavefunction. */
  if (!get_traj_frame(data, data->initatoms, data->numatoms)) {
    return FALSE;
  }

  return TRUE;
}

/******************************************************
 *
 * Reads the info printed after the geometry search
 * has finished or whatever analysis was done in a 
 * single point run, e.g. ESP charges, Hessian, etc.
 * Rewinds to the beginning of the search when done,
 * because we read this part at in the initial phase
 * and might have to look for additional timesteps
 * later.
 *
 ******************************************************/
static int get_final_properties(gamessdata *data) {
  long filepos;
  filepos = ftell(data->file);

  /* Go to end of trajectory */
  fseek(data->file, data->end_of_traj, SEEK_SET);


  if (get_esp_charges(data)) {
    printf("gamessplugin) ESP charges found!\n");
  }

  if (data->runtype == GRADIENT) {
    /* get_singlepoint_gradient(data); */
  }


  if (data->runtype == HESSIAN) {
    /* try reading the hessian matrix in internal and
     * cartesian coordinates as well as the internal
     * coordinates together with their associated
     * force constants */
    if (!get_int_coords(data)) {
      printf("gamessplugin) \n");
      printf("gamessplugin) Could not determine the internal \n");
      printf("gamessplugin) coordinate and Hessian info!! \n");
      printf("gamessplugin) \n");
    }
    
    if (!get_cart_hessian(data)) {
      printf("gamessplugin) \n");
      printf("gamessplugin) Could not determine the cartesian \n");
      printf("gamessplugin) Hessian matrix!! \n");
      printf("gamessplugin) \n");
    }

    /* read the wavenumbers, intensities of the normal modes 
     * as well as the modes themselves */
    if (!get_normal_modes(data)) {
      printf("gamessplugin) \n");
      printf("gamessplugin) Could not scan the normal modes,\n");
      printf("gamessplugin) \n");
    }
  }

  /* Read localized orbitals if there are any */
  read_localized_orbitals(data);


  fseek(data->file, filepos, SEEK_SET);
  return TRUE; 
}


static int read_localized_orbitals(gamessdata *data) {
  int i;
  qm_timestep_t *ts;
  qm_wavefunction_t *wavef;

  /* Move past the listing of the canonical MOs */
  goto_keyline(data->file, "ENERGY COMPONENTS", NULL);
  eatline(data->file, 1);

  ts = data->qm_timestep + data->num_frames-1;

  for (i=0; i<2; i++) {
    wavef = add_wavefunction(ts);

    if (get_wavefunction(data, ts, wavef) == FALSE ||
        (wavef->type!=MOLFILE_WAVE_BOYS &&
         wavef->type!=MOLFILE_WAVE_PIPEK &&
         wavef->type!=MOLFILE_WAVE_RUEDEN)) {
      del_wavefunction(ts);
      return FALSE;
    }
    else {
      char typestr[64];
      if (wavef->spin==SPIN_ALPHA) {
        strcpy(typestr, "alpha");
      }
      else if (wavef->spin==SPIN_BETA) {
        strcpy(typestr, "beta");
      }
      printf("gamessplugin) Localized orbitals (%s) found for timestep %d\n",
             typestr, data->num_frames-1);
    }
  }

  return TRUE;
}


static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
                           int *numatoms) {
  int i = 0;
  int growarray = 0;

  if (*numatoms<0) {
    *atoms = (qm_atom_t*)calloc(1, sizeof(qm_atom_t));
    growarray = 1;
  }

  /* Read in the coordinates until an empty line is reached.
   * We expect 5 entries per line */
  while (1) {
    char buffer[BUFSIZ];
    char atname[BUFSIZ];
    float atomicnum;
    float x,y,z;
    int n;
    qm_atom_t *atm;

    GET_LINE(buffer, file);
    n = sscanf(buffer,"%s %f %f %f %f",atname,&atomicnum,&x,&y,&z);

    if (n!=5) break;

    if (growarray && i>0) {
      *atoms = (qm_atom_t*)realloc(*atoms, (i+1)*sizeof(qm_atom_t));
    }
    atm = (*atoms)+i;

    strncpy(atm->type, atname, sizeof(atm->type));
    atm->atomicnum = floor(atomicnum+0.5); /* nuclear charge */
    
    /* if coordinates are in Bohr convert them to Angstrom */
    if (unit==BOHR) {
      x *= BOHR_TO_ANGS;
      y *= BOHR_TO_ANGS;
      z *= BOHR_TO_ANGS;
    }
    
    atm->x = x;
    atm->y = y;
    atm->z = z; 
    i++;
  }

  /* If file is broken off in the middle of the coordinate block 
   * we cannot use this frame. */
  if (*numatoms>=0 && *numatoms!=i) return FALSE;

  (*numatoms) = i;
  return TRUE;
}


/********************************************************
 *
 * Read basis set and orbital statistics such as
 * # of shells, # of A/B orbitals, # of electrons, 
 * multiplicity and total charge
 *
 ********************************************************/
static int get_basis_stats(gamessdata *data) {

  char buffer[BUFSIZ];
  char word[7][BUFSIZ];
  int i;

  buffer[0] = '\0';
  for (i=0; i<7; i++) word[i][0] = '\0';

  /* look for the orbital/charge/... info section */
  pass_keyline(data->file, "TOTAL NUMBER OF BASIS", NULL);


  /* go ahead reading the info */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%*s %*s %*s %*s %*s %*s %*s %d",
         &(data->wavef_size));

  /* read the number of electrons */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%*s %*s %*s %*s %d",
         &(data->num_electrons));

  /* read the charge of the molecule */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%*s %*s %*s %*s %d",
         &(data->totalcharge));

  /* read the multiplicity of the molecule */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%*s %*s %*s %d",
         &(data->multiplicity));

  /* read number of A orbitals */
  GET_LINE(buffer, data->file);

  /* note the different number of items per line for A/B orbitals
   * due to "(ALPHA)" and "(BETA )" !! */
  sscanf(buffer,"%*s %*s %*s %*s %*s %*s %d",
         &(data->num_occupied_A)); 
    
  /* read number of B orbitals */
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%*s %*s %*s %*s %*s %*s %*s %d",
         &(data->num_occupied_B)); 


  printf("gamessplugin) Number of Electrons: %d \n",
      data->num_electrons);

  printf("gamessplugin) Charge of Molecule : %d \n",
      data->totalcharge);

  printf("gamessplugin) Multiplicity of Wavefunction: %d \n",
      data->multiplicity);

  printf("gamessplugin) Number of occupied A / B orbitals: %d / %d \n",\
      data->num_occupied_A, data->num_occupied_B);

  printf("gamessplugin) Number of gaussian basis functions: %d \n",\
      data->wavef_size);

 
  return TRUE;
}



/*******************************************************
 *
 * Reads in the $GUESS options.
 *
 * ******************************************************/
static int get_guess_options(gamessdata *data)
{
  char word[2][BUFSIZ];
  char buffer[BUFSIZ];
  int i = 0;
  long filepos;
  filepos = ftell(data->file);

  /* initialize buffers */
  buffer[0] = '\0';
  for (i=0; i<2; i++) word[i][0] = '\0';

  /* parse for GUESS field */
  if (!pass_keyline(data->file, "GUESS OPTIONS", NULL)) {
    printf("gamessplugin) No GUESS OPTIONS found!\n");
    fseek(data->file, filepos, SEEK_SET);
    return FALSE;
  }

  /* next line contains all we need */
  eatline(data->file, 1);
  GET_LINE(buffer, data->file);
  sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);

  /* store it */
  strncpy(data->guess,(&word[1][0])+1,sizeof(data->guess));

  printf("gamessplugin) Run was performed with GUESS = %s \n",
        data->guess);

 /* Since this block occurs in the middle of first frame
  * we need to rewind. */
  fseek(data->file, filepos, SEEK_SET);

  return TRUE;
}



/*******************************************************
 *
 * this function reads in the basis set data 
 *
 * ******************************************************/
int get_basis(gamessdata *data) {

  char buffer[BUFSIZ];
  char word[4][BUFSIZ];
  int i = 0; 
  int success = 0;
  int numread, numshells;
  shell_t *shell;
  long filepos;

  if (!strcmp(data->gbasis, "MNDO") ||
      !strcmp(data->gbasis, "AM1")  ||
      !strcmp(data->gbasis, "PM3")) {
    /* Semiempirical methods are based on STOs.
     * The only parameter we need for orbital rendering
     * are the exponents zeta for S, P, D,... shells for
     * each atom. Since GAMESS doesn't print these values
     * we skip reading the basis set but and hardcode the
     * parameters in tables in VMD. */
    return MOLFILE_ERROR;
  }

  /* initialize buffers */
  buffer[0] = '\0';
  for (i=0; i<3; i++) word[i][0] = '\0';
  
  /* Search for "ATOMIC BASIS SET" line */
  /* XXX what is a good stopstring? */
  if (!pass_keyline(data->file, "ATOMIC BASIS SET", NULL))
    printf("gamessplugin) No basis set found!\n");


  /* skip the next 6 lines */
  for(i=0; i<5; i++) { eatline(data->file, 1); }

  /* Allocate space for the basis for all atoms */
  /* When the molecule is symmetric the actual number atoms with
   * a basis set could be smaller */
  data->basis_set = (basis_atom_t*)calloc(data->numatoms, sizeof(basis_atom_t));


  i = 0; /* basis atom counter */

  do {
    prim_t *prim = NULL;
    char shellsymm;
    int numprim = 0;
    int icoeff = 0;
    filepos = ftell(data->file);
    GET_LINE(buffer, data->file);
      
    /* Count the number of relevant words in the line. */
    numread = sscanf(buffer,"%s %s %s %s",&word[0][0], &word[1][0],
           &word[2][0], &word[3][0]);

    switch (numread) {
      case 1:
        /* Next atom */
        strcpy(data->basis_set[i].name, &word[0][0]);

        /* skip initial blank line */
        eatline(data->file, 1);

        /* read the basis set for the current atom */
        shell = (shell_t*)calloc(1, sizeof(shell_t)); 
        numshells = 0;

        do {
          filepos = ftell(data->file);
          numprim = read_shell_primitives(data, &prim, &shellsymm, icoeff);

          if (numprim>0) {
            /* make sure we have eiter S, L, P, D, F or G shells */
            if ( (shellsymm!='S' && shellsymm!='L' && shellsymm!='P' && 
                  shellsymm!='D' && shellsymm!='F' && shellsymm!='G') ) {
              printf("gamessplugin) WARNING ... %c shells are not supported \n", shellsymm);
            }
            
            /* create new shell */
            if (numshells) {
              shell = (shell_t*)realloc(shell, (numshells+1)*sizeof(shell_t));
            }
            shell[numshells].numprims = numprim;
            shell[numshells].symmetry = shellsymm_int(shellsymm);
            shell[numshells].prim = prim;
            data->num_basis_funcs += numprim;

            /* We split L-shells into one S and one P-shell.
             * I.e. for L-shells we have to go back read the shell again
             * this time using the second contraction coefficients.
             * We use L and M instead of S and P for the shell symmetry
             * in order to be able to distinguish SP-type shells. */
            if (shellsymm=='L' && !icoeff) {
              fseek(data->file, filepos, SEEK_SET);
              icoeff++;
            } else if (shellsymm=='L' && icoeff) {
              shell[numshells].symmetry = SP_P_SHELL;
              icoeff = 0;
            }

            numshells++;
          }
        } while (numprim);

        /* store shells in atom */
        data->basis_set[i].numshells = numshells;
        data->basis_set[i].shell = shell;

        /* store the total number of basis functions */
        data->num_shells += numshells;
        i++;

        /* go back one line so that we can read the name of the
         * next atom */
        fseek(data->file, filepos, SEEK_SET);

        break;

      case 4:
        /* this is the very end of the basis set */
        if (!strcmp(&word[0][0],"TOTAL")  &&
            !strcmp(&word[1][0],"NUMBER") && 
            !strcmp(&word[2][0],"OF")     &&
            !strcmp(&word[3][0],"BASIS")) {
          success = 1;
          /* go back one line so that get_basis_stats()
             can use this line as a keystring. */
          fseek(data->file, filepos, SEEK_SET);
        }
        break;
    }

  } while (!success);


  printf("gamessplugin) Parsed %d uncontracted basis functions for %d atoms.\n",
         data->num_basis_funcs, i);

  data->num_basis_atoms = i;


  /* allocate and populate flat arrays needed for molfileplugin */
  return fill_basis_arrays(data);
}


/**************************************************
 *
 * Convert shell symmetry type from char to int.
 *
 ************************************************ */
static int shellsymm_int(char symm) {
  int shell_symmetry;

  switch (symm) {
    case 'L':
      shell_symmetry = SP_S_SHELL;
      break;
    case 'M':
      shell_symmetry = SP_P_SHELL;
      break;
    case 'S':
      shell_symmetry = S_SHELL;
      break;
    case 'P':
      shell_symmetry = P_SHELL;
      break;
    case 'D':
      shell_symmetry = D_SHELL;
      break;
    case 'F':
      shell_symmetry = F_SHELL;
      break;
    case 'G':
      shell_symmetry = G_SHELL;
      break;
    default:
      shell_symmetry = UNK_SHELL;
      break;
  }

  return shell_symmetry;
}



/******************************************************
 *
 * Populate the flat arrays containing the basis
 * set data.
 *
 ******************************************************/
static int fill_basis_arrays(gamessdata *data) {
  int i, j, k;
  int shellcount = 0;
  int primcount = 0;
  int found = 0;

  float *basis;
  int *num_shells_per_atom;
  int *num_prim_per_shell;
  int *shell_symmetry;
  int *atomicnum_per_basisatom;

  /* Count the total number of primitives which
   * determines the size of the basis array. */
  for(i=0; i<data->num_basis_atoms; i++) {
    for (j=0; j<data->basis_set[i].numshells; j++) {
      primcount += data->basis_set[i].shell[j].numprims;
    }
  }

  /* reserve space for pointer to array containing basis
   * info, i.e. contraction coeficients and expansion 
   * coefficients; need 2 entries per basis function, i.e.
   * exponent and contraction coefficient; also,
   * allocate space for the array holding the orbital symmetry
   * information per primitive Gaussian.
   * Finally, initialize the arrays holding the number of 
   * shells per atom and the number of primitives per shell*/
  basis = (float *)calloc(2*primcount,sizeof(float));

  /* make sure memory was allocated properly */
  if (basis == NULL) {
    PRINTERR;
    return MOLFILE_ERROR;
  }

  shell_symmetry = (int *)calloc(data->num_shells, sizeof(int));
  
  /* make sure memory was allocated properly */
  if (shell_symmetry == NULL) {
    PRINTERR; 
    return MOLFILE_ERROR;
  }

  num_shells_per_atom = (int *)calloc(data->num_basis_atoms, sizeof(int));

  /* make sure memory was allocated properly */
  if (num_shells_per_atom == NULL) {
    PRINTERR; 
    return MOLFILE_ERROR;
  }

  num_prim_per_shell = (int *)calloc(data->num_shells, sizeof(int));

  /* make sure memory was allocated properly */
  if (num_prim_per_shell == NULL) {
    PRINTERR;
    return MOLFILE_ERROR;
  }

  atomicnum_per_basisatom = (int *)calloc(data->num_basis_atoms, sizeof(int));

  /* make sure memory was allocated properly */
  if (atomicnum_per_basisatom == NULL) {
    PRINTERR;
    return MOLFILE_ERROR;
  }


  /* store pointers in struct gamessdata */
  data->basis = basis;
  data->shell_symmetry = shell_symmetry;
  data->num_shells_per_atom = num_shells_per_atom;
  data->num_prim_per_shell = num_prim_per_shell;
  data->atomicnum_per_basisatom = atomicnum_per_basisatom;

  primcount = 0;
  for (i=0; i<data->num_basis_atoms; i++) {
    /* assign atomic number */
    for(j=0; j<data->numatoms; j++) {
      if (!strcmp(data->initatoms[j].type, data->basis_set[i].name)) {
        found = 1;
        break;
      }
    }
    if (!found) {
      printf("gamessplugin) ERROR: Couldn't find atomic number for basis set atom %s",
             data->initatoms[j].type);
      return MOLFILE_ERROR;
    }
    data->basis_set[i].atomicnum = data->initatoms[j].atomicnum;
    atomicnum_per_basisatom[i] = data->initatoms[j].atomicnum;

    num_shells_per_atom[i] = data->basis_set[i].numshells;

    for (j=0; j<data->basis_set[i].numshells; j++) {
      shell_symmetry[shellcount] = data->basis_set[i].shell[j].symmetry;
      num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;

      for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
        basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
        basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
        primcount++;
      }
      shellcount++;
    }
  } 

  return TRUE;
}


/******************************************************
 *
 * read all primitives for the current shell
 *
 ******************************************************/
static int read_shell_primitives(gamessdata *data, prim_t **prim, char *shellsymm,
                                 int icoeff) {
  char buffer[BUFSIZ];
  float exponent = 0.0; 
  float contract[2] = {0.0, 0.0};
  int shell, success;
  int primcounter = 0;
  (*prim) = (prim_t*)calloc(1, sizeof(prim_t));

  do {
    GET_LINE(buffer, data->file);
      success = sscanf(buffer,"%i %c %*s %f %f %f", &shell,
                       shellsymm,
                       &exponent, &contract[0], &contract[1]); 

    /* store in basis array and increase the counter */ 
    switch (success) {
      case 4:
        if (primcounter) {
          *prim = (prim_t*)realloc(*prim, (primcounter+1)*sizeof(prim_t));
        }

        /* store exponent */
        (*prim)[primcounter].exponent = exponent;
          
        /* store coefficient */
        (*prim)[primcounter].contraction_coeff = contract[0];
        
        primcounter++;
        break;

      case 5:
        if (primcounter) {
          *prim = (prim_t*)realloc(*prim, (primcounter+1)*sizeof(prim_t));
        }

        /* store exponent */
        (*prim)[primcounter].exponent = exponent;
          
        /* store coefficient */
        (*prim)[primcounter].contraction_coeff = contract[icoeff];
        
        primcounter++;
        break;

      case -1:
        /* otherwise it's an empty line which represents the end of the shell */
        break;

      case 1:
        /* the user had given the next atom a numeric name */
        success = -1;
        break;
    }

  } while(success>0);

  if (!primcounter) free(*prim);

  return primcounter;
}



/******************************************************
 *
 * this function extracts the trajectory information
 * from the output file
 *
 * *****************************************************/
static int get_traj_frame(gamessdata *data, qm_atom_t *atoms,
                          int natoms) {
  qm_timestep_t *cur_ts;
  qm_wavefunction_t *wavef_alpha;
  char buffer[BUFSIZ];
  char word[BUFSIZ];
  buffer[0] = '\0';
  word[0]   = '\0';

  printf("gamessplugin) Timestep %i:\n", data->num_frames_read);
  printf("gamessplugin) ============\n");

  fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);


  if (data->runtype==OPTIMIZE || data->runtype==SADPOINT) {

    pass_keyline(data->file, "COORDINATES OF ALL ATOMS", NULL);
    eatline(data->file, 2);

    if (!get_coordinates(data->file, &data->initatoms, ANGSTROM, &natoms)) {
      printf("gamessplugin) Couldn't find coordinates for timestep %i\n", data->num_frames_read);
    }

  }
  else if (data->runtype==SURFACE) {
    if (!pass_keyline(data->file, "---- SURFACE MAPPING GEOMETRY", NULL)) {
      return FALSE;
    }
  }

 
  /* get a convenient pointer to the current qm timestep */
  cur_ts = data->qm_timestep + data->num_frames_read;

  /* read the SCF energies */
  if (get_scfdata(data, cur_ts) == FALSE) {
    printf("gamessplugin) Couldn't find SCF iterations for timestep %i\n", data->num_frames_read);
  }

  /* Allocate memory for new wavefunction */
  wavef_alpha = add_wavefunction(cur_ts);

  /* Try to read wavefunction and orbital energies */
  if (get_wavefunction(data, cur_ts, wavef_alpha) == FALSE) {
    /* Free the last wavefunction again. */
    del_wavefunction(cur_ts);
    printf("gamessplugin) No canonical wavefunction present for timestep %i\n", data->num_frames_read);

  } else {
    char typestr[64];
    strcpy(typestr, wavef_alpha->info);
    if (data->scftype==UHF) {
      strcat(typestr, ", alpha");
    }
    printf("gamessplugin) Wavefunction (%s) found for timestep %i\n", typestr, data->num_frames_read);

    if (data->scftype==UHF || data->scftype==GVB || data->scftype==MCSCF) {
      /* Try to read second wavefunction
       * (spin beta or GI orbitals or MCSCF optimized orbs) */
      qm_wavefunction_t *wavef_beta;
      wavef_beta = add_wavefunction(cur_ts);

      if (get_wavefunction(data, cur_ts, wavef_beta) == FALSE) {
        /* Free the last wavefunction again. */
        del_wavefunction(cur_ts);

        printf("gamessplugin) No beta wavefunction present for timestep %i\n", data->num_frames_read);
      } else {
        strcpy(typestr, wavef_beta->info);
        if (data->scftype==UHF) {
          strcat(typestr, ", beta");
        }
        printf("gamessplugin) Wavefunction (%s)  found for timestep %i\n", typestr, data->num_frames_read);
      }
    }
  }

  if (get_population(data, cur_ts)) {
    printf("gamessplugin) Mulliken charges found\n");
  }

  if (data->citype!=NONE) {
    if (pass_keyline(data->file, "CI DENSITY MATRIX AND NATURAL ORBITALS",
                       "GRADIENT (HARTREE/BOHR)")) {
      int i, numstates=0;
      qm_wavefunction_t *wave_ci;
      goto_keyline(data->file, "NUMBER OF STATES", NULL);
      GET_LINE(buffer, data->file);
      trimleft(buffer);
      sscanf(buffer, " NUMBER OF STATES = %d", &numstates);
      printf("gamessplugin) Number of CI states = %d\n", numstates);

      for (i=0; i<numstates; i++) {
        float cienergy = 0.f;
        goto_keyline(data->file, "CI EIGENSTATE", NULL);
        GET_LINE(buffer, data->file);
        sscanf(buffer,"%*s %*s %*d %*s %*s %*s %f", &cienergy);
        printf("gamessplugin) CI energy[%d] = %f\n", i, cienergy);

        wave_ci = add_wavefunction(cur_ts);

        if (get_wavefunction(data, cur_ts, wave_ci) == FALSE) {
          del_wavefunction(cur_ts);      
        }
        else {
          int j, canon =-1;
          wave_ci->exci = i;
          printf("gamessplugin) Found CI natural orbitals for timestep %d\n", data->num_frames_read);

          for (j=0; j<cur_ts->numwave; j++) {
            if (cur_ts->wave[j].type==MOLFILE_WAVE_CANON &&
                cur_ts->wave[j].spin==wave_ci->spin) {
              canon = j;
              break;
            }
          }
          if (canon>=0) {
            /* Replace existing canonical wavefunction for this step */
            replace_wavefunction(cur_ts, canon);
          }
        }
      }
    }
  }


  if (get_gradient(data, cur_ts) == FALSE) {
    printf("gamessplugin) No energy gradient present for timestep %i\n", data->num_frames_read);
  } else {
    printf("gamessplugin) Energy gradient found for timestep %i\n", data->num_frames_read);
  }

  /* If this is the last frame of the trajectory and the file
   * wasn't truncated and the program didn't terminate
   * abnormally then read the final wavefunction. */
  if ((data->runtype == OPTIMIZE || data->runtype == SADPOINT) &&
      (data->num_frames_read+1 == data->num_frames &&
       (data->opt_status == MOLFILE_QM_STATUS_UNKNOWN || 
        data->opt_status == MOLFILE_QM_OPT_CONVERGED ||
        data->opt_status == MOLFILE_QM_OPT_NOT_CONV))) {
    qm_wavefunction_t *wave_final;

    /* We need to jump over the end of the trajectory because 
     * this is also the keystring for get_wavefunction() to
     * bail out. */
    if (data->opt_status == MOLFILE_QM_OPT_CONVERGED || 
        data->opt_status == MOLFILE_QM_OPT_NOT_CONV) {
      fseek(data->file, data->end_of_traj, SEEK_SET);
    }

    /* Try to read final wavefunction and orbital energies
     * Any wavefunction previously stored for the final
     * timestep will be overwritten. */
    wave_final = add_wavefunction(cur_ts);

    if (get_wavefunction(data, cur_ts, wave_final) == FALSE) {
      printf("gamessplugin) No final wavefunction present for timestep %d\n", data->num_frames_read);
      /* Free the the pointer an the contents of the last wfn */
      del_wavefunction(cur_ts);

    } else {
      char typestr[MOLFILE_BUFSIZ];

      /* if there exists a canonical wavefunction of the same spin
       * we'll replace it */
      if (cur_ts->numwave>1 && wave_final->type==MOLFILE_WAVE_CANON) {
        int i, found =-1;
        for (i=0; i<cur_ts->numwave; i++) {
          if (cur_ts->wave[i].type==wave_final->type &&
              cur_ts->wave[i].spin==wave_final->spin &&
              cur_ts->wave[i].exci==wave_final->exci &&
              !strncmp(cur_ts->wave[i].info, wave_final->info, MOLFILE_BUFSIZ)) {
            found = i;
            break;
          }
        }
        if (found>=0) {
          /* If the new wavefunction has more orbitals we 
           * replace the old one for this step. */
          if (wave_final->num_orbitals > 
              cur_ts->wave[found].num_orbitals) {
            /* Replace existing wavefunction for this step */
            replace_wavefunction(cur_ts, found);
          } else {
            /* Delete last wavefunction again */
            del_wavefunction(cur_ts);
          }
          wave_final = &cur_ts->wave[cur_ts->numwave-1];
        }
      }


      strcpy(typestr, wave_final->info);
      if (data->scftype==UHF) {
        strcat(typestr, ", alpha");
      }

      printf("gamessplugin) Final wavefunction (%s) found for timestep %d\n",
             typestr, data->num_frames_read);
    }

  }

  data->num_frames_read++;

  return TRUE;
}


/* Analyze the trajectory.
 * Read the parameters comtrolling geometry search and
 * find the end of the trajectory, couinting the frames
 * on the way. Store the filepointer for the beginning of
 * each frame in *filepos_array. */
static int analyze_traj(gamessdata *data) {
  char buffer[BUFSIZ], nserch[BUFSIZ];
  char *line;
  long filepos;
  filepos = ftell(data->file);

  data->filepos_array = (long* )calloc(1, sizeof(long ));


  if (data->runtype==OPTIMIZE || data->runtype==SADPOINT) {
    pass_keyline(data->file,
                   "PARAMETERS CONTROLLING GEOMETRY SEARCH", NULL);
    eatline(data->file, 2);

    GET_LINE(buffer, data->file);
    sscanf(buffer, "NSTEP  = %i", &data->max_opt_steps);
    eatline(data->file, 3);
    GET_LINE(buffer, data->file);
    sscanf(buffer, "OPTTOL = %f", &data->opt_tol);

    /* The $STATP options are followed by the coordinates 
     * but we can skip them here because we rewind after
     * get_guess_options() and try to read them in
     * get_traj_frame(). */
  }
  else if (data->runtype==SURFACE) {
    if (pass_keyline(data->file,
                        "POTENTIAL SURFACE MAP INPUT", NULL)) {
      
      int coord1[2];
      int mplevel1=-1, nstep1;
      float origin1, disp1;
      char runtype1[BUFSIZ];
      char scftype1[BUFSIZ];
      char dfttyp1[BUFSIZ];
      char *tmp;
        
      eatline(data->file, 1);

      GET_LINE(buffer, data->file);
      /* int n=sscanf(buffer, "JOB 1 IS RUNTYP=%s SCFTYP=%s CITYP=%*s",
         runtype1, scftype1); */
      tmp = strstr(buffer, "RUNTYP=") + 7;
      sscanf(tmp, "%s", runtype1);
      tmp = strstr(buffer, "SCFTYP=") + 7;
      sscanf(tmp, "%s", scftype1);
      printf("gamessplugin) JOB 1 IS RUNTYP=%s SCFTYP=%s\n", runtype1, scftype1);

      GET_LINE(buffer, data->file);
      sscanf(buffer, "MPLEVL= %i", &mplevel1);
      tmp = strstr(buffer, "DFTTYP=") + 7;
      sscanf(tmp, "%s", dfttyp1);
      printf("gamessplugin) MPLEVL= %i DFTTYP=%s\n", mplevel1, dfttyp1);

      GET_LINE(buffer, data->file);
      sscanf(buffer, "COORD 1 LYING ALONG ATOM PAIR %i %i",
             coord1, coord1+1);
      GET_LINE(buffer, data->file);
      tmp = strstr(buffer, "ORIGIN=") + 7;
      sscanf(tmp, "%f", &origin1);
      tmp = strstr(buffer, "DISPLACEMENT=") + 13;
      sscanf(tmp, "%f", &disp1);
      tmp = strstr(buffer, "AND") + 3;
      sscanf(tmp, "%i STEPS.", &nstep1);
      printf("gamessplugin) origin=%f, displacement=%f nstep=%i\n", origin1, disp1, nstep1);
    }
  }
  else {
    /* We have just one frame */
    data->num_frames = 1;
    pass_keyline(data->file, "1 ELECTRON INTEGRALS",
                 "ENERGY COMPONENTS");
    data->filepos_array[0] = ftell(data->file);

    /* Check wether SCF has converged */
    if (pass_keyline(data->file, "SCF IS UNCONVERGED, TOO MANY ITERATIONS",
                     "ENERGY COMPONENTS")==FOUND) {
      printf("gamessplugin) SCF IS UNCONVERGED, TOO MANY ITERATIONS\n");
      data->opt_status = MOLFILE_QM_SCF_NOT_CONV;
    } else {
      data->opt_status = MOLFILE_QM_OPT_CONVERGED;
      fseek(data->file, data->filepos_array[0], SEEK_SET);
    }

    pass_keyline(data->file, "ENERGY COMPONENTS", NULL);
    data->end_of_traj = ftell(data->file);

    /* Allocate memory for the frame */
    data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
    memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
    
    return TRUE;
  }

  printf("gamessplugin) Analyzing trajectory...\n");
  data->opt_status = MOLFILE_QM_STATUS_UNKNOWN;

  while (1) {
    if (!fgets(buffer, sizeof(buffer), data->file)) break;
    line = trimleft(buffer);

    /* at this point we have to distinguish between
     * pre="27 JUN 2005 (R2)" and "27 JUN 2005 (R2)"
     * versions since the output format for geometry
     * optimizations has changed */
    if (data->version==1) {
      strcpy(nserch, "1NSERCH=");
    }
    else if (data->version==2) {
      strcpy(nserch, "BEGINNING GEOMETRY SEARCH POINT NSERCH=");
    }

    if (strstr(line, nserch) ||
        strstr(line, "---- SURFACE MAPPING GEOMETRY")) {
      printf("gamessplugin) %s", line);

      if (data->num_frames > 0) {
        data->filepos_array = (long*)realloc(data->filepos_array,
                                (data->num_frames+1)*sizeof(long));
      }
      data->filepos_array[data->num_frames] = ftell(data->file);

      /* Make sure that we have at least a complete coordinate
         block in order to consider this a new frame. */
      if (pass_keyline(data->file, "COORDINATES OF",
               "BEGINNING GEOMETRY SEARCH POINT NSERCH=") == FOUND) {
        if (pass_keyline(data->file, "INTERNUCLEAR DISTANCES",
                           "1 ELECTRON INTEGRALS") ||
            pass_keyline(data->file, "1 ELECTRON INTEGRALS",
               "BEGINNING GEOMETRY SEARCH POINT NSERCH=")) {
          data->num_frames++;
        }
      }
    }
    else if (strstr(line, "***** EQUILIBRIUM GEOMETRY LOCATED") ||
             strstr(line, "... DONE WITH POTENTIAL SURFACE SCAN")) {
      printf("gamessplugin) ==== End of trajectory. ====\n");
      data->opt_status = MOLFILE_QM_OPT_CONVERGED;
      break;
    }
    else if (strstr(line, "***** FAILURE TO LOCATE STATIONARY POINT,")) {
      printf("gamessplugin) %s\n", line);
      if (strstr(strchr(line, ','), "SCF HAS NOT CONVERGED")) {
        data->opt_status = MOLFILE_QM_SCF_NOT_CONV;
        break;
      }
      else if (strstr(strchr(line, ','), "TOO MANY STEPS TAKEN")) {
        data->opt_status = MOLFILE_QM_OPT_NOT_CONV;
        break;
      }
    }
  }

  data->end_of_traj = ftell(data->file);
  fseek(data->file, filepos, SEEK_SET);

  if (data->opt_status == MOLFILE_QM_STATUS_UNKNOWN) {
    /* We didn't find any of the regular key strings,
     * the run was most likely broken off and we have an
     * incomplete file. */
    data->opt_status = MOLFILE_QM_FILE_TRUNCATED;
  }


  /* Allocate memory for all frames */
  data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
                                              sizeof(qm_timestep_t));
  memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));


  if (data->opt_status == MOLFILE_QM_SCF_NOT_CONV ||
      data->opt_status == MOLFILE_QM_FILE_TRUNCATED) {
    return FALSE;  
  }

  return TRUE;
}


/***************************************************************
 *
 * Read the number of scf iterations and the scf energies
 * for the current timestep. 
 * Assumes that the file pointer is somewhere before this:
 * ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE    DIIS ERROR
 * 1  0  0      -39.7266993475   -39.7266993475   0.000000118   0.000000000
 * 2  1  0      -39.7266991566     0.0000001909   0.000000032   0.000000000
 * ...
 * then it reads the block up to the next blank line.
 * The second argument is a pointer to the qm timestep you want to
 * store the data in. Memory for the scfenergies will be
 * allocated.
 *
 ***************************************************************/
static int get_scfdata(gamessdata *data, qm_timestep_t *ts) {
  char buffer[BUFSIZ];
  char word[3][BUFSIZ];
  long filepos;
  int i;
  int numread, numiter=0, dum;
  char *line;

  buffer[0] = '\0';
  for (i=0; i<3; ++i) word[i][0] = '\0';

  /* look for SCF iteration energies */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s", &word[0][0], &word[1][0],
         &word[2][0]);
  } while (!(!strcmp(&word[0][0],"ITER") &&
             (!strcmp(&word[1][0],"EX") || !strcmp(&word[1][0],"TOTAL"))));

  
  /* store current file position since we first have to count iterations */
  filepos = ftell(data->file);

  /* read until the next blank line count the iterations */
  do {
    GET_LINE(buffer, data->file);
    line = trimleft(buffer);
    numread = sscanf(line,"%i %i %*i %*f", &dum, &dum);
    if (numread==2) numiter++;
  } while (strlen(line));

  printf("gamessplugin) %i SCF iterations\n", numiter);

  /* go back and read energies */
  fseek(data->file, filepos, SEEK_SET);
  

  /* allocate memory for scfenergy array */
  ts->scfenergies = (double *)calloc(numiter,sizeof(double));
  
  i=0;
  do {
    GET_LINE(buffer, data->file);
    line = trimleft(buffer);
    numread = sscanf(line,"%i %i %*i %*f", &dum, &dum);
    if (numread==2) {
      sscanf(buffer,"%*i %*i %*i %lf", ts->scfenergies+i);
      i++;
    }
  } while (strlen(line));

#if 0
  for (i=0; i<numiter; i++) {
    printf("scfenergies[%i] = %f\n", i, ts->scfenergies[i]);
  }
#endif

  ts->num_scfiter = numiter;
  
  return TRUE;
}


/*********************************************************
 *
 * this function reads the actual wavefunction, which is
 * punched at the end of the log file
 *
 **********************************************************/
static int get_wavefunction(gamessdata *data, qm_timestep_t *ts, qm_wavefunction_t *wf)
{
  float *orb_energies, orben[5];
  float *wave_coeff;
  char buffer[BUFSIZ];
  char word[6][BUFSIZ];
  int num_orbitals = 0;
  int i = 0, num_values = 0;
  long filepos;
  char *line;
  int have_orben = 0;
  int n[5];

  buffer[0] = '\0';
  for (i=0; i<6; i++) word[i][0] = '\0';

  if (wf == NULL) {
    PRINTERR;         
    return FALSE;
  }

  wf->type = MOLFILE_WAVE_UNKNOWN;
  wf->spin = SPIN_ALPHA;
  wf->exci = 0;
  strncpy(wf->info, "unknown", MOLFILE_BUFSIZ);

  /*
   * Scan for something like this:

          ------------------
          MOLECULAR ORBITALS     <<--- sometimes EIGENVECTORS
          ------------------

                      1          2          3          4          5
                  -11.0297    -0.9121    -0.5205    -0.5205    -0.5205  <<-- orbital energies
                     A          A          A          A          A   
    1  C  1  S    0.991925   0.221431   0.000006  -0.000001   0.000002
    2  C  1  S    0.038356  -0.627585  -0.000021   0.000003  -0.000006
    3  C  1  X    0.000000  -0.000004   0.338169  -0.030481  -0.460283
     ...

                     6          7          8          9
                    0.7192     0.7192     0.7193     0.7611
                     A          A          A          A   
    1  C  1  S    0.000028   0.000012   0.000092   0.252320
    2  C  1  S   -0.000183  -0.000077  -0.000594  -1.632834
    3  C  1  X   -0.890147   0.062618   0.654017  -0.000154
      ...


     ----------------------------------------------------------------
     PROPERTY VALUES FOR THE RHF   SELF-CONSISTENT FIELD WAVEFUNCTION
     ----------------------------------------------------------------
  * 
  */

  /* Remember position in order to go back if no wave function was found */
  filepos = ftell(data->file);

  do {
    GET_LINE(buffer, data->file);

    line = trimleft(buffer);
    if      (strstr(line, "----- ALPHA SET -----")) {
      wf->type = MOLFILE_WAVE_CANON;
      strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
      pass_keyline(data->file, "EIGENVECTORS", NULL);
    }
    else if (strstr(line, "----- BETA SET -----")) {
      wf->type = MOLFILE_WAVE_CANON;
      wf->spin = SPIN_BETA;
      strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
      pass_keyline(data->file, "EIGENVECTORS", NULL);
    }
    else if (strstr(line, "***** BETA ORBITAL LOCALIZATION *****")) {
      wf->spin = SPIN_BETA;
    }
    else if (strstr(line, "EIGENVECTORS")==line) {
      wf->type = MOLFILE_WAVE_CANON;
      strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
    }
    else if (strstr(line, "MOLECULAR ORBITALS") &&
             !strstr(line, "LOCALIZED")) {
      wf->type = MOLFILE_WAVE_CANON;
      strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
    }
    else if (strstr(line, "LOCALIZED ORBITALS ARE")) {
      if (strstr(line, "BOYS"))   {
        wf->type = MOLFILE_WAVE_BOYS;
        strncpy(wf->info, "Boys localized", MOLFILE_BUFSIZ);
      }
      else if (strstr(line, "RUEDEN")) {
        wf->type = MOLFILE_WAVE_RUEDEN;
        strncpy(wf->info, "Ruedenberg localized", MOLFILE_BUFSIZ);
      }
      else if (strstr(line, "PIPEK-MEZEY")) {
        wf->type = MOLFILE_WAVE_PIPEK;
        strncpy(wf->info, "Pipek-Mezey localized", MOLFILE_BUFSIZ);
      }
    }
    else if (strstr(line, "GI ORBITALS")) {
      wf->type = MOLFILE_WAVE_GEMINAL;
      strncpy(wf->info, "GVB geminal pairs", MOLFILE_BUFSIZ);
    }
    else if (strstr(line, "MCSCF NATURAL ORBITALS")) {
      wf->type = MOLFILE_WAVE_MCSCFNAT;
      strncpy(wf->info, "MCSCF natural orbitals", MOLFILE_BUFSIZ);
    }
    else if (strstr(line, "MCSCF OPTIMIZED ORBITALS")) {
      wf->type = MOLFILE_WAVE_MCSCFOPT;
      strncpy(wf->info, "MCSCF optimized orbitals", MOLFILE_BUFSIZ);
    }
    else if (strstr(line, "NATURAL ORBITALS IN ATOMIC")==line) {
      wf->type = MOLFILE_WAVE_CINATUR;
      strncpy(wf->info, "CI natural orbitals", MOLFILE_BUFSIZ);
    }
  } while(wf->type==MOLFILE_WAVE_UNKNOWN &&
          !strstr(line, "ENERGY COMPONENTS") &&
          !strstr(line, "***** EQUILIBRIUM GEOMETRY LOCATED") &&
          !strstr(line, "**** THE GEOMETRY SEARCH IS NOT CONVERGED!"));


  /* If we reach the last line of the rhf section without finding 
   * one of the keywords marking the beginning of a wavefunction
   * table then we return.*/
  if (wf->type==MOLFILE_WAVE_UNKNOWN) {
    fseek(data->file, filepos, SEEK_SET);
    return FALSE;
  }

  /* Reserve space for arrays storing wavefunction and orbital
   * energies. we reserve the max. space (num_orbitals==wavef_size)
   * for now and realloc later if, we have less orbitals. */
  wave_coeff = (float *)calloc(data->wavef_size*data->wavef_size,
                                  sizeof(float)); 

  if (wave_coeff == NULL) {
    PRINTERR;         
    return FALSE;
  }
  
  orb_energies = (float *)calloc(data->wavef_size, sizeof(float));

  if (orb_energies == NULL) {
    free(wave_coeff);
    PRINTERR; 
    return FALSE;
  }

  /* store the pointers */
  wf->wave_coeffs  = wave_coeff;
  wf->orb_energies = orb_energies;

  /* skip the next line which here is typically "-------" */
  eatline(data->file, 1);

  while (1) {
    int over=0;
    float coeff[5];

    if (wf->type == MOLFILE_WAVE_GEMINAL) {
      /* Skip over "PAIR x" header line */
      pass_keyline(data->file, "PAIR ", NULL);
    }

    eatwhitelines(data->file);

    filepos = ftell(data->file);
    GET_LINE(buffer, data->file);
    num_values = sscanf(buffer, "%d %d %d %d %d",
                          &n[0], &n[1], &n[2], &n[3], &n[4]);

    /* If we didn't find this line then coefficient table
     * is complete. */
    if (!num_values) {
      fseek(data->file, filepos, SEEK_SET);
      break;
    }

    eatwhitelines(data->file);

    /* Read first line of orbital energies */
    filepos = ftell(data->file);
    GET_LINE(buffer, data->file);
    have_orben = sscanf(buffer,"%f %f %f %f %f", &orben[0],
                        &orben[1], &orben[2], &orben[3], &orben[4]);
    /* Make sure this is not the first line containing coeffs */
    i = sscanf(buffer, " 1 %*s 1 %*s %f %f %f %f %f",
               &coeff[0], &coeff[1], &coeff[2], &coeff[3], &coeff[4]);
    if (i==num_values) have_orben = 0;

    if (have_orben) {
      /* store the orbital energies in the appropriate arrays 
       * read them until we encounter an empty string */
      for(i=0; i<num_values; i++) {
        *(orb_energies+i) = orben[i];
      }
      
      /* increase orbital energy pointer */
      orb_energies = orb_energies+5;
    }      
    else {
      /* No orbital energies present, go back one line */
      fseek(data->file, filepos, SEEK_SET);
    }

    num_orbitals += num_values;

    /* Find first line containing coefficients */
    filepos = ftell(data->file);
    while (fgets(buffer, sizeof(buffer), data->file)) {
      trimleft(buffer);
      if (strstr(line, "ENERGY COMPONENTS") ||
          strstr(line, "---") ||
          strstr(line, "...")) {
        over = 1; break;
      }
      i = sscanf(buffer, " 1 %*s 1 %*s %f %f %f %f %f",
             &coeff[0], &coeff[1], &coeff[2], &coeff[3], &coeff[4]);
      if (i==num_values) break;
      filepos = ftell(data->file);
    }
    fseek(data->file, filepos, SEEK_SET);

    if (over) break;


    /* Read the wave function coefficient block for up to 5
     * orbitals per line. */
    read_coeff_block(data->file, data->wavef_size,
                     wave_coeff, data->angular_momentum);


    /* move wavefunction pointer to start of next five orbitals */
    if (wf->type == MOLFILE_WAVE_GEMINAL) {
      wave_coeff = wave_coeff + 2*data->wavef_size;
    } else {
      wave_coeff = wave_coeff + 5*data->wavef_size;
    }
  }


  if (!num_orbitals) return FALSE;

  /* resize the array to the actual number of read orbitals */
  if (data->wavef_size!=num_orbitals) {
    wf->orb_energies = (float *)realloc(wf->orb_energies,
          num_orbitals*sizeof(float));

    wf->wave_coeffs  = (float *)realloc(wf->wave_coeffs, data->wavef_size*
                              num_orbitals*sizeof(float)); 
  }


  /* store the number of orbitals read in */
  wf->num_orbitals  = num_orbitals;
  wf->have_energies = TRUE;
  wf->have_occup    = FALSE;

  printf("gamessplugin) Number of orbitals scanned: %d \n",\
             wf->num_orbitals);

  return TRUE;
}


/* Read the wave function coefficient block for up to 5
 * orbitals per line:
 *  1  C  1  S    0.989835   0.155361   0.000000  -0.214258   0.000000
 *  2  C  1  S    0.046228  -0.548915   0.000000   0.645267   0.000000
 *  3  C  1  X    0.000000   0.000000   1.030974   0.000000   0.000000
 */
static int read_coeff_block(FILE *file, int wavef_size,
                            float *wave_coeff, int *angular_momentum) {
  int i, j;
  int truncated = 0;
  char buffer[BUFSIZ];
  /* read in the wavefunction */
  /* XXX This will fail for truncated files */
  for (i=0; i<wavef_size; i++) {
    int xexp=0, yexp=0, zexp=0;
    char symm[BUFSIZ];
    float coeff[5];
    int num_values = 0;
    
    GET_LINE(buffer, file);
    
    /* read in the wavefunction coefficients for 5
     * orbitals at a time line by line */
    num_values = sscanf(buffer,"%*5i%*4s%*2i%4s %f %f %f %f %f", 
                        symm, &coeff[0], &coeff[1], &coeff[2],
                        &coeff[3], &coeff[4]);
    
    if (num_values==0) {
      /* The file must have been truncated! */
      truncated = 1;
      break;
    }
    
    for (j=0; j<strlen(symm); j++) {
      switch (symm[j]) {
      case 'X':
        xexp++;
        break;
      case 'Y':
        yexp++;
        break;
      case 'Z':
        zexp++;
        break;
      }
    }
    angular_momentum[3*i  ] = xexp;
    angular_momentum[3*i+1] = yexp;
    angular_momentum[3*i+2] = zexp;
    
    /* each orbital has data->wavef_size entries, 
     * hence we have to use this number as offset when storing 
     * them in groups of five */
    for (j=0 ; j<num_values-1; j++) {
      wave_coeff[j*wavef_size+i] = coeff[j];
    }
  }
  
  if (truncated) return 0;
  
  return 1;
}

/* Read the population analysis section.
 * Currently we parse only the Mulliken charges
 * but we might want to add support for populations
 * and for Lowdin analysis. */
static int get_population(gamessdata *data, qm_timestep_t *ts) {
  int i;
  char buffer[BUFSIZ];
  data->have_mulliken = FALSE;

  if (pass_keyline(data->file,
                     "TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS",
                     "NSERCH=") != FOUND) {
    return FALSE;
  }

  /* Read Mulliken charges if present */
  ts->mulliken_charges = 
    (double *)calloc(data->numatoms, sizeof(double));

  if (!ts->mulliken_charges) {
    PRINTERR; 
    return FALSE;
  }
  
  eatline(data->file, 1);
  
  for (i=0; i<data->numatoms; i++) {
    int n;
    float mullpop, mullcharge, lowpop, lowcharge;
    GET_LINE(buffer, data->file);
    n = sscanf(buffer,"%*i %*s %f %f %f %f",
                   &mullpop, &mullcharge, &lowpop, &lowcharge);
    if (n!=4) return FALSE;
    ts->mulliken_charges[i] = mullcharge;
  }

  if (i!=data->numatoms) return FALSE;

  data->have_mulliken = TRUE;
  return TRUE;
}


/* Read ESP charges.
 * XXX Right now we don't distinguish between different type of
 * ESP-style charges (CHELPG, CONNOLLY, GEODESIC). 
 * This could be solved by reading in the PTSEL keyword in
 * the $PDC group. */
static int get_esp_charges(gamessdata *data) {
  int i;
  char buffer[BUFSIZ];
  data->have_esp = FALSE;

  if (pass_keyline(data->file,
           "ATOM                CHARGE    E.S.D.",
           "...... END OF PROPERTY EVALUATION ") != FOUND) {
    return FALSE;
  }


  /* Read ESP charges if present */
  data->esp_charges = 
    (double *)calloc(data->numatoms, sizeof(double));

  if (data->esp_charges == NULL) {
    PRINTERR; 
    return FALSE;
  }

  eatline(data->file, 1);

  for (i=0; i<data->numatoms; i++) {
    int n;
    double charge;
    GET_LINE(buffer, data->file);
    n = sscanf(buffer,"%*s %lf ", &charge);
    if (n!=2) return FALSE; /* XXX n!=1 ?? */
    data->esp_charges[i] = charge;
  }

  if (i!=data->numatoms) return FALSE;

  data->have_esp = TRUE;
  return TRUE;
}


static int get_gradient(gamessdata *data, qm_timestep_t *ts) {
  char buffer[BUFSIZ];
  float dx, dy, dz;
  long filepos;
  int i=0;
  int numread;

  buffer[0] = '\0';

  /* remember position in order to go back if no forces were found */
  filepos = ftell(data->file);

  /* look for GRADIENT section */
  if (goto_keystring2(data->file, "GRADIENT (HARTREE",
                "***** EQUILIBRIUM GEOMETRY LOCATED", 
                " BEGINNING GEOMETRY SEARCH") != FOUND) {
    fseek(data->file, filepos, SEEK_SET);
    return FALSE;
  }

  eatline(data->file, 3);

  ts->gradient = (float *)calloc(3*data->numatoms, sizeof(float));

  if (ts->gradient == NULL) {
    PRINTERR;         
    fseek(data->file, filepos, SEEK_SET);
    return FALSE;
  }

  /* read the gradient table */
  do {
    GET_LINE(buffer, data->file);
    numread = sscanf(buffer, "%i %*s %*f %f %f %f", &i, &dx, &dy, &dz);
    if (numread==4) {
      ts->gradient[3*(i-1)  ] = dx;
      ts->gradient[3*(i-1)+1] = dy;
      ts->gradient[3*(i-1)+2] = dz;
    }

  } while(numread==4);

  if (i!=data->numatoms) {
    printf("gamessplugin) Number of gradients != number of atoms!\n");
    fseek(data->file, filepos, SEEK_SET);
    return FALSE;
  }

  return TRUE;
}



/***********************************************************
 *
 * this function reads in the wavenumbers and intensities of
 * the normal modes
 *
 * *********************************************************/
static int get_normal_modes(gamessdata *data) {

  char word[4][BUFSIZ];
  char buffer[BUFSIZ];
  int i = 0, k = 0, j = 0;
  int remaining_columns;
  double entry[6]; 
  char separator;
  char *item;
  int counter;


  /* initialize arrays */
  buffer[0] = '\0';
  memset(entry, 0, sizeof(entry));
  for (i=0; i<4; i++) word[i][0] = '\0';

    
  /* reserve memory for dynamically allocated data
   * arrays */
  item = (char *)calloc(BUFSIZ,sizeof(char));

  if (item==NULL){
    PRINTERR;
    return FALSE;
  }


  data->wavenumbers = 
    (float *)calloc(data->numatoms*3,sizeof(float));

  if (data->wavenumbers==NULL) {
    PRINTERR;
    return FALSE;
  }


  data->intensities = 
    (float *)calloc(data->numatoms*3,sizeof(float));

  if (data->intensities==NULL) {
    PRINTERR;
    return FALSE;
  }


  data->normal_modes = 
    (float *)calloc((data->numatoms*3)*(data->numatoms*3),
                 sizeof(float));

  if (data->normal_modes==NULL) {
    PRINTERR;
    return FALSE;
  }


  /* look for FREQUENCYs and IR INTENSITIES */
  for (i=0; i<data->numatoms*3/5; i++) {
    do {
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s ",&word[0][0]);

    } while(strcmp(&word[0][0],"FREQUENCY:"));


    /* scan the frequencies; 
     * this requires some care since there might
     * be imaginary modes present which leads to
     * an additional char per imaginary mode per
     * line */

    /* check for imaginary modes */
    if (strchr(buffer,'I') != NULL) {
      separator = ' ';
      counter = 0;

      /* read all line elements into individual strings */

      /* skip first entry "FREQUENCY" */
      item = strtok(buffer,&separator);
      
      /* start going through the string */
      while ((item = strtok(NULL,&separator)) != NULL) {
        /* check if item is 'I'; if yes, mark previous mode
         * as imaginary; otherwise save the mode */
        if (*item=='I') data->nimag++;

        else {
          /* save only the first 5 modes - there NEVER should
           * be more in any case, but just to make sure
           * we don't overrun the array */
          if (counter<5) {
            *(data->wavenumbers+(i*5)+counter) = atof(item);
            counter++;
          }
        }
      }
    }

    /* no imaginary mode, reading is straightforward */
    else {
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0],
             &entry[0],&entry[1],&entry[2],&entry[3],&entry[4]); 
    
      for (k=0; k<5; k++) {
        *(data->wavenumbers+(i*5)+k) = entry[k]; 
      }
    }

    eatline(data->file, 1);

    /* next line contains the IR INTENSITIES */
    GET_LINE(buffer, data->file);

    /* scan the IR INTENSITIES */
    sscanf(buffer,"%s %s %lf %lf %lf %lf %lf",&word[0][0],&word[1][0],
           &entry[0],&entry[1],&entry[2],&entry[3],&entry[4]);
 
    for (k=0; k<5; k++) {
      *(data->intensities+(i*5)+k) = entry[k]; 
    }

    eatline(data->file, 1);

    /* read the following five modes */
    for (k=0; k<data->numatoms; k++) {
      /* x */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %s %s %lf %lf %lf %lf %lf",&word[0][0], 
             &word[1][0], &word[2][0], &entry[0], &entry[1], &entry[2],
             &entry[3], &entry[4]);

      for (j=0; j<5; j++) {
        *(data->normal_modes+(3*k)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }

      /* y */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0], &entry[0],
             &entry[1],&entry[2], &entry[3],&entry[4]);

      for (j=0; j<5; j++) {
        *(data->normal_modes+(3*k+1)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }

      /* z */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0], &entry[0],
             &entry[1], &entry[2], &entry[3],&entry[4]);

      for (j=0; j<5; j++) {
        *(data->normal_modes+(3*k+2)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }
    }
  }


  /* read the remaining columns */
  if ( (remaining_columns=(data->numatoms*3)%5) != 0 ) { 

    /* move to next set of modes */
    do {
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s ",&word[0][0]);

    } while(strcmp(&word[0][0],"FREQUENCY:"));


    /* scan the frequencies; 
     * this requires some care since there might
     * be imaginary modes present which leads to
     * an additional char per imaginary mode per
     * line */

    /* check for imaginary modes */
    if ( strchr(buffer,'I') != NULL ) {
      /* initialize */
      separator = ' ';
      counter = 0;

      /* read all line elements into individual strings */

      /* skip first entry "FREQUENCY" */
      item = strtok(buffer,&separator);

      
      /* start going through the string */
      while ( (item = strtok(NULL,&separator)) != NULL ) 
      {
        /* check if item is 'I'; if yes, mark previous mode
         * as imaginary; otherwise save the mode */
        if ( *item == 'I' ) data->nimag++;

        else {
          /* save only the first remaining columns modes - 
           * there NEVER should
           * be more in any case, but just to make sure
           * we don't overrun the array */
          if (counter<remaining_columns) {
            *(data->wavenumbers+(i*5)+counter) = atof(item);
            counter++;
          }
        }
      } 
    }

    /* no imaginary mode, reading is straightforward */
    else {
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0],
             &entry[0],&entry[1],&entry[2],&entry[3],&entry[4]); 
 
      for (k=0; k<remaining_columns; k++) {
        *(data->wavenumbers+(i*5)+k) = entry[k]; 
      }
    }

    eatline(data->file, 1);
    
    /* next line contains the IR INTENSITIES */
    GET_LINE(buffer, data->file);

    /* scan the IR INTENSITIES */
    sscanf(buffer,"%s %s %lf %lf %lf %lf %lf",&word[0][0],&word[1][0],
           &entry[0],&entry[1],&entry[2],&entry[3],&entry[4]);
 
    for (k=0; k<remaining_columns; k++) {
      *(data->intensities+(i*5)+k) = entry[k];
    }

    eatline(data->file, 1);

    /* read the following five modes */
    for (k=0; k<data->numatoms; k++) {
      /* x */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %s %s %lf %lf %lf %lf %lf",&word[0][0], 
             &word[1][0], &word[2][0], &entry[0], &entry[1], &entry[2],
             &entry[3], &entry[4]);

      for (j=0; j<remaining_columns; j++) {
        *(data->normal_modes+(3*k)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }
      
      /* y */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0], &entry[0],
             &entry[1],&entry[2], &entry[3],&entry[4]);
      
      for (j=0; j<remaining_columns; j++) {
        *(data->normal_modes+(3*k+1)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }

      /* z */
      GET_LINE(buffer, data->file);
      sscanf(buffer,"%s %lf %lf %lf %lf %lf",&word[0][0], &entry[0],
             &entry[1], &entry[2], &entry[3],&entry[4]);
      
      for (j=0; j<remaining_columns; j++) {
        *(data->normal_modes+(3*k+2)+((i*5+j)*3*data->numatoms)) = 
          entry[j];
      }
    }
  }

  data->have_normal_modes = TRUE;

  /* release memory that is not needed any more */
  free(item);

  /* print brief message */
  printf("gamessplugin) Successfully scanned normal modes \n");

  return TRUE;
}



/***********************************************************
 *
 * this function reads in the cartesian hessian matrix 
 *
 * *********************************************************/
static int get_cart_hessian(gamessdata *data)
{
  char word[4][BUFSIZ];
  char buffer[BUFSIZ];
  char dummy; 
  int i,j,k;
  float entry[6]; 

  buffer[0] = '\0';
  memset(entry, 0, sizeof(entry));
  for (i=0; i<4; i++) word[i][0] = '\0';


  /* at this point we need to rewind the file, since
   * in case that there is no internal Hessian stuff the
   * previous call to get_int_coords scanned the file
   * until EOF */
  rewind(data->file);


  /* look for CARTESIAN FORCE CONSTANT MATRIX */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %s %s",&word[0][0],&word[1][0],
           &word[2][0],&word[3][0]);

  } while(strcmp(&word[0][0],"CARTESIAN") || 
          strcmp(&word[1][0],"FORCE") ||
          strcmp(&word[2][0],"CONSTANT") ||
          strcmp(&word[3][0],"MATRIX"));


  /* skip next 5 lines */
  for (i=0; i<5; i++) eatline(data->file, 1);


  /* reserve memory for array; 
   * NOTE: this is a lower triangular matrix, but for now
   * we save it in an square matrix of dim(3Nx3N) to 
   * facilitate element access */
  data->carthessian = 
    (double *)calloc((data->numatoms*3)*(data->numatoms*3),
                 sizeof(double));

  
  /* make sure memory was allocated properly */
  if (data->carthessian == NULL) {
    PRINTERR;
    return FALSE;
  }


  /* start scanning; the cartesian hessian matrix is a lower
   * triangular matrix, organized in rows of 6 */

  /* read blocks with complete rows of 6 */
  for (i=0; i<(int)(data->numatoms/2); i++) {
    for (j=0; j<(data->numatoms*3)-(i*6); j++) {
      GET_LINE(buffer, data->file);
 
      if (j%3==0) {
        sscanf(buffer,"%s %s %c %f %f %f %f %f %f",
               &word[0][0],&word[1][0],&dummy,&entry[0],&entry[1],
               &entry[2],&entry[3],&entry[4],&entry[5]);
      }
      else {
        sscanf(buffer,"%1s %f %f %f %f %f %f",
               &dummy,&entry[0],&entry[1],&entry[2],&entry[3],&entry[4],
               &entry[5]);
      }


      /* save entries (lower triangular matrix) in a 
       * square matrix */
      for (k=0; k<=(j<5 ? j : 5); k++) {
        *(data->carthessian+((j+(i*6))*3*data->numatoms)+
          (k+(i*6))) = entry[k];
      }
    }

    /* skip the three line separating the matrix entries */
    eatline(data->file, 4);
  }

  printf("gamessplugin) Scanned Hessian in CARTESIAN coordinates\n");

  data->have_cart_hessian = TRUE;

  return TRUE;
}
  
  
  
/***********************************************************
 *
 * this function reads in the hessian in internal coordinates
 * as well as the internal coordinates 
 *
 * *********************************************************/
static int get_int_coords(gamessdata *data) {

  char word[5][BUFSIZ];
  char buffer[BUFSIZ];
  long position;
  int first, second, third, fourth;
  double value;
  double hess[5];
  int i = 0, j = 0, k = 0, l = 0;
  int n, dummy, remaining_blocks;

  buffer[0] = '\0';
  memset(hess, 0, sizeof(hess));
  for (i=0; i<5; i++) word[i][0] = '\0';

  /* look for list of INTERNAL COORDINATES */
  do {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s ",&word[0][0],&word[1][0]);

  } while(strcmp(&word[0][0],"INTERNAL") || 
          strcmp(&word[1][0],"COORDINATES"));

  
  /* skip next 5 lines */
  for (i=0; i<5; i++) eatline(data->file, 1);

  /* remember current position so we can jump back */
  position = ftell(data->file);

  /* scan the next line */
  GET_LINE(buffer, data->file);
  n = sscanf(buffer,"%s %s", &word[0][0], &word[1][0]); 

  /* read line by line */
  while (n!=-1) {
    /* start counting the number of internal coordinates */
    data->nintcoords++;

    /* count the number of bonds, angles, dihedrals */
    if (!strcmp(&word[1][0],"STRETCH")) {
      data->nbonds++;
    }
    else if (!strcmp(&word[1][0],"BEND")) {
      data->nangles++;
    }
    else if (!strcmp(&word[1][0],"TORSION")) {
      data->ndiheds++;
    }
    else if (!strcmp(&word[1][0],"PLA.BEND")) {
      data->nimprops++;
    }

    /* scan next line */
    GET_LINE(buffer, data->file);
    n = sscanf(buffer,"%s %s", &word[0][0], &word[1][0]); 
  }

  /* now that we know the number of bonds, angles, etc.
   * we can read and store the internal coordinates */
  fseek(data->file,position,SEEK_SET);


  /* reserve memory for the arrays storing the internal
   * coordinates and their values */
  data->bonds = (int *)calloc(2*data->nbonds,sizeof(int));
  data->angles = (int *)calloc(3*data->nangles,sizeof(int));
  data->dihedrals = (int *)calloc(4*data->ndiheds,sizeof(int));
  data->impropers = (int *)calloc(4*data->nimprops,sizeof(int));
  data->internal_coordinates = (double *)calloc(data->nintcoords,
      sizeof(double));


  /* check if we have sufficient memory available */
  if ( (data->bonds == NULL) || 
       (data->angles == NULL) ||
       (data->dihedrals == NULL) || 
       (data->internal_coordinates == NULL)) 
  {
    PRINTERR; 
    return FALSE;
  }


  /* now start going through the internal coordinates
   * and save them in the appropriate arrays; here
   * I drop all safety check since we went through
   * this part of the file already once and should
   * be good */
 
  /* scan the STRETCHES */
  for (i=0; i<data->nbonds; i++) {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %d %d %lf", &word[0][0], &word[1][0], 
           &first, &second, &value);

    *(data->bonds+2*i) = first;
    *(data->bonds+2*i+1) = second;
    *(data->internal_coordinates+i) = value;
  }

  /* scan the BENDS */
  for (j=0; j<data->nangles; j++) {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %d %d %d %lf", &word[0][0], &word[1][0], 
           &first, &second, &third, &value);

    *(data->angles+3*j) = first;
    *(data->angles+3*j+1) = second;
    *(data->angles+3*j+2) = third;
    *(data->internal_coordinates+i+j) = value;
  }

  /* scan the TORSIONS */
  for (k=0; k<data->ndiheds; k++) {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %d %d %d %d %lf", &word[0][0], &word[1][0],
           &first, &second, &third, &fourth, &value);

    *(data->dihedrals+4*k) = first;
    *(data->dihedrals+4*k+1) = second;
    *(data->dihedrals+4*k+2) = third;
    *(data->dihedrals+4*k+3) = fourth;
    *(data->internal_coordinates+i+j+k) = value;
  }

  /* scan the IMPROPERS */
  for (l=0; l<data->nimprops; l++) {
    GET_LINE(buffer, data->file);
    sscanf(buffer,"%s %s %d %d %d %d %lf", &word[0][0], &word[1][0],
           &first, &second, &third, &fourth, &value);

    *(data->impropers+4*l) = first;
    *(data->impropers+4*l+1) = second;
    *(data->impropers+4*l+2) = third;
    *(data->impropers+4*l+3) = fourth;
    *(data->internal_coordinates+i+j+k+l) = value;
  }

  printf("gamessplugin) Scanned %d INTERNAL coordinates \n",
         data->nintcoords);
  printf("gamessplugin)    %d BONDS \n",data->nbonds);
  printf("gamessplugin)    %d ANGLES \n",data->nangles);
  printf("gamessplugin)    %d DIHEDRALS \n",data->ndiheds);
  printf("gamessplugin)    %d IMPROPERS \n",data->nimprops);


  /* next read in the hessian in internal coordinates;
   * we would expect the matrix immediately after the
   * internal coordinates in the output files;
   * we check this first */
  eatline(data->file, 1);


  GET_LINE(buffer, data->file);
  sscanf(buffer,"%s %s %s %s %s", &word[0][0], &word[1][0], 
      &word[2][0], &word[3][0], &word[4][0]);

  if ( strcmp(&word[0][0],"HESSIAN") || 
       strcmp(&word[1][0],"MATRIX") ||
       strcmp(&word[3][0],"INTERNAL") || 
       strcmp(&word[4][0],"COORDINATES")) 
  {
    /* apparently we are out of luck - no Hessian in internal
     * coordinates -- GOOD BYE :) */
    return FALSE;
  }
 

  /* skip to the hessian arrays */
  while (sscanf(buffer,"%d %lf %lf %lf %lf %lf", 
                &dummy, &hess[0], &hess[1], &hess[2], 
                &hess[3], &hess[4]) != 6) {
    GET_LINE(buffer, data->file);
  }

  
  /* reserve memory for inthessian array */
  data->inthessian = 
    (double *)calloc((data->nintcoords)*(data->nintcoords),
                 sizeof(double));


  /* make sure memory was allocated properly */
  if (data->inthessian == NULL) {
    PRINTERR;
    return FALSE;
  }


  /* start scanning; GAMESS organized the output of the
   * internal HESSIAN in rows of 5 */

  /* read blocks with complete rows of 5 */
  for (i=0; i<(int)(data->nintcoords/5); i++) {
    for (j=0; j<data->nintcoords; j++) {
      sscanf(buffer,"%d %lf %lf %lf %lf %lf", &dummy, &hess[0], 
             &hess[1], &hess[2], &hess[3], &hess[4]);

      /* save entries */
      for (k=0; k<5; k++) { 
        *(data->inthessian+(j*data->nintcoords)+(i*5)+k) = hess[k];
      }

      /* next line */
      GET_LINE(buffer, data->file);
    }

    /* skip the two lines separating the matrix entries 
     * and scan next line */
    eatline(data->file, 2);

    GET_LINE(buffer, data->file);
  }

  
  /* read the remaining block with less then 5 rows
   * if present */
  remaining_blocks = data->nintcoords%5;
  
  if (remaining_blocks!=0) {
    for (j=0; j<data->nintcoords; j++) {
      sscanf(buffer,"%d %lf %lf %lf %lf %lf", &dummy, &hess[0], 
             &hess[1], &hess[2], &hess[3], &hess[4]);

      for (k=0; k<remaining_blocks; k++) { 
        *(data->inthessian+(j*data->nintcoords)+(i*5)+k) = hess[k];
      }

      GET_LINE(buffer, data->file);
    }
  }

  printf("gamessplugin) Scanned Hessian in INTERNAL coordinates\n");

  /* finally, dump the diagonal elements of the hessian into the
   * force constant arrays, after converting the units 
   * appropriately;
   * BONDS are in HARTREE/BOHR**2
   * ANGLES,DIHEDRALS,IMPROPERS are in HARTREE/RADIAN**2 */
  
  /* allocate dynamic arrays */
  data->bond_force_const = 
    (double *)calloc(data->nbonds,sizeof(double));
 
  if (data->bond_force_const==NULL) {
    PRINTERR;
    return FALSE;
  }


  data->angle_force_const =
    (double *)calloc(data->nangles,sizeof(double));

  if (data->angle_force_const==NULL) {
    PRINTERR;
    return FALSE;
  }


  data->dihedral_force_const =
    (double *)calloc(data->ndiheds,sizeof(double));

  if (data->dihedral_force_const==NULL) {
    PRINTERR;
    return FALSE;
  }


  data->improper_force_const =
    (double *)calloc(data->nimprops,sizeof(double));

  if (data->improper_force_const==NULL) {
    PRINTERR;
    return FALSE;
  }

  /* scan the bonds */
  for (i=0; i<data->nbonds; i++) {
    *(data->bond_force_const + i) = 
      *(data->inthessian+(i*data->nintcoords)+i) 
      * HARTREE_TO_KCAL / BOHR_TO_ANGS / BOHR_TO_ANGS;

    printf("%3d (BOND) %2d - %2d : %f (CHARMM) %f \n",i, 
           *(data->bonds+2*i), *(data->bonds+2*i+1),
           *(data->bond_force_const +i),
           *(data->bond_force_const +i)*0.5);
  }
  
  /* scan the angles */
  for (j=i; j<i+(data->nangles); j++) {
    *(data->angle_force_const + (j-i)) = 
      *(data->inthessian+(j*data->nintcoords)+j) 
      * HARTREE_TO_KCAL;
    
    printf("%3d (ANGLE) %2d - %2d - %2d : %f (CHARMM) %f \n",j,
           *(data->angles+3*(j-i)), *(data->angles+3*(j-i)+1), 
           *(data->angles+3*(j-i)+2), 
           *(data->angle_force_const + (j-i)),
           *(data->angle_force_const + (j-i))*0.5);
  }

  /* scan the dihedrals */
  for (k=j; k<j+(data->ndiheds); k++) {
    *(data->dihedral_force_const + (k-j)) = 
      *(data->inthessian+(k*data->nintcoords)+k)
      * HARTREE_TO_KCAL;
    
    printf("%3d (DIHEDRAL) %2d - %2d - %2d - %2d : %f \n",k,
           *(data->dihedrals+4*(k-j)), *(data->dihedrals+4*(k-j)+1),
           *(data->dihedrals+4*(k-j)+2), *(data->dihedrals+4*(k-j)+3),
           *(data->dihedral_force_const + (k-j))); 
  }

  /* scan the impropers */
  for (l=k; l<k+(data->nimprops); l++) {
    *(data->improper_force_const + (l-k)) = 
      *(data->inthessian+(l*data->nintcoords)+l)
      * HARTREE_TO_KCAL;
    
    printf("%3d (IMPROPERS) %2d - %2d - %2d - %2d : %f \n",l,
           *(data->impropers+4*(l-k)), *(data->impropers+4*(l-k)+1),
           *(data->impropers+4*(l-k)+2), *(data->impropers+4*(l-k)+3),
           *(data->improper_force_const + (l-k)));
  }


  data->have_internals = TRUE;
  return TRUE;
}


#if 0

/************************************************************
 *
 * this function animates a given normal mode by means of
 * generating mod_num_frames frames away from the equilibrium
 * structure in a direction given by the hessiane 
 *
 ************************************************************/
static int animate_normal_mode(gamessdata *data, int mode)
{
  mode_data *animated_mode = data->animated_mode;
  float *normal_modes = data->normal_modes;
  float scale = animated_mode->mode_scaling;
  int i = 0, k = 0; 
  int l = 0, m = 0;
  int natoms = data->numatoms;
  int num_frames = animated_mode->mode_num_frames;

  /* first sweep to max of interval */
  for ( k = 0; k < num_frames+1; ++k)
  {
    for ( i = 0; i < natoms; ++i)
    {
      *(animated_mode->mode_frames+(k*natoms*3)+(3*i)) = 
        (data->initatoms+i)->x * (1+( k*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i)))));

      *(animated_mode->mode_frames+(k*natoms*3)+(3*i+1)) = 
        (data->initatoms+i)->y * (1+( k*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+1)))));

      *(animated_mode->mode_frames+(k*natoms*3)+(3*i+2)) = 
        (data->initatoms+i)->z * (1+( k*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
    }
  }


  /* second sweep all the way back to min of interval */
  for ( l = 0; l < 2*num_frames+1; ++l)
  {
    for ( i = 0; i < natoms; ++i)
    {
      *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i)) = 
        (data->initatoms+i)->x * (1+((int)(num_frames-l)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i)))));

      *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i+1)) = 
        (data->initatoms+i)->y * (1+((int)(num_frames-l)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+1)))));

      *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i+2)) = 
        (data->initatoms+i)->z * (1+((int)(num_frames-l)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
    }
  }


  /* third sweep back to the native starting structure */
  for ( m = 0; m < num_frames+1; ++m)
  {
    for ( i = 0; i < natoms; ++i)
    {
      *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i)) = 
        (data->initatoms+i)->x * (1+((int)(m-num_frames)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i)))));

      *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i+1)) = 
        (data->initatoms+i)->y * (1+((int)(m-num_frames)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+1)))));

      *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i+2)) = 
        (data->initatoms+i)->z * (1+((int)(m-num_frames)*scale * 
        (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
    }
  }

  printf("gamessplugin) Successfully animated mode %d \n", mode);

  return TRUE;
}
#endif


/********************************************************
 *
 * this function removes trailing spaces/newlines off
 * a character array
 *
 * FIXME: This function is dangerous; if the string is 
 * screwed up and doesn't contain any char that the while
 * loop is looking for we get a nice segfault :(
 ********************************************************/
static char* chop_string_all(char* the_string)
{
  int i = 0;

  while (the_string[i]!='\n' && the_string[i]!=' ' && 
         the_string[i]!='\0') {
    i++;
  }

  the_string[i] = '\0';

  return the_string;
}


/********************************************************
 *
 * this function returns a pointer to the first non-whitespace
 * character in a string.
 * The c-string must be null-terminated.
 *
 ********************************************************/
static char* trimleft(char* the_string)
{
  char *new_string = the_string;
  while ( (*new_string=='\n' || *new_string==' ' || *new_string=='\t') && 
        (*new_string != '\0'))
  {
    new_string++;
  }

  return new_string;
}

static char* trimright(char* s)
{
  int i;
  for (i=strlen(s)-1; i>=0; i--) {
    if (!isspace(s[i])) break;
  }
  s[i+1] = '\0';

  return s;
}

/* Advances the file pointer until the first appearance
 * of a line beginning with the given keystring. Leading
 * whitespace in the lines are ignored, the keystring 
 * should not begin with whitespace otherwise no match
 * will be found.
 * If stopstring is encountered before the keystring 
 * the file is rewound to the position where the search
 * started. If stopstring is NULL then the search stops
 * at EOF. */
static int goto_keyline(FILE *file, const char *keystring,
        const char *stopstring) {
  char buffer[BUFSIZ];
  char *line;
  int found = 0;
  long filepos, curline;
  filepos = ftell(file);

  do {
    curline = ftell(file);
    if (!fgets(buffer, sizeof(buffer), file)) {
      fseek(file, filepos, SEEK_SET);
      return 0;
    }
    line = trimleft(buffer);
    if (strstr(line, keystring)) {
      found = 1;
      fseek(file, curline, SEEK_SET);
      break;
    }
  } while (!stopstring || !strstr(line, stopstring));
    
  if (!found) {
    fseek(file, filepos, SEEK_SET);
    return STOPPED;
  }

  return FOUND;
}

/* places file pointer AFTER the line containing the string */
static int pass_keyline(FILE *file, const char *keystring,
        const char *stopstring) {
  char buffer[BUFSIZ];
  char *line;
  int found = 0;
  long filepos;
  filepos = ftell(file);

  do {
    if (!fgets(buffer, sizeof(buffer), file)) {
      fseek(file, filepos, SEEK_SET);
      return 0;
    }
    line = trimleft(buffer);
    if (strstr(line, keystring)) {
      found = 1;
      break;
    }
  } while (!stopstring || !strstr(line, stopstring));
    
  if (!found) {
    fseek(file, filepos, SEEK_SET);
    return STOPPED;
  }

  return FOUND;
}

static int goto_keystring2(FILE *file, const char *keystring,
        const char *stopstring1, const char *stopstring2) {
  char buffer[BUFSIZ];
  char *line;
  int found = 0;
  long filepos;
  filepos = ftell(file);

  do {
    if (!fgets(buffer, sizeof(buffer), file)) break;
    line = trimleft(buffer);
    if (strstr(line, keystring)) {
      found = 1;
      break;
    }
  } while (!stopstring1 || !strstr(line, stopstring1) || 
           !stopstring2 || !strstr(line, stopstring2));
    
  if (!found) {
    fseek(file, filepos, SEEK_SET);
    return 0;
  }

  return 1;
}

static void eatwhitelines(FILE *file) {
  char buffer[BUFSIZ];
  long filepos;
  filepos = ftell(file);
  while (fgets(buffer, sizeof(buffer), file)) {
    if (strlen(trimright(buffer))) {
      fseek(file, filepos, SEEK_SET);
      break;
    }
    filepos = ftell(file);
  }
}



static void whereami(FILE *file) {
  char buffer[BUFSIZ];
  char *line;
  long filepos;
  filepos = ftell(file);
  do {
    if (!fgets(buffer, sizeof(buffer), file)) {
      if (feof(file)) printf("HERE) EOF\n");
      else printf("HERE) ????\n");
      return;
    }
    line = trimleft(buffer);
  } while (!strlen(line));

  printf("HERE) %s\n", buffer);
  fseek(file, filepos, SEEK_SET);
}

/********************************************************
 *
 * this function removes trailing newlines off
 * a character array
 *
 ********************************************************/
static char* chop_string_nl(char* the_string)
{
  int i = 0;

  while (the_string[i]!='\n' && the_string[i]!='\0') i++;

  the_string[i] = '\0';

  return the_string;
}



/*************************************************************
 *
 * plugin registration 
 *
 **************************************************************/
static molfile_plugin_t plugin;

VMDPLUGIN_API int VMDPLUGIN_init(void) {
  memset(&plugin, 0, sizeof(molfile_plugin_t));
  plugin.abiversion = vmdplugin_ABIVERSION;
  plugin.type = MOLFILE_PLUGIN_TYPE;
  plugin.name = "gamess";
  plugin.prettyname = "GAMESS";
  plugin.author = "Markus Dittrich, Jan Saam";
  plugin.majorv = 0;
  plugin.minorv = 11;
  plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
  plugin.filename_extension = "log";
  plugin.open_file_read = open_gamess_read;
  plugin.read_structure = read_gamess_structure;
  plugin.close_file_read = close_gamess_read;

  plugin.read_qm_metadata = read_gamess_metadata;
  plugin.read_qm_rundata  = read_gamess_rundata;

#if vmdplugin_ABIVERSION > 11
  plugin.read_timestep_metadata    = read_timestep_metadata;
  plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
  plugin.read_timestep = read_timestep;
#endif

  return VMDPLUGIN_SUCCESS;
}

VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
  (*cb)(v, (vmdplugin_t *)&plugin);
  return VMDPLUGIN_SUCCESS;
}

VMDPLUGIN_API int VMDPLUGIN_fini(void) {
  return VMDPLUGIN_SUCCESS;
}