---

gamessplugin.c

Go to the documentation of this file.

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


/* *******************************************************
 *
 *          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 "gamessplugin.h"
 

/* 
 * pre-processor macro to activate grid code 
 */
#ifdef GRID_ACTIVE
#define GRID 1
#else
#define GRID 0
#endif

/*
 * Error reporting macro for use in DEBUG mode
 */
#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 ERR_FALSE(x) if ( x == NULL ) return FALSE;


/***************************************************************
 *
 * subroutine doing the initial reading of the GAMESS
 * input file
 *
 * *************************************************************/
static void *open_gamess_read(const char *filename, 
                  const char *filetype, int *natoms) {

  FILE *fd;
  gamessdata *data;
  char mytime[BUFSIZ];

  /* initialize array */
  mytime[0] = '\0';

  
  /* open the input file */
  fd = fopen(filename, "rb");
 
  /* make sure we were able to load the file properly*/
  if (!fd) 
  {
    PRINTERR;
    return NULL;
  }

  /* obtain the current time */
  get_time(mytime);

  /* print out a short warning message -- we're very alpha */
  printf("\n");
  printf("     *************************************************\n");
  printf("     ***       GAMESSPLUGIN for VMD                ***\n");
  printf("     *** %s      ***\n", mytime);
  printf("     *** email bugs to vmd@ks.uiuc.edu             ***\n");
  printf("     *** v0.4.2 02/20/2007 (c)     Markus Dittrich ***\n");
  printf("     *************************************************\n");
  printf("\n");

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

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


  /**********************************
   * initialize some variables
   *********************************/

  /* volumetric */
  data->have_volumetric = 0; 

  /* runtyp */
  data->runtyp = 0;

  /* initialize the number of atomic shells */
  data->num_shells = 0;
  
  /* have_trajectory flag to 1 */
  data->have_trajectory = 1; 

  /* reset trajectory counter */
  data->num_traj_points = 0;

  /* presence of wavefunction */
  data->got_wavefunction = FALSE;

  /* orbital index of HOMO */
  data->homo_index = 0;

  /* flag indicating the presence of
   * internal coordinates */
  data->have_internals = FALSE;

  /* flag indicating the presence of
   * the cartesian Hessian */
  data->have_cart_hessian = FALSE;

  /* initialize the number of imaginary 
   * modes for HESSIAN type runs */
  data->nimag = 0;

  /* initialize internal coordinate counters */
  data->nintcoords = 0;
  data->nbonds = 0;
  data->nangles = 0;
  data->ndiheds = 0;
  data->nimprops = 0;

  /* initialize the GAMESS version and the PC Gamess flag*/
  data->version = 0;
  data->have_pcgamess = 0;

  /* initialize the basis set info */
  data->ngauss = 0;
  data->npfunc = 0;
  data->ndfunc = 0;
  data->nffunc = 0;
  data->diffs  = 0;
  data->diffsp = 0;

  /* initialize the number of scfenergies */
  data->num_scfenergies = 0;

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

  
  /***********************************
   * done intializing
   **********************************/


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


  /* 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 "static" information from the log file */    
    if ( parse_gamess_log_static(data,natoms) == FALSE ) 
      return NULL;
  }
  else 
  {
    return NULL;
  }


  /* done with the gamess plugin */
  rewind(fd);
  printf("\n");


  /* copy pointer also in bogus tcl data pointer */
  /* tcl_pointer = data; */

  return data;
}


/************************************************************
 * 
 * subroutine reading in the structure of the molecule
 * in the GAMESS output file
 *
 *************************************************************/
static int read_gamess_structure(void *mydata, int *optflags, 
                      molfile_atom_t *atoms) 
{
  gamessdata *data = (gamessdata *)mydata;
  gamess_temp *temp_atom;
  molfile_atom_t *atom;
  unsigned int i = 0;
 
  *optflags = MOLFILE_NOOPTIONS; /* no optional data */

  /* 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 temp arrays */
  temp_atom = data->temporary;

  for(i=0;i<data->numatoms;i++)
  {
    atom=atoms+i;
    strncpy(atom->name,temp_atom->type,sizeof(atom->name)); 
    strncpy(atom->type,temp_atom->type,sizeof(atom->type));
    strncpy(atom->resname,"",sizeof(atom->resname)); 
    atom->resid = 1;
    atom->charge = temp_atom->charge;
    atom->chain[0] = '\0';
    atom->segid[0] = '\0';
   

    temp_atom++;
  }
 
  return MOLFILE_SUCCESS; 
}


/***********************************************************
 *
 * this function reads in the information for the next
 * timestep
 *
 ***********************************************************/
