---

Gromacs.h

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
 ***************************************************************************/
/***************************************************************************
 * RCS INFORMATION:
 *      $RCSfile: Gromacs.h,v $
 *      $Author: johns $       $Locker:  $             $State: Exp $
 *      $Revision: 1.24 $       $Date: 2006/08/12 22:30:52 $
 ***************************************************************************/

/*
 * GROMACS file format reader for VMD
 *
 * This code provides a high level I/O library for reading
 * and writing the following file formats:
 *      gro     GROMACS format or trajectory
 *      g96     GROMOS-96 format or trajectory
 *      trj     Trajectory - x, v and f (binary, full precision)
 *      trr     Trajectory - x, v and f (binary, full precision, portable)
 *      xtc     Trajectory - x only (compressed, portable, any precision)
 *      top
 * Currently supported: gro trj trr g96 [xtc]
 *
 * TODO list
 *   o  velocities are ignored because VMD doesn't use them, but some other 
 *      program might ...
 *   o  gro_rec() assumes positions in .gro files are nanometers and
 *      converts to angstroms, whereas they really could be any unit
 */

#ifndef GROMACS_H
#define GROMACS_H

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#if defined(_AIX)
#include <strings.h>
#endif

#include "endianswap.h"

#if defined(WIN32) || defined(WIN64)
#define strcasecmp stricmp
#endif

#ifndef M_PI_2
#define M_PI_2 1.57079632679489661922
#endif

// Error codes for mdio_errno
#define MDIO_SUCCESS            0
#define MDIO_BADFORMAT          1
#define MDIO_EOF                2
#define MDIO_BADPARAMS          3
#define MDIO_IOERROR            4
#define MDIO_BADPRECISION       5
#define MDIO_BADMALLOC          6
#define MDIO_CANTOPEN           7
#define MDIO_BADEXTENSION       8
#define MDIO_UNKNOWNFMT         9
#define MDIO_CANTCLOSE          10
#define MDIO_WRONGFORMAT        11
#define MDIO_SIZEERROR          12
#define MDIO_UNKNOWNERROR       1000

#define MDIO_READ       0
#define MDIO_WRITE      1

#define MDIO_MAX_ERRVAL         11

// Format extensions
const char *mdio_fmtexts[] = {
    "",
    ".gro",
    ".trr",
    ".g96",
    ".trj",
    ".xtc",
    NULL
};


static int mdio_errcode;        // Last error code

#define TRX_MAGIC       1993    // Magic number for .trX files
#define XTC_MAGIC       1995    // Magic number for .xtc files
#define MAX_GRO_LINE    500     // Maximum line length of .gro files
#define MAX_G96_LINE    500     // Maximum line length of .g96 files
#define MAX_TRX_TITLE   80      // Maximum length of a title in .trX
#define MAX_MDIO_TITLE  80      // Maximum supported title length
#define ANGS_PER_NM     10      // Unit conversion factor
#define ANGS2_PER_NM2   100     // Unit conversion factor


// All the supported file types and their respective extensions
#define MDFMT_GRO               1
#define MDFMT_TRR               2
#define MDFMT_G96               3
#define MDFMT_TRJ               4
#define MDFMT_XTC               5


// A structure to hold .trX file format header information. This
// is an optional member of the md_file structure that is used
// when .trX files are being dealt with.
typedef struct {
        int version;            // File version number
        char title[MAX_TRX_TITLE + 1];  // File title
        int ir_size;
        int e_size;
        int box_size;
        int vir_size;
        int pres_size;
        int top_size;
        int sym_size;
        int x_size;             // Positions of atoms
        int v_size;             // Velocities of atoms
        int f_size;
        int natoms;             // Number of atoms in the system
        int step;
        int nre;
        float t;
        float lambda;
} trx_hdr;


// A generic i/o structure that contains information about the
// file itself and the input/output state
typedef struct {
        FILE *  f;      // Pointer to the file
        int     fmt;    // The file format
        int     prec;   // Real number precision
        int     rev;    // Reverse endiannism?
        trx_hdr * trx;  // Trx files require a great deal more
                        // header data to be stored.
} md_file;


// A format-independent structure to hold header data from files
typedef struct {
        char title[MAX_MDIO_TITLE + 1];
        int natoms;
        float timeval;
} md_header;


// A format-independent structure to hold unit cell data
typedef struct {
  float A, B, C, alpha, beta, gamma;
} md_box;


// Timestep information
typedef struct {
        float *pos;     // Position array (3 * natoms)
        //float *vel;   // Velocity array ** (VMD doesn't use this) **
        //float *f;     // Force array ** (VMD doesn't use this) **
        //float *box;   // Computational box ** (VMD doesn't use this) **
        int natoms;     // Number of atoms
        int step;       // Simulation step
        float time;     // Time of simulation
  md_box *box;
} md_ts;


