/***************************************************************************
 *cr
 *cr            (C) Copyright 1995-2016 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.34 $       $Date: 2016/11/28 05:01:53 $
 ***************************************************************************/

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

/*! \fn static inline bool host_is_little_endian(void)
 * detect endiannes of host machine. returns true on little endian machines. */
static inline int host_is_little_endian(void) 
{
  const union { unsigned 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);

        // comment line?
        if (buf[0] == '#') return mdio_readline(mf,buf,n,strip);

	// 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 && (!(*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);

  return 0; // Done!
}


// 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];
  char xposc[12], yposc[12], zposc[12];
  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%8c%8c%8c", 
             ma->resid, ma->resname, ma->atomname, atomnum, 
             xposc, yposc, zposc);

  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;
  xposc[8] = 0;
  yposc[8] = 0;
  zposc[8] = 0;
 
  if ((sscanf(xposc, "%f", &ma->pos[0]) != 1) ||
      (sscanf(yposc, "%f", &ma->pos[1]) != 1) ||
      (sscanf(zposc, "%f", &ma->pos[2]) != 1)) {
    return mdio_seterror(MDIO_BADFORMAT);
  }

  // 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];
  char xposc[12], yposc[12], zposc[12];

  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%8c%8c%8c", xposc, yposc, zposc);
    if (n != 3) 
      return mdio_seterror(MDIO_BADFORMAT);

    if ((sscanf(xposc, "%f", &ts->pos[coord    ]) != 1) ||
        (sscanf(yposc, "%f", &ts->pos[coord + 1]) != 1) ||
        (sscanf(zposc, "%f", &ts->pos[coord + 2]) != 1)) {
      return mdio_seterror(MDIO_BADFORMAT);
    }

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

    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);
		return 0;
	}
}


// Reads in a timestep frame from the .trX file and returns the
// data in a timestep structure. Returns NULL on error.
static int trx_timestep(md_file *mf, md_ts *ts) {
	int i;
  float x[3], y[3], z[3];
	trx_hdr *hdr;

	if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
	if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
		return mdio_seterror(MDIO_WRONGFORMAT);

	// Read the header
	if (trx_header(mf) < 0) return -1;

	// We need some data from the trX header
	hdr = mf->trx;
	if (!hdr) return mdio_seterror(MDIO_BADPARAMS);

	if (hdr->box_size) { // XXX need to check value of box_size!!
		if (trx_rvector(mf, x) < 0) return -1;
		if (trx_rvector(mf, y) < 0) return -1;
		if (trx_rvector(mf, z) < 0) 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->box);
      ts->box = NULL;
      return -1;
    }
	}

	if (hdr->vir_size) {
		if (trx_rvector(mf, NULL) < 0) return -1;
		if (trx_rvector(mf, NULL) < 0) return -1;
		if (trx_rvector(mf, NULL) < 0) return -1;
	}

        if (hdr->pres_size) {
                if (trx_rvector(mf, NULL) < 0) return -1;
                if (trx_rvector(mf, NULL) < 0) return -1;
                if (trx_rvector(mf, NULL) < 0) return -1;
        }

        if (hdr->x_size) {
                ts->pos = (float *) malloc(sizeof(float) * 3 * hdr->natoms);
                if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);

		ts->natoms = hdr->natoms;

		for (i = 0; i < hdr->natoms; i++) {
			if (trx_rvector(mf, &ts->pos[i * 3]) < 0) {
				mdio_tsfree(ts, 1);
				return -1;
			}
			ts->pos[i * 3] *= ANGS_PER_NM;
			ts->pos[i * 3 + 1] *= ANGS_PER_NM;
			ts->pos[i * 3 + 2] *= ANGS_PER_NM;
		}
        }

        if (hdr->v_size) {
		for (i = 0; i < hdr->natoms; i++) {
			if (trx_rvector(mf, NULL) < 0) {
				mdio_tsfree(ts, 1);
				return -1;
			}
		}
        }

        if (hdr->f_size) {
		for (i = 0; i < hdr->natoms; i++) {
			if (trx_rvector(mf, NULL) < 0) {
				mdio_tsfree(ts, 1);
				return -1;
			}
		}
        }

	return mdio_seterror(MDIO_SUCCESS);
}