static int read_next_timestep(void *mydata, int natoms, 
    molfile_timestep_t *ts) 
{
  gamessdata *data = (gamessdata *)mydata;
  gamess_temp *temp_atom;
  unsigned int i = 0;

#ifdef ANIMATE_MODE
  mode_data *animated_mode = data->animated_mode;
#endif
  
  /* Now, if we are not dealing with a trajectory (a single 
   * * point conformation is considered one element of a trajectory
   */
  if (data->have_trajectory == 0) return MOLFILE_ERROR;
 

  /* in the case of RUNTYP=ENERGY read data from temp
   * arrays and that's it;
   * in the case of RUNTYP=OPTIMIZE and SADPOINT just start with
   * the first 1NSERCH entry */
  if ( data->runtyp == ENERGY ) 
  {
     /* initialize pointer for temporary arrays */
     temp_atom = data->temporary; 
 

     /* now copy the initial coordinates from the temporary
      * arrays into the VMD ones */
     for(i=0;i<natoms;i++)  {

       ts->coords[3*i  ] = temp_atom->x;
       ts->coords[3*i+1] = temp_atom->y;
       ts->coords[3*i+2] = temp_atom->z; 

       temp_atom++;
     }    

     /* Since we are only dealing with a single point which has
      * just been read we have to make sure to let VMD know 
      * and also register the number of trajectory points as 1 */
     data->have_trajectory = 0;

     return MOLFILE_SUCCESS;
  }

#ifndef ANIMATE_MODE
  else if ( data->runtyp == HESSIAN ) 
  {
     /* initialize pointer for temporary arrays */
     temp_atom = data->temporary; 
 

     /* now copy the initial coordinates from the temporary
      * arrays into the VMD ones */
     for(i=0;i<natoms;i++)  {

       ts->coords[3*i  ] = temp_atom->x;
       ts->coords[3*i+1] = temp_atom->y;
       ts->coords[3*i+2] = temp_atom->z; 

       temp_atom++;
     }    

     /* Since we are only dealing with a single point which has
      * just been read we have to make sure to let VMD know 
      * and also register the number of trajectory points as 1 */
     data->have_trajectory = 0;

     return MOLFILE_SUCCESS;
  }
#endif

#ifdef ANIMATE_MODE
  else if ( data->runtyp == HESSIAN ) 
  {
    /* read until we run out of frames for the current mode */
    if ( animated_mode->current_mode_frame < 
              4*animated_mode->mode_num_frames+3 )
    {
      /* now copy the initial coordinates from the temporary
        * arrays into the VMD ones */
      for(i=0; i< natoms; i++)  
      {
        ts->coords[3*i  ] = *(animated_mode->mode_frames +
            (animated_mode->current_mode_frame*3*natoms)+3*i);

        ts->coords[3*i+1] = *(animated_mode->mode_frames +
            (animated_mode->current_mode_frame*3*natoms)+3*i+1);

        ts->coords[3*i+2] = *(animated_mode->mode_frames +
            (animated_mode->current_mode_frame*3*natoms)+3*i+2);
      }    

      /* increase the trajectory and mode counter */ 
      (animated_mode->current_mode_frame)++;
    }
    else 
    {
      return MOLFILE_ERROR;
    }

    return MOLFILE_SUCCESS;
  }
#endif

  else if ( data->runtyp == OPTIMIZE || data->runtyp == SADPOINT ) 
  {
    /* call the routine scanning the output file for the
     * trajectory */
    if ( get_trajectory(data,ts,natoms) == FALSE) {

      /* before we return let's check if we have read
       * any trajectory points whatsoever; if not
       * this could mean that the file might be
       * truncated, and in this case we dump the initial
       * coordinate info, such that the user has something
       * to look at; */
      if ( data->num_traj_points == 0 ) {

        /* initialize pointer for temporary arrays */
        temp_atom = data->temporary; 
 
        
        /* now copy the initial coordinates from the temporary
        * arrays into the VMD ones */
        for(i=0;i<natoms;i++)  {

          ts->coords[3*i  ] = temp_atom->x;
          ts->coords[3*i+1] = temp_atom->y;
          ts->coords[3*i+2] = temp_atom->z; 

          temp_atom++; }


        /* make sure we don't continue reading the
         * trajectory */
        data->have_trajectory = 0;


        /* that's all we can do */
        return MOLFILE_SUCCESS;
      }


      return MOLFILE_ERROR;
    }


    /* increase the trajectory counter */
    (data->num_traj_points)++;
   

    /* success, it seems :) */
    return MOLFILE_SUCCESS;
  }


  /* all other cases should have been rejected allready
   * in open_gamess_read;
   * but just to make sure we return MOLFILE_ERROR 
   * by default */
  return MOLFILE_ERROR;
}


/*********************************************************
 *
 * this subroutine makes the volumetric orbital metadata
 * available to VMD;
 *
 ********************************************************/