// Atom information
typedef struct {
        char resid[7];          // Residue index number
        char resname[7];        // Residue name
        int atomnum;            // Atom index number
        char atomname[7];       // Atom name
        float pos[3];           // Position array (3 * natoms)
        //float vel[3]; // Velocity array ** (VMD doesn't use this) **
} md_atom;


// Open a molecular dynamics file. The second parameter specifies
// the format of the file. If it is zero, the format is determined
// from the file extension. the third argument (if given) decides
// whether to read (==0) or to write (!= 0).
// using a default argument set to read for backward compatibility.
static md_file *mdio_open(const char *, const int, const int=MDIO_READ);

// Closes a molecular dynamics file.
static int mdio_close(md_file *);


// Format-independent file I/O routines
static int mdio_header(md_file *, md_header *);
static int mdio_timestep(md_file *, md_ts *);


// .gro file functions
static int gro_header(md_file *, char *, int, float *, int *, int = 1);
static int gro_rec(md_file *, md_atom *);
static int gro_timestep(md_file *, md_ts *);


// .trX file functions
static int trx_header(md_file *, int = 0);
static int trx_int(md_file *, int *);
static int trx_real(md_file *, float *);

static int trx_rvector(md_file *, float *);
static int trx_string(md_file *, char *, int);
static int trx_timestep(md_file *, md_ts *);

// .g96 file functions
static int g96_header(md_file *, char *, int, float *);
static int g96_timestep(md_file *, md_ts *);
static int g96_rec(md_file *, md_atom *);
static int g96_countatoms(md_file *);


// .xtc file functions
static int xtc_int(md_file *, int *);
static int xtc_float(md_file *, float *);
/* 
static int xtc_receivebits(int *, int);
static void xtc_receiveints(int *, int, int, const unsigned *, int *);
*/
static int xtc_timestep(md_file *, md_ts *);
static int xtc_3dfcoord(md_file *, float *, int *, float *);


// Error reporting functions
static int mdio_errno(void);
static const char *mdio_errmsg(int);
static int mdio_seterror(int);


// Miscellaneous functions
static int strip_white(char *);
static int mdio_readline(md_file *, char *, int, int = 1);
static int mdio_tsfree(md_ts *, int = 0);
static int mdio_readbox(md_box *, float *, float *, float *);



static int xtc_receivebits(int *, int);

// Error descriptions for mdio_errno
static const char *mdio_errdescs[] = {
        "no error",
        "file does not match format",
        "unexpected end-of-file reached",
        "function called with bad parameters",
        "file i/o error",
        "unsupported precision",
        "out of memory",
        "cannot open file",
        "bad file extension",
        "unknown file format",
        "cannot close file",
        "wrong file format for this function",
        "binary i/o error: sizeof(int) != 4",
        NULL
};

static inline int host_is_little_endian(void) 
{
  const union { char c[4]; unsigned int i; } 
  fixed = { { 0x10 , 0x20 , 0x40 , 0x80 } };
  const unsigned int i = 0x80402010U;
        
  if (fixed.i == i) {
    return 1;
  }
  return 0;
}



// Open a molecular dynamics file. The second parameter specifies
// the format of the file. If it is zero, the format is determined
// from the file extension.
md_file *mdio_open(const char *fn, const int fmt, const int rw) {
        md_file *mf;

        if (!fn) {
                mdio_seterror(MDIO_BADPARAMS);
                return NULL;
        }

        // Allocate memory
        mf = (md_file *) malloc(sizeof(md_file));
        if (!mf) {
                mdio_seterror(MDIO_BADMALLOC);
                return NULL;
        }

        // Zero out the structure
        memset(mf, 0, sizeof(md_file));

        // Determine the file type from the extension
        if (!fmt) {
                char *p;
                int n;

                // Seek to the extension part of the filename
                for (p = (char *) &fn[strlen(fn) - 1]; *p != '.' && p > fn; p--);
                if (p == fn) {
                        free(mf);
                        mdio_seterror(MDIO_BADEXTENSION);
                        return NULL;
                }

                // Check the extension against known extensions
                for (n = 1; mdio_fmtexts[n]; n++)
                        if (!strcasecmp(p, mdio_fmtexts[n])) break;

                // If !mdio_fmtexts[n], we failed (unknown ext)
                if (!mdio_fmtexts[n]) {
                        free(mf);
                        mdio_seterror(MDIO_UNKNOWNFMT);
                        return NULL;
                }

                // All set
                mf->fmt = n;
        }
        else {
                mf->fmt = fmt;
        }

        // Differentiate between binary and ascii files. Also,
        // .trX files need a header information structure allocated.
        switch (mf->fmt) {
    case MDFMT_GRO:
        case MDFMT_G96: /* fallthrough */
        if (rw) 
            mf->f = fopen(fn, "wt");
        else
            mf->f = fopen(fn, "rt");

                break;
        case MDFMT_TRR:
        case MDFMT_TRJ: /* fallthrough */
                // Allocate the trx header data struct
                mf->trx = (trx_hdr *) malloc(sizeof(trx_hdr));
                if (!mf->trx) {
                        free(mf);
                        mdio_seterror(MDIO_BADMALLOC);
                        return NULL;
                }
                memset(mf->trx, 0, sizeof(trx_hdr));
        case MDFMT_XTC:  /* fallthrough */
                // Finally, open the file
        if (rw)
            mf->f = fopen(fn, "wb");
        else
            mf->f = fopen(fn, "rb");

                break;
        default:
                free(mf);
                mdio_seterror(MDIO_UNKNOWNFMT);
                return NULL;
        }

        // Check for opening error
        if (!mf->f) {
                if (mf->trx) free(mf->trx);
                free(mf);
                mdio_seterror(MDIO_CANTOPEN);
                return NULL;
        }

        // File is opened, we're all set!
        mdio_seterror(MDIO_SUCCESS);
        return mf;
}