// writes an int in big endian. Returns GMX_SUCCESS
// on success or a negative number on error.
static int put_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 (mf->rev) swap4_aligned(&y, 1);
      if (fwrite(&y, 4, 1, mf->f) != 1)
    return mdio_seterror(MDIO_IOERROR);

  return mdio_seterror(MDIO_SUCCESS);
}

// writes a real in big-endian. Returns GMX_SUCCESS
// on success or a negative number on error.
static int put_trx_real(md_file *mf, float y) {
      if (!mf) return mdio_seterror(MDIO_BADPARAMS);

      if (mf->rev) swap4_aligned(&y, 1);
      if (fwrite(&y, 4, 1, mf->f) != 1)
        return mdio_seterror(MDIO_IOERROR);

      return mdio_seterror(MDIO_SUCCESS);
}


// writes an xdr encoded string. Returns GMX_SUCCESS
// on success or a negative number on error.
static int put_trx_string(md_file *mf, const char *s) {
	if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
        
        // write: size of object, string length, string data
        size_t len = strlen(s);
        if ( put_trx_int(mf, len+1)
             || put_trx_int(mf, len)
             || (fwrite(s, len, 1, mf->f) != 1))
          return mdio_seterror(MDIO_IOERROR);

	return mdio_seterror(MDIO_SUCCESS);
}


// xtc_int() - reads an integer from an xtc file
static int xtc_int(md_file *mf, int *i) {
	unsigned char c[4];

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

	if (fread(c, 1, 4, mf->f) != 4) {
		if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
		else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
		else return mdio_seterror(MDIO_UNKNOWNERROR);
	}

	if (i) *i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
	return mdio_seterror(MDIO_SUCCESS);
}


// xtc_float() - reads a float from an xtc file
static int xtc_float(md_file *mf, float *f) {
	unsigned char c[4];
	int i;

	if (!mf) return mdio_seterror(MDIO_BADPARAMS);

	if (fread(c, 1, 4, mf->f) != 4) {
		if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
		else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
		else return mdio_seterror(MDIO_UNKNOWNERROR);
	}

	if (f) {
		// By reading the number in as an integer and then
		// copying it to a floating point number we can
		// ensure proper endianness
		i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
		memcpy(f, &i, 4);
	}
	return mdio_seterror(MDIO_SUCCESS);
}


// xtc_data() - reads a specific amount of data from an xtc
// file using the xdr format.
static int xtc_data(md_file *mf, char *buf, int len) {
	if (!mf || len < 1) return mdio_seterror(MDIO_BADPARAMS);
  size_t slen = (size_t)len;
	if (buf) {
		if (fread(buf, 1, slen, mf->f) != slen) {
			if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
			if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
			else return mdio_seterror(MDIO_UNKNOWNERROR);
		}
		if (len % 4) {
			if (fseek(mf->f, 4 - (len % 4), SEEK_CUR)) {
				if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
				if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
				else return mdio_seterror(MDIO_UNKNOWNERROR);
			}
		}
	}
	else {
		int newlen;
		newlen = len;
		if (len % 4) newlen += (4 - (len % 4));
		if (fseek(mf->f, newlen, SEEK_CUR)) {
			if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
			if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
			else return mdio_seterror(MDIO_UNKNOWNERROR);
		}
	}
	return len;
}