static int read_orbital_metadata( void *mydata, int *nsets, 
                          molfile_volumetric_t **metadata) 
{
    gamessdata *data = (gamessdata *)mydata;

    /* first of all we have to test the presence of
     * volumetric data; if not set number of volumetric
     * data sets to zero and return */

    if ( data->have_volumetric != 1 ) 
    {
      *nsets = 0;
      return MOLFILE_SUCCESS;
    }

    /* we only have one set of volumetric data */
    *nsets = 1; 

    /* all the necessary data structures have been filled
     * in the function get_system_dimensions and we only
     * need to copy the pointer. */
    *metadata = data->vol;

    return MOLFILE_SUCCESS;
}


/**********************************************************
 *
 * this subroutine makes the volumetric data available 
 * to VMD
 *
 *********************************************************/
static int read_orbital_data( void *mydata, int set, 
    float *datablock, float *colorblock) 
{
  gamessdata *data = (gamessdata *)mydata;
  gamess_temp *temp_data;
  molfile_volumetric_t *vol_metadata;
  unsigned int i = 0;
  int numxyz;

  /* get pointer of temporary data storage */
  temp_data = data->temporary; 

  /* get pointer to volumetric metadata */
  vol_metadata = data->vol;

  /* move array containing volumetric over to datablock 
   * the array contains xsize*ysize*zsize elements */
  numxyz = ( vol_metadata->xsize * vol_metadata->ysize * 
             vol_metadata->zsize );

  for ( i = 0 ; i < numxyz ; ++i )
  {
    *(datablock + i) = *(data->orbital_grid+i);
  }

  return MOLFILE_SUCCESS;
}


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

  gamessdata *data = (gamessdata *)mydata;
  fclose(data->file);

  /* free memory */
  free(data->temporary);
  free(data->basis);
  free(data->system_dimensions);
  free(data->system_center); 
  free(data->orbital_grid);
  free(data->basis_counter);
  free(data->wave_function);
  free(data->orbital_symmetry);
  free(data->atomic_shells);
  free(data->shell_primitives);
  free(data->orbital_energy);
  free(data->atomic_number);
  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->vol);
  free(data->scfenergies);
  free(data->wavenumbers);
  free(data->intensities);
  free(data->normal_modes);
  
  if ( data->runtyp == HESSIAN )
  {
#ifdef ANIMATE_MODE
    free(data->animated_mode->mode_frames);
    free(data->animated_mode);
#endif
  }

  free(data);
}

#if vmdplugin_ABIVERSION > 9

/*****************************************************
 *
 * 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 *gamess_metadata)
{
  gamessdata *data = (gamessdata *)mydata;


  /* only reserve memory for cartesian/internal 
   * hessian if present */
  if ( data->runtyp == HESSIAN )
  {
    gamess_metadata->ncart = (3*data->numatoms);
    gamess_metadata->nimag = data->nimag;             
   
    /* possibly resever memory for internal hessian */
    if ( data->have_internals )
    {
      gamess_metadata->nintcoords = data->nintcoords; 
    }
    else
    {
      gamess_metadata->nintcoords = 0;
    }
  }
  else
  {
    gamess_metadata->ncart = 0;
    gamess_metadata->nimag = 0;
    gamess_metadata->nintcoords = 0;
  }

  
  /* fill in orbital data if present */
  gamess_metadata->num_basis_funcs = data->num_basis_funcs;
  gamess_metadata->orbital_counter = data->orbital_counter;  


  /* fill in the trajectory information */
  gamess_metadata->have_trajectory = data->have_trajectory;
  gamess_metadata->num_traj_points = data->num_traj_points;

  return MOLFILE_SUCCESS;
}



/******************************************************
 * 
 * file the VMD internal QM data structures with their
 * respective values
 *
 ******************************************************/