// Closes a molecular dynamics file.
static int mdio_close(md_file *mf) {
        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        if (fclose(mf->f) == EOF) return mdio_seterror(MDIO_CANTCLOSE);

        // Free the dynamically allocated memory
        if (mf->trx) free(mf->trx);
        free(mf);
        return mdio_seterror(MDIO_SUCCESS);
}


// Returns the last error code reported by any of the mdio functions
static int mdio_errno(void) {
        return mdio_errcode;
}


// Returns a textual message regarding an mdio error code
static const char *mdio_errmsg(int n) {
        if (n < 0 || n > MDIO_MAX_ERRVAL) return (char *) "unknown error";
        else return mdio_errdescs[n];
}


// Sets the error code and returns an appropriate return value
// for the calling function to return to its parent
static int mdio_seterror(int code) {
        mdio_errcode = code;
        return code ? -1 : 0;
}


// Reads a line from the text file, strips leading/trailing whitespace
// and newline, checks for errors, and returns the number of characters
// in the string on success or -1 on error.
static int mdio_readline(md_file *mf, char *buf, int n, int strip) {
        if (!buf || n < 1 || !mf) return mdio_seterror(MDIO_BADPARAMS);

        // Read the line
        fgets(buf, n, mf->f);

        // End of file reached?
        if (feof(mf->f)) return mdio_seterror(MDIO_EOF);

        // File I/O error?
        if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);

        // Strip whitespace
        if (strip) strip_white(buf);

        return strlen(buf);
}


// Strips leading and trailing whitespace from a string. Tabs,
// spaces, newlines and carriage returns are stripped. Example:
// "\n   hello\t \r" becomes "hello".
static int strip_white(char *buf) {
        int i, j, k;

        // Protect against NULL pointer
        if (!buf) return -1;
        if (!strlen(buf)) return -1;

        // Kill trailing whitespace first
        for (i = strlen(buf) - 1;
             buf[i] == ' ' || buf[i] == '\t' ||
             buf[i] == '\n' || buf[i] == '\r';
             i--)
                buf[i] = 0;

        // Skip past leading whitespace
        for (i = 0; buf[i] == ' ' || buf[i] == '\t' ||
             buf[i] == '\n' || buf[i] == '\r'; i++);
        if (i) {
                k = 0;
                for (j = i; buf[j]; j++)
                        buf[k++] = buf[j];
                buf[k] = 0;
        }

        return strlen(buf);
}


// Frees the memory allocated in a ts structure. The holderror
// parameter defaults to zero. Programs that are calling this
// function because of an error reported by another function should
// set holderror so that mdio_tsfree() does not overwrite the error
// code with mdio_seterror().
static int mdio_tsfree(md_ts *ts, int holderror) {
        if (!ts) {
                if (holderror) return -1;
                else return mdio_seterror(MDIO_BADPARAMS);
        }

        if (ts->pos && ts->natoms > 0) free(ts->pos);

  if (ts->box) free(ts->box);

        if (holderror) return -1;
        else return mdio_seterror(MDIO_SUCCESS);
}


// Converts box basis vectors to A, B, C, alpha, beta, and gamma.  
// Stores values in md_box struct, which should be allocated before calling
// this function.
static int mdio_readbox(md_box *box, float *x, float *y, float *z) {
  float A, B, C;

  if (!box) {
    return mdio_seterror(MDIO_BADPARAMS);
  }

  // A, B, C are the lengths of the x, y, z vectors, respectively
  A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] ) * ANGS_PER_NM;
  B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] ) * ANGS_PER_NM;
  C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] ) * ANGS_PER_NM;
  if ((A<=0) || (B<=0) || (C<=0)) {
    /* Use zero-length box size and set angles to 90. */
    box->A = box->B = box->C = 0;
    box->alpha = box->beta = box->gamma = 90;
  } else {
    box->A = A;
    box->B = B;
    box->C = C;
  
    // gamma, beta, alpha are the angles between the x & y, x & z, y & z
    // vectors, respectively
    box->gamma = acos( (x[0]*y[0]+x[1]*y[1]+x[2]*y[2])*ANGS2_PER_NM2/(A*B) ) * 90.0/M_PI_2;
    box->beta = acos( (x[0]*z[0]+x[1]*z[1]+x[2]*z[2])*ANGS2_PER_NM2/(A*C) ) * 90.0/M_PI_2;
    box->alpha = acos( (y[0]*z[0]+y[1]*z[1]+y[2]*z[2])*ANGS2_PER_NM2/(B*C) ) * 90.0/M_PI_2; 
  }
  return mdio_seterror(MDIO_SUCCESS);
}


