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: johns $       $Locker:  $             $State: Exp $
00012  *      $Revision: 1.36 $       $Date: 2020/01/09 16:21:28 $
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 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 
00216 // .g96 file functions
00217 static int g96_header(md_file *, char *, int, float *);
00218 static int g96_timestep(md_file *, md_ts *);
00219 static int g96_rec(md_file *, md_atom *);
00220 static int g96_countatoms(md_file *);
00221 
00222 
00223 // .xtc file functions
00224 static int xtc_int(md_file *, int *);
00225 static int xtc_float(md_file *, float *);
00226 /* 
00227 static int xtc_receivebits(int *, int);
00228 static void xtc_receiveints(int *, int, int, const unsigned *, int *);
00229 */
00230 static int xtc_timestep(md_file *, md_ts *);
00231 static int xtc_3dfcoord(md_file *, float *, int *, float *);
00232 
00233 
00234 // Error reporting functions
00235 static int mdio_errno(void);
00236 static const char *mdio_errmsg(int);
00237 static int mdio_seterror(int);
00238 
00239 
00240 // Miscellaneous functions
00241 static int strip_white(char *);
00242 static int mdio_readline(md_file *, char *, int, int = 1);
00243 static int mdio_tsfree(md_ts *, int = 0);
00244 static int mdio_readbox(md_box *, float *, float *, float *);
00245 
00246 
00247 
00248 static int xtc_receivebits(int *, int);
00249 
00250 // Error descriptions for mdio_errno
00251 static const char *mdio_errdescs[] = {
00252         "no error",
00253         "file does not match format",
00254         "unexpected end-of-file reached",
00255         "function called with bad parameters",
00256         "file i/o error",
00257         "unsupported precision",
00258         "out of memory",
00259         "cannot open file",
00260         "bad file extension",
00261         "unknown file format",
00262         "cannot close file",
00263         "wrong file format for this function",
00264         "binary i/o error: sizeof(int) != 4",
00265         NULL
00266 };
00267 
00270 static inline int host_is_little_endian(void) 
00271 {
00272   const union { unsigned char c[4]; unsigned int i; } 
00273   fixed = { { 0x10 , 0x20 , 0x40 , 0x80 } };
00274   const unsigned int i = 0x80402010U;
00275         
00276   if (fixed.i == i) {
00277     return 1;
00278   }
00279   return 0;
00280 }
00281 
00282 
00283 
00284 // Open a molecular dynamics file. The second parameter specifies
00285 // the format of the file. If it is zero, the format is determined
00286 // from the file extension.
00287 md_file *mdio_open(const char *fn, const int fmt, const int rw) {
00288         md_file *mf;
00289 
00290         if (!fn) {
00291                 mdio_seterror(MDIO_BADPARAMS);
00292                 return NULL;
00293         }
00294 
00295         // Allocate memory
00296         mf = (md_file *) malloc(sizeof(md_file));
00297         if (!mf) {
00298                 mdio_seterror(MDIO_BADMALLOC);
00299                 return NULL;
00300         }
00301 
00302         // Zero out the structure
00303         memset(mf, 0, sizeof(md_file));
00304 
00305         // Determine the file type from the extension
00306         if (!fmt) {
00307                 char *p;
00308                 int n;
00309 
00310                 // Seek to the extension part of the filename
00311                 for (p = (char *) &fn[strlen(fn) - 1]; *p != '.' && p > fn; p--);
00312                 if (p == fn) {
00313                         free(mf);
00314                         mdio_seterror(MDIO_BADEXTENSION);
00315                         return NULL;
00316                 }
00317 
00318                 // Check the extension against known extensions
00319                 for (n = 1; mdio_fmtexts[n]; n++)
00320                         if (!strcasecmp(p, mdio_fmtexts[n])) break;
00321 
00322                 // If !mdio_fmtexts[n], we failed (unknown ext)
00323                 if (!mdio_fmtexts[n]) {
00324                         free(mf);
00325                         mdio_seterror(MDIO_UNKNOWNFMT);
00326                         return NULL;
00327                 }
00328 
00329                 // All set
00330                 mf->fmt = n;
00331         }
00332         else {
00333                 mf->fmt = fmt;
00334         }
00335 
00336         // Differentiate between binary and ascii files. Also,
00337         // .trX files need a header information structure allocated.
00338         switch (mf->fmt) {
00339     case MDFMT_GRO:
00340         case MDFMT_G96: /* fallthrough */
00341         if (rw) 
00342             mf->f = fopen(fn, "wt");
00343         else
00344             mf->f = fopen(fn, "rt");
00345 
00346                 break;
00347         case MDFMT_TRR:
00348         case MDFMT_TRJ: /* fallthrough */
00349                 // Allocate the trx header data struct
00350                 mf->trx = (trx_hdr *) malloc(sizeof(trx_hdr));
00351                 if (!mf->trx) {
00352                         free(mf);
00353                         mdio_seterror(MDIO_BADMALLOC);
00354                         return NULL;
00355                 }
00356                 memset(mf->trx, 0, sizeof(trx_hdr));
00357         case MDFMT_XTC:  /* fallthrough */
00358                 // Finally, open the file
00359         if (rw)
00360             mf->f = fopen(fn, "wb");
00361         else
00362             mf->f = fopen(fn, "rb");
00363 
00364                 break;
00365         default:
00366                 free(mf);
00367                 mdio_seterror(MDIO_UNKNOWNFMT);
00368                 return NULL;
00369         }
00370 
00371         // Check for opening error
00372         if (!mf->f) {
00373                 if (mf->trx) free(mf->trx);
00374                 free(mf);
00375                 mdio_seterror(MDIO_CANTOPEN);
00376                 return NULL;
00377         }
00378 
00379         // File is opened, we're all set!
00380         mdio_seterror(MDIO_SUCCESS);
00381         return mf;
00382 }
00383 
00384 
00385 // Closes a molecular dynamics file.
00386 static int mdio_close(md_file *mf) {
00387         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00388 
00389         if (fclose(mf->f) == EOF) return mdio_seterror(MDIO_CANTCLOSE);
00390 
00391         // Free the dynamically allocated memory
00392         if (mf->trx) free(mf->trx);
00393         free(mf);
00394         return mdio_seterror(MDIO_SUCCESS);
00395 }
00396 
00397 
00398 // Returns the last error code reported by any of the mdio functions
00399 static int mdio_errno(void) {
00400         return mdio_errcode;
00401 }
00402 
00403 
00404 // Returns a textual message regarding an mdio error code
00405 static const char *mdio_errmsg(int n) {
00406         if (n < 0 || n > MDIO_MAX_ERRVAL) return (char *) "unknown error";
00407         else return mdio_errdescs[n];
00408 }
00409 
00410 
00411 // Sets the error code and returns an appropriate return value
00412 // for the calling function to return to its parent
00413 static int mdio_seterror(int code) {
00414         mdio_errcode = code;
00415         return code ? -1 : 0;
00416 }
00417 
00418 
00419 // Reads a line from the text file, strips leading/trailing whitespace
00420 // and newline, checks for errors, and returns the number of characters
00421 // in the string on success or -1 on error.
00422 static int mdio_readline(md_file *mf, char *buf, int n, int strip) {
00423         if (!buf || n < 1 || !mf) return mdio_seterror(MDIO_BADPARAMS);
00424 
00425         // Read the line
00426         fgets(buf, n, mf->f);
00427 
00428         // End of file reached?
00429         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
00430 
00431         // File I/O error?
00432         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
00433 
00434         // comment line?
00435         if (buf[0] == '#') return mdio_readline(mf,buf,n,strip);
00436 
00437         // Strip whitespace
00438         if (strip) strip_white(buf);
00439 
00440         return strlen(buf);
00441 }
00442 
00443 
00444 // Strips leading and trailing whitespace from a string. Tabs,
00445 // spaces, newlines and carriage returns are stripped. Example:
00446 // "\n   hello\t \r" becomes "hello".
00447 static int strip_white(char *buf) {
00448         int i, j, k;
00449 
00450         // Protect against NULL pointer
00451         if (!buf) return -1;
00452         if (!strlen(buf)) return -1;
00453 
00454         // Kill trailing whitespace first
00455         for (i = strlen(buf) - 1;
00456              buf[i] == ' ' || buf[i] == '\t' ||
00457              buf[i] == '\n' || buf[i] == '\r';
00458              i--)
00459                 buf[i] = 0;
00460 
00461         // Skip past leading whitespace
00462         for (i = 0; buf[i] == ' ' || buf[i] == '\t' ||
00463              buf[i] == '\n' || buf[i] == '\r'; i++);
00464         if (i) {
00465                 k = 0;
00466                 for (j = i; buf[j]; j++)
00467                         buf[k++] = buf[j];
00468                 buf[k] = 0;
00469         }
00470 
00471         return strlen(buf);
00472 }
00473 
00474 
00475 // Frees the memory allocated in a ts structure. The holderror
00476 // parameter defaults to zero. Programs that are calling this
00477 // function because of an error reported by another function should
00478 // set holderror so that mdio_tsfree() does not overwrite the error
00479 // code with mdio_seterror().
00480 static int mdio_tsfree(md_ts *ts, int holderror) {
00481         if (!ts) {
00482                 if (holderror) return -1;
00483                 else return mdio_seterror(MDIO_BADPARAMS);
00484         }
00485 
00486         if (ts->pos && ts->natoms > 0) free(ts->pos);
00487 
00488   if (ts->box) free(ts->box);
00489 
00490         if (holderror) return -1;
00491         else return mdio_seterror(MDIO_SUCCESS);
00492 }
00493 
00494 
00495 // Converts box basis vectors to A, B, C, alpha, beta, and gamma.  
00496 // Stores values in md_box struct, which should be allocated before calling
00497 // this function.
00498 static int mdio_readbox(md_box *box, float *x, float *y, float *z) {
00499   float A, B, C;
00500 
00501   if (!box) {
00502     return mdio_seterror(MDIO_BADPARAMS);
00503   }
00504 
00505   // A, B, C are the lengths of the x, y, z vectors, respectively
00506   A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] ) * ANGS_PER_NM;
00507   B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] ) * ANGS_PER_NM;
00508   C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] ) * ANGS_PER_NM;
00509   if ((A<=0) || (B<=0) || (C<=0)) {
00510     /* Use zero-length box size and set angles to 90. */
00511     box->A = box->B = box->C = 0;
00512     box->alpha = box->beta = box->gamma = 90;
00513   } else {
00514     box->A = A;
00515     box->B = B;
00516     box->C = C;
00517   
00518     // gamma, beta, alpha are the angles between the x & y, x & z, y & z
00519     // vectors, respectively
00520     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;
00521     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;
00522     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; 
00523   }
00524   return mdio_seterror(MDIO_SUCCESS);
00525 }
00526 
00527 
00528 // Reads the header of a file (format independent)
00529 static int mdio_header(md_file *mf, md_header *mdh) {
00530         int n;
00531         if (!mf || !mdh) return mdio_seterror(MDIO_BADPARAMS);
00532         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00533 
00534         switch (mf->fmt) {
00535         case MDFMT_GRO:
00536                 if (gro_header(mf, mdh->title, MAX_MDIO_TITLE,
00537                 &mdh->timeval, &mdh->natoms, 1) < 0)
00538                         return -1;
00539                 return 0;
00540 
00541         case MDFMT_TRR: 
00542         case MDFMT_TRJ: /* fallthrough */
00543                 if (trx_header(mf, 1) < 0) return -1;
00544                 mdh->natoms = mf->trx->natoms;
00545                 mdh->timeval = (float) mf->trx->t;
00546                 strncpy(mdh->title, mf->trx->title, MAX_MDIO_TITLE);
00547                 return 0;
00548 
00549         case MDFMT_G96:
00550                 if (g96_header(mf, mdh->title, MAX_MDIO_TITLE,
00551                 &mdh->timeval) < 0) return -1;
00552                 mdh->natoms = -1;
00553                 return 0;
00554 
00555         case MDFMT_XTC:
00556                 memset(mdh, 0, sizeof(md_header));
00557                 // Check magic number
00558                 if (xtc_int(mf, &n) < 0) return -1;
00559                 if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
00560 
00561                 // Get number of atoms
00562                 if (xtc_int(mf, &n) < 0) return -1;
00563                 mdh->natoms = n;
00564                 rewind(mf->f);
00565                 return 0;
00566 
00567         default:
00568                 return mdio_seterror(MDIO_UNKNOWNFMT);
00569         }
00570 }
00571 
00572 
00573 // Reads in a timestep from a file (format independent)
00574 static int mdio_timestep(md_file *mf, md_ts *ts) {
00575         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00576         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00577 
00578         switch (mf->fmt) {
00579         case MDFMT_GRO:
00580                 return gro_timestep(mf, ts);
00581 
00582         case MDFMT_TRR:
00583         case MDFMT_TRJ: /* fallthrough */
00584                 return trx_timestep(mf, ts);
00585 
00586         case MDFMT_G96:
00587                 return g96_timestep(mf, ts);
00588 
00589         case MDFMT_XTC:
00590                 return xtc_timestep(mf, ts);
00591 
00592         default:
00593                 return mdio_seterror(MDIO_UNKNOWNFMT);
00594         }
00595 }
00596 
00597 
00598 
00599 static int g96_header(md_file *mf, char *title, int titlelen, float *timeval) {
00600         char buf[MAX_G96_LINE + 1];
00601         char *p;
00602 
00603         // Check parameters
00604         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00605 
00606         // The header consists of blocks. The title block
00607         // is mandatory, and a TIMESTEP block is optional.
00608         // Example:
00609         //
00610         // TITLE
00611         // Generated by trjconv :  t=  90.00000
00612         // more title info
00613         // .
00614         // .
00615         // .
00616         // END
00617         // .
00618         // .
00619         // .
00620 
00621         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00622         if (strcasecmp(buf, "TITLE")) return mdio_seterror(MDIO_BADFORMAT);
00623 
00624         // Read in the title itself
00625         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00626 
00627         // The timevalue can be included in the title string
00628         // after a "t=" prefix.
00629         if ((p = (char *) strstr(buf, "t="))) {
00630                 char *q = p;
00631                 *(q--) = 0;
00632 
00633                 // Skip the `t=' and strip whitespace from
00634                 // the resulting strings
00635                 p += 2;
00636                 strip_white(p);
00637                 strip_white(buf);
00638 
00639                 // Grab the timevalue from the title string
00640                 if (timeval) *timeval = (float) atof(p);
00641         }
00642         else {
00643                 // No timevalue - just copy the string and strip
00644                 // any leading/trailing whitespace
00645                 if (timeval) *timeval = 0;
00646                 strip_white(buf);
00647         }
00648 
00649         // Copy the title string
00650         if (title && titlelen) strncpy(title, buf, titlelen);
00651 
00652         // Now ignore subsequent title lines and get the END string
00653         while (strcasecmp(buf, "END"))
00654                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00655 
00656         // Done!
00657         return mdio_seterror(MDIO_SUCCESS);
00658 }
00659 
00660 
00661 // Used to determine the number of atoms in a g96 file, because
00662 // VMD needs to know this for some reason.
00663 static int g96_countatoms(md_file *mf) {
00664         char buf[MAX_G96_LINE + 1];
00665         int natoms;
00666         int n;
00667         long fpos;
00668         float lastf;
00669 
00670         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00671 
00672         fpos = ftell(mf->f);
00673 
00674         natoms = 0;
00675         for (;;) {
00676                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00677                         break;
00678                 n = sscanf(buf, "%*6c%*6c%*6c%*6c %*f %*f %f", &lastf);
00679                 if (n == 1) natoms++;
00680                 else {
00681                         strip_white(buf);
00682                         if (!strcasecmp(buf, "END")) break;
00683                 }
00684         }
00685 
00686         fseek(mf->f, fpos, SEEK_SET);
00687         return natoms;
00688 }
00689 
00690 
00691 // Reads an atom line from the G96 file
00692 static int g96_rec(md_file *mf, md_atom *ma) {
00693         char buf[MAX_G96_LINE + 1];
00694         char atomnum[7];
00695         int n;
00696 
00697         // Check parameters
00698         if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);
00699 
00700         // Read in a line, assuming it is an atom line
00701         do {
00702                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0) return -1;
00703         } while (buf[0] == '#' || strlen(buf) == 0);
00704 
00705         n = sscanf(buf, "%6c%6c%6c%6c %f %f %f",
00706                 ma->resid, ma->resname, ma->atomname, atomnum,
00707                 &ma->pos[0], &ma->pos[1], &ma->pos[2]);
00708         if (n == 7) {
00709                 atomnum[6] = 0;
00710                 ma->resid[6] = 0;
00711                 ma->resname[6] = 0;
00712                 ma->atomname[6] = 0;
00713 
00714                 strip_white(atomnum);
00715                 strip_white(ma->resid);
00716                 strip_white(ma->resname);
00717                 strip_white(ma->atomname);
00718 
00719                 ma->atomnum = atoi(atomnum);
00720 
00721                 ma->pos[0] *= ANGS_PER_NM;
00722                 ma->pos[1] *= ANGS_PER_NM;
00723                 ma->pos[2] *= ANGS_PER_NM;
00724 
00725                 return 0;
00726         }
00727 
00728         return mdio_seterror(MDIO_BADFORMAT);
00729 }
00730 
00731 
00732 // Reads a timestep from a G96 file and stores the data in
00733 // the generic md_ts structure. Returns 0 on success or a
00734 // negative number on error and sets mdio_errcode.
00735 static int g96_timestep(md_file *mf, md_ts *ts) {
00736         char            buf[MAX_G96_LINE + 1];
00737         char            stripbuf[MAX_G96_LINE + 1];
00738         float           pos[3], x[3], y[3], z[3], *currAtom;
00739         long            fpos;
00740         int             n, i, boxItems;
00741 
00742         // Check parameters
00743         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00744 
00745   // Allocate data space for the timestep, using the number of atoms
00746   // determined by open_g96_read().
00747         ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
00748         if (!ts->pos) {
00749                 return mdio_seterror(MDIO_BADMALLOC);
00750         }
00751   currAtom = ts->pos;
00752 
00753         // The timesteps follow the header in a fixed block
00754         // format:
00755         //
00756         // TIMESTEP
00757         //         <step number> <time value>
00758         // END
00759         // POSITIONRED
00760         //     <x float> <y float> <z float>
00761         //     .         .         .
00762         //     .         .         .
00763         //     .         .         .
00764         // END
00765         // VELOCITYRED
00766         //     <x float> <y float> <z float>
00767         //     .         .         .
00768         //     .         .         .
00769         //     .         .         .
00770         // END
00771         // BOX
00772         //     <x float> <y float> <z float>
00773         // END
00774         //
00775         // -----
00776         //
00777         // The TIMESTEP, VELOCITY and BOX blocks are optional.
00778         // Floats are written in 15.9 precision.
00779         //
00780         // Reference: GROMACS 2.0 user manual
00781         //            http://rugmd4.chem.rug.nl/~gmx/online2.0/g96.html
00782 
00783         // First, look for an (optional) title block and skip it
00784         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00785 
00786   if (!strcasecmp(buf, "TITLE")) {
00787     // skip over the text until we reach 'END'
00788     while (strcasecmp(buf, "END")) {
00789       if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00790     }
00791 
00792     // Read in the next line
00793     if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00794   }
00795 
00796         // Next, look for a timestep block
00797         if (!strcasecmp(buf, "TIMESTEP")) {
00798                 // Read in the value line
00799                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00800 
00801                 // Extract the time value and the timestep index
00802                 n = sscanf(buf, "%d %f", &ts->step, &ts->time);
00803                 if (n != 2) return mdio_seterror(MDIO_BADFORMAT);
00804 
00805                 // Read the "END" line
00806                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00807                 if (strcasecmp(buf, "END"))
00808                         return mdio_seterror(MDIO_BADFORMAT);
00809 
00810                 // Read in the next line
00811                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00812         }
00813         else {
00814                 // No timestep specified -- set to zero
00815                 ts->step = 0;
00816                 ts->time = 0;
00817         }
00818 
00819         // At this point a POSITION or POSITIONRED block
00820         // is REQUIRED by the format
00821         if (!strcasecmp(buf, "POSITIONRED")) {
00822 
00823     // So now we read in some atoms
00824     i = 0;
00825                 while (i < ts->natoms) {
00826                         // Read in an atom
00827                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00828                                 return -1;
00829  
00830       // We shouldn't reach the end yet
00831       if (!strcasecmp(buf, "END"))
00832         return mdio_seterror(MDIO_BADFORMAT);
00833 
00834                         // Get the x,y,z coordinates
00835                         n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
00836       
00837       // Ignore improperly formatted lines
00838                         if (n == 3) {
00839                                 pos[0] *= ANGS_PER_NM;
00840                                 pos[1] *= ANGS_PER_NM;
00841                                 pos[2] *= ANGS_PER_NM;
00842 
00843                                 // Copy the atom data into the array
00844                                 memcpy(currAtom, pos, sizeof(float) * 3);
00845         currAtom += 3;
00846         i++;
00847                         }
00848                 }
00849         }
00850         else if (!strcasecmp(buf, "POSITION") || !strcasecmp(buf, "REFPOSITION")) {
00851                 /*
00852                 char resnum[7];
00853                 char resname[7];
00854                 char atomname[7];
00855                 char atomnum[7];
00856                 */
00857 
00858                 // So now we read in some atoms
00859     i = 0;
00860                 while (i < ts->natoms) {
00861                         // Read in the first line
00862                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00863                                 return -1;
00864  
00865       // We shouldn't reach the end yet
00866       strcpy(stripbuf, buf);
00867       strip_white(stripbuf); 
00868       if (!strcasecmp(stripbuf, "END"))
00869         return mdio_seterror(MDIO_BADFORMAT);
00870 
00871                         // Get the x,y,z coordinates and name data
00872                         n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
00873                                 &pos[0], &pos[1], &pos[2]);
00874 
00875       // Ignore improperly formatted lines
00876                         if (n == 3) {
00877                                 pos[0] *= ANGS_PER_NM;
00878                                 pos[1] *= ANGS_PER_NM;
00879                                 pos[2] *= ANGS_PER_NM;
00880 
00881                                 // Copy the atom data into the linked list item
00882                                 memcpy(currAtom, pos, sizeof(float) * 3);
00883                                 currAtom += 3;
00884         i++;
00885                         }
00886                 }
00887         }
00888         else {
00889                 return mdio_seterror(MDIO_BADFORMAT);
00890         }
00891 
00892   // Read the END keyword
00893   if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00894     return -1;
00895   if (strcasecmp(buf, "END"))
00896     return mdio_seterror(MDIO_BADFORMAT);
00897 
00898         // ... another problem: there may or may not be a VELOCITY
00899         // block or a BOX block, so we need to read one line beyond
00900         // the POSITION block to determine this. If neither VEL. nor
00901         // BOX are present we've read a line too far and infringed
00902         // on the next timestep, so we need to keep track of the
00903         // position now for a possible fseek() later to backtrack.
00904         fpos = ftell(mf->f);
00905 
00906         // Now we must read in the velocities and the box, if present
00907         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00908     // It's okay if we end the file here; any other errors need to be
00909     // reported.
00910     if (mdio_errcode == MDIO_EOF) 
00911       return mdio_seterror(MDIO_SUCCESS);
00912     else 
00913       return -1;
00914   }
00915 
00916         // Is there a velocity block present ?
00917         if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
00918                 // Ignore all the coordinates - VMD doesn't use them
00919                 for (;;) {
00920                         if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00921                                 return -1;
00922                         if (!strcasecmp(buf, "END")) break;
00923                 }
00924 
00925                 // Again, record our position because we may need
00926                 // to fseek here later if we read too far.
00927                 fpos = ftell(mf->f);
00928 
00929                 // Go ahead and read the next line.
00930                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00931         }
00932 
00933         // Is there a box present ?
00934         if (!strcasecmp(buf, "BOX")) {
00935                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00936     boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
00937                &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
00938     if (boxItems == 3) {
00939       x[1] = x[2] = 0;
00940       y[0] = y[2] = 0;
00941       z[0] = z[1] = 0;
00942     }
00943     else if (boxItems != 9) 
00944       return mdio_seterror(MDIO_BADFORMAT);
00945 
00946     // Allocate the box and convert the vectors.
00947     ts->box = (md_box *) malloc(sizeof(md_box));
00948     if (mdio_readbox(ts->box, x, y, z) < 0) {
00949       free(ts->box);
00950       ts->box = NULL;
00951       return mdio_seterror(MDIO_BADFORMAT);
00952     }
00953 
00954                 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00955       free(ts->box);
00956       ts->box = NULL;
00957       return -1;
00958     }
00959                 if (strcasecmp(buf, "END")) {
00960       free(ts->box);
00961       ts->box = NULL;
00962                         return mdio_seterror(MDIO_BADFORMAT);
00963     }
00964         }
00965         else {
00966                 // We have read too far, so fseek back to the
00967                 // last known safe position so we don't return
00968                 // with the file pointer set infringing on the
00969                 // next timestep data.
00970                 fseek(mf->f, fpos, SEEK_SET);
00971         }
00972 
00973         // We're done!
00974         return mdio_seterror(MDIO_SUCCESS);
00975 }
00976 
00977 
00978 // Attempts to read header data from a GROMACS structure file
00979 // The GROMACS header format is as follows (fixed, 2 lines ASCII):
00980 // <title> [ n= <timevalue> ]
00981 //     <num atoms>
00982 static int gro_header(md_file *mf, char *title, int titlelen, float *timeval,
00983                int *natoms, int rewind) {
00984   char buf[MAX_GRO_LINE + 1];
00985   long fpos;
00986   char *p;
00987 
00988   // Check parameters
00989   if (!mf)
00990     return mdio_seterror(MDIO_BADPARAMS);
00991 
00992   // Get the current file position for rewinding later
00993   fpos = ftell(mf->f);
00994 
00995   // The header consists of 2 lines - get the first line
00996   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
00997 
00998   // The timevalue can be included in the title string
00999   // after a "t=" prefix.
01000   if ((p = (char *) strstr(buf, "t="))) {
01001     char *q = p;
01002     *(q--) = 0;
01003 
01004     // Skip the `t=' and strip whitespace from
01005     // the resulting strings
01006     p += 2;
01007     strip_white(p);
01008     strip_white(buf);
01009 
01010     // Grab the timevalue from the title string
01011     if (timeval) *timeval = (float) atof(p);
01012   } else {
01013     // No timevalue - just copy the string
01014     if (timeval) *timeval = 0;
01015   }
01016 
01017   // Copy the title string
01018   if (title && titlelen) strncpy(title, buf, titlelen);
01019 
01020   // Get the second line and grab the number of atoms
01021   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
01022 
01023   // Store the number of atoms
01024   if (natoms && (!(*natoms = atoi(buf))))
01025     return mdio_seterror(MDIO_BADFORMAT);
01026 
01027   // Now we rewind the file so that subsequent calls to
01028   // gro_timestep() will succeed. gro_timestep() requires
01029   // the header to be at the current file pointer.
01030   if (rewind)
01031     fseek(mf->f, fpos, SEEK_SET);
01032 
01033   return 0; // Done!
01034 }
01035 
01036 
01037 // Reads one atom record from a GROMACS file. Returns GMX_SUCCESS
01038 // on success or a negative number on error.
01039 //
01040 // Record format (one line, fixed):
01041 //    rrrrrRRRRRaaaaaAAAAA <x pos> <y pos> <z pos> <x vel> <y vel> <z vel>
01042 //
01043 //    r = residue number
01044 //    R = residue name
01045 //    a = atom name
01046 //    A = atom number
01047 //
01048 static int gro_rec(md_file *mf, md_atom *ma) {
01049   char buf[MAX_GRO_LINE + 1];
01050   char atomnum[6];
01051   char xposc[12], yposc[12], zposc[12];
01052   int n;
01053 
01054   if (!mf)
01055     return mdio_seterror(MDIO_BADPARAMS);
01056 
01057   do {
01058     if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0)
01059       return -1;
01060   } while (buf[0] == '#' || !strlen(buf));
01061 
01062   // Read in the fields
01063   n = sscanf(buf, "%5c%5c%5c%5c%8c%8c%8c", 
01064              ma->resid, ma->resname, ma->atomname, atomnum, 
01065              xposc, yposc, zposc);
01066 
01067   if (n != 7)
01068     return mdio_seterror(MDIO_BADFORMAT);
01069 
01070   // Null terminate the strings
01071   ma->resname[5] = 0;
01072   ma->atomname[5] = 0;
01073   ma->resid[5] = 0;
01074   atomnum[5] = 0;
01075   xposc[8] = 0;
01076   yposc[8] = 0;
01077   zposc[8] = 0;
01078  
01079   if ((sscanf(xposc, "%f", &ma->pos[0]) != 1) ||
01080       (sscanf(yposc, "%f", &ma->pos[1]) != 1) ||
01081       (sscanf(zposc, "%f", &ma->pos[2]) != 1)) {
01082     return mdio_seterror(MDIO_BADFORMAT);
01083   }
01084 
01085   // Convert strings to numbers
01086   strip_white(atomnum);
01087   ma->atomnum = atoi(atomnum);
01088 
01089   // Convert nanometers to angstroms
01090   ma->pos[0] *= ANGS_PER_NM;
01091   ma->pos[1] *= ANGS_PER_NM;
01092   ma->pos[2] *= ANGS_PER_NM;
01093 
01094   // Strip leading and trailing whitespace
01095   strip_white(ma->atomname);
01096   strip_white(ma->resname);
01097   strip_white(ma->resid);
01098 
01099   return 0;
01100 }
01101 
01102 
01103 // Reads in a timestep from a .gro file. Ignores the data
01104 // not needed for a timestep, so is a little faster than
01105 // calling gro_rec() for each atom. Also reads in the
01106 // header block.
01107 //
01108 static int gro_timestep(md_file *mf, md_ts *ts) {
01109         char buf[MAX_GRO_LINE + 1];
01110         long coord;
01111         int i, n, boxItems;
01112   float x[3], y[3], z[3];
01113   char xposc[12], yposc[12], zposc[12];
01114 
01115   if (!mf || !ts) 
01116     return mdio_seterror(MDIO_BADPARAMS);
01117 
01118   if (gro_header(mf, NULL, 0, &ts->time, &ts->natoms, 0) < 0)
01119     return -1;
01120 
01121   ts->pos = (float *) malloc(3 * sizeof(float) * ts->natoms);
01122   if (!ts->pos)
01123     return mdio_seterror(MDIO_BADMALLOC);
01124 
01125   coord = 0;
01126   for (i = 0; i < ts->natoms; i++) {
01127     if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01128       free(ts->pos);
01129       return -1;
01130     }
01131         
01132     n = sscanf(buf, "%*5c%*5c%*5c%*5c%8c%8c%8c", xposc, yposc, zposc);
01133     if (n != 3) 
01134       return mdio_seterror(MDIO_BADFORMAT);
01135 
01136     if ((sscanf(xposc, "%f", &ts->pos[coord    ]) != 1) ||
01137         (sscanf(yposc, "%f", &ts->pos[coord + 1]) != 1) ||
01138         (sscanf(zposc, "%f", &ts->pos[coord + 2]) != 1)) {
01139       return mdio_seterror(MDIO_BADFORMAT);
01140     }
01141 
01142     ts->pos[coord    ] *= ANGS_PER_NM;
01143     ts->pos[coord + 1] *= ANGS_PER_NM;
01144     ts->pos[coord + 2] *= ANGS_PER_NM;
01145 
01146     coord += 3;
01147   }
01148 
01149   // Read the box, stored as three vectors representing its edges
01150   if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01151     free(ts->pos);
01152     return -1;
01153   }
01154 
01155   boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f", 
01156              &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
01157 
01158   // File may only include three scalars for the box information -- if
01159   // that's the case, the box is orthoganal.
01160   if (boxItems == 3) {
01161     x[1] = x[2] = 0;
01162     y[0] = y[2] = 0;
01163     z[0] = z[1] = 0;
01164   } else if (boxItems != 9) {
01165     free(ts->pos);
01166     return -1;
01167   }
01168 
01169   // Allocate the box and convert the vectors.
01170   ts->box = (md_box *) malloc(sizeof(md_box));
01171   if (mdio_readbox(ts->box, x, y, z) < 0) {
01172     free(ts->pos);
01173     free(ts->box);
01174     ts->box = NULL;
01175     return -1;
01176   }
01177 
01178   return 0;
01179 }
01180 
01181 
01182 // Attempts to read header data from a .trX trajectory file
01183 //
01184 // The .trX header format is as follows:
01185 //
01186 //      4 bytes         - magic number (0x07C9)
01187 //      ...
01188 //
01189 static int trx_header(md_file *mf, int rewind) {
01190         int magic;
01191         trx_hdr *hdr;
01192         long fpos;
01193 
01194         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01195 
01196         // In case we need to rewind
01197         fpos = ftell(mf->f);
01198 
01199         // We need to store some data to the trX header data
01200         // structure inside the md_file structure
01201         hdr = mf->trx;
01202         if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);
01203 
01204         // Read the magic number
01205         if (trx_int(mf, &magic) < 0) return -1;
01206         if (magic != TRX_MAGIC) {
01207                 // Try reverse endianism
01208                 swap4_aligned(&magic, 1);
01209                 if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01210 
01211                 // Enable byte swapping (actually works, too!)
01212                 mf->rev = 1;
01213         }
01214 
01215         // Read the version number. 
01216         // XXX. this is not the version number, but the storage size
01217         // of the following XDR encoded string.
01218         // the 'title' string is in fact the version identifier.
01219         // since VMD does not use any of that, it does no harm,
01220         // but is should still be fixed occasionally. AK 2005/01/08.
01221 
01222         if(mf->fmt!=MDFMT_TRJ) {
01223                 // It appears that TRJ files either don't contain a version
01224                 // number or don't have a length-delimiter on the string,
01225                 // whereas TRR files do contain both.  Thus, with TRJ, we just
01226                 // assume that the version number is the string length and 
01227                 // just hope for the best. -- WLD 2006/07/09
01228                 if (trx_int(mf, &hdr->version) < 0) return -1;
01229         }
01230 
01231         // Read in the title string
01232         if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
01233                 return -1;
01234 
01235         // Read in some size data
01236         if (trx_int(mf, &hdr->ir_size) < 0) return -1;
01237         if (trx_int(mf, &hdr->e_size) < 0) return -1;
01238         if (trx_int(mf, &hdr->box_size) < 0) return -1;
01239         if (trx_int(mf, &hdr->vir_size) < 0) return -1;
01240         if (trx_int(mf, &hdr->pres_size) < 0) return -1;
01241         if (trx_int(mf, &hdr->top_size) < 0) return -1;
01242         if (trx_int(mf, &hdr->sym_size) < 0) return -1;
01243         if (trx_int(mf, &hdr->x_size) < 0) return -1;
01244         if (trx_int(mf, &hdr->v_size) < 0) return -1;
01245         if (trx_int(mf, &hdr->f_size) < 0) return -1;
01246         if (trx_int(mf, &hdr->natoms) < 0) return -1;
01247         if (trx_int(mf, &hdr->step) < 0) return -1;
01248         if (trx_int(mf, &hdr->nre) < 0) return -1;
01249 
01250         // Make sure there are atoms...
01251         if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);
01252 
01253         // Try to determine precision (float? double?)
01254         if (hdr->x_size) mf->prec = hdr->x_size / (hdr->natoms * 3);
01255         else if (hdr->v_size) mf->prec = hdr->v_size / (hdr->natoms * 3);
01256         else if (hdr->f_size) mf->prec = hdr->f_size / (hdr->natoms * 3);
01257         else return mdio_seterror(MDIO_BADPRECISION);
01258 
01259         if (mf->prec != sizeof(float) && mf->prec != sizeof(double)) {
01260                 // We have no data types this size! The
01261                 // file must've been generated on another
01262                 // platform
01263                 return mdio_seterror(MDIO_BADPRECISION);
01264         }
01265 
01266         // Read in timestep and lambda
01267         if (trx_real(mf, &hdr->t) < 0) return -1;
01268         if (trx_real(mf, &hdr->lambda) < 0) return -1;
01269 
01270         // Rewind if necessary
01271         if (rewind) fseek(mf->f, fpos, SEEK_SET);
01272 
01273         return 0;
01274 }
01275 
01276 
01277 // Reads in an integer and stores it in y. Returns GMX_SUCCESS
01278 // on success or a negative number on error.
01279 static int trx_int(md_file *mf, int *y) {
01280         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01281 
01282         // sanity check.
01283         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01284 
01285         if (y) {
01286                 if (fread(y, 4, 1, mf->f) != 1)
01287                         return mdio_seterror(MDIO_IOERROR);
01288                 if (mf->rev) swap4_aligned(y, 1);
01289         }
01290         else if (fseek(mf->f, 4, SEEK_CUR) != 0)
01291                 return mdio_seterror(MDIO_IOERROR);
01292 
01293         return mdio_seterror(MDIO_SUCCESS);
01294 }
01295 
01296 
01297 // Reads in either a float or a double, depending on the
01298 // precision, and stores that number in y. Returns
01299 // GMX_SUCCESS on success or a negative number on error.
01300 static int trx_real(md_file *mf, float *y) {
01301         double x;
01302 
01303         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01304 
01305         switch (mf->prec) {
01306                 case sizeof(float):
01307                         if (!y) {
01308                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01309                                         return mdio_seterror(MDIO_IOERROR);
01310                         } else {
01311                                 if (fread(y, mf->prec, 1, mf->f) != 1)
01312                                         return mdio_seterror(MDIO_IOERROR);
01313                                 if (mf->rev) swap4_aligned(y, 1);
01314                         }
01315                         return mdio_seterror(MDIO_SUCCESS);
01316 
01317                 case sizeof(double):
01318                         if (!y) {
01319                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01320                                         return mdio_seterror(MDIO_IOERROR);
01321                         } else {
01322                                 if (fread(&x, mf->prec, 1, mf->f) != 1)
01323                                         return mdio_seterror(MDIO_IOERROR);
01324                                 if (mf->rev) swap8_aligned(&x, 1);
01325                                 *y = (float) x;
01326                         }
01327                         return mdio_seterror(MDIO_SUCCESS);
01328 
01329                 default:
01330                         return mdio_seterror(MDIO_BADPRECISION);
01331         }
01332 
01333 }
01334 
01335 
01336 // Reads in a real-valued vector (taking precision into account).
01337 // Stores the vector in vec, and returns GMX_SUCCESS on success
01338 // or a negative number on error.
01339 static int trx_rvector(md_file *mf, float *vec) {
01340         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01341 
01342         if (!vec) {
01343                 if (trx_real(mf, NULL) < 0) return -1;
01344                 if (trx_real(mf, NULL) < 0) return -1;
01345                 if (trx_real(mf, NULL) < 0) return -1;
01346                 return mdio_seterror(MDIO_SUCCESS);
01347         } else {
01348                 if (trx_real(mf, &vec[0]) < 0) return -1;
01349                 if (trx_real(mf, &vec[1]) < 0) return -1;
01350                 if (trx_real(mf, &vec[2]) < 0) return -1;
01351                 return mdio_seterror(MDIO_SUCCESS);
01352         }
01353 }
01354 
01355 
01356 // Reads in a string by first reading an integer containing the
01357 // string's length, then reading in the string itself and storing
01358 // it in str. If the length is greater than max, it is truncated
01359 // and the rest of the string is skipped in the file. Returns the
01360 // length of the string on success or a negative number on error.
01361 static int trx_string(md_file *mf, char *str, int max) {
01362         int size;
01363   size_t ssize;
01364 
01365         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01366 
01367         if (trx_int(mf, &size) < 0) return -1;
01368   ssize = (size_t)size;
01369 
01370         if (str && size <= max) {
01371                 if (fread(str, 1, size, mf->f) != ssize)
01372                         return mdio_seterror(MDIO_IOERROR);
01373                 str[size] = 0;
01374                 return size;
01375         } else if (str) {
01376                 if (fread(str, 1, max, mf->f) != ssize)
01377                         return mdio_seterror(MDIO_IOERROR);
01378                 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01379                         return mdio_seterror(MDIO_IOERROR);
01380                 str[max] = 0;
01381                 return max;
01382         } else {
01383                 if (fseek(mf->f, size, SEEK_CUR) != 0)
01384                         return mdio_seterror(MDIO_IOERROR);
01385                 return 0;
01386         }
01387 }
01388 
01389 
01390 // Reads in a timestep frame from the .trX file and returns the
01391 // data in a timestep structure. Returns NULL on error.
01392 static int trx_timestep(md_file *mf, md_ts *ts) {
01393   int i;
01394   float x[3], y[3], z[3];
01395   trx_hdr *hdr;
01396 
01397   if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01398   if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
01399     return mdio_seterror(MDIO_WRONGFORMAT);
01400 
01401   // Read the header
01402   if (trx_header(mf) < 0) return -1;
01403 
01404   // We need some data from the trX header
01405   hdr = mf->trx;
01406   if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
01407 
01408   if (hdr->box_size) { // XXX need to check value of box_size!!
01409     if (trx_rvector(mf, x) < 0) return -1;
01410     if (trx_rvector(mf, y) < 0) return -1;
01411     if (trx_rvector(mf, z) < 0) return -1;
01412 
01413     // Allocate the box and convert the vectors.
01414     ts->box = (md_box *) malloc(sizeof(md_box));
01415     if (mdio_readbox(ts->box, x, y, z) < 0) {
01416       free(ts->box);
01417       ts->box = NULL;
01418       return -1;
01419     }
01420   }
01421 
01422   if (hdr->vir_size) {
01423     if (trx_rvector(mf, NULL) < 0) return -1;
01424     if (trx_rvector(mf, NULL) < 0) return -1;
01425     if (trx_rvector(mf, NULL) < 0) return -1;
01426   }
01427 
01428   if (hdr->pres_size) {
01429     if (trx_rvector(mf, NULL) < 0) return -1;
01430     if (trx_rvector(mf, NULL) < 0) return -1;
01431     if (trx_rvector(mf, NULL) < 0) return -1;
01432   }
01433 
01434   if (hdr->x_size) {
01435     ts->pos = (float *) malloc(sizeof(float) * 3 * hdr->natoms);
01436     if (!ts->pos) 
01437       return mdio_seterror(MDIO_BADMALLOC);
01438 
01439     ts->natoms = hdr->natoms;
01440     for (i = 0; i < hdr->natoms; i++) {
01441       if (trx_rvector(mf, &ts->pos[i * 3]) < 0) {
01442         mdio_tsfree(ts, 1);
01443         return -1;
01444       }
01445 
01446       ts->pos[i * 3    ] *= ANGS_PER_NM;
01447       ts->pos[i * 3 + 1] *= ANGS_PER_NM;
01448       ts->pos[i * 3 + 2] *= ANGS_PER_NM;
01449     }
01450   }
01451 
01452   if (hdr->v_size) {
01453     for (i = 0; i < hdr->natoms; i++) {
01454       if (trx_rvector(mf, NULL) < 0) {
01455         mdio_tsfree(ts, 1);
01456         return -1;
01457       }
01458     }
01459   }
01460 
01461   if (hdr->f_size) {
01462     for (i = 0; i < hdr->natoms; i++) {
01463       if (trx_rvector(mf, NULL) < 0) {
01464         mdio_tsfree(ts, 1);
01465         return -1;
01466       }
01467     }
01468   }
01469 
01470   return mdio_seterror(MDIO_SUCCESS);
01471 }
01472 
01473 
01474 // writes an int in big endian. Returns GMX_SUCCESS
01475 // on success or a negative number on error.
01476 static int put_trx_int(md_file *mf, int y) {
01477       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01478 
01479       // sanity check.
01480       if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01481 
01482       if (mf->rev) swap4_aligned(&y, 1);
01483       if (fwrite(&y, 4, 1, mf->f) != 1)
01484     return mdio_seterror(MDIO_IOERROR);
01485 
01486   return mdio_seterror(MDIO_SUCCESS);
01487 }
01488 
01489 // writes a real in big-endian. Returns GMX_SUCCESS
01490 // on success or a negative number on error.
01491 static int put_trx_real(md_file *mf, float y) {
01492       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01493 
01494       if (mf->rev) swap4_aligned(&y, 1);
01495       if (fwrite(&y, 4, 1, mf->f) != 1)
01496         return mdio_seterror(MDIO_IOERROR);
01497 
01498       return mdio_seterror(MDIO_SUCCESS);
01499 }
01500 
01501 
01502 // writes an xdr encoded string. Returns GMX_SUCCESS
01503 // on success or a negative number on error.
01504 static int put_trx_string(md_file *mf, const char *s) {
01505         if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
01506         
01507         // write: size of object, string length, string data
01508         size_t len = strlen(s);
01509         if ( put_trx_int(mf, len+1)
01510              || put_trx_int(mf, len)
01511              || (fwrite(s, len, 1, mf->f) != 1))
01512           return mdio_seterror(MDIO_IOERROR);
01513 
01514         return mdio_seterror(MDIO_SUCCESS);
01515 }
01516 
01517 
01518 // xtc_int() - reads an integer from an xtc file
01519 static int xtc_int(md_file *mf, int *i) {
01520         unsigned char c[4];
01521 
01522         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01523         // sanity check.
01524         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01525 
01526         if (fread(c, 1, 4, mf->f) != 4) {
01527                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01528                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01529                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01530         }
01531 
01532         if (i) *i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01533         return mdio_seterror(MDIO_SUCCESS);
01534 }
01535 
01536 
01537 // xtc_float() - reads a float from an xtc file
01538 static int xtc_float(md_file *mf, float *f) {
01539         unsigned char c[4];
01540         int i;
01541 
01542         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01543 
01544         if (fread(c, 1, 4, mf->f) != 4) {
01545                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01546                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01547                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01548         }
01549 
01550         if (f) {
01551                 // By reading the number in as an integer and then
01552                 // copying it to a floating point number we can
01553                 // ensure proper endianness
01554                 i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01555                 memcpy(f, &i, 4);
01556         }
01557         return mdio_seterror(MDIO_SUCCESS);
01558 }
01559 
01560 
01561 // xtc_data() - reads a specific amount of data from an xtc
01562 // file using the xdr format.
01563 static int xtc_data(md_file *mf, char *buf, int len) {
01564         if (!mf || len < 1) return mdio_seterror(MDIO_BADPARAMS);
01565   size_t slen = (size_t)len;
01566         if (buf) {
01567                 if (fread(buf, 1, slen, mf->f) != slen) {
01568                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01569                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01570                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01571                 }
01572                 if (len % 4) {
01573                         if (fseek(mf->f, 4 - (len % 4), SEEK_CUR)) {
01574                                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01575                                 if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01576                                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01577                         }
01578                 }
01579         }
01580         else {
01581                 int newlen;
01582                 newlen = len;
01583                 if (len % 4) newlen += (4 - (len % 4));
01584                 if (fseek(mf->f, newlen, SEEK_CUR)) {
01585                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01586                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01587                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01588                 }
01589         }
01590         return len;
01591 }
01592 
01593 
01594 // xtc_timestep() - reads a timestep from an .xtc file.
01595 static int xtc_timestep(md_file *mf, md_ts *ts) {
01596         int n;
01597         float f, x[3], y[3], z[3];
01598 
01599         int size = 0; // explicitly initialized to zero.
01600         float precision;
01601 
01602         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01603         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
01604         if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);
01605 
01606         // Check magic number
01607         if (xtc_int(mf, &n) < 0) return -1;
01608         if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01609 
01610         // Get number of atoms
01611         if (xtc_int(mf, &n) < 0) return -1;
01612         ts->natoms = n;
01613 
01614         // Get the simulation step
01615         if (xtc_int(mf, &n) < 0) return -1;
01616         ts->step = n;
01617 
01618         // Get the time value
01619         if (xtc_float(mf, &f) < 0) return -1;
01620         ts->time = f;
01621 
01622         // Read the basis vectors of the box
01623   if ( (xtc_float(mf, &x[0]) < 0) ||
01624        (xtc_float(mf, &x[1]) < 0) ||
01625        (xtc_float(mf, &x[2]) < 0) ||
01626        (xtc_float(mf, &y[0]) < 0) ||
01627        (xtc_float(mf, &y[1]) < 0) ||
01628        (xtc_float(mf, &y[2]) < 0) ||
01629        (xtc_float(mf, &z[0]) < 0) ||
01630        (xtc_float(mf, &z[1]) < 0) ||
01631        (xtc_float(mf, &z[2]) < 0) )
01632     return -1;
01633   // Allocate the box and convert the vectors.
01634   ts->box = (md_box *) malloc(sizeof(md_box));
01635   if (mdio_readbox(ts->box, x, y, z) < 0) {
01636     free(ts->box);
01637     ts->box = NULL;
01638     return -1;
01639   }
01640 
01641         ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
01642         if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
01643         n = xtc_3dfcoord(mf, ts->pos, &size, &precision);
01644         if (n < 0) return -1;
01645 
01646         /* Now we're left with the job of scaling... */
01647         for (n = 0; n < ts->natoms * 3; n++)
01648                 ts->pos[n] *= ANGS_PER_NM;
01649 
01650         return mdio_seterror(MDIO_SUCCESS);
01651 }
01652 
01653 
01655 // This algorithm is an implementation of the 3dfcoord algorithm
01656 // written by Frans van Hoesel (hoesel@chem.rug.nl) as part of the
01657 // Europort project in 1995.
01659 
01660 // integer table used in decompression
01661 static int xtc_magicints[] = {
01662         0, 0, 0, 0, 0, 0, 0, 0, 0,8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
01663         80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
01664         1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, 16384,
01665         20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031, 131072,
01666         165140, 208063, 262144, 330280, 416127, 524287, 660561, 832255,
01667         1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
01668         5284491, 6658042, 8388607, 10568983, 13316085, 16777216 };
01669 
01670 #define FIRSTIDX 9
01671 /* note that magicints[FIRSTIDX-1] == 0 */
01672 #define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))
01673 
01674 
01675 // returns the number of bits in the binary expansion of
01676 // the given integer.
01677 static int xtc_sizeofint(int size) {
01678         unsigned int num = 1;
01679   unsigned int ssize = (unsigned int)size;
01680         int nbits = 0;
01681 
01682         while (ssize >= num && nbits < 32) {
01683                 nbits++;
01684                 num <<= 1;
01685         }
01686         return nbits;
01687 }
01688 
01689 // calculates the number of bits a set of integers, when compressed,
01690 // will take up.
01691 static int xtc_sizeofints(int nints, unsigned int *sizes) {
01692         int i;
01693   unsigned int num;
01694         unsigned int nbytes, nbits, bytes[32], bytecnt, tmp;
01695         nbytes = 1;
01696         bytes[0] = 1;
01697         nbits = 0;
01698         for (i=0; i < nints; i++) {     
01699                 tmp = 0;
01700                 for (bytecnt = 0; bytecnt < nbytes; bytecnt++) {
01701                         tmp = bytes[bytecnt] * sizes[i] + tmp;
01702                         bytes[bytecnt] = tmp & 0xff;
01703                         tmp >>= 8;
01704                 }
01705                 while (tmp != 0) {
01706                         bytes[bytecnt++] = tmp & 0xff;
01707                         tmp >>= 8;
01708                 }
01709                 nbytes = bytecnt;
01710         }
01711         num = 1;
01712         nbytes--;
01713         while (bytes[nbytes] >= num) {
01714                 nbits++;
01715                 num *= 2;
01716         }
01717         return nbits + nbytes * 8;
01718 }
01719 
01720 // reads bits from a buffer.    
01721 static int xtc_receivebits(int *buf, int nbits) {
01722         int cnt, num; 
01723         unsigned int lastbits, lastbyte;
01724         unsigned char * cbuf;
01725         int mask = (1 << nbits) -1;
01726 
01727         cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
01728         cnt = buf[0];
01729         lastbits = (unsigned int) buf[1];
01730         lastbyte = (unsigned int) buf[2];
01731 
01732         num = 0;
01733         while (nbits >= 8) {
01734                 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
01735                 num |=  (lastbyte >> lastbits) << (nbits - 8);
01736                 nbits -=8;
01737         }
01738         if (nbits > 0) {
01739                 if (lastbits < (unsigned int)nbits) {
01740                         lastbits += 8;
01741                         lastbyte = (lastbyte << 8) | cbuf[cnt++];
01742                 }
01743                 lastbits -= nbits;
01744                 num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
01745         }
01746         num &= mask;
01747         buf[0] = cnt;
01748         buf[1] = lastbits;
01749         buf[2] = lastbyte;
01750         return num; 
01751 }
01752 
01753 // decompresses small integers from the buffer
01754 // sizes parameter has to be non-zero to prevent divide-by-zero
01755 static void xtc_receiveints(int *buf, const int nints, int nbits,
01756                         unsigned int *sizes, int *nums) {
01757         int bytes[32];
01758         int i, j, nbytes, p, num;
01759 
01760         bytes[1] = bytes[2] = bytes[3] = 0;
01761         nbytes = 0;
01762         while (nbits > 8) {
01763                 bytes[nbytes++] = xtc_receivebits(buf, 8);
01764                 nbits -= 8;
01765         }
01766         if (nbits > 0) {
01767                 bytes[nbytes++] = xtc_receivebits(buf, nbits);
01768         }
01769         for (i = nints-1; i > 0; i--) {
01770                 num = 0;
01771                 for (j = nbytes-1; j >=0; j--) {
01772                         num = (num << 8) | bytes[j];
01773                         p = num / sizes[i];
01774                         bytes[j] = p;
01775                         num = num - p * sizes[i];
01776                 }
01777                 nums[i] = num;
01778         }
01779         nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
01780 }
01781 
01782 // function that actually reads compressed coordinates    
01783 static int xtc_3dfcoord(md_file *mf, float *fp, int *size, float *precision) {
01784         static int *ip = NULL;
01785         static int oldsize;
01786         static int *buf;
01787 
01788         int minint[3], maxint[3], *lip;
01789         int smallidx;
01790         unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
01791         int flag, k;
01792         int small, smaller, i, is_smaller, run;
01793         float *lfp;
01794         int tmp, *thiscoord,  prevcoord[3];
01795 
01796         int bufsize, lsize;
01797         unsigned int bitsize;
01798         float inv_precision;
01799 
01800         /* avoid uninitialized data compiler warnings */
01801         bitsizeint[0] = 0;
01802         bitsizeint[1] = 0;
01803         bitsizeint[2] = 0;
01804 
01805         if (xtc_int(mf, &lsize) < 0) return -1;
01806 
01807         if (*size != 0 && lsize != *size) return mdio_seterror(MDIO_BADFORMAT);
01808 
01809         *size = lsize;
01810         size3 = *size * 3;
01811         if (*size <= 9) {
01812                 for (i = 0; i < *size; i++) {
01813                         if (xtc_float(mf, fp + (3 * i)) < 0) return -1;
01814                         if (xtc_float(mf, fp + (3 * i) + 1) < 0) return -1;
01815                         if (xtc_float(mf, fp + (3 * i) + 2) < 0) return -1;
01816                 }
01817                 return *size;
01818         }
01819         xtc_float(mf, precision);
01820         if (ip == NULL) {
01821                 ip = (int *)malloc(size3 * sizeof(*ip));
01822                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
01823                 bufsize = (int) (size3 * 1.2);
01824                 buf = (int *)malloc(bufsize * sizeof(*buf));
01825                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
01826                 oldsize = *size;
01827         } else if (*size > oldsize) {
01828                 ip = (int *)realloc(ip, size3 * sizeof(*ip));
01829                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
01830                 bufsize = (int) (size3 * 1.2);
01831                 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
01832                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
01833                 oldsize = *size;
01834         }
01835         buf[0] = buf[1] = buf[2] = 0;
01836 
01837         xtc_int(mf, &(minint[0]));
01838         xtc_int(mf, &(minint[1]));
01839         xtc_int(mf, &(minint[2]));
01840 
01841         xtc_int(mf, &(maxint[0]));
01842         xtc_int(mf, &(maxint[1]));
01843         xtc_int(mf, &(maxint[2]));
01844                 
01845         sizeint[0] = maxint[0] - minint[0]+1;
01846         sizeint[1] = maxint[1] - minint[1]+1;
01847         sizeint[2] = maxint[2] - minint[2]+1;
01848         
01849         /* check if one of the sizes is to big to be multiplied */
01850         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
01851                 bitsizeint[0] = xtc_sizeofint(sizeint[0]);
01852                 bitsizeint[1] = xtc_sizeofint(sizeint[1]);
01853                 bitsizeint[2] = xtc_sizeofint(sizeint[2]);
01854                 bitsize = 0; /* flag the use of large sizes */
01855         } else {
01856                 bitsize = xtc_sizeofints(3, sizeint);
01857         }
01858 
01859         xtc_int(mf, &smallidx);
01860         smaller = xtc_magicints[FIRSTIDX > smallidx - 1 ? FIRSTIDX : smallidx - 1] / 2;
01861         small = xtc_magicints[smallidx] / 2;
01862         sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx];
01863 
01864         /* check for zero values that would yield corrupted data */
01865         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
01866                 printf("XTC corrupted, sizesmall==0 (case 1)\n");
01867                 return -1;
01868         }
01869 
01870 
01871         /* buf[0] holds the length in bytes */
01872         if (xtc_int(mf, &(buf[0])) < 0) return -1;
01873 
01874         if (xtc_data(mf, (char *) &buf[3], (int) buf[0]) < 0) return -1;
01875 
01876         buf[0] = buf[1] = buf[2] = 0;
01877 
01878         lfp = fp;
01879         inv_precision = 1.0f / (*precision);
01880         run = 0;
01881         i = 0;
01882         lip = ip;
01883         while (i < lsize) {
01884                 thiscoord = (int *)(lip) + i * 3;
01885 
01886                 if (bitsize == 0) {
01887                         thiscoord[0] = xtc_receivebits(buf, bitsizeint[0]);
01888                         thiscoord[1] = xtc_receivebits(buf, bitsizeint[1]);
01889                         thiscoord[2] = xtc_receivebits(buf, bitsizeint[2]);
01890                 } else {
01891                         xtc_receiveints(buf, 3, bitsize, sizeint, thiscoord);
01892                 }
01893 
01894                 i++;
01895                 thiscoord[0] += minint[0];
01896                 thiscoord[1] += minint[1];
01897                 thiscoord[2] += minint[2];
01898 
01899                 prevcoord[0] = thiscoord[0];
01900                 prevcoord[1] = thiscoord[1];
01901                 prevcoord[2] = thiscoord[2];
01902  
01903 
01904                 flag = xtc_receivebits(buf, 1);
01905                 is_smaller = 0;
01906                 if (flag == 1) {
01907                         run = xtc_receivebits(buf, 5);
01908                         is_smaller = run % 3;
01909                         run -= is_smaller;
01910                         is_smaller--;
01911                 }
01912                 if (run > 0) {
01913                         thiscoord += 3;
01914                         for (k = 0; k < run; k+=3) {
01915                                 xtc_receiveints(buf, 3, smallidx, sizesmall, thiscoord);
01916                                 i++;
01917                                 thiscoord[0] += prevcoord[0] - small;
01918                                 thiscoord[1] += prevcoord[1] - small;
01919                                 thiscoord[2] += prevcoord[2] - small;
01920                                 if (k == 0) {
01921                                         /* interchange first with second atom for better
01922                                          * compression of water molecules
01923                                          */
01924                                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
01925                                         prevcoord[0] = tmp;
01926                                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
01927                                         prevcoord[1] = tmp;
01928                                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
01929                                         prevcoord[2] = tmp;
01930                                         *lfp++ = prevcoord[0] * inv_precision;
01931                                         *lfp++ = prevcoord[1] * inv_precision;
01932                                         *lfp++ = prevcoord[2] * inv_precision;
01933 
01934                                         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
01935                                                 printf("XTC corrupted, sizesmall==0 (case 2)\n");
01936                                                 return -1;
01937                                         }
01938 
01939                                 } else {
01940                                         prevcoord[0] = thiscoord[0];
01941                                         prevcoord[1] = thiscoord[1];
01942                                         prevcoord[2] = thiscoord[2];
01943                                 }
01944                                 *lfp++ = thiscoord[0] * inv_precision;
01945                                 *lfp++ = thiscoord[1] * inv_precision;
01946                                 *lfp++ = thiscoord[2] * inv_precision;
01947                         }
01948                 } else {
01949                         *lfp++ = thiscoord[0] * inv_precision;
01950                         *lfp++ = thiscoord[1] * inv_precision;
01951                         *lfp++ = thiscoord[2] * inv_precision;          
01952                 }
01953                 smallidx += is_smaller;
01954                 if (is_smaller < 0) {
01955                         small = smaller;
01956                         if (smallidx > FIRSTIDX) {
01957                                 smaller = xtc_magicints[smallidx - 1] /2;
01958                         } else {
01959                                 smaller = 0;
01960                         }
01961                 } else if (is_smaller > 0) {
01962                         smaller = small;
01963                         small = xtc_magicints[smallidx] / 2;
01964                 }
01965                 sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx] ;
01966         }
01967         return 1;
01968 }
01969 
01970 
01971 
01972 #endif

Generated on Tue Jan 21 02:56:13 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002