static int read_gamess_rundata(void *mydata, 
    molfile_qm_t *my_qm_data)
{
  gamessdata *data = (gamessdata *)mydata;
  unsigned int i;
  unsigned int ncart;
  molfile_qm_hessian_t hessian_data = my_qm_data->hess;
  molfile_qm_orbital_t orbital_data = my_qm_data->orb;
  molfile_qm_sysinfo_t sys_data = my_qm_data->run;
  /* molfile_qm_timestep_t *timestep_data = my_qm_data->timestep; */


  /* fill in molfile_qm_hessian_t if HESSIAN information
   * is available */
  if ( data->runtyp == HESSIAN)
  {
    /* initialize sizes */
    ncart = 3*data->numatoms;

    /* hessian/normalmodes in cartesian coordinates */ 
    for ( i = 0; i < ncart*ncart; ++i)
    {
      *(hessian_data.carthessian+i) = *(data->carthessian+i);
      *(hessian_data.normalmodes+i) = *(data->normal_modes+i);
    }

    /* wavenumbers/intensities */
    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);
    }

    /* internal hessian if present */
    if ( data->have_internals )
    {
      for ( i = 0; i < (data->nintcoords)*(data->nintcoords); ++i)
      {
        *(hessian_data.inthessian+i) = *(data->inthessian+i);
      }
    }
  }

  
  /* fill in molfile_qm_sysinfo_t */
  sys_data.runtyp = data->runtyp;
  sys_data.nproc = data->nproc;
  sys_data.num_gauss_basis_funcs = data->num_gauss_basis_funcs;
  sys_data.num_electrons = data->num_electrons;
  sys_data.totalcharge = data->totalcharge;
  sys_data.multiplicity = data->multiplicity;
  sys_data.num_orbitals_A = data->num_orbitals_A;
  sys_data.num_orbitals_B = data->num_orbitals_B;
  sys_data.scftyp = data->scftyp;
  sys_data.ngauss = data->ngauss; 
  sys_data.npfunc = data->npfunc;
  sys_data.ndfunc = data->ndfunc;
  sys_data.nffunc = data->nffunc;
  sys_data.diffs = data->diffs;
  sys_data.diffsp = data->diffsp;

  for ( i = 0; i < data->numatoms; ++i)
  {
    *(sys_data.nuc_charge+i) = (data->temporary+i)->charge;
  }

  strncpy(sys_data.runtyp_string, data->runtyp_string, 
          sizeof(sys_data.runtyp_string));

  strncpy(sys_data.scftyp_string, data->scftyp_string,
          sizeof(sys_data.scftyp_string));

  strncpy(sys_data.basis_string, data->basis_string,
          sizeof(sys_data.basis_string));
  
  strncpy(sys_data.memory, data->memory, sizeof(sys_data.memory)); 
  strncpy(sys_data.gbasis, data->gbasis, sizeof(sys_data.gbasis));
  strncpy(sys_data.runtitle, data->runtitle, sizeof(sys_data.runtitle));
  strncpy(sys_data.guess, data->guess, sizeof(sys_data.guess));
  strncpy(sys_data.geometry, data->geometry, sizeof(sys_data.geometry));
  strncpy(sys_data.version_string, data->version_string,
          sizeof(sys_data.version_string));


  /* fill in molfile_qm_orbital_t */
  for ( i = 0; i < data->numatoms; ++i)
  {
    *(orbital_data.basis_counter+i) = *(data->basis_counter+i);
    *(orbital_data.atomic_shells+i) = *(data->atomic_shells+i);
  }

  for ( i = 0; i < data->num_shells; ++i)
  {
    *(orbital_data.shell_primitives+i) = *(data->shell_primitives+i);
  }

  for ( i = 0; i < 2*data->num_basis_funcs; ++i)
  {
    *(orbital_data.basis+i) = *(data->basis+i);
  }

  for ( i = 0; i < data->num_basis_funcs; ++i)
  {
    *(orbital_data.orbital_symmetry+i) = *(data->orbital_symmetry+i);
  }

 
  /* NOTE: all the quantities in molfile_qm_timestep_t are
   * currently only read for the last frame and are therefore
   * not available for all timesteps;
   * ideally we don't want to the size of molfile_qm_timestep_t
   * to be equal to num_traj_points, but rather leave it up
   * to the user to select a single point, a certain number
   * of points, or all of them. Hence, num_traj_points should
   * eventually be replaced by num_available_frames */

  return MOLFILE_SUCCESS;
}

#endif

/*************************************************************
 *
 * registration stuff 
 *
 **************************************************************/
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";
  plugin.majorv = 0;
  plugin.minorv = 6;
  plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
  plugin.filename_extension = "log";
  plugin.open_file_read = open_gamess_read;
  plugin.read_structure = read_gamess_structure;
  plugin.read_next_timestep = read_next_timestep;
  plugin.close_file_read = close_gamess_read;
  plugin.read_volumetric_metadata = read_orbital_metadata;
  plugin.read_volumetric_data = read_orbital_data;

#if vmdplugin_ABIVERSION > 9
  plugin.read_qm_metadata = read_gamess_metadata;
  plugin.read_qm_rundata = read_gamess_rundata;
#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;
}




/********************************************************
 *
 * THE NEXT SECTION CONTAINs THE GAMESS OUTPUT FILE
 * SPECIFIC SUBROUTINES
 *
 ********************************************************/





/********************************************************
 *
 * this routine is the main gamess log file
 * parser responsible for static, i.e. 
 * non-trajectory information 
 *
 ********************************************************/