// Reads the header of a file (format independent)
static int mdio_header(md_file *mf, md_header *mdh) {
        int n;
        if (!mf || !mdh) return mdio_seterror(MDIO_BADPARAMS);
        if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);

        switch (mf->fmt) {
        case MDFMT_GRO:
                if (gro_header(mf, mdh->title, MAX_MDIO_TITLE,
                &mdh->timeval, &mdh->natoms, 1) < 0)
                        return -1;
                return 0;

        case MDFMT_TRR: 
        case MDFMT_TRJ: /* fallthrough */
                if (trx_header(mf, 1) < 0) return -1;
                mdh->natoms = mf->trx->natoms;
                mdh->timeval = (float) mf->trx->t;
                strncpy(mdh->title, mf->trx->title, MAX_MDIO_TITLE);
                return 0;

        case MDFMT_G96:
                if (g96_header(mf, mdh->title, MAX_MDIO_TITLE,
                &mdh->timeval) < 0) return -1;
                mdh->natoms = -1;
                return 0;

        case MDFMT_XTC:
                memset(mdh, 0, sizeof(md_header));
                // Check magic number
                if (xtc_int(mf, &n) < 0) return -1;
                if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);

                // Get number of atoms
                if (xtc_int(mf, &n) < 0) return -1;
                mdh->natoms = n;
                rewind(mf->f);
                return 0;

        default:
                return mdio_seterror(MDIO_UNKNOWNFMT);
        }
}


// Reads in a timestep from a file (format independent)
static int mdio_timestep(md_file *mf, md_ts *ts) {
        if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
        if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);

        switch (mf->fmt) {
        case MDFMT_GRO:
                return gro_timestep(mf, ts);

        case MDFMT_TRR:
        case MDFMT_TRJ: /* fallthrough */
                return trx_timestep(mf, ts);

        case MDFMT_G96:
                return g96_timestep(mf, ts);

        case MDFMT_XTC:
                return xtc_timestep(mf, ts);

        default:
                return mdio_seterror(MDIO_UNKNOWNFMT);
        }
}



static int g96_header(md_file *mf, char *title, int titlelen, float *timeval) {
        char buf[MAX_G96_LINE + 1];
        char *p;

        // Check parameters
        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        // The header consists of blocks. The title block
        // is mandatory, and a TIMESTEP block is optional.
        // Example:
        //
        // TITLE
        // Generated by trjconv :  t=  90.00000
        // more title info
        // .
        // .
        // .
        // END
        // .
        // .
        // .

        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
        if (strcasecmp(buf, "TITLE")) return mdio_seterror(MDIO_BADFORMAT);

        // Read in the title itself
        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;

        // The timevalue can be included in the title string
        // after a "t=" prefix.
        if ((p = (char *) strstr(buf, "t="))) {
                char *q = p;
                *(q--) = 0;

                // Skip the `t=' and strip whitespace from
                // the resulting strings
                p += 2;
                strip_white(p);
                strip_white(buf);

                // Grab the timevalue from the title string
                if (timeval) *timeval = (float) atof(p);
        }
        else {
                // No timevalue - just copy the string and strip
                // any leading/trailing whitespace
                if (timeval) *timeval = 0;
                strip_white(buf);
        }

        // Copy the title string
        if (title && titlelen) strncpy(title, buf, titlelen);

        // Now ignore subsequent title lines and get the END string
        while (strcasecmp(buf, "END"))
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;

        // Done!
        return mdio_seterror(MDIO_SUCCESS);
}


// Used to determine the number of atoms in a g96 file, because
// VMD needs to know this for some reason.
static int g96_countatoms(md_file *mf) {
        char buf[MAX_G96_LINE + 1];
        int natoms;
        int n;
        long fpos;
        float lastf;

        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        fpos = ftell(mf->f);

        natoms = 0;
        for (;;) {
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
                        break;
                n = sscanf(buf, "%*6c%*6c%*6c%*6c %*f %*f %f", &lastf);
                if (n == 1) natoms++;
                else {
                        strip_white(buf);
                        if (!strcasecmp(buf, "END")) break;
                }
        }

        fseek(mf->f, fpos, SEEK_SET);
        return natoms;
}