// xtc_timestep() - reads a timestep from an .xtc file.
static int xtc_timestep(md_file *mf, md_ts *ts) {
	int n;
	float f, x[3], y[3], z[3];

	int size = 0; // explicitly initialized to zero.
	float precision;

	if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
	if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
	if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);

	// 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;
	ts->natoms = n;

	// Get the simulation step
	if (xtc_int(mf, &n) < 0) return -1;
	ts->step = n;

	// Get the time value
	if (xtc_float(mf, &f) < 0) return -1;
	ts->time = f;

	// Read the basis vectors of the box
  if ( (xtc_float(mf, &x[0]) < 0) ||
       (xtc_float(mf, &x[1]) < 0) ||
       (xtc_float(mf, &x[2]) < 0) ||
       (xtc_float(mf, &y[0]) < 0) ||
       (xtc_float(mf, &y[1]) < 0) ||
       (xtc_float(mf, &y[2]) < 0) ||
       (xtc_float(mf, &z[0]) < 0) ||
       (xtc_float(mf, &z[1]) < 0) ||
       (xtc_float(mf, &z[2]) < 0) )
    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->box);
    ts->box = NULL;
    return -1;
  }

	ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
	if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
	n = xtc_3dfcoord(mf, ts->pos, &size, &precision);
	if (n < 0) return -1;

	/* Now we're left with the job of scaling... */
	for (n = 0; n < ts->natoms * 3; n++)
		ts->pos[n] *= ANGS_PER_NM;

	return mdio_seterror(MDIO_SUCCESS);
}


///////////////////////////////////////////////////////////////////////
// This algorithm is an implementation of the 3dfcoord algorithm
// written by Frans van Hoesel (hoesel@chem.rug.nl) as part of the
// Europort project in 1995.
///////////////////////////////////////////////////////////////////////

// integer table used in decompression
static int xtc_magicints[] = {
	0, 0, 0, 0, 0, 0, 0, 0, 0,8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
	80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
	1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, 16384,
	20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031, 131072,
	165140, 208063, 262144, 330280, 416127, 524287, 660561, 832255,
	1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
	5284491, 6658042, 8388607, 10568983, 13316085, 16777216 };

#define FIRSTIDX 9
/* note that magicints[FIRSTIDX-1] == 0 */
#define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))


// returns the number of bits in the binary expansion of
// the given integer.
static int xtc_sizeofint(int size) {
	unsigned int num = 1;
  unsigned int ssize = (unsigned int)size;
	int nbits = 0;

	while (ssize >= num && nbits < 32) {
		nbits++;
		num <<= 1;
	}
	return nbits;
}

// calculates the number of bits a set of integers, when compressed,
// will take up.
static int xtc_sizeofints(int nints, unsigned int *sizes) {
	int i;
  unsigned int num;
	unsigned int nbytes, nbits, bytes[32], bytecnt, tmp;
	nbytes = 1;
	bytes[0] = 1;
	nbits = 0;
	for (i=0; i < nints; i++) {	
		tmp = 0;
		for (bytecnt = 0; bytecnt < nbytes; bytecnt++) {
			tmp = bytes[bytecnt] * sizes[i] + tmp;
			bytes[bytecnt] = tmp & 0xff;
			tmp >>= 8;
		}
		while (tmp != 0) {
			bytes[bytecnt++] = tmp & 0xff;
			tmp >>= 8;
		}
		nbytes = bytecnt;
	}
	num = 1;
	nbytes--;
	while (bytes[nbytes] >= num) {
		nbits++;
		num *= 2;
	}
	return nbits + nbytes * 8;
}

// reads bits from a buffer.    
static int xtc_receivebits(int *buf, int nbits) {
	int cnt, num; 
	unsigned int lastbits, lastbyte;
	unsigned char * cbuf;
	int mask = (1 << nbits) -1;

	cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
	cnt = buf[0];
	lastbits = (unsigned int) buf[1];
	lastbyte = (unsigned int) buf[2];

	num = 0;
	while (nbits >= 8) {
		lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
		num |=  (lastbyte >> lastbits) << (nbits - 8);
		nbits -=8;
	}
	if (nbits > 0) {
		if (lastbits < (unsigned int)nbits) {
			lastbits += 8;
			lastbyte = (lastbyte << 8) | cbuf[cnt++];
		}
		lastbits -= nbits;
		num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
	}
	num &= mask;
	buf[0] = cnt;
	buf[1] = lastbits;
	buf[2] = lastbyte;
	return num; 
}