static int parse_gamess_log_static(void *mydata, int *natoms) 
{
  gamessdata *data = (gamessdata *)mydata;

  /* determine the number of procs and the amount
   * of memory requested */
  if ( get_proc_mem(data) == FALSE ) return FALSE;


  /* determine the GBASIS used for the run */
  if ( get_gbasis(data) == FALSE ) return FALSE;


  /* read the run title */
  if ( get_runtitle(data) == FALSE ) return FALSE;


  /* read the contrl group; here we determine the
   * RUNTYP and bomb out if we don't support
   * it; also read in some controlflow related
   * info as well as the COORD type */
  if ( check_contrl(data) == FALSE ) return FALSE;


  /* now call routine get_initial_info to obtain 
   * the number of atoms, the atom types, and the
   * coordinates, which will be stored in gamess_temp */
  if ( get_initial_info(data) == FALSE ) return FALSE;


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


  /* read in the number of orbitals, charge, 
   * multiplicity, ... */
  if ( get_num_orbitals(data) == FALSE ) return FALSE;


  /* check if the wavefunction and orbital stuff
   * is suported for current GBASIS; if not stop
   * here and return data pointer */
  if ( have_supported_gbasis(data) == FALSE) return TRUE;


  /* read in the guess options */
  if ( get_guess(data) == FALSE ) return FALSE;

  
  /* read in the basis set information */
  if ( get_basis(data) == FALSE ) return FALSE; 


  /* read in the final wavefunction; do this only for 
   * RUNTYP=ENERGY for now. For the other runs I have to
   * put in some more infrastructure to make sure to
   * correlate each of the possible many geometries
   * during e.g. an optimization run with the correstpoding
   * wavefunction 
   * for now we also have to restrict reading the wavefunctions
   * for singlets only, since otherwise the punched wavefunction
   * output differs;
   * At least the call to generate orbital grid should
   * not be done automagically but rather after a TCL call
   * only. 
   * NOTE: For PC Gamess we skip this step for now;
   * more work to be done here. */
  if ( data->runtyp == ENERGY 
         && data->num_orbitals_A == data->num_orbitals_B 
         && GRID
         && data->have_pcgamess != 1 )
  {
    /* read the wavefunction */
    if ( get_wavefunction(data) == FALSE) return FALSE;  


    /* figure out which one of the orbitals is the HOMO */
    if ( find_homo(data) == FALSE) return FALSE; 

 
    /* generate the grid containing the volumetric data
     * for an orbital */
    if( orbital_grid_driver(data) == FALSE ) return FALSE;  

  }

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

  return TRUE;
}

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

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

  /* QMData *data = (QMData *)mydata; */
  printf(" %10s WORDS OF MEMORY AVAILABLE\n", data->memory);
  printf("\n");
  printf("     BASIS OPTIONS\n");
  printf("     -------------\n");
  printf("     GBASIS=%8s     IGAUSS=%8i      POLAR=%8s\n", data->gbasis, data->ngauss, "NONE");
  printf("     NDFUNC=%8i", data->ndfunc);
  printf("     NFFUNC=%8i", data->nffunc);
  printf("     DIFFSP=%8s\n", TORF(data->diffsp));
  printf("     NPFUNC=%8i", data->npfunc);
  printf("      DIFFS=%8s\n", TORF(data->diffs));
  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 %6.1f", data->temporary[i].type, data->temporary[i].charge);
    
    printf("%17.10f", ANGS_TO_BOHR*data->temporary[i].x);
    printf("%20.10f", ANGS_TO_BOHR*data->temporary[i].y);
    printf("%20.10f\n", ANGS_TO_BOHR*data->temporary[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");

  for (i=0; i<data->numatoms; i++) {
    printf("%-8s\n\n", data->temporary[i].type);
    printf("\n");
    printf("nshells=%i\n", data->atomic_shells[i]);
    for (j=0; j<data->atomic_shells[i]; j++) {
      printf("nprim=%i\n", data->shell_primitives[j]);
      for (k=0; k<data->shell_primitives[j]; k++) {
        printf("%6i   %c %7i %22f", j, data->orbital_symmetry[count], count+1, data->basis[2*count]);
        if (data->orbital_symmetry[count] == 'S') {
          printf("%22f", data->basis[2*count+1]);
        }
        printf("\n");
        count++;
      }
      printf("\n");
    }
  }
  printf("\n");
  printf(" TOTAL NUMBER OF BASIS SET SHELLS             =%5i\n", data->num_shells);
  printf(" NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =%5i\n", data->num_gauss_basis_funcs);
  printf(" NUMBER OF ELECTRONS                          =%5i\n", data->num_electrons);
  printf(" CHARGE OF MOLECULE                           =%5i\n", data->totalcharge);
  printf(" SPIN MULTIPLICITY                            =%5i\n", data->multiplicity);
  printf(" NUMBER OF OCCUPIED ORBITALS (ALPHA)          =%5i\n", data->num_orbitals_A);
  printf(" NUMBER OF OCCUPIED ORBITALS (BETA )          =%5i\n", data->num_orbitals_B);
  printf(" TOTAL NUMBER OF ATOMS                        =%5i\n", data->numatoms);
  printf("\n");
}


/******************************************************
 *
 * this function extracts the trajectory information
 * from the output file
 *
 * *****************************************************/
static int get_trajectory(void *mydata, molfile_timestep_t *ts,
                   int natoms) 
{
  gamessdata *data = (gamessdata *)mydata;
  char buffer[BUFSIZ];
  char word[4][BUFSIZ];
  char type[8];
  unsigned int i = 0;
  float charge,x,y,z;
  double scfenergy = 0.0;

  /* zero out the word arrays */
  buffer[0] = '\0';
  type[0] = '\0';
  for ( i = 0; i < 4; ++i) word[i][0] = '\0';


  /* search for the first coordinate frame during optimization */
  ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )


  /* 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)
  {
     /* scan for the next optimized geometry punched in the 
      * logfile */
     do 
     {
       sscanf(buffer,"%s %s",&word[0][0],&word[1][0]); 
       ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )

     } while (strcmp(&word[0][0],"1NSERCH="));
  }
  else if ( data->version == 2 )
  {
     /* scan for the next optimized geometry punched in the 
      * logfile */
     do 
     {
       sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],
            &word[2][0]); 
       ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )

     } while (strcmp(&word[0][0],"BEGINNING") || 
              strcmp(&word[1][0],"GEOMETRY")  ||
              strcmp(&word[2][0],"SEARCH"));
  }


  /* next we need to check that we're reading the 
   * coordinates of all atoms not only the symmetry
   * unique ones in case the point group is not C1 */
  ERR_FALSE( fgets(buffer, sizeof(buffer), data->file) )
  sscanf(buffer,"%s %s %s", &word[0][0], &word[1][0],
      &word[2][0]);

  /* if symmetry unique atoms follow we have to advance
   * some more */
  if ( !strcmp(&word[2][0],"SYMMETRY") )
  {
    ERR_FALSE( fgets(buffer, sizeof(buffer), data->file) )
   
    do 
    {
      sscanf(buffer,"%s",&word[0][0]); 
      ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )

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

    /* skip next line */
    eatline(data->file);
  }
  else
  {
    /* skip the next two lines */
    eatline(data->file);
    eatline(data->file);
  }
  

  /* read the actual coordinates */
  for(i=0;i<natoms;i++) 
  {
    ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )
    sscanf(buffer,"%s %f %f %f %f",type,&charge,&x,&y,&z);  

    /* store the coordinates */
    ts->coords[3*i  ] = x; 
    ts->coords[3*i+1] = y;
    ts->coords[3*i+2] = z; 
  }     


  /* now we look for the SCF energy of this trajectory 
   * point */
   do 
   {
    ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )
    sscanf(buffer,"%s %s %s %lf",&word[0][0],&word[1][0],
        &word[2][0],&scfenergy); 
 
    } while (strcmp(&word[2][0],"ENERGY=") && 
             strcmp(&word[1][0],"ENERGY=")); 
  
  /* allocate more memory for scfenergy array */
  data->scfenergies = (double *)realloc(data->scfenergies,
        (data->num_traj_points+1)*sizeof(double));


  /* append current scfenergy to array and increase counter */
  *(data->scfenergies+data->num_traj_points) = scfenergy;
  data->num_scfenergies++;


  /* done with this trajectory point */ 
  return TRUE;
}