// Reads an atom line from the G96 file
static int g96_rec(md_file *mf, md_atom *ma) {
        char buf[MAX_G96_LINE + 1];
        char atomnum[7];
        int n;

        // Check parameters
        if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);

        // Read in a line, assuming it is an atom line
        do {
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0) return -1;
        } while (buf[0] == '#' || strlen(buf) == 0);

        n = sscanf(buf, "%6c%6c%6c%6c %f %f %f",
                ma->resid, ma->resname, ma->atomname, atomnum,
                &ma->pos[0], &ma->pos[1], &ma->pos[2]);
        if (n == 7) {
                atomnum[6] = 0;
                ma->resid[6] = 0;
                ma->resname[6] = 0;
                ma->atomname[6] = 0;

                strip_white(atomnum);
                strip_white(ma->resid);
                strip_white(ma->resname);
                strip_white(ma->atomname);

                ma->atomnum = atoi(atomnum);

                ma->pos[0] *= ANGS_PER_NM;
                ma->pos[1] *= ANGS_PER_NM;
                ma->pos[2] *= ANGS_PER_NM;

                return 0;
        }

        return mdio_seterror(MDIO_BADFORMAT);
}


// Reads a timestep from a G96 file and stores the data in
// the generic md_ts structure. Returns 0 on success or a
// negative number on error and sets mdio_errcode.
static int g96_timestep(md_file *mf, md_ts *ts) {
        char            buf[MAX_G96_LINE + 1];
        char            stripbuf[MAX_G96_LINE + 1];
        float           pos[3], x[3], y[3], z[3], *currAtom;
        long            fpos;
        int             n, i, boxItems;

        // Check parameters
        if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);

  // Allocate data space for the timestep, using the number of atoms
  // determined by open_g96_read().
        ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
        if (!ts->pos) {
                return mdio_seterror(MDIO_BADMALLOC);
        }
  currAtom = ts->pos;

        // The timesteps follow the header in a fixed block
        // format:
        //
        // TIMESTEP
        //         <step number> <time value>
        // END
        // POSITIONRED
        //     <x float> <y float> <z float>
        //     .         .         .
        //     .         .         .
        //     .         .         .
        // END
        // VELOCITYRED
        //     <x float> <y float> <z float>
        //     .         .         .
        //     .         .         .
        //     .         .         .
        // END
        // BOX
        //     <x float> <y float> <z float>
        // END
        //
        // -----
        //
        // The TIMESTEP, VELOCITY and BOX blocks are optional.
        // Floats are written in 15.9 precision.
        //
        // Reference: GROMACS 2.0 user manual
        //            http://rugmd4.chem.rug.nl/~gmx/online2.0/g96.html

        // First, look for an (optional) title block and skip it
        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;

  if (!strcasecmp(buf, "TITLE")) {
    // skip over the text until we reach 'END'
    while (strcasecmp(buf, "END")) {
      if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
    }

    // Read in the next line
    if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
  }

        // Next, look for a timestep block
        if (!strcasecmp(buf, "TIMESTEP")) {
                // Read in the value line
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;

                // Extract the time value and the timestep index
                n = sscanf(buf, "%d %f", &ts->step, &ts->time);
                if (n != 2) return mdio_seterror(MDIO_BADFORMAT);

                // Read the "END" line
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
                if (strcasecmp(buf, "END"))
                        return mdio_seterror(MDIO_BADFORMAT);

                // Read in the next line
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
        }
        else {
                // No timestep specified -- set to zero
                ts->step = 0;
                ts->time = 0;
        }

        // At this point a POSITION or POSITIONRED block
        // is REQUIRED by the format
        if (!strcasecmp(buf, "POSITIONRED")) {

    // So now we read in some atoms
    i = 0;
                while (i < ts->natoms) {
                        // Read in an atom
                        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
                                return -1;
 
      // We shouldn't reach the end yet
      if (!strcasecmp(buf, "END"))
        return mdio_seterror(MDIO_BADFORMAT);

                        // Get the x,y,z coordinates
                        n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
      
      // Ignore improperly formatted lines
                        if (n == 3) {
                                pos[0] *= ANGS_PER_NM;
                                pos[1] *= ANGS_PER_NM;
                                pos[2] *= ANGS_PER_NM;

                                // Copy the atom data into the array
                                memcpy(currAtom, pos, sizeof(float) * 3);
        currAtom += 3;
        i++;
                        }
                }
        }
        else if (!strcasecmp(buf, "POSITION") || !strcasecmp(buf, "REFPOSITION")) {
                /*
                char resnum[7];
                char resname[7];
                char atomname[7];
                char atomnum[7];
                */

                // So now we read in some atoms
    i = 0;
                while (i < ts->natoms) {
                        // Read in the first line
                        if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
                                return -1;
 
      // We shouldn't reach the end yet
      strcpy(stripbuf, buf);
      strip_white(stripbuf); 
      if (!strcasecmp(stripbuf, "END"))
        return mdio_seterror(MDIO_BADFORMAT);

                        // Get the x,y,z coordinates and name data
                        n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
                                &pos[0], &pos[1], &pos[2]);

      // Ignore improperly formatted lines
                        if (n == 3) {
                                pos[0] *= ANGS_PER_NM;
                                pos[1] *= ANGS_PER_NM;
                                pos[2] *= ANGS_PER_NM;

                                // Copy the atom data into the linked list item
                                memcpy(currAtom, pos, sizeof(float) * 3);
                                currAtom += 3;
        i++;
                        }
                }
        }
        else {
                return mdio_seterror(MDIO_BADFORMAT);
        }

  // Read the END keyword
  if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
    return -1;
  if (strcasecmp(buf, "END"))
    return mdio_seterror(MDIO_BADFORMAT);

        // ... another problem: there may or may not be a VELOCITY
        // block or a BOX block, so we need to read one line beyond
        // the POSITION block to determine this. If neither VEL. nor
        // BOX are present we've read a line too far and infringed
        // on the next timestep, so we need to keep track of the
        // position now for a possible fseek() later to backtrack.
        fpos = ftell(mf->f);

        // Now we must read in the velocities and the box, if present
        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
    // It's okay if we end the file here; any other errors need to be
    // reported.
    if (mdio_errcode == MDIO_EOF) 
      return mdio_seterror(MDIO_SUCCESS);
    else 
      return -1;
  }

        // Is there a velocity block present ?
        if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
                // Ignore all the coordinates - VMD doesn't use them
                for (;;) {
                        if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
                                return -1;
                        if (!strcasecmp(buf, "END")) break;
                }

                // Again, record our position because we may need
                // to fseek here later if we read too far.
                fpos = ftell(mf->f);

                // Go ahead and read the next line.
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
        }

        // Is there a box present ?
        if (!strcasecmp(buf, "BOX")) {
                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
    boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
               &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
    if (boxItems == 3) {
      x[1] = x[2] = 0;
      y[0] = y[2] = 0;
      z[0] = z[1] = 0;
    }
    else if (boxItems != 9) 
      return mdio_seterror(MDIO_BADFORMAT);

    // Allocate the box and convert the vectors.
    ts->box = (md_box *) malloc(sizeof(md_box));
    if (mdio_readbox(ts->box, x, y, z) < 0) {
      free(ts->box);
      ts->box = NULL;
      return mdio_seterror(MDIO_BADFORMAT);
    }

                if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
      free(ts->box);
      ts->box = NULL;
      return -1;
    }
                if (strcasecmp(buf, "END")) {
      free(ts->box);
      ts->box = NULL;
                        return mdio_seterror(MDIO_BADFORMAT);
    }
        }
        else {
                // We have read too far, so fseek back to the
                // last known safe position so we don't return
                // with the file pointer set infringing on the
                // next timestep data.
                fseek(mf->f, fpos, SEEK_SET);
        }

        // We're done!
        return mdio_seterror(MDIO_SUCCESS);
}


