Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

Gromacs.h

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 /***************************************************************************
00009  * RCS INFORMATION:
00010  *      $RCSfile: Gromacs.h,v $
00011  *      $Author: dgomes $       $Locker:  $             $State: Exp $
00012  *      $Revision: 1.38 $       $Date: 2024/05/16 14:03:31 $
00013  ***************************************************************************/
00014 
00015 /*
00016  * GROMACS file format reader for VMD
00017  *
00018  * This code provides a high level I/O library for reading
00019  * and writing the following file formats:
00020  *      gro     GROMACS format or trajectory
00021  *      g96     GROMOS-96 format or trajectory
00022  *      trj     Trajectory - x, v and f (binary, full precision)
00023  *      trr     Trajectory - x, v and f (binary, full precision, portable)
00024  *      xtc     Trajectory - x only (compressed, portable, any precision)
00025  *      top
00026  * Currently supported: gro trj trr g96 [xtc]
00027  *
00028  * TODO list
00029  *   o  velocities are ignored because VMD doesn't use them, but some other 
00030  *      program might ...
00031  *   o  gro_rec() assumes positions in .gro files are nanometers and
00032  *      converts to angstroms, whereas they really could be any unit
00033  */
00034 
00035 #ifndef GROMACS_H
00036 #define GROMACS_H
00037 
00038 #include <math.h>
00039 #include <stdio.h>
00040 #include <stdlib.h>
00041 #include <string.h>
00042 #include <ctype.h>
00043 
00044 #if defined(_AIX)
00045 #include <strings.h>
00046 #endif
00047 
00048 #include "endianswap.h"
00049 
00050 #if defined(WIN32) || defined(WIN64)
00051 #define strcasecmp stricmp
00052 #endif
00053 
00054 #ifndef M_PI_2
00055 #define M_PI_2 1.57079632679489661922
00056 #endif
00057 
00058 // Error codes for mdio_errno
00059 #define MDIO_SUCCESS            0
00060 #define MDIO_BADFORMAT          1
00061 #define MDIO_EOF                2
00062 #define MDIO_BADPARAMS          3
00063 #define MDIO_IOERROR            4
00064 #define MDIO_BADPRECISION       5
00065 #define MDIO_BADMALLOC          6
00066 #define MDIO_CANTOPEN           7
00067 #define MDIO_BADEXTENSION       8
00068 #define MDIO_UNKNOWNFMT         9
00069 #define MDIO_CANTCLOSE          10
00070 #define MDIO_WRONGFORMAT        11
00071 #define MDIO_SIZEERROR          12
00072 #define MDIO_UNKNOWNERROR       1000
00073 
00074 #define MDIO_READ       0
00075 #define MDIO_WRITE      1
00076 
00077 #define MDIO_MAX_ERRVAL         11
00078 
00079 // Format extensions
00080 static const char *mdio_fmtexts[] = {
00081     "",
00082     ".gro",
00083     ".trr",
00084     ".g96",
00085     ".trj",
00086     ".xtc",
00087     NULL
00088 };
00089 
00090 
00091 static int mdio_errcode;        // Last error code
00092 
00093 #define TRX_MAGIC       1993    // Magic number for .trX files
00094 #define XTC_MAGIC       1995    // Magic number for .xtc files
00095 #define MAX_GRO_LINE    500     // Maximum line length of .gro files
00096 #define MAX_G96_LINE    500     // Maximum line length of .g96 files
00097 #define MAX_TRX_TITLE   80      // Maximum length of a title in .trX
00098 #define MAX_MDIO_TITLE  80      // Maximum supported title length
00099 #define ANGS_PER_NM     10      // Unit conversion factor
00100 #define ANGS2_PER_NM2   100     // Unit conversion factor
00101 
00102 
00103 // All the supported file types and their respective extensions
00104 #define MDFMT_GRO               1
00105 #define MDFMT_TRR               2
00106 #define MDFMT_G96               3
00107 #define MDFMT_TRJ               4
00108 #define MDFMT_XTC               5
00109 
00110 
00111 // A structure to hold .trX file format header information. This
00112 // is an optional member of the md_file structure that is used
00113 // when .trX files are being dealt with.
00114 typedef struct {
00115         int version;            // File version number
00116         char title[MAX_TRX_TITLE + 1];  // File title
00117         int ir_size;
00118         int e_size;
00119         int box_size;
00120         int vir_size;
00121         int pres_size;
00122         int top_size;
00123         int sym_size;
00124         int x_size;             // Positions of atoms
00125         int v_size;             // Velocities of atoms
00126         int f_size;
00127         int natoms;             // Number of atoms in the system
00128         int step;
00129         int nre;
00130         float t;
00131         float lambda;
00132 } trx_hdr;
00133 
00134 
00135 // A generic i/o structure that contains information about the
00136 // file itself and the input/output state
00137 typedef struct {
00138         FILE *  f;      // Pointer to the file
00139         int     fmt;    // The file format
00140         int     prec;   // Real number precision
00141         int     rev;    // Reverse endiannism?
00142         trx_hdr * trx;  // Trx files require a great deal more
00143                         // header data to be stored.
00144 } md_file;
00145 
00146 
00147 // A format-independent structure to hold header data from files
00148 typedef struct {
00149         char title[MAX_MDIO_TITLE + 1];
00150         int natoms;
00151         float timeval;
00152 } md_header;
00153 
00154 
00155 // A format-independent structure to hold unit cell data
00156 typedef struct {
00157   float A, B, C, alpha, beta, gamma;
00158 } md_box;
00159 
00160 
00161 // Timestep information
00162 typedef struct {
00163         float *pos;     // Position array (3 * natoms)
00164         //float *vel;   // Velocity array ** (VMD doesn't use this) **
00165         //float *f;     // Force array ** (VMD doesn't use this) **
00166         //float *box;   // Computational box ** (VMD doesn't use this) **
00167         int natoms;     // Number of atoms
00168         int step;       // Simulation step
00169         float time;     // Time of simulation
00170   md_box *box;
00171 } md_ts;
00172 
00173 
00174 // Atom information
00175 typedef struct {
00176         char resid[7];          // Residue index number
00177         char resname[7];        // Residue name
00178         int atomnum;            // Atom index number
00179         char atomname[7];       // Atom name
00180         float pos[3];           // Position array (3 * natoms)
00181         //float vel[3]; // Velocity array ** (VMD doesn't use this) **
00182 } md_atom;
00183 
00184 
00185 // Open a molecular dynamics file. The second parameter specifies
00186 // the format of the file. If it is zero, the format is determined
00187 // from the file extension. the third argument (if given) decides
00188 // whether to read (==0) or to write (!= 0).
00189 // using a default argument set to read for backward compatibility.
00190 static md_file *mdio_open(const char *, const int, const int=MDIO_READ);
00191 
00192 // Closes a molecular dynamics file.
00193 static int mdio_close(md_file *);
00194 
00195 
00196 // Format-independent file I/O routines
00197 static int mdio_header(md_file *, md_header *);
00198 static int mdio_timestep(md_file *, md_ts *);
00199 
00200 
00201 // .gro file functions
00202 static int gro_header(md_file *, char *, int, float *, int *, int = 1);
00203 static int gro_rec(md_file *, md_atom *);
00204 static int gro_timestep(md_file *, md_ts *);
00205 
00206 
00207 // .trX file functions
00208 static int trx_header(md_file *, int = 0);
00209 static int trx_int(md_file *, int *);
00210 static int trx_real(md_file *, float *);
00211 
00212 static int trx_rvector(md_file *, float *);
00213 static int trx_string(md_file *, char *, int);
00214 static int trx_timestep(md_file *, md_ts *);
00215 static int trx_timeskip(md_file *, md_ts *);
00216 
00217 // .g96 file functions
00218 static int g96_header(md_file *, char *, int, float *);
00219 static int g96_timestep(md_file *, md_ts *);
00220 static int g96_rec(md_file *, md_atom *);
00221 static int g96_countatoms(md_file *);
00222 
00223 
00224 // .xtc file functions
00225 static int xtc_int(md_file *, int *);
00226 static int xtc_float(md_file *, float *);
00227 /* 
00228 static int xtc_receivebits(int *, int);
00229 static void xtc_receiveints(int *, int, int, const unsigned *, int *);
00230 */
00231 static int xtc_timestep(md_file *, md_ts *);
00232 static int xtc_timeskip(md_file *, md_ts *);
00233 static int xtc_3dfcoord(md_file *, float *, int *, float *);
00234 
00235 
00236 // Error reporting functions
00237 static int mdio_errno(void);
00238 static const char *mdio_errmsg(int);
00239 static int mdio_seterror(int);
00240 
00241 
00242 // Miscellaneous functions
00243 static int strip_white(char *);
00244 static int mdio_readline(md_file *, char *, int, int = 1);
00245 static int mdio_tsfree(md_ts *, int = 0);
00246 static int mdio_readbox(md_box *, float *, float *, float *);
00247 
00248 
00249 
00250 static int xtc_receivebits(int *, int);
00251 
00252 // Error descriptions for mdio_errno
00253 static const char *mdio_errdescs[] = {
00254         "no error",
00255         "file does not match format",
00256         "unexpected end-of-file reached",
00257         "function called with bad parameters",
00258         "file i/o error",
00259         "unsupported precision",
00260         "out of memory",
00261         "cannot open file",
00262         "bad file extension",
00263         "unknown file format",
00264         "cannot close file",
00265         "wrong file format for this function",
00266         "binary i/o error: sizeof(int) != 4",
00267         NULL
00268 };
00269 
00272 static inline int host_is_little_endian(void) 
00273 {
00274   const union { unsigned char c[4]; unsigned int i; } 
00275   fixed = { { 0x10 , 0x20 , 0x40 , 0x80 } };
00276   const unsigned int i = 0x80402010U;
00277         
00278   if (fixed.i == i) {
00279     return 1;
00280   }
00281   return 0;
00282 }
00283 
00284 
00285 
00286 // Open a molecular dynamics file. The second parameter specifies
00287 // the format of the file. If it is zero, the format is determined
00288 // from the file extension.
00289 md_file *mdio_open(const char *fn, const int fmt, const int rw) {
00290         md_file *mf;
00291 
00292         if (!fn) {
00293                 mdio_seterror(MDIO_BADPARAMS);
00294                 return NULL;
00295         }
00296 
00297         // Allocate memory
00298         mf = (md_file *) malloc(sizeof(md_file));
00299         if (!mf) {
00300                 mdio_seterror(MDIO_BADMALLOC);
00301                 return NULL;
00302         }
00303 
00304         // Zero out the structure
00305         memset(mf, 0, sizeof(md_file));
00306 
00307         // Determine the file type from the extension
00308         if (!fmt) {
00309                 char *p;
00310                 int n;
00311 
00312                 // Seek to the extension part of the filename
00313                 for (p = (char *) &fn[strlen(fn) - 1]; *p != '.' && p > fn; p--);
00314                 if (p == fn) {
00315                         free(mf);
00316                         mdio_seterror(MDIO_BADEXTENSION);
00317                         return NULL;
00318                 }
00319 
00320                 // Check the extension against known extensions
00321                 for (n = 1; mdio_fmtexts[n]; n++)
00322                         if (!strcasecmp(p, mdio_fmtexts[n])) break;
00323 
00324                 // If !mdio_fmtexts[n], we failed (unknown ext)
00325                 if (!mdio_fmtexts[n]) {
00326                         free(mf);
00327                         mdio_seterror(MDIO_UNKNOWNFMT);
00328                         return NULL;
00329                 }
00330 
00331                 // All set
00332                 mf->fmt = n;
00333         }
00334         else {
00335                 mf->fmt = fmt;
00336         }
00337 
00338         // Differentiate between binary and ascii files. Also,
00339         // .trX files need a header information structure allocated.
00340         switch (mf->fmt) {
00341     case MDFMT_GRO:
00342         case MDFMT_G96: /* fallthrough */
00343         if (rw) 
00344             mf->f = fopen(fn, "wt");
00345         else
00346             mf->f = fopen(fn, "rt");
00347 
00348                 break;
00349         case MDFMT_TRR:
00350         case MDFMT_TRJ: /* fallthrough */
00351                 // Allocate the trx header data struct
00352                 mf->trx = (trx_hdr *) malloc(sizeof(trx_hdr));
00353                 if (!mf->trx) {
00354                         free(mf);
00355                         mdio_seterror(MDIO_BADMALLOC);
00356                         return NULL;
00357                 }
00358                 memset(mf->trx, 0, sizeof(trx_hdr));
00359         case MDFMT_XTC:  /* fallthrough */
00360                 // Finally, open the file
00361         if (rw)
00362             mf->f = fopen(fn, "wb");
00363         else
00364             mf->f = fopen(fn, "rb");
00365 
00366                 break;
00367         default:
00368                 free(mf);
00369                 mdio_seterror(MDIO_UNKNOWNFMT);
00370                 return NULL;
00371         }
00372 
00373         // Check for opening error
00374         if (!mf->f) {
00375                 if (mf->trx) free(mf->trx);
00376                 free(mf);
00377                 mdio_seterror(MDIO_CANTOPEN);
00378                 return NULL;
00379         }
00380 
00381         // File is opened, we're all set!
00382         mdio_seterror(MDIO_SUCCESS);
00383         return mf;
00384 }
00385 
00386 
00387 // Closes a molecular dynamics file.
00388 static int mdio_close(md_file *mf) {
00389         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00390 
00391         if (fclose(mf->f) == EOF) return mdio_seterror(MDIO_CANTCLOSE);
00392 
00393         // Free the dynamically allocated memory
00394         if (mf->trx) free(mf->trx);
00395         free(mf);
00396         return mdio_seterror(MDIO_SUCCESS);
00397 }
00398 
00399 
00400 // Returns the last error code reported by any of the mdio functions
00401 static int mdio_errno(void) {
00402         return mdio_errcode;
00403 }
00404 
00405 
00406 // Returns a textual message regarding an mdio error code
00407 static const char *mdio_errmsg(int n) {
00408         if (n < 0 || n > MDIO_MAX_ERRVAL) return (char *) "unknown error";
00409         else return mdio_errdescs[n];
00410 }
00411 
00412 
00413 // Sets the error code and returns an appropriate return value
00414 // for the calling function to return to its parent
00415 static int mdio_seterror(int code) {
00416         mdio_errcode = code;
00417         return code ? -1 : 0;
00418 }
00419 
00420 
00421 // Reads a line from the text file, strips leading/trailing whitespace
00422 // and newline, checks for errors, and returns the number of characters
00423 // in the string on success or -1 on error.
00424 static int mdio_readline(md_file *mf, char *buf, int n, int strip) {
00425         if (!buf || n < 1 || !mf) return mdio_seterror(MDIO_BADPARAMS);
00426 
00427         // Read the line
00428         fgets(buf, n, mf->f);
00429 
00430         // End of file reached?
00431         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
00432 
00433         // File I/O error?
00434         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
00435 
00436         // comment line?
00437         if (buf[0] == '#') return mdio_readline(mf,buf,n,strip);
00438 
00439         // Strip whitespace
00440         if (strip) strip_white(buf);
00441 
00442         return strlen(buf);
00443 }
00444 
00445 
00446 // Strips leading and trailing whitespace from a string. Tabs,
00447 // spaces, newlines and carriage returns are stripped. Example:
00448 // "\n   hello\t \r" becomes "hello".
00449 static int strip_white(char *buf) {
00450         int i, j, k;
00451 
00452         // Protect against NULL pointer
00453         if (!buf) return -1;
00454         if (!strlen(buf)) return -1;
00455 
00456         // Kill trailing whitespace first
00457         for (i = strlen(buf) - 1;
00458              buf[i] == ' ' || buf[i] == '\t' ||
00459              buf[i] == '\n' || buf[i] == '\r';
00460              i--)
00461                 buf[i] = 0;
00462 
00463         // Skip past leading whitespace
00464         for (i = 0; buf[i] == ' ' || buf[i] == '\t' ||
00465              buf[i] == '\n' || buf[i] == '\r'; i++);
00466         if (i) {
00467                 k = 0;
00468                 for (j = i; buf[j]; j++)
00469                         buf[k++] = buf[j];
00470                 buf[k] = 0;
00471         }
00472 
00473         return strlen(buf);
00474 }
00475 
00476 
00477 // Frees the memory allocated in a ts structure. The holderror
00478 // parameter defaults to zero. Programs that are calling this
00479 // function because of an error reported by another function should
00480 // set holderror so that mdio_tsfree() does not overwrite the error
00481 // code with mdio_seterror().
00482 static int mdio_tsfree(md_ts *ts, int holderror) {
00483         if (!ts) {
00484                 if (holderror) return -1;
00485                 else return mdio_seterror(MDIO_BADPARAMS);
00486         }
00487 
00488         if (ts->pos && ts->natoms > 0) free(ts->pos);
00489 
00490   if (ts->box) free(ts->box);
00491 
00492         if (holderror) return -1;
00493         else return mdio_seterror(MDIO_SUCCESS);
00494 }
00495 
00496 
00497 // Converts box basis vectors to A, B, C, alpha, beta, and gamma.  
00498 // Stores values in md_box struct, which should be allocated before calling
00499 // this function.
00500 static int mdio_readbox(md_box *box, float *x, float *y, float *z) {
00501   float A, B, C;
00502 
00503   if (!box) {
00504     return mdio_seterror(MDIO_BADPARAMS);
00505   }
00506 
00507   // A, B, C are the lengths of the x, y, z vectors, respectively
00508   A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] ) * ANGS_PER_NM;
00509   B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] ) * ANGS_PER_NM;
00510   C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] ) * ANGS_PER_NM;
00511   if ((A<=0) || (B<=0) || (C<=0)) {
00512     /* Use zero-length box size and set angles to 90. */
00513     box->A = box->B = box->C = 0;
00514     box->alpha = box->beta = box->gamma = 90;
00515   } else {
00516     box->A = A;
00517     box->B = B;
00518     box->C = C;
00519   
00520     // gamma, beta, alpha are the angles between the x & y, x & z, y & z
00521     // vectors, respectively
00522     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;
00523     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;
00524     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; 
00525   }
00526   return mdio_seterror(MDIO_SUCCESS);
00527 }
00528 
00529 
00530 // Reads the header of a file (format independent)
00531 static int mdio_header(md_file *mf, md_header *mdh) {
00532         int n;
00533         if (!mf || !mdh) return mdio_seterror(MDIO_BADPARAMS);
00534         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00535 
00536         switch (mf->fmt) {
00537         case MDFMT_GRO:
00538                 if (gro_header(mf, mdh->title, MAX_MDIO_TITLE,
00539                 &mdh->timeval, &mdh->natoms, 1) < 0)
00540                         return -1;
00541                 return 0;
00542 
00543         case MDFMT_TRR: 
00544         case MDFMT_TRJ: /* fallthrough */
00545                 if (trx_header(mf, 1) < 0) return -1;
00546                 mdh->natoms = mf->trx->natoms;
00547                 mdh->timeval = (float) mf->trx->t;
00548                 strncpy(mdh->title, mf->trx->title, MAX_MDIO_TITLE);
00549                 return 0;
00550 
00551         case MDFMT_G96:
00552                 if (g96_header(mf, mdh->title, MAX_MDIO_TITLE,
00553                 &mdh->timeval) < 0) return -1;
00554                 mdh->natoms = -1;
00555                 return 0;
00556 
00557         case MDFMT_XTC:
00558                 memset(mdh, 0, sizeof(md_header));
00559                 // Check magic number
00560                 if (xtc_int(mf, &n) < 0) return -1;
00561                 if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
00562 
00563                 // Get number of atoms
00564                 if (xtc_int(mf, &n) < 0) return -1;
00565                 mdh->natoms = n;
00566                 rewind(mf->f);
00567                 return 0;
00568 
00569         default:
00570                 return mdio_seterror(MDIO_UNKNOWNFMT);
00571         }
00572 }
00573 
00574 
00575 // Reads in a timestep from a file (format independent)
00576 static int mdio_timestep(md_file *mf, md_ts *ts) {
00577         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00578         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00579 
00580         switch (mf->fmt) {
00581         case MDFMT_GRO:
00582                 return gro_timestep(mf, ts);
00583 
00584         case MDFMT_TRR:
00585         case MDFMT_TRJ: /* fallthrough */
00586                 return trx_timestep(mf, ts);
00587 
00588         case MDFMT_G96:
00589                 return g96_timestep(mf, ts);
00590 
00591         case MDFMT_XTC:
00592                 return xtc_timestep(mf, ts);
00593 
00594         default:
00595                 return mdio_seterror(MDIO_UNKNOWNFMT);
00596         }
00597 }
00598 
00599 
00600 // Reads in a timestep from a file (format independent)
00601 static int mdio_timeskip(md_file *mf, md_ts *ts) {
00602         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00603         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00604 
00605         switch (mf->fmt) {
00606         case MDFMT_GRO:
00607                 return gro_timestep(mf, ts);
00608 
00609         case MDFMT_TRR:
00610         case MDFMT_TRJ: /* fallthrough */
00611                 return trx_timeskip(mf, ts);
00612 
00613         case MDFMT_G96:
00614                 return g96_timestep(mf, ts);
00615 
00616         case MDFMT_XTC:
00617                 return xtc_timeskip(mf, ts);
00618 
00619         default:
00620                 return mdio_seterror(MDIO_UNKNOWNFMT);
00621         }
00622 }
00623 
00624 static int g96_header(md_file *mf, char *title, int titlelen, float *timeval) {
00625         char buf[MAX_G96_LINE + 1];
00626         char *p;
00627 
00628         // Check parameters
00629         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00630 
00631         // The header consists of blocks. The title block
00632         // is mandatory, and a TIMESTEP block is optional.
00633         // Example:
00634         //
00635         // TITLE
00636         // Generated by trjconv :  t=  90.00000
00637         // more title info
00638         // .
00639         // .
00640         // .
00641         // END
00642         // .
00643         // .
00644         // .
00645 
00646         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00647         if (strcasecmp(buf, "TITLE")) return mdio_seterror(MDIO_BADFORMAT);
00648 
00649         // Read in the title itself
00650         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00651 
00652         // The timevalue can be included in the title string
00653         // after a "t=" prefix.
00654         if ((p = (char *) strstr(buf, "t="))) {
00655                 char *q = p;
00656                 *(q--) = 0;
00657 
00658                 // Skip the `t=' and strip whitespace from
00659                 // the resulting strings
00660                 p += 2;
00661                 strip_white(p);
00662                 strip_white(buf);
00663 
00664                 // Grab the timevalue from the title string
00665                 if (timeval) *timeval = (float) atof(p);
00666         }
00667         else {
00668                 // No timevalue - just copy the string and strip
00669                 // any leading/trailing whitespace
00670                 if (timeval) *timeval = 0;
00671                 strip_white(buf);
00672         }
00673 
00674         // Copy the title string
00675         if (title && titlelen) strncpy(title, buf, titlelen);
00676 
00677         // Now ignore subsequent title lines and get the END string
00678         while (strcasecmp(buf, "END"))
00679                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00680 
00681         // Done!
00682         return mdio_seterror(MDIO_SUCCESS);
00683 }
00684 
00685 
00686 // Used to determine the number of atoms in a g96 file, because
00687 // VMD needs to know this for some reason.
00688 static int g96_countatoms(md_file *mf) {
00689         char buf[MAX_G96_LINE + 1];
00690         int natoms;
00691         int n;
00692         long fpos;
00693         float lastf;
00694 
00695         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00696 
00697         fpos = ftell(mf->f);
00698 
00699         natoms = 0;
00700         for (;;) {
00701                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00702                         break;
00703                 n = sscanf(buf, "%*6c%*6c%*6c%*6c %*f %*f %f", &lastf);
00704                 if (n == 1) natoms++;
00705                 else {
00706                         strip_white(buf);
00707                         if (!strcasecmp(buf, "END")) break;
00708                 }
00709         }
00710 
00711         fseek(mf->f, fpos, SEEK_SET);
00712         return natoms;
00713 }
00714 
00715 
00716 // Reads an atom line from the G96 file
00717 static int g96_rec(md_file *mf, md_atom *ma) {
00718         char buf[MAX_G96_LINE + 1];
00719         char atomnum[7];
00720         int n;
00721 
00722         // Check parameters
00723         if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);
00724 
00725         // Read in a line, assuming it is an atom line
00726         do {
00727                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0) return -1;
00728         } while (buf[0] == '#' || strlen(buf) == 0);
00729 
00730         n = sscanf(buf, "%6c%6c%6c%6c %f %f %f",
00731                 ma->resid, ma->resname, ma->atomname, atomnum,
00732                 &ma->pos[0], &ma->pos[1], &ma->pos[2]);
00733         if (n == 7) {
00734                 atomnum[6] = 0;
00735                 ma->resid[6] = 0;
00736                 ma->resname[6] = 0;
00737                 ma->atomname[6] = 0;
00738 
00739                 strip_white(atomnum);
00740                 strip_white(ma->resid);
00741                 strip_white(ma->resname);
00742                 strip_white(ma->atomname);
00743 
00744                 ma->atomnum = atoi(atomnum);
00745 
00746                 ma->pos[0] *= ANGS_PER_NM;
00747                 ma->pos[1] *= ANGS_PER_NM;
00748                 ma->pos[2] *= ANGS_PER_NM;
00749 
00750                 return 0;
00751         }
00752 
00753         return mdio_seterror(MDIO_BADFORMAT);
00754 }
00755 
00756 
00757 // Reads a timestep from a G96 file and stores the data in
00758 // the generic md_ts structure. Returns 0 on success or a
00759 // negative number on error and sets mdio_errcode.
00760 static int g96_timestep(md_file *mf, md_ts *ts) {
00761         char            buf[MAX_G96_LINE + 1];
00762         char            stripbuf[MAX_G96_LINE + 1];
00763         float           pos[3], x[3], y[3], z[3], *currAtom;
00764         long            fpos;
00765         int             n, i, boxItems;
00766 
00767         // Check parameters
00768         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00769 
00770   // Allocate data space for the timestep, using the number of atoms
00771   // determined by open_g96_read().
00772         ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
00773         if (!ts->pos) {
00774                 return mdio_seterror(MDIO_BADMALLOC);
00775         }
00776   currAtom = ts->pos;
00777 
00778         // The timesteps follow the header in a fixed block
00779         // format:
00780         //
00781         // TIMESTEP
00782         //         <step number> <time value>
00783         // END
00784         // POSITIONRED
00785         //     <x float> <y float> <z float>
00786         //     .         .         .
00787         //     .         .         .
00788         //     .         .         .
00789         // END
00790         // VELOCITYRED
00791         //     <x float> <y float> <z float>
00792         //     .         .         .
00793         //     .         .         .
00794         //     .         .         .
00795         // END
00796         // BOX
00797         //     <x float> <y float> <z float>
00798         // END
00799         //
00800         // -----
00801         //
00802         // The TIMESTEP, VELOCITY and BOX blocks are optional.
00803         // Floats are written in 15.9 precision.
00804         //
00805         // Reference: GROMACS 2.0 user manual
00806         //            http://rugmd4.chem.rug.nl/~gmx/online2.0/g96.html
00807 
00808         // First, look for an (optional) title block and skip it
00809         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00810 
00811   if (!strcasecmp(buf, "TITLE")) {
00812     // skip over the text until we reach 'END'
00813     while (strcasecmp(buf, "END")) {
00814       if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00815     }
00816 
00817     // Read in the next line
00818     if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00819   }
00820 
00821         // Next, look for a timestep block
00822         if (!strcasecmp(buf, "TIMESTEP")) {
00823                 // Read in the value line
00824                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00825 
00826                 // Extract the time value and the timestep index
00827                 n = sscanf(buf, "%d %f", &ts->step, &ts->time);
00828                 if (n != 2) return mdio_seterror(MDIO_BADFORMAT);
00829 
00830                 // Read the "END" line
00831                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00832                 if (strcasecmp(buf, "END"))
00833                         return mdio_seterror(MDIO_BADFORMAT);
00834 
00835                 // Read in the next line
00836                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00837         }
00838         else {
00839                 // No timestep specified -- set to zero
00840                 ts->step = 0;
00841                 ts->time = 0;
00842         }
00843 
00844         // At this point a POSITION or POSITIONRED block
00845         // is REQUIRED by the format
00846         if (!strcasecmp(buf, "POSITIONRED")) {
00847 
00848     // So now we read in some atoms
00849     i = 0;
00850                 while (i < ts->natoms) {
00851                         // Read in an atom
00852                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00853                                 return -1;
00854  
00855       // We shouldn't reach the end yet
00856       if (!strcasecmp(buf, "END"))
00857         return mdio_seterror(MDIO_BADFORMAT);
00858 
00859                         // Get the x,y,z coordinates
00860                         n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
00861       
00862       // Ignore improperly formatted lines
00863                         if (n == 3) {
00864                                 pos[0] *= ANGS_PER_NM;
00865                                 pos[1] *= ANGS_PER_NM;
00866                                 pos[2] *= ANGS_PER_NM;
00867 
00868                                 // Copy the atom data into the array
00869                                 memcpy(currAtom, pos, sizeof(float) * 3);
00870         currAtom += 3;
00871         i++;
00872                         }
00873                 }
00874         }
00875         else if (!strcasecmp(buf, "POSITION") || !strcasecmp(buf, "REFPOSITION")) {
00876                 /*
00877                 char resnum[7];
00878                 char resname[7];
00879                 char atomname[7];
00880                 char atomnum[7];
00881                 */
00882 
00883                 // So now we read in some atoms
00884     i = 0;
00885                 while (i < ts->natoms) {
00886                         // Read in the first line
00887                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00888                                 return -1;
00889  
00890       // We shouldn't reach the end yet
00891       strcpy(stripbuf, buf);
00892       strip_white(stripbuf); 
00893       if (!strcasecmp(stripbuf, "END"))
00894         return mdio_seterror(MDIO_BADFORMAT);
00895 
00896                         // Get the x,y,z coordinates and name data
00897                         n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
00898                                 &pos[0], &pos[1], &pos[2]);
00899 
00900       // Ignore improperly formatted lines
00901                         if (n == 3) {
00902                                 pos[0] *= ANGS_PER_NM;
00903                                 pos[1] *= ANGS_PER_NM;
00904                                 pos[2] *= ANGS_PER_NM;
00905 
00906                                 // Copy the atom data into the linked list item
00907                                 memcpy(currAtom, pos, sizeof(float) * 3);
00908                                 currAtom += 3;
00909         i++;
00910                         }
00911                 }
00912         }
00913         else {
00914                 return mdio_seterror(MDIO_BADFORMAT);
00915         }
00916 
00917   // Read the END keyword
00918   if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00919     return -1;
00920   if (strcasecmp(buf, "END"))
00921     return mdio_seterror(MDIO_BADFORMAT);
00922 
00923         // ... another problem: there may or may not be a VELOCITY
00924         // block or a BOX block, so we need to read one line beyond
00925         // the POSITION block to determine this. If neither VEL. nor
00926         // BOX are present we've read a line too far and infringed
00927         // on the next timestep, so we need to keep track of the
00928         // position now for a possible fseek() later to backtrack.
00929         fpos = ftell(mf->f);
00930 
00931         // Now we must read in the velocities and the box, if present
00932         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00933     // It's okay if we end the file here; any other errors need to be
00934     // reported.
00935     if (mdio_errcode == MDIO_EOF) 
00936       return mdio_seterror(MDIO_SUCCESS);
00937     else 
00938       return -1;
00939   }
00940 
00941         // Is there a velocity block present ?
00942         if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
00943                 // Ignore all the coordinates - VMD doesn't use them
00944                 for (;;) {
00945                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00946                                 return -1;
00947                         if (!strcasecmp(buf, "END")) break;
00948                 }
00949 
00950                 // Again, record our position because we may need
00951                 // to fseek here later if we read too far.
00952                 fpos = ftell(mf->f);
00953 
00954                 // Go ahead and read the next line.
00955                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00956         }
00957 
00958         // Is there a box present ?
00959         if (!strcasecmp(buf, "BOX")) {
00960                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00961     boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
00962                &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
00963     if (boxItems == 3) {
00964       x[1] = x[2] = 0;
00965       y[0] = y[2] = 0;
00966       z[0] = z[1] = 0;
00967     }
00968     else if (boxItems != 9) 
00969       return mdio_seterror(MDIO_BADFORMAT);
00970 
00971     // Allocate the box and convert the vectors.
00972     ts->box = (md_box *) malloc(sizeof(md_box));
00973     if (mdio_readbox(ts->box, x, y, z) < 0) {
00974       free(ts->box);
00975       ts->box = NULL;
00976       return mdio_seterror(MDIO_BADFORMAT);
00977     }
00978 
00979                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00980       free(ts->box);
00981       ts->box = NULL;
00982       return -1;
00983     }
00984                 if (strcasecmp(buf, "END")) {
00985       free(ts->box);
00986       ts->box = NULL;
00987                         return mdio_seterror(MDIO_BADFORMAT);
00988     }
00989         }
00990         else {
00991                 // We have read too far, so fseek back to the
00992                 // last known safe position so we don't return
00993                 // with the file pointer set infringing on the
00994                 // next timestep data.
00995                 fseek(mf->f, fpos, SEEK_SET);
00996         }
00997 
00998         // We're done!
00999         return mdio_seterror(MDIO_SUCCESS);
01000 }
01001 
01002 
01003 // Attempts to read header data from a GROMACS structure file
01004 // The GROMACS header format is as follows (fixed, 2 lines ASCII):
01005 // <title> [ n= <timevalue> ]
01006 //     <num atoms>
01007 static int gro_header(md_file *mf, char *title, int titlelen, float *timeval,
01008                int *natoms, int rewind) {
01009   char buf[MAX_GRO_LINE + 1];
01010   long fpos;
01011   char *p;
01012 
01013   // Check parameters
01014   if (!mf)
01015     return mdio_seterror(MDIO_BADPARAMS);
01016 
01017   // Get the current file position for rewinding later
01018   fpos = ftell(mf->f);
01019 
01020   // The header consists of 2 lines - get the first line
01021   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
01022 
01023   // The timevalue can be included in the title string
01024   // after a "t=" prefix.
01025   if ((p = (char *) strstr(buf, "t="))) {
01026     char *q = p;
01027     *(q--) = 0;
01028 
01029     // Skip the `t=' and strip whitespace from
01030     // the resulting strings
01031     p += 2;
01032     strip_white(p);
01033     strip_white(buf);
01034 
01035     // Grab the timevalue from the title string
01036     if (timeval) *timeval = (float) atof(p);
01037   } else {
01038     // No timevalue - just copy the string
01039     if (timeval) *timeval = 0;
01040   }
01041 
01042   // Copy the title string
01043   if (title && titlelen) strncpy(title, buf, titlelen);
01044 
01045   // Get the second line and grab the number of atoms
01046   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
01047 
01048   // Store the number of atoms
01049   if (natoms && (!(*natoms = atoi(buf))))
01050     return mdio_seterror(MDIO_BADFORMAT);
01051 
01052   // Now we rewind the file so that subsequent calls to
01053   // gro_timestep() will succeed. gro_timestep() requires
01054   // the header to be at the current file pointer.
01055   if (rewind)
01056     fseek(mf->f, fpos, SEEK_SET);
01057 
01058   return 0; // Done!
01059 }
01060 
01061 
01062 // Reads one atom record from a GROMACS file. Returns GMX_SUCCESS
01063 // on success or a negative number on error.
01064 //
01065 // Record format (one line, fixed):
01066 //    rrrrrRRRRRaaaaaAAAAA <x pos> <y pos> <z pos> <x vel> <y vel> <z vel>
01067 //
01068 //    r = residue number
01069 //    R = residue name
01070 //    a = atom name
01071 //    A = atom number
01072 //
01073 static int gro_rec(md_file *mf, md_atom *ma) {
01074   char buf[MAX_GRO_LINE + 1];
01075   char atomnum[6];
01076   char xposc[12], yposc[12], zposc[12];
01077   int n;
01078 
01079   if (!mf)
01080     return mdio_seterror(MDIO_BADPARAMS);
01081 
01082   do {
01083     if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0)
01084       return -1;
01085   } while (buf[0] == '#' || !strlen(buf));
01086 
01087   // Read in the fields
01088   n = sscanf(buf, "%5c%5c%5c%5c%8c%8c%8c", 
01089              ma->resid, ma->resname, ma->atomname, atomnum, 
01090              xposc, yposc, zposc);
01091 
01092   if (n != 7)
01093     return mdio_seterror(MDIO_BADFORMAT);
01094 
01095   // Null terminate the strings
01096   ma->resname[5] = 0;
01097   ma->atomname[5] = 0;
01098   ma->resid[5] = 0;
01099   atomnum[5] = 0;
01100   xposc[8] = 0;
01101   yposc[8] = 0;
01102   zposc[8] = 0;
01103  
01104   if ((sscanf(xposc, "%f", &ma->pos[0]) != 1) ||
01105       (sscanf(yposc, "%f", &ma->pos[1]) != 1) ||
01106       (sscanf(zposc, "%f", &ma->pos[2]) != 1)) {
01107     return mdio_seterror(MDIO_BADFORMAT);
01108   }
01109 
01110   // Convert strings to numbers
01111   strip_white(atomnum);
01112   ma->atomnum = atoi(atomnum);
01113 
01114   // Convert nanometers to angstroms
01115   ma->pos[0] *= ANGS_PER_NM;
01116   ma->pos[1] *= ANGS_PER_NM;
01117   ma->pos[2] *= ANGS_PER_NM;
01118 
01119   // Strip leading and trailing whitespace
01120   strip_white(ma->atomname);
01121   strip_white(ma->resname);
01122   strip_white(ma->resid);
01123 
01124   return 0;
01125 }
01126 
01127 
01128 // Reads in a timestep from a .gro file. Ignores the data
01129 // not needed for a timestep, so is a little faster than
01130 // calling gro_rec() for each atom. Also reads in the
01131 // header block.
01132 //
01133 static int gro_timestep(md_file *mf, md_ts *ts) {
01134         char buf[MAX_GRO_LINE + 1];
01135         long coord;
01136         int i, n, boxItems;
01137   float x[3], y[3], z[3];
01138   char xposc[12], yposc[12], zposc[12];
01139 
01140   if (!mf || !ts) 
01141     return mdio_seterror(MDIO_BADPARAMS);
01142 
01143   if (gro_header(mf, NULL, 0, &ts->time, &ts->natoms, 0) < 0)
01144     return -1;
01145 
01146   ts->pos = (float *) malloc(3 * sizeof(float) * ts->natoms);
01147   if (!ts->pos)
01148     return mdio_seterror(MDIO_BADMALLOC);
01149 
01150   coord = 0;
01151   for (i = 0; i < ts->natoms; i++) {
01152     if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01153       free(ts->pos);
01154       return -1;
01155     }
01156         
01157     n = sscanf(buf, "%*5c%*5c%*5c%*5c%8c%8c%8c", xposc, yposc, zposc);
01158     if (n != 3) 
01159       return mdio_seterror(MDIO_BADFORMAT);
01160 
01161     if ((sscanf(xposc, "%f", &ts->pos[coord    ]) != 1) ||
01162         (sscanf(yposc, "%f", &ts->pos[coord + 1]) != 1) ||
01163         (sscanf(zposc, "%f", &ts->pos[coord + 2]) != 1)) {
01164       return mdio_seterror(MDIO_BADFORMAT);
01165     }
01166 
01167     ts->pos[coord    ] *= ANGS_PER_NM;
01168     ts->pos[coord + 1] *= ANGS_PER_NM;
01169     ts->pos[coord + 2] *= ANGS_PER_NM;
01170 
01171     coord += 3;
01172   }
01173 
01174   // Read the box, stored as three vectors representing its edges
01175   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01176     free(ts->pos);
01177     return -1;
01178   }
01179 
01180   boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
01181              &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
01182 
01183   // File may only include three scalars for the box information -- if
01184   // that's the case, the box is orthoganal.
01185   if (boxItems == 3) {
01186     x[1] = x[2] = 0;
01187     y[0] = y[2] = 0;
01188     z[0] = z[1] = 0;
01189   } else if (boxItems != 9) {
01190     free(ts->pos);
01191     return -1;
01192   }
01193 
01194   // Allocate the box and convert the vectors.
01195   ts->box = (md_box *) malloc(sizeof(md_box));
01196   if (mdio_readbox(ts->box, x, y, z) < 0) {
01197     free(ts->pos);
01198     free(ts->box);
01199     ts->box = NULL;
01200     return -1;
01201   }
01202 
01203   return 0;
01204 }
01205 
01206 
01207 // Attempts to read header data from a .trX trajectory file
01208 //
01209 // The .trX header format is as follows:
01210 //
01211 //      4 bytes         - magic number (0x07C9)
01212 //      ...
01213 //
01214 static int trx_header(md_file *mf, int rewind) {
01215         int magic;
01216         trx_hdr *hdr;
01217         long fpos;
01218 
01219         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01220 
01221         // In case we need to rewind
01222         fpos = ftell(mf->f);
01223 
01224         // We need to store some data to the trX header data
01225         // structure inside the md_file structure
01226         hdr = mf->trx;
01227         if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);
01228 
01229         // Read the magic number
01230         if (trx_int(mf, &magic) < 0) return -1;
01231         if (magic != TRX_MAGIC) {
01232                 // Try reverse endianism
01233                 swap4_aligned(&magic, 1);
01234                 if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01235 
01236                 // Enable byte swapping (actually works, too!)
01237                 mf->rev = 1;
01238         }
01239 
01240         // Read the version number. 
01241         // XXX. this is not the version number, but the storage size
01242         // of the following XDR encoded string.
01243         // the 'title' string is in fact the version identifier.
01244         // since VMD does not use any of that, it does no harm,
01245         // but is should still be fixed occasionally. AK 2005/01/08.
01246 
01247         if(mf->fmt!=MDFMT_TRJ) {
01248                 // It appears that TRJ files either don't contain a version
01249                 // number or don't have a length-delimiter on the string,
01250                 // whereas TRR files do contain both.  Thus, with TRJ, we just
01251                 // assume that the version number is the string length and 
01252                 // just hope for the best. -- WLD 2006/07/09
01253                 if (trx_int(mf, &hdr->version) < 0) return -1;
01254         }
01255 
01256         // Read in the title string
01257         if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
01258                 return -1;
01259 
01260         // Read in some size data
01261         if (trx_int(mf, &hdr->ir_size) < 0) return -1;
01262         if (trx_int(mf, &hdr->e_size) < 0) return -1;
01263         if (trx_int(mf, &hdr->box_size) < 0) return -1;
01264         if (trx_int(mf, &hdr->vir_size) < 0) return -1;
01265         if (trx_int(mf, &hdr->pres_size) < 0) return -1;
01266         if (trx_int(mf, &hdr->top_size) < 0) return -1;
01267         if (trx_int(mf, &hdr->sym_size) < 0) return -1;
01268         if (trx_int(mf, &hdr->x_size) < 0) return -1;
01269         if (trx_int(mf, &hdr->v_size) < 0) return -1;
01270         if (trx_int(mf, &hdr->f_size) < 0) return -1;
01271         if (trx_int(mf, &hdr->natoms) < 0) return -1;
01272         if (trx_int(mf, &hdr->step) < 0) return -1;
01273         if (trx_int(mf, &hdr->nre) < 0) return -1;
01274 
01275         // Make sure there are atoms...
01276         if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);
01277 
01278         // Try to determine precision (float? double?)
01279         if (hdr->x_size) mf->prec = hdr->x_size / (hdr->natoms * 3);
01280         else if (hdr->v_size) mf->prec = hdr->v_size / (hdr->natoms * 3);
01281         else if (hdr->f_size) mf->prec = hdr->f_size / (hdr->natoms * 3);
01282         else return mdio_seterror(MDIO_BADPRECISION);
01283 
01284         if (mf->prec != sizeof(float) && mf->prec != sizeof(double)) {
01285                 // We have no data types this size! The
01286                 // file must've been generated on another
01287                 // platform
01288                 return mdio_seterror(MDIO_BADPRECISION);
01289         }
01290 
01291         // Read in timestep and lambda
01292         if (trx_real(mf, &hdr->t) < 0) return -1;
01293         if (trx_real(mf, &hdr->lambda) < 0) return -1;
01294 
01295         // Rewind if necessary
01296         if (rewind) fseek(mf->f, fpos, SEEK_SET);
01297 
01298         return 0;
01299 }
01300 
01301 
01302 // Reads in an integer and stores it in y. Returns GMX_SUCCESS
01303 // on success or a negative number on error.
01304 static int trx_int(md_file *mf, int *y) {
01305         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01306 
01307         // sanity check.
01308         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01309 
01310         if (y) {
01311                 if (fread(y, 4, 1, mf->f) != 1)
01312                         return mdio_seterror(MDIO_IOERROR);
01313                 if (mf->rev) swap4_aligned(y, 1);
01314         }
01315         else if (fseek(mf->f, 4, SEEK_CUR) != 0)
01316                 return mdio_seterror(MDIO_IOERROR);
01317 
01318         return mdio_seterror(MDIO_SUCCESS);
01319 }
01320 
01321 static int trx_long(md_file *mf, long *y) {
01322         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01323 
01324         // sanity check.
01325         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01326 
01327         if (y) {
01328                 if (fread(y, 8, 1, mf->f) != 1)
01329                         return mdio_seterror(MDIO_IOERROR);
01330                 if (mf->rev) swap8_aligned(y, 1);
01331         }
01332         else if (fseek(mf->f, 4, SEEK_CUR) != 0)
01333                 return mdio_seterror(MDIO_IOERROR);
01334 
01335         return mdio_seterror(MDIO_SUCCESS);
01336 }
01337 
01338 
01339 // Reads in either a float or a double, depending on the
01340 // precision, and stores that number in y. Returns
01341 // GMX_SUCCESS on success or a negative number on error.
01342 static int trx_real(md_file *mf, float *y) {
01343         double x;
01344 
01345         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01346 
01347         switch (mf->prec) {
01348                 case sizeof(float):
01349                         if (!y) {
01350                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01351                                         return mdio_seterror(MDIO_IOERROR);
01352                         } else {
01353                                 if (fread(y, mf->prec, 1, mf->f) != 1)
01354                                         return mdio_seterror(MDIO_IOERROR);
01355                                 if (mf->rev) swap4_aligned(y, 1);
01356                         }
01357                         return mdio_seterror(MDIO_SUCCESS);
01358 
01359                 case sizeof(double):
01360                         if (!y) {
01361                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01362                                         return mdio_seterror(MDIO_IOERROR);
01363                         } else {
01364                                 if (fread(&x, mf->prec, 1, mf->f) != 1)
01365                                         return mdio_seterror(MDIO_IOERROR);
01366                                 if (mf->rev) swap8_aligned(&x, 1);
01367                                 *y = (float) x;
01368                         }
01369                         return mdio_seterror(MDIO_SUCCESS);
01370 
01371                 default:
01372                         return mdio_seterror(MDIO_BADPRECISION);
01373         }
01374 
01375 }
01376 
01377 static int trx_double(md_file *mf, double *y) {
01378         double x;
01379 
01380         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01381 
01382         if (!y) {
01383                 if (fseek(mf->f, 8, SEEK_CUR) != 0)
01384                         return mdio_seterror(MDIO_IOERROR);
01385         } else {
01386                 if (fread(&x, 8, 1, mf->f) != 1)
01387                         return mdio_seterror(MDIO_IOERROR);
01388                 if (mf->rev) swap8_aligned(&x, 1);
01389                 *y = (float) x;
01390         }
01391         return mdio_seterror(MDIO_SUCCESS);
01392 
01393 }
01394 
01395 
01396 // Reads in a real-valued vector (taking precision into account).
01397 // Stores the vector in vec, and returns GMX_SUCCESS on success
01398 // or a negative number on error.
01399 static int trx_rvector(md_file *mf, float *vec) {
01400         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01401 
01402         if (!vec) {
01403                 if (trx_real(mf, NULL) < 0) return -1;
01404                 if (trx_real(mf, NULL) < 0) return -1;
01405                 if (trx_real(mf, NULL) < 0) return -1;
01406                 return mdio_seterror(MDIO_SUCCESS);
01407         } else {
01408                 if (trx_real(mf, &vec[0]) < 0) return -1;
01409                 if (trx_real(mf, &vec[1]) < 0) return -1;
01410                 if (trx_real(mf, &vec[2]) < 0) return -1;
01411                 return mdio_seterror(MDIO_SUCCESS);
01412         }
01413 }
01414 
01415 static int tpr_rvector (md_file *mf, float *arr, int len) {
01416     int i;
01417     for (i = 0; i < len; i++) {
01418         if (trx_real(mf, &(arr[i]))) return -1;
01419     }
01420     return mdio_seterror(MDIO_SUCCESS);
01421 }
01422 
01423 static int tpr_ivector (md_file *mf, int *arr, int len) {
01424     int i;
01425     for (i = 0; i < len; i++) {
01426         if (trx_int(mf, &(arr[i]))) return -1;
01427     }
01428     return mdio_seterror(MDIO_SUCCESS);
01429 }
01430 
01431 // Reads in a string by first reading an integer containing the
01432 // string's length, then reading in the string itself and storing
01433 // it in str. If the length is greater than max, it is truncated
01434 // and the rest of the string is skipped in the file. Returns the
01435 // length of the string on success or a negative number on error.
01436 static int trx_string(md_file *mf, char *str, int max) {
01437         int size;
01438   size_t ssize;
01439 
01440         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01441 
01442         if (trx_int(mf, &size) < 0) return -1;
01443   ssize = (size_t)size;
01444 
01445         if (str && size <= max) {
01446                 if (fread(str, 1, size, mf->f) != ssize)
01447                         return mdio_seterror(MDIO_IOERROR);
01448                 str[size] = 0;
01449                 return size;
01450         } else if (str) {
01451                 if (fread(str, 1, max, mf->f) != ssize)
01452                         return mdio_seterror(MDIO_IOERROR);
01453                 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01454                         return mdio_seterror(MDIO_IOERROR);
01455                 str[max] = 0;
01456                 return max;
01457         } else {
01458                 if (fseek(mf->f, size, SEEK_CUR) != 0)
01459                         return mdio_seterror(MDIO_IOERROR);
01460                 return 0;
01461         }
01462 }
01463 
01464 
01465 // Reads in a string by first reading an integer containing the
01466 // string's length, then reading in the string itself and storing
01467 // it in str. If the length is greater than max, it is truncated
01468 // and the rest of the string is skipped in the file. Returns the
01469 // length of the string on success or a negative number on error.
01470 // Unlike trx_string, this will read in a 4-byte aligned way.
01471 static int tpr_string(md_file *mf, char *str, int max) {
01472         int size;
01473   size_t ssize;
01474 
01475         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01476 
01477         if (trx_int(mf, &size) < 0) return -1;
01478         if (size % 4) {
01479                 size += 4 - (size % 4);
01480         }
01481   ssize = (size_t)size;
01482 
01483         if (str && size <= max) {
01484                 if (fread(str, 1, size, mf->f) != ssize)
01485                         return mdio_seterror(MDIO_IOERROR);
01486                 str[size] = 0;
01487                 return size;
01488         } else if (str) {
01489                 if (fread(str, 1, max, mf->f) != ssize)
01490                         return mdio_seterror(MDIO_IOERROR);
01491                 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01492                         return mdio_seterror(MDIO_IOERROR);
01493                 str[max] = 0;
01494                 return max;
01495         } else {
01496                 if (fseek(mf->f, size, SEEK_CUR) != 0)
01497                         return mdio_seterror(MDIO_IOERROR);
01498                 return 0;
01499         }
01500 }
01501 
01502 // Reads in a timestep frame from the .trX file and returns the
01503 // data in a timestep structure. Returns NULL on error.
01504 static int trx_timestep(md_file *mf, md_ts *ts) {
01505   int i;
01506   float x[3], y[3], z[3];
01507   trx_hdr *hdr;
01508 
01509   if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01510   if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
01511     return mdio_seterror(MDIO_WRONGFORMAT);
01512 
01513   // Read the header
01514   if (trx_header(mf) < 0) return -1;
01515 
01516   // We need some data from the trX header
01517   hdr = mf->trx;
01518   if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
01519 
01520   if (hdr->box_size) { // XXX need to check value of box_size!!
01521     if (trx_rvector(mf, x) < 0) return -1;
01522     if (trx_rvector(mf, y) < 0) return -1;
01523     if (trx_rvector(mf, z) < 0) return -1;
01524 
01525     // Allocate the box and convert the vectors.
01526     ts->box = (md_box *) malloc(sizeof(md_box));
01527     if (mdio_readbox(ts->box, x, y, z) < 0) {
01528       free(ts->box);
01529       ts->box = NULL;
01530       return -1;
01531     }
01532   }
01533 
01534   if (hdr->vir_size) {
01535     if (trx_rvector(mf, NULL) < 0) return -1;
01536     if (trx_rvector(mf, NULL) < 0) return -1;
01537     if (trx_rvector(mf, NULL) < 0) return -1;
01538   }
01539 
01540   if (hdr->pres_size) {
01541     if (trx_rvector(mf, NULL) < 0) return -1;
01542     if (trx_rvector(mf, NULL) < 0) return -1;
01543     if (trx_rvector(mf, NULL) < 0) return -1;
01544   }
01545 
01546   if (hdr->x_size) {
01547     ts->pos = (float *) malloc(sizeof(float) * 3 * hdr->natoms);
01548     if (!ts->pos) 
01549       return mdio_seterror(MDIO_BADMALLOC);
01550 
01551     ts->natoms = hdr->natoms;
01552     for (i = 0; i < hdr->natoms; i++) {
01553       if (trx_rvector(mf, &ts->pos[i * 3]) < 0) {
01554         mdio_tsfree(ts, 1);
01555         return -1;
01556       }
01557 
01558       ts->pos[i * 3    ] *= ANGS_PER_NM;
01559       ts->pos[i * 3 + 1] *= ANGS_PER_NM;
01560       ts->pos[i * 3 + 2] *= ANGS_PER_NM;
01561     }
01562   }
01563 
01564   if (hdr->v_size) {
01565     for (i = 0; i < hdr->natoms; i++) {
01566       if (trx_rvector(mf, NULL) < 0) {
01567         mdio_tsfree(ts, 1);
01568         return -1;
01569       }
01570     }
01571   }
01572 
01573   if (hdr->f_size) {
01574     for (i = 0; i < hdr->natoms; i++) {
01575       if (trx_rvector(mf, NULL) < 0) {
01576         mdio_tsfree(ts, 1);
01577         return -1;
01578       }
01579     }
01580   }
01581 
01582   return mdio_seterror(MDIO_SUCCESS);
01583 }
01584 
01585 // Reads in a timestep header from the .trX file and skips the rest.
01586 static int trx_timeskip(md_file *mf, md_ts *ts) {
01587   int i;
01588   float x[3], y[3], z[3];
01589   trx_hdr *hdr;
01590 
01591   if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01592   if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
01593     return mdio_seterror(MDIO_WRONGFORMAT);
01594 
01595   // Read the header
01596   if (trx_header(mf) < 0) return -1;
01597 
01598   // We need some data from the trX header
01599   hdr = mf->trx;
01600   if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
01601 
01602   if (hdr->box_size) { // XXX need to check value of box_size!!
01603         fseek(mf->f, mf->prec * 9, SEEK_CUR);
01604   }
01605 
01606   if (hdr->vir_size) {
01607         fseek(mf->f, mf->prec * 9, SEEK_CUR);
01608   }
01609 
01610   if (hdr->pres_size) {
01611         fseek(mf->f, mf->prec * 9, SEEK_CUR);
01612   }
01613 
01614   //Almost all the I/O will come from the coordinates/velocities... Which we can just skip.
01615   if (hdr->x_size) {
01616         fseek(mf->f, mf->prec * hdr->natoms * 3, SEEK_CUR);
01617   }
01618 
01619   if (hdr->v_size) {
01620     fseek(mf->f, mf->prec * hdr->natoms * 3, SEEK_CUR);
01621   }
01622 
01623   if (hdr->f_size) {
01624     fseek(mf->f, mf->prec * hdr->natoms * 3, SEEK_CUR);
01625   }
01626 
01627   return mdio_seterror(MDIO_SUCCESS);
01628 }
01629 
01630 
01631 // writes an int in big endian. Returns GMX_SUCCESS
01632 // on success or a negative number on error.
01633 static int put_trx_int(md_file *mf, int y) {
01634       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01635 
01636       // sanity check.
01637       if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01638 
01639       if (mf->rev) swap4_aligned(&y, 1);
01640       if (fwrite(&y, 4, 1, mf->f) != 1)
01641     return mdio_seterror(MDIO_IOERROR);
01642 
01643   return mdio_seterror(MDIO_SUCCESS);
01644 }
01645 
01646 // writes a real in big-endian. Returns GMX_SUCCESS
01647 // on success or a negative number on error.
01648 static int put_trx_real(md_file *mf, float y) {
01649       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01650 
01651       if (mf->rev) swap4_aligned(&y, 1);
01652       if (fwrite(&y, 4, 1, mf->f) != 1)
01653         return mdio_seterror(MDIO_IOERROR);
01654 
01655       return mdio_seterror(MDIO_SUCCESS);
01656 }
01657 
01658 
01659 // writes an xdr encoded string. Returns GMX_SUCCESS
01660 // on success or a negative number on error.
01661 static int put_trx_string(md_file *mf, const char *s) {
01662         if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
01663         
01664         // write: size of object, string length, string data
01665         size_t len = strlen(s);
01666         if ( put_trx_int(mf, len+1)
01667              || put_trx_int(mf, len)
01668              || (fwrite(s, len, 1, mf->f) != 1))
01669           return mdio_seterror(MDIO_IOERROR);
01670 
01671         return mdio_seterror(MDIO_SUCCESS);
01672 }
01673 
01674 
01675 // xtc_int() - reads an integer from an xtc file
01676 static int xtc_int(md_file *mf, int *i) {
01677         unsigned char c[4];
01678 
01679         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01680         // sanity check.
01681         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01682 
01683         if (fread(c, 1, 4, mf->f) != 4) {
01684                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01685                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01686                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01687         }
01688 
01689         if (i) *i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01690         return mdio_seterror(MDIO_SUCCESS);
01691 }
01692 
01693 
01694 // xtc_float() - reads a float from an xtc file
01695 static int xtc_float(md_file *mf, float *f) {
01696         unsigned char c[4];
01697         int i;
01698 
01699         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01700 
01701         if (fread(c, 1, 4, mf->f) != 4) {
01702                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01703                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01704                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01705         }
01706 
01707         if (f) {
01708                 // By reading the number in as an integer and then
01709                 // copying it to a floating point number we can
01710                 // ensure proper endianness
01711                 i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01712                 memcpy(f, &i, 4);
01713         }
01714         return mdio_seterror(MDIO_SUCCESS);
01715 }
01716 
01717 
01718 // xtc_data() - reads a specific amount of data from an xtc
01719 // file using the xdr format.
01720 static int xtc_data(md_file *mf, char *buf, int len) {
01721         if (!mf || len < 1) return mdio_seterror(MDIO_BADPARAMS);
01722   size_t slen = (size_t)len;
01723         if (buf) {
01724                 if (fread(buf, 1, slen, mf->f) != slen) {
01725                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01726                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01727                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01728                 }
01729                 if (len % 4) {
01730                         if (fseek(mf->f, 4 - (len % 4), SEEK_CUR)) {
01731                                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01732                                 if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01733                                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01734                         }
01735                 }
01736         }
01737         else {
01738                 int newlen;
01739                 newlen = len;
01740                 if (len % 4) newlen += (4 - (len % 4));
01741                 if (fseek(mf->f, newlen, SEEK_CUR)) {
01742                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01743                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01744                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01745                 }
01746         }
01747         return len;
01748 }
01749 
01750 
01751 // xtc_timestep() - reads a timestep from an .xtc file.
01752 static int xtc_timestep(md_file *mf, md_ts *ts) {
01753         int n;
01754         float f, x[3], y[3], z[3];
01755 
01756         int size = 0; // explicitly initialized to zero.
01757         float precision;
01758 
01759         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01760         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
01761         if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);
01762 
01763         // Check magic number
01764         if (xtc_int(mf, &n) < 0) return -1;
01765         if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01766 
01767         // Get number of atoms
01768         if (xtc_int(mf, &n) < 0) return -1;
01769         ts->natoms = n;
01770 
01771         // Get the simulation step
01772         if (xtc_int(mf, &n) < 0) return -1;
01773         ts->step = n;
01774 
01775         // Get the time value
01776         if (xtc_float(mf, &f) < 0) return -1;
01777         ts->time = f;
01778 
01779         // Read the basis vectors of the box
01780   if ( (xtc_float(mf, &x[0]) < 0) ||
01781        (xtc_float(mf, &x[1]) < 0) ||
01782        (xtc_float(mf, &x[2]) < 0) ||
01783        (xtc_float(mf, &y[0]) < 0) ||
01784        (xtc_float(mf, &y[1]) < 0) ||
01785        (xtc_float(mf, &y[2]) < 0) ||
01786        (xtc_float(mf, &z[0]) < 0) ||
01787        (xtc_float(mf, &z[1]) < 0) ||
01788        (xtc_float(mf, &z[2]) < 0) )
01789     return -1;
01790   // Allocate the box and convert the vectors.
01791   ts->box = (md_box *) malloc(sizeof(md_box));
01792   if (mdio_readbox(ts->box, x, y, z) < 0) {
01793     free(ts->box);
01794     ts->box = NULL;
01795     return -1;
01796   }
01797 
01798         ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
01799         if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
01800         n = xtc_3dfcoord(mf, ts->pos, &size, &precision);
01801         if (n < 0) return -1;
01802 
01803         /* Now we're left with the job of scaling... */
01804         for (n = 0; n < ts->natoms * 3; n++)
01805                 ts->pos[n] *= ANGS_PER_NM;
01806 
01807         return mdio_seterror(MDIO_SUCCESS);
01808 }
01809 //This method returns 0 if there is no error, but its not the header
01810 //1 if it IS the header
01811 //-1 if there is a I/O error.
01812 //This is translated from libxdrf.cpp from GROMACS.
01813 static int xtc_at_header_start(md_file *mf, int natoms, int* timestep)
01814 {
01815     int       i_inp[3];
01816     float     f_inp[10];
01817     int       i;
01818     long off;
01819 
01820 
01821     if ((off = ftell(mf->f)) < 0)
01822     {
01823         return -1;
01824     }
01825     /* read magic natoms and timestep */
01826     for (i = 0; i < 3; i++)
01827     {
01828         if (xtc_int(mf, &(i_inp[i])) < 0)
01829         {
01830                 fseek(mf->f, off, SEEK_SET);
01831                 return -1;
01832       }
01833     }
01834     /* quick return */
01835     if (i_inp[0] != XTC_MAGIC)
01836     {
01837                 //Advance the file pointer to the next int.
01838         if (fseek(mf->f, off + 4, SEEK_SET))
01839         {
01840             return -1;
01841         }
01842         return 0;
01843     }
01844     /* read time and box */
01845     for (i = 0; i < 10; i++)
01846     {
01847         if (xtc_float(mf, &(f_inp[i])) < 0)
01848         {
01849             fseek(mf->f, off, SEEK_SET);
01850             return -1;
01851         }
01852     }
01853     /* Make a rigourous check to see if we are in the beggining of a header
01854        Hopefully there are no ambiguous cases */
01855     /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
01856        right triangle and that the first element must be nonzero unless the entire matrix is zero
01857      */
01858     if (i_inp[1] == natoms
01859         && ((f_inp[1] != 0 && f_inp[6] == 0) || (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
01860     {
01861         if (fseek(mf->f, off, SEEK_SET))
01862         {
01863             return -1;
01864         }
01865         *timestep = i_inp[2];
01866         return 1;
01867     }
01868     //Advance the file pointer to the next int.
01869     if (fseek(mf->f, off + 4, SEEK_SET))
01870     {
01871         return -1;
01872     }
01873     return 0;
01874 }
01875 
01876 // xtc_timeskip() - reads a timestep from an .xtc file, and just ignores it.
01877 static int xtc_timeskip(md_file *mf, md_ts *ts) {
01878         int n, i;
01879         float f, x[3], y[3], z[3];
01880         long fpos, epos;
01881         long low;
01882         int size = 0; // explicitly initialized to zero.
01883         float precision;
01884 
01885         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01886         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
01887         if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);
01888 
01889         // Check magic number
01890         if (xtc_int(mf, &n) < 0) return -1;
01891         if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01892 
01893         // Get number of atoms
01894         if (xtc_int(mf, &n) < 0) return -1;
01895         ts->natoms = n;
01896 
01897         // Get the simulation step
01898         if (xtc_int(mf, &n) < 0) return -1;
01899         ts->step = n;
01900         low = (long(1.1 * ts->natoms * 3) / 4) * 4; //Low bound estimate for where the next header will be found.
01901         //The 1.1 is safe but empirical. Typically, you find the next magic int somewhere in the 1.15-1.25 range.
01902         fpos = ftell(mf->f); // Current position
01903         fseek(mf->f, 0, SEEK_END);
01904         epos = ftell(mf->f); // End position
01905         //In case we would overrun the end of the file, we just need to return end of file.
01906         if (low > (epos - fpos)) {
01907                 return mdio_seterror(MDIO_EOF);
01908         }
01909         fseek(mf->f, fpos + low, SEEK_SET);
01910         for (i = 0; i < low; i++) {
01911                 switch(xtc_at_header_start(mf, ts->natoms, &n)){
01912                         case -1:
01913                                 return mdio_seterror(MDIO_IOERROR);
01914                         case 1:
01915                                 return mdio_seterror(MDIO_SUCCESS);
01916                         case 0:
01917                                 break;
01918                 }
01919         }
01920         //Something went very wrong if we made it here.
01921         return mdio_seterror(MDIO_UNKNOWNERROR);
01922 }
01923 
01925 // This algorithm is an implementation of the 3dfcoord algorithm
01926 // written by Frans van Hoesel (hoesel@chem.rug.nl) as part of the
01927 // Europort project in 1995.
01929 
01930 // integer table used in decompression
01931 static int xtc_magicints[] = {
01932         0, 0, 0, 0, 0, 0, 0, 0, 0,8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
01933         80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
01934         1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, 16384,
01935         20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031, 131072,
01936         165140, 208063, 262144, 330280, 416127, 524287, 660561, 832255,
01937         1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
01938         5284491, 6658042, 8388607, 10568983, 13316085, 16777216 };
01939 
01940 #define FIRSTIDX 9
01941 /* note that magicints[FIRSTIDX-1] == 0 */
01942 #define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))
01943 
01944 
01945 // returns the number of bits in the binary expansion of
01946 // the given integer.
01947 static int xtc_sizeofint(int size) {
01948         unsigned int num = 1;
01949   unsigned int ssize = (unsigned int)size;
01950         int nbits = 0;
01951 
01952         while (ssize >= num && nbits < 32) {
01953                 nbits++;
01954                 num <<= 1;
01955         }
01956         return nbits;
01957 }
01958 
01959 // calculates the number of bits a set of integers, when compressed,
01960 // will take up.
01961 static int xtc_sizeofints(int nints, unsigned int *sizes) {
01962         int i;
01963   unsigned int num;
01964         unsigned int nbytes, nbits, bytes[32], bytecnt, tmp;
01965         nbytes = 1;
01966         bytes[0] = 1;
01967         nbits = 0;
01968         for (i=0; i < nints; i++) {     
01969                 tmp = 0;
01970                 for (bytecnt = 0; bytecnt < nbytes; bytecnt++) {
01971                         tmp = bytes[bytecnt] * sizes[i] + tmp;
01972                         bytes[bytecnt] = tmp & 0xff;
01973                         tmp >>= 8;
01974                 }
01975                 while (tmp != 0) {
01976                         bytes[bytecnt++] = tmp & 0xff;
01977                         tmp >>= 8;
01978                 }
01979                 nbytes = bytecnt;
01980         }
01981         num = 1;
01982         nbytes--;
01983         while (bytes[nbytes] >= num) {
01984                 nbits++;
01985                 num *= 2;
01986         }
01987         return nbits + nbytes * 8;
01988 }
01989 
01990 // reads bits from a buffer.    
01991 static int xtc_receivebits(int *buf, int nbits) {
01992         int cnt, num; 
01993         unsigned int lastbits, lastbyte;
01994         unsigned char * cbuf;
01995         int mask = (1 << nbits) -1;
01996 
01997         cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
01998         cnt = buf[0];
01999         lastbits = (unsigned int) buf[1];
02000         lastbyte = (unsigned int) buf[2];
02001 
02002         num = 0;
02003         while (nbits >= 8) {
02004                 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
02005                 num |=  (lastbyte >> lastbits) << (nbits - 8);
02006                 nbits -=8;
02007         }
02008         if (nbits > 0) {
02009                 if (lastbits < (unsigned int)nbits) {
02010                         lastbits += 8;
02011                         lastbyte = (lastbyte << 8) | cbuf[cnt++];
02012                 }
02013                 lastbits -= nbits;
02014                 num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
02015         }
02016         num &= mask;
02017         buf[0] = cnt;
02018         buf[1] = lastbits;
02019         buf[2] = lastbyte;
02020         return num; 
02021 }
02022 
02023 // decompresses small integers from the buffer
02024 // sizes parameter has to be non-zero to prevent divide-by-zero
02025 static void xtc_receiveints(int *buf, const int nints, int nbits,
02026                         unsigned int *sizes, int *nums) {
02027         int bytes[32];
02028         int i, j, nbytes, p, num;
02029 
02030         bytes[1] = bytes[2] = bytes[3] = 0;
02031         nbytes = 0;
02032         while (nbits > 8) {
02033                 bytes[nbytes++] = xtc_receivebits(buf, 8);
02034                 nbits -= 8;
02035         }
02036         if (nbits > 0) {
02037                 bytes[nbytes++] = xtc_receivebits(buf, nbits);
02038         }
02039         for (i = nints-1; i > 0; i--) {
02040                 num = 0;
02041                 for (j = nbytes-1; j >=0; j--) {
02042                         num = (num << 8) | bytes[j];
02043                         p = num / sizes[i];
02044                         bytes[j] = p;
02045                         num = num - p * sizes[i];
02046                 }
02047                 nums[i] = num;
02048         }
02049         nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
02050 }
02051 
02052 // function that actually reads compressed coordinates    
02053 static int xtc_3dfcoord(md_file *mf, float *fp, int *size, float *precision) {
02054         static int *ip = NULL;
02055         static int oldsize;
02056         static int *buf;
02057 
02058         int minint[3], maxint[3], *lip;
02059         int smallidx;
02060         unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
02061         int flag, k;
02062         int small, smaller, i, is_smaller, run;
02063         float *lfp;
02064         int tmp, *thiscoord,  prevcoord[3];
02065 
02066         int bufsize, lsize;
02067         unsigned int bitsize;
02068         float inv_precision;
02069 
02070         /* avoid uninitialized data compiler warnings */
02071         bitsizeint[0] = 0;
02072         bitsizeint[1] = 0;
02073         bitsizeint[2] = 0;
02074 
02075         if (xtc_int(mf, &lsize) < 0) return -1;
02076 
02077         if (*size != 0 && lsize != *size) return mdio_seterror(MDIO_BADFORMAT);
02078 
02079         *size = lsize;
02080         size3 = *size * 3;
02081         if (*size <= 9) {
02082                 for (i = 0; i < *size; i++) {
02083                         if (xtc_float(mf, fp + (3 * i)) < 0) return -1;
02084                         if (xtc_float(mf, fp + (3 * i) + 1) < 0) return -1;
02085                         if (xtc_float(mf, fp + (3 * i) + 2) < 0) return -1;
02086                 }
02087                 return *size;
02088         }
02089         xtc_float(mf, precision);
02090         if (ip == NULL) {
02091                 ip = (int *)malloc(size3 * sizeof(*ip));
02092                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
02093                 bufsize = (int) (size3 * 1.2);
02094                 buf = (int *)malloc(bufsize * sizeof(*buf));
02095                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
02096                 oldsize = *size;
02097         } else if (*size > oldsize) {
02098                 ip = (int *)realloc(ip, size3 * sizeof(*ip));
02099                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
02100                 bufsize = (int) (size3 * 1.2);
02101                 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
02102                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
02103                 oldsize = *size;
02104         }
02105         buf[0] = buf[1] = buf[2] = 0;
02106 
02107         xtc_int(mf, &(minint[0]));
02108         xtc_int(mf, &(minint[1]));
02109         xtc_int(mf, &(minint[2]));
02110 
02111         xtc_int(mf, &(maxint[0]));
02112         xtc_int(mf, &(maxint[1]));
02113         xtc_int(mf, &(maxint[2]));
02114                 
02115         sizeint[0] = maxint[0] - minint[0]+1;
02116         sizeint[1] = maxint[1] - minint[1]+1;
02117         sizeint[2] = maxint[2] - minint[2]+1;
02118         
02119         /* check if one of the sizes is to big to be multiplied */
02120         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
02121                 bitsizeint[0] = xtc_sizeofint(sizeint[0]);
02122                 bitsizeint[1] = xtc_sizeofint(sizeint[1]);
02123                 bitsizeint[2] = xtc_sizeofint(sizeint[2]);
02124                 bitsize = 0; /* flag the use of large sizes */
02125         } else {
02126                 bitsize = xtc_sizeofints(3, sizeint);
02127         }
02128 
02129         xtc_int(mf, &smallidx);
02130         smaller = xtc_magicints[FIRSTIDX > smallidx - 1 ? FIRSTIDX : smallidx - 1] / 2;
02131         small = xtc_magicints[smallidx] / 2;
02132         sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx];
02133 
02134         /* check for zero values that would yield corrupted data */
02135         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
02136                 printf("XTC corrupted, sizesmall==0 (case 1)\n");
02137                 return -1;
02138         }
02139 
02140 
02141         /* buf[0] holds the length in bytes */
02142         if (xtc_int(mf, &(buf[0])) < 0) return -1;
02143 
02144         if (xtc_data(mf, (char *) &buf[3], (int) buf[0]) < 0) return -1;
02145 
02146         buf[0] = buf[1] = buf[2] = 0;
02147 
02148         lfp = fp;
02149         inv_precision = 1.0f / (*precision);
02150         run = 0;
02151         i = 0;
02152         lip = ip;
02153         while (i < lsize) {
02154                 thiscoord = (int *)(lip) + i * 3;
02155 
02156                 if (bitsize == 0) {
02157                         thiscoord[0] = xtc_receivebits(buf, bitsizeint[0]);
02158                         thiscoord[1] = xtc_receivebits(buf, bitsizeint[1]);
02159                         thiscoord[2] = xtc_receivebits(buf, bitsizeint[2]);
02160                 } else {
02161                         xtc_receiveints(buf, 3, bitsize, sizeint, thiscoord);
02162                 }
02163 
02164                 i++;
02165                 thiscoord[0] += minint[0];
02166                 thiscoord[1] += minint[1];
02167                 thiscoord[2] += minint[2];
02168 
02169                 prevcoord[0] = thiscoord[0];
02170                 prevcoord[1] = thiscoord[1];
02171                 prevcoord[2] = thiscoord[2];
02172  
02173 
02174                 flag = xtc_receivebits(buf, 1);
02175                 is_smaller = 0;
02176                 if (flag == 1) {
02177                         run = xtc_receivebits(buf, 5);
02178                         is_smaller = run % 3;
02179                         run -= is_smaller;
02180                         is_smaller--;
02181                 }
02182                 if (run > 0) {
02183                         thiscoord += 3;
02184                         for (k = 0; k < run; k+=3) {
02185                                 xtc_receiveints(buf, 3, smallidx, sizesmall, thiscoord);
02186                                 i++;
02187                                 thiscoord[0] += prevcoord[0] - small;
02188                                 thiscoord[1] += prevcoord[1] - small;
02189                                 thiscoord[2] += prevcoord[2] - small;
02190                                 if (k == 0) {
02191                                         /* interchange first with second atom for better
02192                                          * compression of water molecules
02193                                          */
02194                                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
02195                                         prevcoord[0] = tmp;
02196                                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
02197                                         prevcoord[1] = tmp;
02198                                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
02199                                         prevcoord[2] = tmp;
02200                                         *lfp++ = prevcoord[0] * inv_precision;
02201                                         *lfp++ = prevcoord[1] * inv_precision;
02202                                         *lfp++ = prevcoord[2] * inv_precision;
02203 
02204                                         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
02205                                                 printf("XTC corrupted, sizesmall==0 (case 2)\n");
02206                                                 return -1;
02207                                         }
02208 
02209                                 } else {
02210                                         prevcoord[0] = thiscoord[0];
02211                                         prevcoord[1] = thiscoord[1];
02212                                         prevcoord[2] = thiscoord[2];
02213                                 }
02214                                 *lfp++ = thiscoord[0] * inv_precision;
02215                                 *lfp++ = thiscoord[1] * inv_precision;
02216                                 *lfp++ = thiscoord[2] * inv_precision;
02217                         }
02218                 } else {
02219                         *lfp++ = thiscoord[0] * inv_precision;
02220                         *lfp++ = thiscoord[1] * inv_precision;
02221                         *lfp++ = thiscoord[2] * inv_precision;          
02222                 }
02223                 smallidx += is_smaller;
02224                 if (is_smaller < 0) {
02225                         small = smaller;
02226                         if (smallidx > FIRSTIDX) {
02227                                 smaller = xtc_magicints[smallidx - 1] /2;
02228                         } else {
02229                                 smaller = 0;
02230                         }
02231                 } else if (is_smaller > 0) {
02232                         smaller = small;
02233                         small = xtc_magicints[smallidx] / 2;
02234                 }
02235                 sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx] ;
02236         }
02237         return 1;
02238 }
02239 
02240 
02241 
02242 #endif

Generated on Wed May 22 03:10:55 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002