/******************************************************
 *
 * this function performs an initial file read
 * to check its validity and determine the number
 * of atoms in the system 
 *
 * *****************************************************/
static int get_initial_info(void *mydata)
{
  gamessdata *data = (gamessdata *)mydata;
  gamess_temp *atoms;
  char buffer[BUFSIZ];
  char word[4][BUFSIZ];
  char atname[BUFSIZ];
  float charge;
  float x,y,z;
  unsigned int numatoms;
  unsigned int bohr;
  unsigned int i,n;
  unsigned int have_normal_modes = TRUE;
  double scfenergy;
  char *status;

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

  /* alocate temporary data structure */
  atoms = (gamess_temp *)calloc(MAXQMATOMS,sizeof(gamess_temp));

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


  /* save pointer */
  data->temporary = atoms;


  /* counts the number of atoms in input file */
  numatoms=0;  

  /* look for the initial coordinate section 
   * NOTE: I also have to check if the coordinates are
   * in Angstrom or Bohr and convert them to A in the
   * latter case */
  rewind(data->file);

  do 
  {
    ERR_FALSE( fgets(buffer,sizeof(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)"));


  /* skip next line */
  eatline(data->file);


  /* now read in the coordinates until an empty line is being
   * encoutered */
  ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )
  n = sscanf(buffer,"%s %f %f %f %f",atname,&charge,&x,&y,&z);

  /* we expect 5 entries per line */
  while( n == 5 ) 
  {
    strncpy(atoms->type,atname,sizeof(atoms->type));
    atoms->charge = charge;


    /* if coordinates are in Bohr convert them to Angstrom here */
    if (bohr) 
    {
      x = BOHR_TO_ANGS * x;
      y = BOHR_TO_ANGS * y;
      z = BOHR_TO_ANGS * z;
    }

    atoms->x = x;
    atoms->y = y;
    atoms->z = z; 


    atoms++;
    numatoms++;

    ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )
    n = sscanf(buffer,"%s %f %f %f %f",atname,&charge,&x,&y,&z);
  }

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



  /* save a copy of atom charges in atomic_number array
   * of gamessdata */

  /* alocate temporary data structure */
  data->atomic_number = (int *)calloc(numatoms,sizeof(int));

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

  
  /* move data over */ 
  for ( n = 0; n < numatoms; n++) 
  {
    *(data->atomic_number + n) = (int)(data->temporary + n)->charge;
  }


  /* at this point we create an empty scfenergy arrar */
  data->scfenergies = (double *)calloc(0,sizeof(double));


  /* of we are dealing with a single point energy run 
   * we read in the final SCF energy and store it;
   * TODO: for now also include HESSIAN runs here even though
   * this is not completely correct in the case of numerical
   * runs*/
  if ( (data->runtyp == ENERGY) || (data->runtyp == HESSIAN) )
  {
    status = fgets(buffer,sizeof(buffer),data->file);
     
    while ( strcmp(&word[0][0],"FINAL") || 
            strcmp(&word[2][0],"ENERGY") ||
            strcmp(&word[3][0],"IS"))
    {
      sscanf(buffer,"%s %s %s %s %lf",&word[0][0],&word[1][0],
             &word[2][0],&word[3][0],&scfenergy);

      if ( fgets(buffer,sizeof(buffer),data->file) == NULL) break;
    }

    data->scfenergies = (double *)realloc(data->scfenergies,
             sizeof(double));
    *(data->scfenergies) = scfenergy;
    data->num_scfenergies++;
  }


  /* next we check for the Mulliken charges; if they are present
   * we set the have_mulliken charge flag, but don't bomb out */

  /* set have_mulliken initially to true */
  data->have_mulliken = TRUE;

  do 
  {
    status = fgets(buffer,sizeof(buffer),data->file);

    /* check for EOF */
    if (status == NULL) {
     
      data->have_mulliken = FALSE;
      break;
    }

    sscanf(buffer,"%s %s %s %s",&word[0][0],&word[1][0],
        &word[2][0],&word[3][0]);

  } while(strcmp(&word[0][0],"TOTAL") || 
          strcmp(&word[1][0],"MULLIKEN") ||
          strcmp(&word[2][0],"AND")   || 
          strcmp(&word[3][0],"LOWDIN") );


  /* if Mulliken Charges are present read them, if not, rewind
   * the file and look for the ESP charges */
  if ( data->have_mulliken == TRUE ) 
  {
    /* reserve memory for dynamically allocated
     * Mulliken charge array */
    data->mulliken_charges = 
      (double *)calloc(numatoms,sizeof(double));

    /* check for success */
    if (data->mulliken_charges == NULL) 
    {
      PRINTERR; 
      return FALSE;
    }


    /* skip next line */
    eatline(data->file);

    for ( n = 0; n < numatoms; ++n) 
    {
      /* make sure we don't get nuked by a bad file */
      ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )

      sscanf(buffer,"%s %s %s %lf ",&word[0][0], &word[1][0], 
          &word[2][0], (data->mulliken_charges+n));
    }
  }

  /* no Mulliken charges found; in this case we don't know where
   * we are in the file; hence rewind file */
  else 
  {
    rewind(data->file); 
  }


  /* next we check for the ESP charges; if they are present
   * we set the have_esp charge flag, but don't bomb out */

  /* set have_esd initially to true */
  data->have_esp = TRUE;

  do 
  {
    status = fgets(buffer,sizeof(buffer),data->file);

    /* check for EOF */
    if (status == NULL) {
     
      data->have_esp = FALSE;
      break;
    }

    sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],&word[2][0]);

  } while(strcmp(&word[0][0],"ATOM") || 
          strcmp(&word[1][0],"CHARGE") ||
          strcmp(&word[2][0],"E.S.D.") );


  /* if ESP Charges are present read them, if not, rewind
   * the file and return */
  if ( data->have_esp == TRUE ) 
  {
    /* reserve memory for dynamically allocated
     * Mulliken charge array */
    data->esp_charges = 
      (double *)calloc(numatoms,sizeof(double));

    /* check for success */
    if (data->esp_charges == NULL) 
    {
      PRINTERR; 
      return FALSE;
    }


    /* skip next line */
    eatline(data->file);

    for ( n = 0; n < numatoms; ++n) {

      /* make sure we don't get nuked by a bad file */
      ERR_FALSE( fgets(buffer,sizeof(buffer),data->file) )
      sscanf(buffer,"%s %lf ",&word[0][0], (data->esp_charges+n));
    }
  }

  /* no ESP charges found; in this case we don't know where
   * we are in the file; hence rewind file */
  else 
  {
    rewind(data->file); 
  }


  /* short message */ 
  printf("gamessplugin> Detected %d atoms in the input file. \n",
          data->numatoms);


  if ( data->runtyp == HESSIAN ) 
  {

#ifdef ANIMATE_MODE
    initialize_animated_mode(data);
#endif

    /* 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) == TRUE ) 
    {
      data->have_internals = TRUE;
    }
    else 
    {
      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) == TRUE ) 
    {
      data->have_cart_hessian = TRUE;
    }
    else 
    {
      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) == FALSE ) 
    {
      printf("gamessplugin> \n");
      printf("gamessplugin> Could not scan the normal modes,\n");
      printf("gamessplugin> \n");

      have_normal_modes = FALSE;
    }

#ifdef ANIMATE_MODE
    /* generate an animated trajectory for mode i */ 
    if (data->have_cart_hessian && have_normal_modes) 
    {
      if ( animate_normal_mode(data, 62) == FALSE )
      {
        printf("gamessplugin> \n");
        printf("gamessplugin> Error generating animated normal mode");
        printf("gamessplugin> \n \n");
      }
    }
#endif

  }

  /* done with this one */
  return TRUE; 
}