// Attempts to read header data from a GROMACS structure file
// The GROMACS header format is as follows (fixed, 2 lines ASCII):
// <title> [ n= <timevalue> ]
//     <num atoms>
static int gro_header(md_file *mf, char *title, int titlelen, float *timeval,
               int *natoms, int rewind) {
        char buf[MAX_GRO_LINE + 1];
        long fpos;
        char *p;

        // Check parameters
        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        // Get the current file position for rewinding later
        fpos = ftell(mf->f);

        // The header consists of 2 lines - get the first line
        if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;

        // The timevalue can be included in the title string
        // after a "t=" prefix.
        if ((p = (char *) strstr(buf, "t="))) {
                char *q = p;
                *(q--) = 0;

                // Skip the `t=' and strip whitespace from
                // the resulting strings
                p += 2;
                strip_white(p);
                strip_white(buf);

                // Grab the timevalue from the title string
                if (timeval) *timeval = (float) atof(p);
        }
        else {
                // No timevalue - just copy the string
                if (timeval) *timeval = 0;
        }

        // Copy the title string
        if (title && titlelen) strncpy(title, buf, titlelen);

        // Get the second line and grab the number of atoms
        if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;

        // Store the number of atoms
        if (natoms) if (!(*natoms = atoi(buf)))
                return mdio_seterror(MDIO_BADFORMAT);

        // Now we rewind the file so that subsequent calls to
        // gro_timestep() will succeed. gro_timestep() requires
        // the header to be at the current file pointer.
        if (rewind) fseek(mf->f, fpos, SEEK_SET);

        // Done!
        return 0;
}


