00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
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;
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
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
00112
00113
00114 typedef struct {
00115 int version;
00116 char title[MAX_TRX_TITLE + 1];
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;
00125 int v_size;
00126 int f_size;
00127 int natoms;
00128 int step;
00129 int nre;
00130 float t;
00131 float lambda;
00132 } trx_hdr;
00133
00134
00135
00136
00137 typedef struct {
00138 FILE * f;
00139 int fmt;
00140 int prec;
00141 int rev;
00142 trx_hdr * trx;
00143
00144 } md_file;
00145
00146
00147
00148 typedef struct {
00149 char title[MAX_MDIO_TITLE + 1];
00150 int natoms;
00151 float timeval;
00152 } md_header;
00153
00154
00155
00156 typedef struct {
00157 float A, B, C, alpha, beta, gamma;
00158 } md_box;
00159
00160
00161
00162 typedef struct {
00163 float *pos;
00164
00165
00166
00167 int natoms;
00168 int step;
00169 float time;
00170 md_box *box;
00171 } md_ts;
00172
00173
00174
00175 typedef struct {
00176 char resid[7];
00177 char resname[7];
00178 int atomnum;
00179 char atomname[7];
00180 float pos[3];
00181
00182 } md_atom;
00183
00184
00185
00186
00187
00188
00189
00190 static md_file *mdio_open(const char *, const int, const int=MDIO_READ);
00191
00192
00193 static int mdio_close(md_file *);
00194
00195
00196
00197 static int mdio_header(md_file *, md_header *);
00198 static int mdio_timestep(md_file *, md_ts *);
00199
00200
00201
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
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
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
00224 static int xtc_int(md_file *, int *);
00225 static int xtc_float(md_file *, float *);
00226
00227
00228
00229
00230 static int xtc_timestep(md_file *, md_ts *);
00231 static int xtc_3dfcoord(md_file *, float *, int *, float *);
00232
00233
00234
00235 static int mdio_errno(void);
00236 static const char *mdio_errmsg(int);
00237 static int mdio_seterror(int);
00238
00239
00240
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
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
00285
00286
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
00296 mf = (md_file *) malloc(sizeof(md_file));
00297 if (!mf) {
00298 mdio_seterror(MDIO_BADMALLOC);
00299 return NULL;
00300 }
00301
00302
00303 memset(mf, 0, sizeof(md_file));
00304
00305
00306 if (!fmt) {
00307 char *p;
00308 int n;
00309
00310
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
00319 for (n = 1; mdio_fmtexts[n]; n++)
00320 if (!strcasecmp(p, mdio_fmtexts[n])) break;
00321
00322
00323 if (!mdio_fmtexts[n]) {
00324 free(mf);
00325 mdio_seterror(MDIO_UNKNOWNFMT);
00326 return NULL;
00327 }
00328
00329
00330 mf->fmt = n;
00331 }
00332 else {
00333 mf->fmt = fmt;
00334 }
00335
00336
00337
00338 switch (mf->fmt) {
00339 case MDFMT_GRO:
00340 case MDFMT_G96:
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:
00349
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:
00358
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
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
00380 mdio_seterror(MDIO_SUCCESS);
00381 return mf;
00382 }
00383
00384
00385
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
00392 if (mf->trx) free(mf->trx);
00393 free(mf);
00394 return mdio_seterror(MDIO_SUCCESS);
00395 }
00396
00397
00398
00399 static int mdio_errno(void) {
00400 return mdio_errcode;
00401 }
00402
00403
00404
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
00412
00413 static int mdio_seterror(int code) {
00414 mdio_errcode = code;
00415 return code ? -1 : 0;
00416 }
00417
00418
00419
00420
00421
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
00426 fgets(buf, n, mf->f);
00427
00428
00429 if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
00430
00431
00432 if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
00433
00434
00435 if (buf[0] == '#') return mdio_readline(mf,buf,n,strip);
00436
00437
00438 if (strip) strip_white(buf);
00439
00440 return strlen(buf);
00441 }
00442
00443
00444
00445
00446
00447 static int strip_white(char *buf) {
00448 int i, j, k;
00449
00450
00451 if (!buf) return -1;
00452 if (!strlen(buf)) return -1;
00453
00454
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
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
00476
00477
00478
00479
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
00496
00497
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
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
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
00519
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
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:
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
00558 if (xtc_int(mf, &n) < 0) return -1;
00559 if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
00560
00561
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
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:
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
00604 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
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
00625 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00626
00627
00628
00629 if ((p = (char *) strstr(buf, "t="))) {
00630 char *q = p;
00631 *(q--) = 0;
00632
00633
00634
00635 p += 2;
00636 strip_white(p);
00637 strip_white(buf);
00638
00639
00640 if (timeval) *timeval = (float) atof(p);
00641 }
00642 else {
00643
00644
00645 if (timeval) *timeval = 0;
00646 strip_white(buf);
00647 }
00648
00649
00650 if (title && titlelen) strncpy(title, buf, titlelen);
00651
00652
00653 while (strcasecmp(buf, "END"))
00654 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00655
00656
00657 return mdio_seterror(MDIO_SUCCESS);
00658 }
00659
00660
00661
00662
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
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
00698 if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);
00699
00700
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
00733
00734
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
00743 if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00744
00745
00746
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
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00785
00786 if (!strcasecmp(buf, "TITLE")) {
00787
00788 while (strcasecmp(buf, "END")) {
00789 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00790 }
00791
00792
00793 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00794 }
00795
00796
00797 if (!strcasecmp(buf, "TIMESTEP")) {
00798
00799 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00800
00801
00802 n = sscanf(buf, "%d %f", &ts->step, &ts->time);
00803 if (n != 2) return mdio_seterror(MDIO_BADFORMAT);
00804
00805
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
00811 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00812 }
00813 else {
00814
00815 ts->step = 0;
00816 ts->time = 0;
00817 }
00818
00819
00820
00821 if (!strcasecmp(buf, "POSITIONRED")) {
00822
00823
00824 i = 0;
00825 while (i < ts->natoms) {
00826
00827 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00828 return -1;
00829
00830
00831 if (!strcasecmp(buf, "END"))
00832 return mdio_seterror(MDIO_BADFORMAT);
00833
00834
00835 n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
00836
00837
00838 if (n == 3) {
00839 pos[0] *= ANGS_PER_NM;
00840 pos[1] *= ANGS_PER_NM;
00841 pos[2] *= ANGS_PER_NM;
00842
00843
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
00853
00854
00855
00856
00857
00858
00859 i = 0;
00860 while (i < ts->natoms) {
00861
00862 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00863 return -1;
00864
00865
00866 strcpy(stripbuf, buf);
00867 strip_white(stripbuf);
00868 if (!strcasecmp(stripbuf, "END"))
00869 return mdio_seterror(MDIO_BADFORMAT);
00870
00871
00872 n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
00873 &pos[0], &pos[1], &pos[2]);
00874
00875
00876 if (n == 3) {
00877 pos[0] *= ANGS_PER_NM;
00878 pos[1] *= ANGS_PER_NM;
00879 pos[2] *= ANGS_PER_NM;
00880
00881
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
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
00899
00900
00901
00902
00903
00904 fpos = ftell(mf->f);
00905
00906
00907 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00908
00909
00910 if (mdio_errcode == MDIO_EOF)
00911 return mdio_seterror(MDIO_SUCCESS);
00912 else
00913 return -1;
00914 }
00915
00916
00917 if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
00918
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
00926
00927 fpos = ftell(mf->f);
00928
00929
00930 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00931 }
00932
00933
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
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
00967
00968
00969
00970 fseek(mf->f, fpos, SEEK_SET);
00971 }
00972
00973
00974 return mdio_seterror(MDIO_SUCCESS);
00975 }
00976
00977
00978
00979
00980
00981
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
00989 if (!mf)
00990 return mdio_seterror(MDIO_BADPARAMS);
00991
00992
00993 fpos = ftell(mf->f);
00994
00995
00996 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
00997
00998
00999
01000 if ((p = (char *) strstr(buf, "t="))) {
01001 char *q = p;
01002 *(q--) = 0;
01003
01004
01005
01006 p += 2;
01007 strip_white(p);
01008 strip_white(buf);
01009
01010
01011 if (timeval) *timeval = (float) atof(p);
01012 } else {
01013
01014 if (timeval) *timeval = 0;
01015 }
01016
01017
01018 if (title && titlelen) strncpy(title, buf, titlelen);
01019
01020
01021 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
01022
01023
01024 if (natoms && (!(*natoms = atoi(buf))))
01025 return mdio_seterror(MDIO_BADFORMAT);
01026
01027
01028
01029
01030 if (rewind)
01031 fseek(mf->f, fpos, SEEK_SET);
01032
01033 return 0;
01034 }
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
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
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
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
01086 strip_white(atomnum);
01087 ma->atomnum = atoi(atomnum);
01088
01089
01090 ma->pos[0] *= ANGS_PER_NM;
01091 ma->pos[1] *= ANGS_PER_NM;
01092 ma->pos[2] *= ANGS_PER_NM;
01093
01094
01095 strip_white(ma->atomname);
01096 strip_white(ma->resname);
01097 strip_white(ma->resid);
01098
01099 return 0;
01100 }
01101
01102
01103
01104
01105
01106
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
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
01159
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
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
01183
01184
01185
01186
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
01197 fpos = ftell(mf->f);
01198
01199
01200
01201 hdr = mf->trx;
01202 if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);
01203
01204
01205 if (trx_int(mf, &magic) < 0) return -1;
01206 if (magic != TRX_MAGIC) {
01207
01208 swap4_aligned(&magic, 1);
01209 if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01210
01211
01212 mf->rev = 1;
01213 }
01214
01215
01216
01217
01218
01219
01220
01221
01222 if(mf->fmt!=MDFMT_TRJ) {
01223
01224
01225
01226
01227
01228 if (trx_int(mf, &hdr->version) < 0) return -1;
01229 }
01230
01231
01232 if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
01233 return -1;
01234
01235
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
01251 if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);
01252
01253
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
01261
01262
01263 return mdio_seterror(MDIO_BADPRECISION);
01264 }
01265
01266
01267 if (trx_real(mf, &hdr->t) < 0) return -1;
01268 if (trx_real(mf, &hdr->lambda) < 0) return -1;
01269
01270
01271 if (rewind) fseek(mf->f, fpos, SEEK_SET);
01272
01273 return 0;
01274 }
01275
01276
01277
01278
01279 static int trx_int(md_file *mf, int *y) {
01280 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01281
01282
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
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
01315
01316
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
01372
01373
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
01407
01408
01409
01410
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
01441
01442
01443
01444
01445
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
01478
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
01489 if (trx_header(mf) < 0) return -1;
01490
01491
01492 hdr = mf->trx;
01493 if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
01494
01495 if (hdr->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
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
01562
01563 static int put_trx_int(md_file *mf, int y) {
01564 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01565
01566
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
01577
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
01590
01591 static int put_trx_string(md_file *mf, const char *s) {
01592 if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
01593
01594
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
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
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
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
01639
01640
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
01649
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
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;
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
01694 if (xtc_int(mf, &n) < 0) return -1;
01695 if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01696
01697
01698 if (xtc_int(mf, &n) < 0) return -1;
01699 ts->natoms = n;
01700
01701
01702 if (xtc_int(mf, &n) < 0) return -1;
01703 ts->step = n;
01704
01705
01706 if (xtc_float(mf, &f) < 0) return -1;
01707 ts->time = f;
01708
01709
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
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
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
01743
01744
01746
01747
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
01759 #define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))
01760
01761
01762
01763
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
01777
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
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
01841
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
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
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
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;
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
01952 if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
01953 printf("XTC corrupted, sizesmall==0 (case 1)\n");
01954 return -1;
01955 }
01956
01957
01958
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
02009
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