/******************************************************
 *
 * this function generates animated frames for normal
 * mode mode_to_animate
 *
 * *****************************************************/
static int initialize_animated_mode(void *mydata)
{
  gamessdata *data = (gamessdata *)mydata;
  mode_data *animated_mode;

  /* allocate memory for a mode_data struct */
  animated_mode = (mode_data *)calloc(1,sizeof(mode_data)); 

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

  /* save pointer in gamessdata struct */
  data->animated_mode = animated_mode;

  /* set the number of frames per animated mode */
  animated_mode->mode_num_frames = 25;

  /* initialize the animated mode tracker */
  animated_mode->current_mode_frame = 0;

  /* set the scaling factor for the animated mode */
  animated_mode->mode_scaling = 0.01;

  /* reserve memory for animated mode trajectory */
  animated_mode->mode_frames = (double *)calloc((4 *
      animated_mode->mode_num_frames+3)*data->numatoms*3, 
      sizeof(double));

  if ( animated_mode->mode_frames == NULL )
  {
    PRINTERR;
    return FALSE;
  }

  return TRUE;
}



/************************************************************
 *
 * 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(void *mydata, unsigned int mode)
{
  gamessdata *data = (gamessdata *)mydata;
  mode_data *animated_mode = data->animated_mode;
  double *normal_modes = data->normal_modes;
  double scale = animated_mode->mode_scaling;
  unsigned int i = 0, k = 0; 
  int l = 0, m = 0;
  unsigned int natoms = data->numatoms;
  unsigned 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->temporary+i)->x * (1+( k*scale * 
          (*(normal_modes+(mode*natoms*3)+(3*i)))));

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

      *(animated_mode->mode_frames+(k*natoms*3)+(3*i+2)) = 
          (data->temporary+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->temporary+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->temporary+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->temporary+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->temporary+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->temporary+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->temporary+i)->z * (1+((int)(m-num_frames)*scale * 
          (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
    }
  }

  /* short message */
  printf("gamessplugin> Successfully animated mode %d \n", mode);

  return TRUE;
}



