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 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 { 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 (strip) strip_white(buf);
00436
00437 return strlen(buf);
00438 }
00439
00440
00441
00442
00443
00444 static int strip_white(char *buf) {
00445 int i, j, k;
00446
00447
00448 if (!buf) return -1;
00449 if (!strlen(buf)) return -1;
00450
00451
00452 for (i = strlen(buf) - 1;
00453 buf[i] == ' ' || buf[i] == '\t' ||
00454 buf[i] == '\n' || buf[i] == '\r';
00455 i--)
00456 buf[i] = 0;
00457
00458
00459 for (i = 0; buf[i] == ' ' || buf[i] == '\t' ||
00460 buf[i] == '\n' || buf[i] == '\r'; i++);
00461 if (i) {
00462 k = 0;
00463 for (j = i; buf[j]; j++)
00464 buf[k++] = buf[j];
00465 buf[k] = 0;
00466 }
00467
00468 return strlen(buf);
00469 }
00470
00471
00472
00473
00474
00475
00476
00477 static int mdio_tsfree(md_ts *ts, int holderror) {
00478 if (!ts) {
00479 if (holderror) return -1;
00480 else return mdio_seterror(MDIO_BADPARAMS);
00481 }
00482
00483 if (ts->pos && ts->natoms > 0) free(ts->pos);
00484
00485 if (ts->box) free(ts->box);
00486
00487 if (holderror) return -1;
00488 else return mdio_seterror(MDIO_SUCCESS);
00489 }
00490
00491
00492
00493
00494
00495 static int mdio_readbox(md_box *box, float *x, float *y, float *z) {
00496 float A, B, C;
00497
00498 if (!box) {
00499 return mdio_seterror(MDIO_BADPARAMS);
00500 }
00501
00502
00503 A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] ) * ANGS_PER_NM;
00504 B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] ) * ANGS_PER_NM;
00505 C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] ) * ANGS_PER_NM;
00506 if ((A<=0) || (B<=0) || (C<=0)) {
00507
00508 box->A = box->B = box->C = 0;
00509 box->alpha = box->beta = box->gamma = 90;
00510 } else {
00511 box->A = A;
00512 box->B = B;
00513 box->C = C;
00514
00515
00516
00517 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;
00518 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;
00519 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;
00520 }
00521 return mdio_seterror(MDIO_SUCCESS);
00522 }
00523
00524
00525
00526 static int mdio_header(md_file *mf, md_header *mdh) {
00527 int n;
00528 if (!mf || !mdh) return mdio_seterror(MDIO_BADPARAMS);
00529 if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00530
00531 switch (mf->fmt) {
00532 case MDFMT_GRO:
00533 if (gro_header(mf, mdh->title, MAX_MDIO_TITLE,
00534 &mdh->timeval, &mdh->natoms, 1) < 0)
00535 return -1;
00536 return 0;
00537
00538 case MDFMT_TRR:
00539 case MDFMT_TRJ:
00540 if (trx_header(mf, 1) < 0) return -1;
00541 mdh->natoms = mf->trx->natoms;
00542 mdh->timeval = (float) mf->trx->t;
00543 strncpy(mdh->title, mf->trx->title, MAX_MDIO_TITLE);
00544 return 0;
00545
00546 case MDFMT_G96:
00547 if (g96_header(mf, mdh->title, MAX_MDIO_TITLE,
00548 &mdh->timeval) < 0) return -1;
00549 mdh->natoms = -1;
00550 return 0;
00551
00552 case MDFMT_XTC:
00553 memset(mdh, 0, sizeof(md_header));
00554
00555 if (xtc_int(mf, &n) < 0) return -1;
00556 if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
00557
00558
00559 if (xtc_int(mf, &n) < 0) return -1;
00560 mdh->natoms = n;
00561 rewind(mf->f);
00562 return 0;
00563
00564 default:
00565 return mdio_seterror(MDIO_UNKNOWNFMT);
00566 }
00567 }
00568
00569
00570
00571 static int mdio_timestep(md_file *mf, md_ts *ts) {
00572 if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00573 if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
00574
00575 switch (mf->fmt) {
00576 case MDFMT_GRO:
00577 return gro_timestep(mf, ts);
00578
00579 case MDFMT_TRR:
00580 case MDFMT_TRJ:
00581 return trx_timestep(mf, ts);
00582
00583 case MDFMT_G96:
00584 return g96_timestep(mf, ts);
00585
00586 case MDFMT_XTC:
00587 return xtc_timestep(mf, ts);
00588
00589 default:
00590 return mdio_seterror(MDIO_UNKNOWNFMT);
00591 }
00592 }
00593
00594
00595
00596 static int g96_header(md_file *mf, char *title, int titlelen, float *timeval) {
00597 char buf[MAX_G96_LINE + 1];
00598 char *p;
00599
00600
00601 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00619 if (strcasecmp(buf, "TITLE")) return mdio_seterror(MDIO_BADFORMAT);
00620
00621
00622 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00623
00624
00625
00626 if ((p = (char *) strstr(buf, "t="))) {
00627 char *q = p;
00628 *(q--) = 0;
00629
00630
00631
00632 p += 2;
00633 strip_white(p);
00634 strip_white(buf);
00635
00636
00637 if (timeval) *timeval = (float) atof(p);
00638 }
00639 else {
00640
00641
00642 if (timeval) *timeval = 0;
00643 strip_white(buf);
00644 }
00645
00646
00647 if (title && titlelen) strncpy(title, buf, titlelen);
00648
00649
00650 while (strcasecmp(buf, "END"))
00651 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00652
00653
00654 return mdio_seterror(MDIO_SUCCESS);
00655 }
00656
00657
00658
00659
00660 static int g96_countatoms(md_file *mf) {
00661 char buf[MAX_G96_LINE + 1];
00662 int natoms;
00663 int n;
00664 long fpos;
00665 float lastf;
00666
00667 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00668
00669 fpos = ftell(mf->f);
00670
00671 natoms = 0;
00672 for (;;) {
00673 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00674 break;
00675 n = sscanf(buf, "%*6c%*6c%*6c%*6c %*f %*f %f", &lastf);
00676 if (n == 1) natoms++;
00677 else {
00678 strip_white(buf);
00679 if (!strcasecmp(buf, "END")) break;
00680 }
00681 }
00682
00683 fseek(mf->f, fpos, SEEK_SET);
00684 return natoms;
00685 }
00686
00687
00688
00689 static int g96_rec(md_file *mf, md_atom *ma) {
00690 char buf[MAX_G96_LINE + 1];
00691 char atomnum[7];
00692 int n;
00693
00694
00695 if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);
00696
00697
00698 do {
00699 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0) return -1;
00700 } while (buf[0] == '#' || strlen(buf) == 0);
00701
00702 n = sscanf(buf, "%6c%6c%6c%6c %f %f %f",
00703 ma->resid, ma->resname, ma->atomname, atomnum,
00704 &ma->pos[0], &ma->pos[1], &ma->pos[2]);
00705 if (n == 7) {
00706 atomnum[6] = 0;
00707 ma->resid[6] = 0;
00708 ma->resname[6] = 0;
00709 ma->atomname[6] = 0;
00710
00711 strip_white(atomnum);
00712 strip_white(ma->resid);
00713 strip_white(ma->resname);
00714 strip_white(ma->atomname);
00715
00716 ma->atomnum = atoi(atomnum);
00717
00718 ma->pos[0] *= ANGS_PER_NM;
00719 ma->pos[1] *= ANGS_PER_NM;
00720 ma->pos[2] *= ANGS_PER_NM;
00721
00722 return 0;
00723 }
00724
00725 return mdio_seterror(MDIO_BADFORMAT);
00726 }
00727
00728
00729
00730
00731
00732 static int g96_timestep(md_file *mf, md_ts *ts) {
00733 char buf[MAX_G96_LINE + 1];
00734 char stripbuf[MAX_G96_LINE + 1];
00735 float pos[3], x[3], y[3], z[3], *currAtom;
00736 long fpos;
00737 int n, i, boxItems;
00738
00739
00740 if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
00741
00742
00743
00744 ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
00745 if (!ts->pos) {
00746 return mdio_seterror(MDIO_BADMALLOC);
00747 }
00748 currAtom = ts->pos;
00749
00750
00751
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 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00782
00783 if (!strcasecmp(buf, "TITLE")) {
00784
00785 while (strcasecmp(buf, "END")) {
00786 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00787 }
00788
00789
00790 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00791 }
00792
00793
00794 if (!strcasecmp(buf, "TIMESTEP")) {
00795
00796 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00797
00798
00799 n = sscanf(buf, "%d %f", &ts->step, &ts->time);
00800 if (n != 2) return mdio_seterror(MDIO_BADFORMAT);
00801
00802
00803 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00804 if (strcasecmp(buf, "END"))
00805 return mdio_seterror(MDIO_BADFORMAT);
00806
00807
00808 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00809 }
00810 else {
00811
00812 ts->step = 0;
00813 ts->time = 0;
00814 }
00815
00816
00817
00818 if (!strcasecmp(buf, "POSITIONRED")) {
00819
00820
00821 i = 0;
00822 while (i < ts->natoms) {
00823
00824 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00825 return -1;
00826
00827
00828 if (!strcasecmp(buf, "END"))
00829 return mdio_seterror(MDIO_BADFORMAT);
00830
00831
00832 n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
00833
00834
00835 if (n == 3) {
00836 pos[0] *= ANGS_PER_NM;
00837 pos[1] *= ANGS_PER_NM;
00838 pos[2] *= ANGS_PER_NM;
00839
00840
00841 memcpy(currAtom, pos, sizeof(float) * 3);
00842 currAtom += 3;
00843 i++;
00844 }
00845 }
00846 }
00847 else if (!strcasecmp(buf, "POSITION") || !strcasecmp(buf, "REFPOSITION")) {
00848
00849
00850
00851
00852
00853
00854
00855
00856 i = 0;
00857 while (i < ts->natoms) {
00858
00859 if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
00860 return -1;
00861
00862
00863 strcpy(stripbuf, buf);
00864 strip_white(stripbuf);
00865 if (!strcasecmp(stripbuf, "END"))
00866 return mdio_seterror(MDIO_BADFORMAT);
00867
00868
00869 n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
00870 &pos[0], &pos[1], &pos[2]);
00871
00872
00873 if (n == 3) {
00874 pos[0] *= ANGS_PER_NM;
00875 pos[1] *= ANGS_PER_NM;
00876 pos[2] *= ANGS_PER_NM;
00877
00878
00879 memcpy(currAtom, pos, sizeof(float) * 3);
00880 currAtom += 3;
00881 i++;
00882 }
00883 }
00884 }
00885 else {
00886 return mdio_seterror(MDIO_BADFORMAT);
00887 }
00888
00889
00890 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00891 return -1;
00892 if (strcasecmp(buf, "END"))
00893 return mdio_seterror(MDIO_BADFORMAT);
00894
00895
00896
00897
00898
00899
00900
00901 fpos = ftell(mf->f);
00902
00903
00904 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00905
00906
00907 if (mdio_errcode == MDIO_EOF)
00908 return mdio_seterror(MDIO_SUCCESS);
00909 else
00910 return -1;
00911 }
00912
00913
00914 if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
00915
00916 for (;;) {
00917 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
00918 return -1;
00919 if (!strcasecmp(buf, "END")) break;
00920 }
00921
00922
00923
00924 fpos = ftell(mf->f);
00925
00926
00927 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00928 }
00929
00930
00931 if (!strcasecmp(buf, "BOX")) {
00932 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
00933 boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f",
00934 &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
00935 if (boxItems == 3) {
00936 x[1] = x[2] = 0;
00937 y[0] = y[2] = 0;
00938 z[0] = z[1] = 0;
00939 }
00940 else if (boxItems != 9)
00941 return mdio_seterror(MDIO_BADFORMAT);
00942
00943
00944 ts->box = (md_box *) malloc(sizeof(md_box));
00945 if (mdio_readbox(ts->box, x, y, z) < 0) {
00946 free(ts->box);
00947 ts->box = NULL;
00948 return mdio_seterror(MDIO_BADFORMAT);
00949 }
00950
00951 if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
00952 free(ts->box);
00953 ts->box = NULL;
00954 return -1;
00955 }
00956 if (strcasecmp(buf, "END")) {
00957 free(ts->box);
00958 ts->box = NULL;
00959 return mdio_seterror(MDIO_BADFORMAT);
00960 }
00961 }
00962 else {
00963
00964
00965
00966
00967 fseek(mf->f, fpos, SEEK_SET);
00968 }
00969
00970
00971 return mdio_seterror(MDIO_SUCCESS);
00972 }
00973
00974
00975
00976
00977
00978
00979 static int gro_header(md_file *mf, char *title, int titlelen, float *timeval,
00980 int *natoms, int rewind) {
00981 char buf[MAX_GRO_LINE + 1];
00982 long fpos;
00983 char *p;
00984
00985
00986 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
00987
00988
00989 fpos = ftell(mf->f);
00990
00991
00992 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
00993
00994
00995
00996 if ((p = (char *) strstr(buf, "t="))) {
00997 char *q = p;
00998 *(q--) = 0;
00999
01000
01001
01002 p += 2;
01003 strip_white(p);
01004 strip_white(buf);
01005
01006
01007 if (timeval) *timeval = (float) atof(p);
01008 }
01009 else {
01010
01011 if (timeval) *timeval = 0;
01012 }
01013
01014
01015 if (title && titlelen) strncpy(title, buf, titlelen);
01016
01017
01018 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
01019
01020
01021 if (natoms) if (!(*natoms = atoi(buf)))
01022 return mdio_seterror(MDIO_BADFORMAT);
01023
01024
01025
01026
01027 if (rewind) fseek(mf->f, fpos, SEEK_SET);
01028
01029
01030 return 0;
01031 }
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045 static int gro_rec(md_file *mf, md_atom *ma) {
01046 char buf[MAX_GRO_LINE + 1];
01047 char atomnum[6];
01048 int n;
01049
01050 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01051
01052 do {
01053 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) return -1;
01054 } while (buf[0] == '#' || !strlen(buf));
01055
01056
01057 n = sscanf(buf, "%5c%5c%5c%5c%f %f %f", ma->resid,
01058 ma->resname, ma->atomname, atomnum, ma->pos,
01059 &ma->pos[1], &ma->pos[2]);
01060 if (n != 7) return mdio_seterror(MDIO_BADFORMAT);
01061
01062
01063 ma->resname[5] = 0;
01064 ma->atomname[5] = 0;
01065 ma->resid[5] = 0;
01066 atomnum[5] = 0;
01067
01068
01069 strip_white(atomnum);
01070 ma->atomnum = atoi(atomnum);
01071
01072
01073 ma->pos[0] *= ANGS_PER_NM;
01074 ma->pos[1] *= ANGS_PER_NM;
01075 ma->pos[2] *= ANGS_PER_NM;
01076
01077
01078 strip_white(ma->atomname);
01079 strip_white(ma->resname);
01080 strip_white(ma->resid);
01081
01082 return 0;
01083 }
01084
01085
01086
01087
01088
01089
01090
01091 static int gro_timestep(md_file *mf, md_ts *ts) {
01092 char buf[MAX_GRO_LINE + 1];
01093 long coord;
01094 int i, n, boxItems;
01095 float x[3], y[3], z[3];
01096
01097 if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
01098
01099 if (gro_header(mf, NULL, 0, &ts->time, &ts->natoms, 0) < 0)
01100 return -1;
01101 ts->pos = (float *) malloc(3 * sizeof(float) * ts->natoms);
01102 if (!ts->pos)
01103 return mdio_seterror(MDIO_BADMALLOC);
01104
01105 coord = 0;
01106 for (i = 0; i < ts->natoms; i++) {
01107 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01108 free(ts->pos);
01109 return -1;
01110 }
01111
01112 n = sscanf(buf, "%*5c%*5c%*5c%*5c%f %f %f",
01113 &ts->pos[coord], &ts->pos[coord + 1],
01114 &ts->pos[coord + 2]);
01115
01116 ts->pos[coord] *= ANGS_PER_NM;
01117 ts->pos[coord + 1] *= ANGS_PER_NM;
01118 ts->pos[coord + 2] *= ANGS_PER_NM;
01119
01120 if (n != 3) return mdio_seterror(MDIO_BADFORMAT);
01121 coord += 3;
01122 }
01123
01124
01125 if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
01126 free(ts->pos);
01127 return -1;
01128 }
01129 boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f",
01130 &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
01131
01132
01133 if (boxItems == 3) {
01134 x[1] = x[2] = 0;
01135 y[0] = y[2] = 0;
01136 z[0] = z[1] = 0;
01137 }
01138 else if (boxItems != 9) {
01139 free(ts->pos);
01140 return -1;
01141 }
01142
01143
01144 ts->box = (md_box *) malloc(sizeof(md_box));
01145 if (mdio_readbox(ts->box, x, y, z) < 0) {
01146 free(ts->pos);
01147 free(ts->box);
01148 ts->box = NULL;
01149 return -1;
01150 }
01151
01152 return 0;
01153 }
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163 static int trx_header(md_file *mf, int rewind) {
01164 int magic;
01165 trx_hdr *hdr;
01166 long fpos;
01167
01168 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01169
01170
01171 fpos = ftell(mf->f);
01172
01173
01174
01175 hdr = mf->trx;
01176 if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);
01177
01178
01179 if (trx_int(mf, &magic) < 0) return -1;
01180 if (magic != TRX_MAGIC) {
01181
01182 swap4_aligned(&magic, 1);
01183 if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
01184
01185
01186 mf->rev = 1;
01187 }
01188
01189
01190
01191
01192
01193
01194
01195
01196 if(mf->fmt!=MDFMT_TRJ) {
01197
01198
01199
01200
01201
01202 if (trx_int(mf, &hdr->version) < 0) return -1;
01203 }
01204
01205
01206 if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
01207 return -1;
01208
01209
01210 if (trx_int(mf, &hdr->ir_size) < 0) return -1;
01211 if (trx_int(mf, &hdr->e_size) < 0) return -1;
01212 if (trx_int(mf, &hdr->box_size) < 0) return -1;
01213 if (trx_int(mf, &hdr->vir_size) < 0) return -1;
01214 if (trx_int(mf, &hdr->pres_size) < 0) return -1;
01215 if (trx_int(mf, &hdr->top_size) < 0) return -1;
01216 if (trx_int(mf, &hdr->sym_size) < 0) return -1;
01217 if (trx_int(mf, &hdr->x_size) < 0) return -1;
01218 if (trx_int(mf, &hdr->v_size) < 0) return -1;
01219 if (trx_int(mf, &hdr->f_size) < 0) return -1;
01220 if (trx_int(mf, &hdr->natoms) < 0) return -1;
01221 if (trx_int(mf, &hdr->step) < 0) return -1;
01222 if (trx_int(mf, &hdr->nre) < 0) return -1;
01223
01224
01225 if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);
01226
01227
01228 if (hdr->x_size) mf->prec = hdr->x_size / (hdr->natoms * 3);
01229 else if (hdr->v_size) mf->prec = hdr->v_size / (hdr->natoms * 3);
01230 else if (hdr->f_size) mf->prec = hdr->f_size / (hdr->natoms * 3);
01231 else return mdio_seterror(MDIO_BADPRECISION);
01232
01233 if (mf->prec != sizeof(float) && mf->prec != sizeof(double)) {
01234
01235
01236
01237 return mdio_seterror(MDIO_BADPRECISION);
01238 }
01239
01240
01241 if (trx_real(mf, &hdr->t) < 0) return -1;
01242 if (trx_real(mf, &hdr->lambda) < 0) return -1;
01243
01244
01245 if (rewind) fseek(mf->f, fpos, SEEK_SET);
01246
01247 return 0;
01248 }
01249
01250
01251
01252
01253 static int trx_int(md_file *mf, int *y) {
01254 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01255
01256
01257 if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
01258
01259 if (y) {
01260 if (fread(y, 4, 1, mf->f) != 1)
01261 return mdio_seterror(MDIO_IOERROR);
01262 if (mf->rev) swap4_aligned(y, 1);
01263 }
01264 else if (fseek(mf->f, 4, SEEK_CUR) != 0)
01265 return mdio_seterror(MDIO_IOERROR);
01266
01267 return mdio_seterror(MDIO_SUCCESS);
01268 }
01269
01270
01271
01272
01273
01274 static int trx_real(md_file *mf, float *y) {
01275 double x;
01276
01277 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01278
01279 switch (mf->prec) {
01280 case sizeof(float):
01281 if (!y) {
01282 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01283 return mdio_seterror(MDIO_IOERROR);
01284 } else {
01285 if (fread(y, mf->prec, 1, mf->f) != 1)
01286 return mdio_seterror(MDIO_IOERROR);
01287 if (mf->rev) swap4_aligned(y, 1);
01288 }
01289 return mdio_seterror(MDIO_SUCCESS);
01290
01291 case sizeof(double):
01292 if (!y) {
01293 if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
01294 return mdio_seterror(MDIO_IOERROR);
01295 } else {
01296 if (fread(&x, mf->prec, 1, mf->f) != 1)
01297 return mdio_seterror(MDIO_IOERROR);
01298 if (mf->rev) swap8_aligned(&x, 1);
01299 *y = (float) x;
01300 }
01301 return mdio_seterror(MDIO_SUCCESS);
01302
01303 default:
01304 return mdio_seterror(MDIO_BADPRECISION);
01305 }
01306
01307 }
01308
01309
01310
01311
01312
01313 static int trx_rvector(md_file *mf, float *vec) {
01314 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01315
01316 if (!vec) {
01317 if (trx_real(mf, NULL) < 0) return -1;
01318 if (trx_real(mf, NULL) < 0) return -1;
01319 if (trx_real(mf, NULL) < 0) return -1;
01320 return mdio_seterror(MDIO_SUCCESS);
01321 } else {
01322 if (trx_real(mf, &vec[0]) < 0) return -1;
01323 if (trx_real(mf, &vec[1]) < 0) return -1;
01324 if (trx_real(mf, &vec[2]) < 0) return -1;
01325 return mdio_seterror(MDIO_SUCCESS);
01326 }
01327 }
01328
01329
01330
01331
01332
01333
01334
01335 static int trx_string(md_file *mf, char *str, int max) {
01336 int size;
01337 size_t ssize;
01338
01339 if (!mf) return mdio_seterror(MDIO_BADPARAMS);
01340
01341 if (trx_int(mf, &size) < 0) return -1;
01342 ssize = (size_t)size;
01343
01344 if (str && size <= max) {
01345 if (fread(str, 1, size, mf->f) != ssize)
01346 return mdio_seterror(MDIO_IOERROR);
01347 str[size] = 0;
01348 return size;
01349 } else if (str) {
01350 if (fread(str, 1, max, mf->f) != ssize)
01351 return mdio_seterror(MDIO_IOERROR);
01352 if (fseek(mf->f, size - max, SEEK_CUR) != 0)
01353 return mdio_seterror(MDIO_IOERROR);
01354 str[max] = 0;
01355 return max;
01356 } else {
01357 if (fseek(mf->f, size, SEEK_CUR) != 0)
01358 return mdio_seterror(MDIO_IOERROR);
01359