// decompresses small integers from the buffer
// sizes parameter has to be non-zero to prevent divide-by-zero
static void xtc_receiveints(int *buf, const int nints, int nbits,
			unsigned int *sizes, int *nums) {
	int bytes[32];
	int i, j, nbytes, p, num;

	bytes[1] = bytes[2] = bytes[3] = 0;
	nbytes = 0;
	while (nbits > 8) {
		bytes[nbytes++] = xtc_receivebits(buf, 8);
		nbits -= 8;
	}
	if (nbits > 0) {
		bytes[nbytes++] = xtc_receivebits(buf, nbits);
	}
	for (i = nints-1; i > 0; i--) {
		num = 0;
		for (j = nbytes-1; j >=0; j--) {
			num = (num << 8) | bytes[j];
			p = num / sizes[i];
			bytes[j] = p;
			num = num - p * sizes[i];
		}
		nums[i] = num;
	}
	nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
}

// function that actually reads and writes compressed coordinates    
static int xtc_3dfcoord(md_file *mf, float *fp, int *size, float *precision) {
	static int *ip = NULL;
	static int oldsize;
	static int *buf;

	int minint[3], maxint[3], *lip;
	int smallidx;
	unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
	int flag, k;
	int small, smaller, i, is_smaller, run;
	float *lfp;
	int tmp, *thiscoord,  prevcoord[3];

	int bufsize, lsize;
	unsigned int bitsize;
	float inv_precision;

        /* avoid uninitialized data compiler warnings */
        bitsizeint[0] = 0;
        bitsizeint[1] = 0;
        bitsizeint[2] = 0;

	if (xtc_int(mf, &lsize) < 0) return -1;

	if (*size != 0 && lsize != *size) return mdio_seterror(MDIO_BADFORMAT);

	*size = lsize;
	size3 = *size * 3;
	if (*size <= 9) {
		for (i = 0; i < *size; i++) {
			if (xtc_float(mf, fp + (3 * i)) < 0) return -1;
			if (xtc_float(mf, fp + (3 * i) + 1) < 0) return -1;
			if (xtc_float(mf, fp + (3 * i) + 2) < 0) return -1;
		}
		return *size;
	}
	xtc_float(mf, precision);
	if (ip == NULL) {
		ip = (int *)malloc(size3 * sizeof(*ip));
		if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
		bufsize = (int) (size3 * 1.2);
		buf = (int *)malloc(bufsize * sizeof(*buf));
		if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
		oldsize = *size;
	} else if (*size > oldsize) {
		ip = (int *)realloc(ip, size3 * sizeof(*ip));
		if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
		bufsize = (int) (size3 * 1.2);
		buf = (int *)realloc(buf, bufsize * sizeof(*buf));
		if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
		oldsize = *size;
	}
	buf[0] = buf[1] = buf[2] = 0;

	xtc_int(mf, &(minint[0]));
	xtc_int(mf, &(minint[1]));
	xtc_int(mf, &(minint[2]));

	xtc_int(mf, &(maxint[0]));
	xtc_int(mf, &(maxint[1]));
	xtc_int(mf, &(maxint[2]));
		
	sizeint[0] = maxint[0] - minint[0]+1;
	sizeint[1] = maxint[1] - minint[1]+1;
	sizeint[2] = maxint[2] - minint[2]+1;
	
	/* check if one of the sizes is to big to be multiplied */
	if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
		bitsizeint[0] = xtc_sizeofint(sizeint[0]);
		bitsizeint[1] = xtc_sizeofint(sizeint[1]);
		bitsizeint[2] = xtc_sizeofint(sizeint[2]);
		bitsize = 0; /* flag the use of large sizes */
	} else {
		bitsize = xtc_sizeofints(3, sizeint);
	}

	xtc_int(mf, &smallidx);
	smaller = xtc_magicints[FIRSTIDX > smallidx - 1 ? FIRSTIDX : smallidx - 1] / 2;
	small = xtc_magicints[smallidx] / 2;
	sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx];

	/* check for zero values that would yield corrupted data */
	if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
		printf("XTC corrupted, sizesmall==0 (case 1)\n");
		return -1;
	}


	/* buf[0] holds the length in bytes */
	if (xtc_int(mf, &(buf[0])) < 0) return -1;

	if (xtc_data(mf, (char *) &buf[3], (int) buf[0]) < 0) return -1;

	buf[0] = buf[1] = buf[2] = 0;

	lfp = fp;
	inv_precision = 1.0f / (*precision);
	run = 0;
	i = 0;
	lip = ip;
	while (i < lsize) {
		thiscoord = (int *)(lip) + i * 3;

		if (bitsize == 0) {
			thiscoord[0] = xtc_receivebits(buf, bitsizeint[0]);
			thiscoord[1] = xtc_receivebits(buf, bitsizeint[1]);
			thiscoord[2] = xtc_receivebits(buf, bitsizeint[2]);
		} else {
			xtc_receiveints(buf, 3, bitsize, sizeint, thiscoord);
		}

		i++;
		thiscoord[0] += minint[0];
		thiscoord[1] += minint[1];
		thiscoord[2] += minint[2];

		prevcoord[0] = thiscoord[0];
		prevcoord[1] = thiscoord[1];
		prevcoord[2] = thiscoord[2];
 

		flag = xtc_receivebits(buf, 1);
		is_smaller = 0;
		if (flag == 1) {
			run = xtc_receivebits(buf, 5);
			is_smaller = run % 3;
			run -= is_smaller;
			is_smaller--;
		}
		if (run > 0) {
			thiscoord += 3;
			for (k = 0; k < run; k+=3) {
				xtc_receiveints(buf, 3, smallidx, sizesmall, thiscoord);
				i++;
				thiscoord[0] += prevcoord[0] - small;
				thiscoord[1] += prevcoord[1] - small;
				thiscoord[2] += prevcoord[2] - small;
				if (k == 0) {
					/* interchange first with second atom for better
					 * compression of water molecules
					 */
					tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
					prevcoord[0] = tmp;
					tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
					prevcoord[1] = tmp;
					tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
					prevcoord[2] = tmp;
					*lfp++ = prevcoord[0] * inv_precision;
					*lfp++ = prevcoord[1] * inv_precision;
					*lfp++ = prevcoord[2] * inv_precision;

					if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
						printf("XTC corrupted, sizesmall==0 (case 2)\n");
						return -1;
					}

				} else {
					prevcoord[0] = thiscoord[0];
					prevcoord[1] = thiscoord[1];
					prevcoord[2] = thiscoord[2];
				}
				*lfp++ = thiscoord[0] * inv_precision;
				*lfp++ = thiscoord[1] * inv_precision;
				*lfp++ = thiscoord[2] * inv_precision;
			}
		} else {
			*lfp++ = thiscoord[0] * inv_precision;
			*lfp++ = thiscoord[1] * inv_precision;
			*lfp++ = thiscoord[2] * inv_precision;		
		}
		smallidx += is_smaller;
		if (is_smaller < 0) {
			small = smaller;
			if (smallidx > FIRSTIDX) {
				smaller = xtc_magicints[smallidx - 1] /2;
			} else {
				smaller = 0;
			}
		} else if (is_smaller > 0) {
			smaller = small;
			small = xtc_magicints[smallidx] / 2;
		}
		sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx] ;
	}
	return 1;
}
#endif