/***********************************************************
 *
 * this function reads in the wavenumbers and intensities of
 * the normal modes
 *
 * *********************************************************/
static int get_normal_modes(void *mydata)
{
  gamessdata *data = (gamessdata *)mydata;
  char word[4][BUFSIZ];
  char buffer[BUFSIZ];
  unsigned int i = 0, k = 0, j = 0;
  unsigned int remaining_columns;
  double entry[6]; 
  char separator;
  char *item;
  unsigned 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 = 
    (double *)calloc(data->numatoms*3,sizeof(double));

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


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

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


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

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


  /* look for FREQUENCYs and IR INTENSITIES */
  for ( i = 0; i < (unsigned int)(data->numatoms*3/5); ++i) 
  {
    do 
    {
      ERR_FALSE( fgets(buffer,sizeof(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 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]); 
    
      /* save 'em */
      for ( k = 0; k < 5; ++k) 
      {
        *(data->wavenumbers+(i*5)+k) = entry[k]; 
      }
    }


    /* skip next line */
    eatline(data->file);


    /* next line contains the IR INTENSITIES */
    ERR_FALSE( fgets(buffer,sizeof(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]);
 

    /* save 'em */
    for ( k = 0; k < 5; ++k) 
    {
      *(data->intensities+(i*5)+k) = entry[k]; 
    }

    /* skip the next line */
    eatline(data->file);

    /* read the following five modes */
    for ( k = 0; k < data->numatoms; ++k)
    {
      /* x */
      ERR_FALSE( fgets(buffer,sizeof(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]);

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


      /* y */
      ERR_FALSE( fgets(buffer,sizeof