// Reads one atom record from a GROMACS file. Returns GMX_SUCCESS
// on success or a negative number on error.
//
// Record format (one line, fixed):
//    rrrrrRRRRRaaaaaAAAAA <x pos> <y pos> <z pos> <x vel> <y vel> <z vel>
//
//    r = residue number
//    R = residue name
//    a = atom name
//    A = atom number
//
static int gro_rec(md_file *mf, md_atom *ma) {
        char    buf[MAX_GRO_LINE + 1];
        char    atomnum[6];
        int     n;

        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        do {
                if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) return -1;
        } while (buf[0] == '#' || !strlen(buf));

        // Read in the fields
        n = sscanf(buf, "%5c%5c%5c%5c%f %f %f", ma->resid,
                ma->resname, ma->atomname, atomnum, ma->pos,
                &ma->pos[1], &ma->pos[2]);
        if (n != 7) return mdio_seterror(MDIO_BADFORMAT);

        // Null terminate the strings
        ma->resname[5] = 0;
        ma->atomname[5] = 0;
        ma->resid[5] = 0;
        atomnum[5] = 0;

        // Convert strings to numbers
        strip_white(atomnum);
        ma->atomnum = atoi(atomnum);

        // Convert nanometers to angstroms
        ma->pos[0] *= ANGS_PER_NM;
        ma->pos[1] *= ANGS_PER_NM;
        ma->pos[2] *= ANGS_PER_NM;

        // Strip leading and trailing whitespace
        strip_white(ma->atomname);
        strip_white(ma->resname);
        strip_white(ma->resid);

        return 0;
}


// Reads in a timestep from a .gro file. Ignores the data
// not needed for a timestep, so is a little faster than
// calling gro_rec() for each atom. Also reads in the
// header block.
//
static int gro_timestep(md_file *mf, md_ts *ts) {
        char buf[MAX_GRO_LINE + 1];
        long coord;
        int i, n, boxItems;
  float x[3], y[3], z[3];

        if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);

        if (gro_header(mf, NULL, 0, &ts->time, &ts->natoms, 0) < 0)
                return -1;
        ts->pos = (float *) malloc(3 * sizeof(float) * ts->natoms);
        if (!ts->pos)
                return mdio_seterror(MDIO_BADMALLOC);

        coord = 0;
        for (i = 0; i < ts->natoms; i++) {
                if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
                        free(ts->pos);
                        return -1;
                }
        
                n = sscanf(buf, "%*5c%*5c%*5c%*5c%f %f %f",
                        &ts->pos[coord], &ts->pos[coord + 1],
                        &ts->pos[coord + 2]);

                ts->pos[coord] *= ANGS_PER_NM;
                ts->pos[coord + 1] *= ANGS_PER_NM;
                ts->pos[coord + 2] *= ANGS_PER_NM;

                if (n != 3) return mdio_seterror(MDIO_BADFORMAT);
                coord += 3;
        }

        // Read the box, stored as three vectors representing its edges
        if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
                free(ts->pos);
                return -1;
        }
  boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
             &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
  // File may only include three scalars for the box information -- if
  // that's the case, the box is orthoganal.
  if (boxItems == 3) {
    x[1] = x[2] = 0;
    y[0] = y[2] = 0;
    z[0] = z[1] = 0;
  }
  else if (boxItems != 9) {
    free(ts->pos);
    return -1;
  }

  // Allocate the box and convert the vectors.
  ts->box = (md_box *) malloc(sizeof(md_box));
  if (mdio_readbox(ts->box, x, y, z) < 0) {
    free(ts->pos);
    free(ts->box);
    ts->box = NULL;
    return -1;
  }

        return 0;
}


// Attempts to read header data from a .trX trajectory file
//
// The .trX header format is as follows:
//
//      4 bytes         - magic number (0x07C9)
//      ...
//
static int trx_header(md_file *mf, int rewind) {
        int magic;
        trx_hdr *hdr;
        long fpos;

        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        // In case we need to rewind
        fpos = ftell(mf->f);

        // We need to store some data to the trX header data
        // structure inside the md_file structure
        hdr = mf->trx;
        if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);

        // Read the magic number
        if (trx_int(mf, &magic) < 0) return -1;
        if (magic != TRX_MAGIC) {
                // Try reverse endianism
                swap4_aligned(&magic, 1);
                if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);

                // Enable byte swapping (actually works, too!)
                mf->rev = 1;
        }

        // Read the version number. 
        // XXX. this is not the version number, but the storage size
        // of the following XDR encoded string.
        // the 'title' string is in fact the version identifier.
        // since VMD does not use any of that, it does no harm,
        // but is should still be fixed occasionally. AK 2005/01/08.

        if(mf->fmt!=MDFMT_TRJ) {
                // It appears that TRJ files either don't contain a version
                // number or don't have a length-delimiter on the string,
                // whereas TRR files do contain both.  Thus, with TRJ, we just
                // assume that the version number is the string length and 
                // just hope for the best. -- WLD 2006/07/09
                if (trx_int(mf, &hdr->version) < 0) return -1;
        }

        // Read in the title string
        if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
                return -1;

        // Read in some size data
        if (trx_int(mf, &hdr->ir_size) < 0) return -1;
        if (trx_int(mf, &hdr->e_size) < 0) return -1;
        if (trx_int(mf, &hdr->box_size) < 0) return -1;
        if (trx_int(mf, &hdr->vir_size) < 0) return -1;
        if (trx_int(mf, &hdr->pres_size) < 0) return -1;
        if (trx_int(mf, &hdr->top_size) < 0) return -1;
        if (trx_int(mf, &hdr->sym_size) < 0) return -1;
        if (trx_int(mf, &hdr->x_size) < 0) return -1;
        if (trx_int(mf, &hdr->v_size) < 0) return -1;
        if (trx_int(mf, &hdr->f_size) < 0) return -1;
        if (trx_int(mf, &hdr->natoms) < 0) return -1;
        if (trx_int(mf, &hdr->step) < 0) return -1;
        if (trx_int(mf, &hdr->nre) < 0) return -1;

        // Make sure there are atoms...
        if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);

        // Try to determine precision (float? double?)
        if (hdr->x_size) mf->prec = hdr->x_size / (hdr->natoms * 3);
        else if (hdr->v_size) mf->prec = hdr->v_size / (hdr->natoms * 3);
        else if (hdr->f_size) mf->prec = hdr->f_size / (hdr->natoms * 3);
        else return mdio_seterror(MDIO_BADPRECISION);

        if (mf->prec != sizeof(float) && mf->prec != sizeof(double)) {
                // We have no data types this size! The
                // file must've been generated on another
                // platform
                return mdio_seterror(MDIO_BADPRECISION);
        }

        // Read in timestep and lambda
        if (trx_real(mf, &hdr->t) < 0) return -1;
        if (trx_real(mf, &hdr->lambda) < 0) return -1;

        // Rewind if necessary
        if (rewind) fseek(mf->f, fpos, SEEK_SET);

        return 0;
}


