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.37 $       $Date: 2023/04/22 06:10:42 $
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 
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 static int trx_long(md_file *mf, long *y) {
01297         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01298 
01299         // sanity check.
01300         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01301 
01302         if (y) {
01303                 if (fread(y, 8, 1, mf->f) != 1)
01304                         return mdio_seterror(MDIO_IOERROR);
01305                 if (mf->rev) swap8_aligned(y, 1);
01306         }
01307         else if (fseek(mf->f, 4, SEEK_CUR) != 0)
01308                 return mdio_seterror(MDIO_IOERROR);
01309 
01310         return mdio_seterror(MDIO_SUCCESS);
01311 }
01312 
01313 
01314 // Reads in either a float or a double, depending on the
01315 // precision, and stores that number in y. Returns
01316 // GMX_SUCCESS on success or a negative number on error.
01317 static int trx_real(md_file *mf, float *y) {
01318         double x;
01319 
01320         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01321 
01322         switch (mf->prec) {
01323                 case sizeof(float):
01324                         if (!y) {
01325                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01326                                         return mdio_seterror(MDIO_IOERROR);
01327                         } else {
01328                                 if (fread(y, mf->prec, 1, mf->f) != 1)
01329                                         return mdio_seterror(MDIO_IOERROR);
01330                                 if (mf->rev) swap4_aligned(y, 1);
01331                         }
01332                         return mdio_seterror(MDIO_SUCCESS);
01333 
01334                 case sizeof(double):
01335                         if (!y) {
01336                                 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01337                                         return mdio_seterror(MDIO_IOERROR);
01338                         } else {
01339                                 if (fread(&x, mf->prec, 1, mf->f) != 1)
01340                                         return mdio_seterror(MDIO_IOERROR);
01341                                 if (mf->rev) swap8_aligned(&x, 1);
01342                                 *y = (float) x;
01343                         }
01344                         return mdio_seterror(MDIO_SUCCESS);
01345 
01346                 default:
01347                         return mdio_seterror(MDIO_BADPRECISION);
01348         }
01349 
01350 }
01351 
01352 static int trx_double(md_file *mf, double *y) {
01353         double x;
01354 
01355         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01356 
01357         if (!y) {
01358                 if (fseek(mf->f, 8, SEEK_CUR) != 0)
01359                         return mdio_seterror(MDIO_IOERROR);
01360         } else {
01361                 if (fread(&x, 8, 1, mf->f) != 1)
01362                         return mdio_seterror(MDIO_IOERROR);
01363                 if (mf->rev) swap8_aligned(&x, 1);
01364                 *y = (float) x;
01365         }
01366         return mdio_seterror(MDIO_SUCCESS);
01367 
01368 }
01369 
01370 
01371 // Reads in a real-valued vector (taking precision into account).
01372 // Stores the vector in vec, and returns GMX_SUCCESS on success
01373 // or a negative number on error.
01374 static int trx_rvector(md_file *mf, float *vec) {
01375         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01376 
01377         if (!vec) {
01378                 if (trx_real(mf, NULL) < 0) return -1;
01379                 if (trx_real(mf, NULL) < 0) return -1;
01380                 if (trx_real(mf, NULL) < 0) return -1;
01381                 return mdio_seterror(MDIO_SUCCESS);
01382         } else {
01383                 if (trx_real(mf, &vec[0]) < 0) return -1;
01384                 if (trx_real(mf, &vec[1]) < 0) return -1;
01385                 if (trx_real(mf, &vec[2]) < 0) return -1;
01386                 return mdio_seterror(MDIO_SUCCESS);
01387         }
01388 }
01389 
01390 static int tpr_rvector (md_file *mf, float *arr, int len) {
01391     int i;
01392     for (i = 0; i < len; i++) {
01393         if (trx_real(mf, &(arr[i]))) return -1;
01394     }
01395     return mdio_seterror(MDIO_SUCCESS);
01396 }
01397 
01398 static int tpr_ivector (md_file *mf, int *arr, int len) {
01399     int i;
01400     for (i = 0; i < len; i++) {
01401         if (trx_int(mf, &(arr[i]))) return -1;
01402     }
01403     return mdio_seterror(MDIO_SUCCESS);
01404 }
01405 
01406 // Reads in a string by first reading an integer containing the
01407 // string's length, then reading in the string itself and storing
01408 // it in str. If the length is greater than max, it is truncated
01409 // and the rest of the string is skipped in the file. Returns the
01410 // length of the string on success or a negative number on error.
01411 static int trx_string(md_file *mf, char *str, int max) {
01412         int size;
01413   size_t ssize;
01414 
01415         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01416 
01417         if (trx_int(mf, &size) < 0) return -1;
01418   ssize = (size_t)size;
01419 
01420         if (str && size <= max) {
01421                 if (fread(str, 1, size, mf->f) != ssize)
01422                         return mdio_seterror(MDIO_IOERROR);
01423                 str[size] = 0;
01424                 return size;
01425         } else if (str) {
01426                 if (fread(str, 1, max, mf->f) != ssize)
01427                         return mdio_seterror(MDIO_IOERROR);
01428                 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01429                         return mdio_seterror(MDIO_IOERROR);
01430                 str[max] = 0;
01431                 return max;
01432         } else {
01433                 if (fseek(mf->f, size, SEEK_CUR) != 0)
01434                         return mdio_seterror(MDIO_IOERROR);
01435                 return 0;
01436         }
01437 }
01438 
01439 
01440 // Reads in a string by first reading an integer containing the
01441 // string's length, then reading in the string itself and storing
01442 // it in str. If the length is greater than max, it is truncated
01443 // and the rest of the string is skipped in the file. Returns the
01444 // length of the string on success or a negative number on error.
01445 // Unlike trx_string, this will read in a 4-byte aligned way.
01446 static int tpr_string(md_file *mf, char *str, int max) {
01447         int size;
01448   size_t ssize;
01449 
01450         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01451 
01452         if (trx_int(mf, &size) < 0) return -1;
01453         if (size % 4) {
01454                 size += 4 - (size % 4);
01455         }
01456   ssize = (size_t)size;
01457 
01458         if (str && size <= max) {
01459                 if (fread(str, 1, size, mf->f) != ssize)
01460                         return mdio_seterror(MDIO_IOERROR);
01461                 str[size] = 0;
01462                 return size;
01463         } else if (str) {
01464                 if (fread(str, 1, max, mf->f) != ssize)
01465                         return mdio_seterror(MDIO_IOERROR);
01466                 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01467                         return mdio_seterror(MDIO_IOERROR);
01468                 str[max] = 0;
01469                 return max;
01470         } else {
01471                 if (fseek(mf->f, size, SEEK_CUR) != 0)
01472                         return mdio_seterror(MDIO_IOERROR);
01473                 return 0;
01474         }
01475 }
01476 
01477 // Reads in a timestep frame from the .trX file and returns the
01478 // data in a timestep structure. Returns NULL on error.
01479 static int trx_timestep(md_file *mf, md_ts *ts) {
01480   int i;
01481   float x[3], y[3], z[3];
01482   trx_hdr *hdr;
01483 
01484   if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01485   if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
01486     return mdio_seterror(MDIO_WRONGFORMAT);
01487 
01488   // Read the header
01489   if (trx_header(mf) < 0) return -1;
01490 
01491   // We need some data from the trX header
01492   hdr = mf->trx;
01493   if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
01494 
01495   if (hdr->box_size) { // XXX need to check value of box_size!!
01496     if (trx_rvector(mf, x) < 0) return -1;
01497     if (trx_rvector(mf, y) < 0) return -1;
01498     if (trx_rvector(mf, z) < 0) return -1;
01499 
01500     // Allocate the box and convert the vectors.
01501     ts->box = (md_box *) malloc(sizeof(md_box));
01502     if (mdio_readbox(ts->box, x, y, z) < 0) {
01503       free(ts->box);
01504       ts->box = NULL;
01505       return -1;
01506     }
01507   }
01508 
01509   if (hdr->vir_size) {
01510     if (trx_rvector(mf, NULL) < 0) return -1;
01511     if (trx_rvector(mf, NULL) < 0) return -1;
01512     if (trx_rvector(mf, NULL) < 0) return -1;
01513   }
01514 
01515   if (hdr->pres_size) {
01516     if (trx_rvector(mf, NULL) < 0) return -1;
01517     if (trx_rvector(mf, NULL) < 0) return -1;
01518     if (trx_rvector(mf, NULL) < 0) return -1;
01519   }
01520 
01521   if (hdr->x_size) {
01522     ts->pos = (float *) malloc(sizeof(float) * 3 * hdr->natoms);
01523     if (!ts->pos) 
01524       return mdio_seterror(MDIO_BADMALLOC);
01525 
01526     ts->natoms = hdr->natoms;
01527     for (i = 0; i < hdr->natoms; i++) {
01528       if (trx_rvector(mf, &ts->pos[i * 3]) < 0) {
01529         mdio_tsfree(ts, 1);
01530         return -1;
01531       }
01532 
01533       ts->pos[i * 3    ] *= ANGS_PER_NM;
01534       ts->pos[i * 3 + 1] *= ANGS_PER_NM;
01535       ts->pos[i * 3 + 2] *= ANGS_PER_NM;
01536     }
01537   }
01538 
01539   if (hdr->v_size) {
01540     for (i = 0; i < hdr->natoms; i++) {
01541       if (trx_rvector(mf, NULL) < 0) {
01542         mdio_tsfree(ts, 1);
01543         return -1;
01544       }
01545     }
01546   }
01547 
01548   if (hdr->f_size) {
01549     for (i = 0; i < hdr->natoms; i++) {
01550       if (trx_rvector(mf, NULL) < 0) {
01551         mdio_tsfree(ts, 1);
01552         return -1;
01553       }
01554     }
01555   }
01556 
01557   return mdio_seterror(MDIO_SUCCESS);
01558 }
01559 
01560 
01561 // writes an int in big endian. Returns GMX_SUCCESS
01562 // on success or a negative number on error.
01563 static int put_trx_int(md_file *mf, int y) {
01564       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01565 
01566       // sanity check.
01567       if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01568 
01569       if (mf->rev) swap4_aligned(&y, 1);
01570       if (fwrite(&y, 4, 1, mf->f) != 1)
01571     return mdio_seterror(MDIO_IOERROR);
01572 
01573   return mdio_seterror(MDIO_SUCCESS);
01574 }
01575 
01576 // writes a real in big-endian. Returns GMX_SUCCESS
01577 // on success or a negative number on error.
01578 static int put_trx_real(md_file *mf, float y) {
01579       if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01580 
01581       if (mf->rev) swap4_aligned(&y, 1);
01582       if (fwrite(&y, 4, 1, mf->f) != 1)
01583         return mdio_seterror(MDIO_IOERROR);
01584 
01585       return mdio_seterror(MDIO_SUCCESS);
01586 }
01587 
01588 
01589 // writes an xdr encoded string. Returns GMX_SUCCESS
01590 // on success or a negative number on error.
01591 static int put_trx_string(md_file *mf, const char *s) {
01592         if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
01593         
01594         // write: size of object, string length, string data
01595         size_t len = strlen(s);
01596         if ( put_trx_int(mf, len+1)
01597              || put_trx_int(mf, len)
01598              || (fwrite(s, len, 1, mf->f) != 1))
01599           return mdio_seterror(MDIO_IOERROR);
01600 
01601         return mdio_seterror(MDIO_SUCCESS);
01602 }
01603 
01604 
01605 // xtc_int() - reads an integer from an xtc file
01606 static int xtc_int(md_file *mf, int *i) {
01607         unsigned char c[4];
01608 
01609         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01610         // sanity check.
01611         if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01612 
01613         if (fread(c, 1, 4, mf->f) != 4) {
01614                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01615                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01616                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01617         }
01618 
01619         if (i) *i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01620         return mdio_seterror(MDIO_SUCCESS);
01621 }
01622 
01623 
01624 // xtc_float() - reads a float from an xtc file
01625 static int xtc_float(md_file *mf, float *f) {
01626         unsigned char c[4];
01627         int i;
01628 
01629         if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01630 
01631         if (fread(c, 1, 4, mf->f) != 4) {
01632                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01633                 else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01634                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01635         }
01636 
01637         if (f) {
01638                 // By reading the number in as an integer and then
01639                 // copying it to a floating point number we can
01640                 // ensure proper endianness
01641                 i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
01642                 memcpy(f, &i, 4);
01643         }
01644         return mdio_seterror(MDIO_SUCCESS);
01645 }
01646 
01647 
01648 // xtc_data() - reads a specific amount of data from an xtc
01649 // file using the xdr format.
01650 static int xtc_data(md_file *mf, char *buf, int len) {
01651         if (!mf || len < 1) return mdio_seterror(MDIO_BADPARAMS);
01652   size_t slen = (size_t)len;
01653         if (buf) {
01654                 if (fread(buf, 1, slen, mf->f) != slen) {
01655                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01656                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01657                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01658                 }
01659                 if (len % 4) {
01660                         if (fseek(mf->f, 4 - (len % 4), SEEK_CUR)) {
01661                                 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01662                                 if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01663                                 else return mdio_seterror(MDIO_UNKNOWNERROR);
01664                         }
01665                 }
01666         }
01667         else {
01668                 int newlen;
01669                 newlen = len;
01670                 if (len % 4) newlen += (4 - (len % 4));
01671                 if (fseek(mf->f, newlen, SEEK_CUR)) {
01672                         if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
01673                         if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
01674                         else return mdio_seterror(MDIO_UNKNOWNERROR);
01675                 }
01676         }
01677         return len;
01678 }
01679 
01680 
01681 // xtc_timestep() - reads a timestep from an .xtc file.
01682 static int xtc_timestep(md_file *mf, md_ts *ts) {
01683         int n;
01684         float f, x[3], y[3], z[3];
01685 
01686         int size = 0; // explicitly initialized to zero.
01687         float precision;
01688 
01689         if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01690         if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
01691         if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);
01692 
01693         // Check magic number
01694         if (xtc_int(mf, &n) < 0) return -1;
01695         if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01696 
01697         // Get number of atoms
01698         if (xtc_int(mf, &n) < 0) return -1;
01699         ts->natoms = n;
01700 
01701         // Get the simulation step
01702         if (xtc_int(mf, &n) < 0) return -1;
01703         ts->step = n;
01704 
01705         // Get the time value
01706         if (xtc_float(mf, &f) < 0) return -1;
01707         ts->time = f;
01708 
01709         // Read the basis vectors of the box
01710   if ( (xtc_float(mf, &x[0]) < 0) ||
01711        (xtc_float(mf, &x[1]) < 0) ||
01712        (xtc_float(mf, &x[2]) < 0) ||
01713        (xtc_float(mf, &y[0]) < 0) ||
01714        (xtc_float(mf, &y[1]) < 0) ||
01715        (xtc_float(mf, &y[2]) < 0) ||
01716        (xtc_float(mf, &z[0]) < 0) ||
01717        (xtc_float(mf, &z[1]) < 0) ||
01718        (xtc_float(mf, &z[2]) < 0) )
01719     return -1;
01720   // Allocate the box and convert the vectors.
01721   ts->box = (md_box *) malloc(sizeof(md_box));
01722   if (mdio_readbox(ts->box, x, y, z) < 0) {
01723     free(ts->box);
01724     ts->box = NULL;
01725     return -1;
01726   }
01727 
01728         ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
01729         if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
01730         n = xtc_3dfcoord(mf, ts->pos, &size, &precision);
01731         if (n < 0) return -1;
01732 
01733         /* Now we're left with the job of scaling... */
01734         for (n = 0; n < ts->natoms * 3; n++)
01735                 ts->pos[n] *= ANGS_PER_NM;
01736 
01737         return mdio_seterror(MDIO_SUCCESS);
01738 }
01739 
01740 
01742 // This algorithm is an implementation of the 3dfcoord algorithm
01743 // written by Frans van Hoesel (hoesel@chem.rug.nl) as part of the
01744 // Europort project in 1995.
01746 
01747 // integer table used in decompression
01748 static int xtc_magicints[] = {
01749         0, 0, 0, 0, 0, 0, 0, 0, 0,8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
01750         80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
01751         1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, 16384,
01752         20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031, 131072,
01753         165140, 208063, 262144, 330280, 416127, 524287, 660561, 832255,
01754         1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
01755         5284491, 6658042, 8388607, 10568983, 13316085, 16777216 };
01756 
01757 #define FIRSTIDX 9
01758 /* note that magicints[FIRSTIDX-1] == 0 */
01759 #define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))
01760 
01761 
01762 // returns the number of bits in the binary expansion of
01763 // the given integer.
01764 static int xtc_sizeofint(int size) {
01765         unsigned int num = 1;
01766   unsigned int ssize = (unsigned int)size;
01767         int nbits = 0;
01768 
01769         while (ssize >= num && nbits < 32) {
01770                 nbits++;
01771                 num <<= 1;
01772         }
01773         return nbits;
01774 }
01775 
01776 // calculates the number of bits a set of integers, when compressed,
01777 // will take up.
01778 static int xtc_sizeofints(int nints, unsigned int *sizes) {
01779         int i;
01780   unsigned int num;
01781         unsigned int nbytes, nbits, bytes[32], bytecnt, tmp;
01782         nbytes = 1;
01783         bytes[0] = 1;
01784         nbits = 0;
01785         for (i=0; i < nints; i++) {     
01786                 tmp = 0;
01787                 for (bytecnt = 0; bytecnt < nbytes; bytecnt++) {
01788                         tmp = bytes[bytecnt] * sizes[i] + tmp;
01789                         bytes[bytecnt] = tmp & 0xff;
01790                         tmp >>= 8;
01791                 }
01792                 while (tmp != 0) {
01793                         bytes[bytecnt++] = tmp & 0xff;
01794                         tmp >>= 8;
01795                 }
01796                 nbytes = bytecnt;
01797         }
01798         num = 1;
01799         nbytes--;
01800         while (bytes[nbytes] >= num) {
01801                 nbits++;
01802                 num *= 2;
01803         }
01804         return nbits + nbytes * 8;
01805 }
01806 
01807 // reads bits from a buffer.    
01808 static int xtc_receivebits(int *buf, int nbits) {
01809         int cnt, num; 
01810         unsigned int lastbits, lastbyte;
01811         unsigned char * cbuf;
01812         int mask = (1 << nbits) -1;
01813 
01814         cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
01815         cnt = buf[0];
01816         lastbits = (unsigned int) buf[1];
01817         lastbyte = (unsigned int) buf[2];
01818 
01819         num = 0;
01820         while (nbits >= 8) {
01821                 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
01822                 num |=  (lastbyte >> lastbits) << (nbits - 8);
01823                 nbits -=8;
01824         }
01825         if (nbits > 0) {
01826                 if (lastbits < (unsigned int)nbits) {
01827                         lastbits += 8;
01828                         lastbyte = (lastbyte << 8) | cbuf[cnt++];
01829                 }
01830                 lastbits -= nbits;
01831                 num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
01832         }
01833         num &= mask;
01834         buf[0] = cnt;
01835         buf[1] = lastbits;
01836         buf[2] = lastbyte;
01837         return num; 
01838 }
01839 
01840 // decompresses small integers from the buffer
01841 // sizes parameter has to be non-zero to prevent divide-by-zero
01842 static void xtc_receiveints(int *buf, const int nints, int nbits,
01843                         unsigned int *sizes, int *nums) {
01844         int bytes[32];
01845         int i, j, nbytes, p, num;
01846 
01847         bytes[1] = bytes[2] = bytes[3] = 0;
01848         nbytes = 0;
01849         while (nbits > 8) {
01850                 bytes[nbytes++] = xtc_receivebits(buf, 8);
01851                 nbits -= 8;
01852         }
01853         if (nbits > 0) {
01854                 bytes[nbytes++] = xtc_receivebits(buf, nbits);
01855         }
01856         for (i = nints-1; i > 0; i--) {
01857                 num = 0;
01858                 for (j = nbytes-1; j >=0; j--) {
01859                         num = (num << 8) | bytes[j];
01860                         p = num / sizes[i];
01861                         bytes[j] = p;
01862                         num = num - p * sizes[i];
01863                 }
01864                 nums[i] = num;
01865         }
01866         nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
01867 }
01868 
01869 // function that actually reads compressed coordinates    
01870 static int xtc_3dfcoord(md_file *mf, float *fp, int *size, float *precision) {
01871         static int *ip = NULL;
01872         static int oldsize;
01873         static int *buf;
01874 
01875         int minint[3], maxint[3], *lip;
01876         int smallidx;
01877         unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
01878         int flag, k;
01879         int small, smaller, i, is_smaller, run;
01880         float *lfp;
01881         int tmp, *thiscoord,  prevcoord[3];
01882 
01883         int bufsize, lsize;
01884         unsigned int bitsize;
01885         float inv_precision;
01886 
01887         /* avoid uninitialized data compiler warnings */
01888         bitsizeint[0] = 0;
01889         bitsizeint[1] = 0;
01890         bitsizeint[2] = 0;
01891 
01892         if (xtc_int(mf, &lsize) < 0) return -1;
01893 
01894         if (*size != 0 && lsize != *size) return mdio_seterror(MDIO_BADFORMAT);
01895 
01896         *size = lsize;
01897         size3 = *size * 3;
01898         if (*size <= 9) {
01899                 for (i = 0; i < *size; i++) {
01900                         if (xtc_float(mf, fp + (3 * i)) < 0) return -1;
01901                         if (xtc_float(mf, fp + (3 * i) + 1) < 0) return -1;
01902                         if (xtc_float(mf, fp + (3 * i) + 2) < 0) return -1;
01903                 }
01904                 return *size;
01905         }
01906         xtc_float(mf, precision);
01907         if (ip == NULL) {
01908                 ip = (int *)malloc(size3 * sizeof(*ip));
01909                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
01910                 bufsize = (int) (size3 * 1.2);
01911                 buf = (int *)malloc(bufsize * sizeof(*buf));
01912                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
01913                 oldsize = *size;
01914         } else if (*size > oldsize) {
01915                 ip = (int *)realloc(ip, size3 * sizeof(*ip));
01916                 if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
01917                 bufsize = (int) (size3 * 1.2);
01918                 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
01919                 if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
01920                 oldsize = *size;
01921         }
01922         buf[0] = buf[1] = buf[2] = 0;
01923 
01924         xtc_int(mf, &(minint[0]));
01925         xtc_int(mf, &(minint[1]));
01926         xtc_int(mf, &(minint[2]));
01927 
01928         xtc_int(mf, &(maxint[0]));
01929         xtc_int(mf, &(maxint[1]));
01930         xtc_int(mf, &(maxint[2]));
01931                 
01932         sizeint[0] = maxint[0] - minint[0]+1;
01933         sizeint[1] = maxint[1] - minint[1]+1;
01934         sizeint[2] = maxint[2] - minint[2]+1;
01935         
01936         /* check if one of the sizes is to big to be multiplied */
01937         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
01938                 bitsizeint[0] = xtc_sizeofint(sizeint[0]);
01939                 bitsizeint[1] = xtc_sizeofint(sizeint[1]);
01940                 bitsizeint[2] = xtc_sizeofint(sizeint[2]);
01941                 bitsize = 0; /* flag the use of large sizes */
01942         } else {
01943                 bitsize = xtc_sizeofints(3, sizeint);
01944         }
01945 
01946         xtc_int(mf, &smallidx);
01947         smaller = xtc_magicints[FIRSTIDX > smallidx - 1 ? FIRSTIDX : smallidx - 1] / 2;
01948         small = xtc_magicints[smallidx] / 2;
01949         sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx];
01950 
01951         /* check for zero values that would yield corrupted data */
01952         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
01953                 printf("XTC corrupted, sizesmall==0 (case 1)\n");
01954                 return -1;
01955         }
01956 
01957 
01958         /* buf[0] holds the length in bytes */
01959         if (xtc_int(mf, &(buf[0])) < 0) return -1;
01960 
01961         if (xtc_data(mf, (char *) &buf[3], (int) buf[0]) < 0) return -1;
01962 
01963         buf[0] = buf[1] = buf[2] = 0;
01964 
01965         lfp = fp;
01966         inv_precision = 1.0f / (*precision);
01967         run = 0;
01968         i = 0;
01969         lip = ip;
01970         while (i < lsize) {
01971                 thiscoord = (int *)(lip) + i * 3;
01972 
01973                 if (bitsize == 0) {
01974                         thiscoord[0] = xtc_receivebits(buf, bitsizeint[0]);
01975                         thiscoord[1] = xtc_receivebits(buf, bitsizeint[1]);
01976                         thiscoord[2] = xtc_receivebits(buf, bitsizeint[2]);
01977                 } else {
01978                         xtc_receiveints(buf, 3, bitsize, sizeint, thiscoord);
01979                 }
01980 
01981                 i++;
01982                 thiscoord[0] += minint[0];
01983                 thiscoord[1] += minint[1];
01984                 thiscoord[2] += minint[2];
01985 
01986                 prevcoord[0] = thiscoord[0];
01987                 prevcoord[1] = thiscoord[1];
01988                 prevcoord[2] = thiscoord[2];
01989  
01990 
01991                 flag = xtc_receivebits(buf, 1);
01992                 is_smaller = 0;
01993                 if (flag == 1) {
01994                         run = xtc_receivebits(buf, 5);
01995                         is_smaller = run % 3;
01996                         run -= is_smaller;
01997                         is_smaller--;
01998                 }
01999                 if (run > 0) {
02000                         thiscoord += 3;
02001                         for (k = 0; k < run; k+=3) {
02002                                 xtc_receiveints(buf, 3, smallidx, sizesmall, thiscoord);
02003                                 i++;
02004                                 thiscoord[0] += prevcoord[0] - small;
02005                                 thiscoord[1] += prevcoord[1] - small;
02006                                 thiscoord[2] += prevcoord[2] - small;
02007                                 if (k == 0) {
02008                                         /* interchange first with second atom for better
02009                                          * compression of water molecules
02010                                          */
02011                                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
02012                                         prevcoord[0] = tmp;
02013                                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
02014                                         prevcoord[1] = tmp;
02015                                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
02016                                         prevcoord[2] = tmp;
02017                                         *lfp++ = prevcoord[0] * inv_precision;
02018                                         *lfp++ = prevcoord[1] * inv_precision;
02019                                         *lfp++ = prevcoord[2] * inv_precision;
02020 
02021                                         if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
02022                                                 printf("XTC corrupted, sizesmall==0 (case 2)\n");
02023                                                 return -1;
02024                                         }
02025 
02026                                 } else {
02027                                         prevcoord[0] = thiscoord[0];
02028                                         prevcoord[1] = thiscoord[1];
02029                                         prevcoord[2] = thiscoord[2];
02030                                 }
02031                                 *lfp++ = thiscoord[0] * inv_precision;
02032                                 *lfp++ = thiscoord[1] * inv_precision;
02033                                 *lfp++ = thiscoord[2] * inv_precision;
02034                         }
02035                 } else {
02036                         *lfp++ = thiscoord[0] * inv_precision;
02037                         *lfp++ = thiscoord[1] * inv_precision;
02038                         *lfp++ = thiscoord[2] * inv_precision;          
02039                 }
02040                 smallidx += is_smaller;
02041                 if (is_smaller < 0) {
02042                         small = smaller;
02043                         if (smallidx > FIRSTIDX) {
02044                                 smaller = xtc_magicints[smallidx - 1] /2;
02045                         } else {
02046                                 smaller = 0;
02047                         }
02048                 } else if (is_smaller > 0) {
02049                         smaller = small;
02050                         small = xtc_magicints[smallidx] / 2;
02051                 }
02052                 sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx] ;
02053         }
02054         return 1;
02055 }
02056 
02057 
02058 
02059 #endif

Generated on Fri Mar 29 03:10:08 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002