// Reads in an integer and stores it in y. Returns GMX_SUCCESS
// on success or a negative number on error.
static int trx_int(md_file *mf, int *y) {
        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        // sanity check.
        if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);

        if (y) {
                if (fread(y, 4, 1, mf->f) != 1)
                        return mdio_seterror(MDIO_IOERROR);
                if (mf->rev) swap4_aligned(y, 1);
        }
        else if (fseek(mf->f, 4, SEEK_CUR) != 0)
                return mdio_seterror(MDIO_IOERROR);

        return mdio_seterror(MDIO_SUCCESS);
}


// Reads in either a float or a double, depending on the
// precision, and stores that number in y. Returns
// GMX_SUCCESS on success or a negative number on error.
static int trx_real(md_file *mf, float *y) {
        double x;

        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        switch (mf->prec) {
                case sizeof(float):
                        if (!y) {
                                if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
                                        return mdio_seterror(MDIO_IOERROR);
                        } else {
                                if (fread(y, mf->prec, 1, mf->f) != 1)
                                        return mdio_seterror(MDIO_IOERROR);
                                if (mf->rev) swap4_aligned(y, 1);
                        }
                        return mdio_seterror(MDIO_SUCCESS);

                case sizeof(double):
                        if (!y) {
                                if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
                                        return mdio_seterror(MDIO_IOERROR);
                        } else {
                                if (fread(&x, mf->prec, 1, mf->f) != 1)
                                        return mdio_seterror(MDIO_IOERROR);
                                if (mf->rev) swap8_aligned(&x, 1);
                                *y = (float) x;
                        }
                        return mdio_seterror(MDIO_SUCCESS);

                default:
                        return mdio_seterror(MDIO_BADPRECISION);
        }

}


// Reads in a real-valued vector (taking precision into account).
// Stores the vector in vec, and returns GMX_SUCCESS on success
// or a negative number on error.
static int trx_rvector(md_file *mf, float *vec) {
        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        if (!vec) {
                if (trx_real(mf, NULL) < 0) return -1;
                if (trx_real(mf, NULL) < 0) return -1;
                if (trx_real(mf, NULL) < 0) return -1;
                return mdio_seterror(MDIO_SUCCESS);
        } else {
                if (trx_real(mf, &vec[0]) < 0) return -1;
                if (trx_real(mf, &vec[1]) < 0) return -1;
                if (trx_real(mf, &vec[2]) < 0) return -1;
                return mdio_seterror(MDIO_SUCCESS);
        }
}


// Reads in a string by first reading an integer containing the
// string's length, then reading in the string itself and storing
// it in str. If the length is greater than max, it is truncated
// and the rest of the string is skipped in the file. Returns the
// length of the string on success or a negative number on error.
static int trx_string(md_file *mf, char *str, int max) {
        int size;
  size_t ssize;

        if (!mf) return mdio_seterror(MDIO_BADPARAMS);

        if (trx_int(mf, &size) < 0) return -1;
  ssize = (size_t)size;

        if (str && size <= max) {
                if (fread(str, 1, size, mf->f) != ssize)
                        return mdio_seterror(MDIO_IOERROR);
                str[size] = 0;
                return size;
        } else if (str) {
                if (fread(str, 1, max, mf->f) != ssize)
                        return mdio_seterror(MDIO_IOERROR);
                if (fseek(mf->f, size - max, SEEK_CUR) != 0)
                        return mdio_seterror(MDIO_IOERROR);
                str[max] = 0;
                return max;
        } else {
                if (fseek(mf->f, size, SEEK_CUR) != 0)
                        return mdio_seterror(MDIO_IOERROR);