src/gmxlib/libxdrf.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <limits.h>
  40 #include <math.h>
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  44 #include "statutil.h"
  45 #include "xdrf.h"
  46
  47
  48
  49
  50 #ifdef HAVE_FSEEKO
  51 #  define gmx_fseek(A,B,C) fseeko(A,B,C)
  52 #  define gmx_ftell(A) ftello(A)
  53 #  define gmx_off_t off_t
  54 #else
  55 #  define gmx_fseek(A,B,C) fseek(A,B,C)
  56 #  define gmx_ftell(A) ftell(A)
  57 #  define gmx_off_t int
  58 #endif
  59
  60
  61 #define MAXID 256
  62 static FILE *xdrfiles[MAXID];
  63 static XDR *xdridptr[MAXID];
  64 static char xdrmodes[MAXID];
  65 static unsigned int cnt;
  66
  67 /* This is just for clarity - it can never be anything but 4! */
  68 #define XDR_INT_SIZE 4
  69
  70 #ifdef GMX_FORTRAN
  71
  72 typedef void (* F77_FUNC(xdrfproc,XDRFPROC))(int *, void *, int *);
  73
  74 int ftocstr(char *ds, int dl, char *ss, int sl)
  75     /* dst, src ptrs */
  76     /* dst max len */
  77     /* src len */
  78 {
  79     char *p;
  80
  81     p = ss + sl;
  82     while ( --p >= ss && *p == ' ' );
  83     sl = p - ss + 1;
  84     dl--;
  85     ds[0] = 0;
  86     if (sl > dl)
  87       return 1;
  88     while (sl--)
  89       (*ds++ = *ss++);
  90     *ds = '\0';
  91     return 0;
  92 }
  93
  94
  95 int ctofstr(char *ds, int dl, char *ss)
  96      /* dest space */
  97      /* max dest length */
  98      /* src string (0-term) */
  99 {
 100     while (dl && *ss) {
 101         *ds++ = *ss++;
 102         dl--;
 103     }
 104     while (dl--)
 105         *ds++ = ' ';
 106     return 0;
 107 }
 108
 109 void
 110 F77_FUNC(xdrfbool,XDRFBOOL)(int *xdrid, int *pb, int *ret)
 111 {
 112         *ret = xdr_bool(xdridptr[*xdrid], pb);
 113         cnt += XDR_INT_SIZE;
 114 }
 115
 116 void
 117 F77_FUNC(xdrfchar,XDRFCHAR)(int *xdrid, char *cp, int *ret)
 118 {
 119         *ret = xdr_char(xdridptr[*xdrid], cp);
 120         cnt += sizeof(char);
 121 }
 122
 123 void
 124 F77_FUNC(xdrfdouble,XDRFDOUBLE)(int *xdrid, double *dp, int *ret)
 125 {
 126         *ret = xdr_double(xdridptr[*xdrid], dp);
 127         cnt += sizeof(double);
 128 }
 129
 130 void
 131 F77_FUNC(xdrffloat,XDRFFLOAT)(int *xdrid, float *fp, int *ret)
 132 {
 133         *ret = xdr_float(xdridptr[*xdrid], fp);
 134         cnt += sizeof(float);
 135 }
 136
 137 void
 138 F77_FUNC(xdrfint,XDRFINT)(int *xdrid, int *ip, int *ret)
 139 {
 140         *ret = xdr_int(xdridptr[*xdrid], ip);
 141         cnt += XDR_INT_SIZE;
 142 }
 143
 144 F77_FUNC(xdrfshort,XDRFSHORT)(int *xdrid, short *sp, int *ret)
 145 {
 146         *ret = xdr_short(xdridptr[*xdrid], sp);
 147         cnt += sizeof(sp);
 148 }
 149
 150 void
 151 F77_FUNC(xdrfuchar,XDRFUCHAR)(int *xdrid, unsigned char *ucp, int *ret)
 152 {
 153         *ret = xdr_u_char(xdridptr[*xdrid], (u_char *)ucp);
 154         cnt += sizeof(char);
 155 }
 156
 157
 158 void
 159 F77_FUNC(xdrfushort,XDRFUSHORT)(int *xdrid, unsigned short *usp, int *ret)
 160 {
 161         *ret = xdr_u_short(xdridptr[*xdrid], (unsigned short *)usp);
 162         cnt += sizeof(unsigned short);
 163 }
 164
 165 void
 166 F77_FUNC(xdrf3dfcoord,XDRF3DFCOORD)(int *xdrid, float *fp, int *size, float *precision, int *ret)
 167 {
 168         *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
 169 }
 170
 171 void
 172 F77_FUNC(xdrfstring,XDRFSTRING)(int *xdrid, char * sp_ptr,
 173                                 int *maxsize, int *ret, int sp_len)
 174 {
 175         char *tsp;
 176
 177         tsp = (char*) malloc((size_t)(((sp_len) + 1) * sizeof(char)));
 178         if (tsp == NULL) {
 179             *ret = -1;
 180             return;
 181         }
 182         if (ftocstr(tsp, *maxsize+1, sp_ptr, sp_len)) {
 183             *ret = -1;
 184             free(tsp);
 185             return;
 186         }
 187         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (unsigned int) *maxsize);
 188         ctofstr( sp_ptr, sp_len , tsp);
 189         cnt += *maxsize;
 190         free(tsp);
 191 }
 192
 193 void
 194 F77_FUNC(xdrfwrapstring,XDRFWRAPSTRING)(int *xdrid, char *sp_ptr,
 195                                         int *ret, int sp_len)
 196 {
 197         char *tsp;
 198         int maxsize;
 199         maxsize = (sp_len) + 1;
 200         tsp = (char*) malloc((size_t)(maxsize * sizeof(char)));
 201         if (tsp == NULL) {
 202             *ret = -1;
 203             return;
 204         }
 205         if (ftocstr(tsp, maxsize, sp_ptr, sp_len)) {
 206             *ret = -1;
 207             free(tsp);
 208             return;
 209         }
 210         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
 211         ctofstr( sp_ptr, sp_len, tsp);
 212         cnt += maxsize;
 213         free(tsp);
 214 }
 215
 216 void
 217 F77_FUNC(xdrfopaque,XDRFOPAQUE)(int *xdrid, caddr_t *cp, int *ccnt, int *ret)
 218 {
 219         *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
 220         cnt += *ccnt;
 221 }
 222
 223 void
 224 F77_FUNC(xdrfsetpos,XDRFSETPOS)(int *xdrid, int *pos, int *ret)
 225 {
 226         *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
 227 }
 228
 229
 230 void
 231 F77_FUNC(xdrf,XDRF)(int *xdrid, int *pos)
 232 {
 233         *pos = xdr_getpos(xdridptr[*xdrid]);
 234 }
 235
 236 void
 237 F77_FUNC(xdrfvector,XDRFVECTOR)(int *xdrid, char *cp, int *size, F77_FUNC(xdrfproc,XDRFPROC) elproc, int *ret)
 238 {
 239         int lcnt;
 240         cnt = 0;
 241         for (lcnt = 0; lcnt < *size; lcnt++) {
 242                 elproc(xdrid, (cp+cnt) , ret);
 243         }
 244 }
 245
 246
 247 void
 248 F77_FUNC(xdrfclose,XDRFCLOSE)(int *xdrid, int *ret)
 249 {
 250         *ret = xdrclose(xdridptr[*xdrid]);
 251         cnt = 0;
 252 }
 253
 254 void
 255 F77_FUNC(xdrfopen,XDRFOPEN)(int *xdrid, char *fp_ptr, char *mode_ptr,
 256                             int *ret, int fp_len, int mode_len)
 257 {
 258         char fname[512];
 259         char fmode[3];
 260
 261         if (ftocstr(fname, sizeof(fname), fp_ptr, fp_len)) {
 262                 *ret = 0;
 263         }
 264         if (ftocstr(fmode, sizeof(fmode), mode_ptr,
 265                         mode_len)) {
 266                 *ret = 0;
 267         }
 268
 269         *xdrid = xdropen(NULL, fname, fmode);
 270         if (*xdrid == 0)
 271                 *ret = 0;
 272         else
 273                 *ret = 1;
 274 }
 275 #endif /* GMX_FORTRAN */
 276
 277 /*___________________________________________________________________________
 278  |
 279  | what follows are the C routines for opening, closing xdr streams
 280  | and the routine to read/write compressed coordinates together
 281  | with some routines to assist in this task (those are marked
 282  | static and cannot be called from user programs)
 283 */
 284 #define MAXABS INT_MAX-2
 285
 286 #ifndef MIN
 287 #define MIN(x,y) ((x) < (y) ? (x):(y))
 288 #endif
 289 #ifndef MAX
 290 #define MAX(x,y) ((x) > (y) ? (x):(y))
 291 #endif
 292 #ifndef SQR
 293 #define SQR(x) ((x)*(x))
 294 #endif
 295 static int magicints[] = {
 296     0, 0, 0, 0, 0, 0, 0, 0, 0,
 297     8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
 298     80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
 299     812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
 300     8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
 301     82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
 302     832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
 303     8388607, 10568983, 13316085, 16777216 };
 304
 305 #define FIRSTIDX 9
 306 /* note that magicints[FIRSTIDX-1] == 0 */
 307 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
 308
 309
 310 /*__________________________________________________________________________
 311  |
 312  | xdropen - open xdr file
 313  |
 314  | This versions differs from xdrstdio_create, because I need to know
 315  | the state of the file (read or write) so I can use xdr3dfcoord
 316  | in eigther read or write mode, and the file descriptor
 317  | so I can close the file (something xdr_destroy doesn't do).
 318  |
 319 */
 320
 321 int xdropen(XDR *xdrs, const char *filename, const char *type) {
 322     static int init_done = 0;
 323     enum xdr_op lmode;
 324     int xdrid;
 325     char newtype[5];
 326
 327     if (init_done == 0) {
 328         for (xdrid = 1; xdrid < MAXID; xdrid++) {
 329             xdridptr[xdrid] = NULL;
 330         }
 331         init_done = 1;
 332     }
 333     xdrid = 1;
 334     while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
 335         xdrid++;
 336     }
 337     if (xdrid == MAXID) {
 338         return 0;
 339     }
 340     if (*type == 'w' || *type == 'W')
 341     {
 342         xdrmodes[xdrid] = 'w';
 343         strcpy(newtype, "wb+");
 344         lmode = XDR_ENCODE;
 345     }
 346     else if (*type == 'a' || *type == 'A')
 347     {
 348         xdrmodes[xdrid] = 'a';
 349         strcpy(newtype, "ab+");
 350         lmode = XDR_ENCODE;
 351     }
 352     else if (strncasecmp(type, "r+", 2) == 0)
 353     {
 354         xdrmodes[xdrid] = 'a';
 355         strcpy(newtype, "rb+");
 356         lmode = XDR_ENCODE;
 357     }
 358     else
 359     {
 360         xdrmodes[xdrid] = 'r';
 361         strcpy(newtype, "rb");
 362         lmode = XDR_DECODE;
 363     }
 364     xdrfiles[xdrid] = fopen(filename, newtype);
 365
 366     if (xdrfiles[xdrid] == NULL) {
 367         xdrs = NULL;
 368         return 0;
 369     }
 370
 371     /* next test isn't usefull in the case of C language
 372      * but is used for the Fortran interface
 373      * (C users are expected to pass the address of an already allocated
 374      * XDR staructure)
 375      */
 376     if (xdrs == NULL) {
 377         xdridptr[xdrid] = (XDR *) malloc((size_t)sizeof(XDR));
 378         xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
 379     } else {
 380         xdridptr[xdrid] = xdrs;
 381         xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
 382     }
 383     return xdrid;
 384 }
 385
 386 /*_________________________________________________________________________
 387  |
 388  | xdrclose - close a xdr file
 389  |
 390  | This will flush the xdr buffers, and destroy the xdr stream.
 391  | It also closes the associated file descriptor (this is *not*
 392  | done by xdr_destroy).
 393  |
 394 */
 395
 396 int xdrclose(XDR *xdrs) {
 397     int xdrid;
 398     int rc = 0;
 399
 400     if (xdrs == NULL) {
 401         fprintf(stderr, "xdrclose: passed a NULL pointer\n");
 402         exit(1);
 403     }
 404     for (xdrid = 1; xdrid < MAXID && rc==0; xdrid++) {
 405         if (xdridptr[xdrid] == xdrs) {
 406
 407             xdr_destroy(xdrs);
 408             rc = fclose(xdrfiles[xdrid]);
 409             xdridptr[xdrid] = NULL;
 410             return !rc; /* xdr routines return 0 when ok */
 411         }
 412     }
 413     fprintf(stderr, "xdrclose: no such open xdr file\n");
 414     exit(1);
 415
 416     /* to make some compilers happy: */
 417     return 0;
 418 }
 419
 420 FILE *
 421 xdr_get_fp(int xdrid)
 422 {
 423         return xdrfiles[xdrid];
 424 }
 425
 426
 427 /*____________________________________________________________________________
 428  |
 429  | sendbits - encode num into buf using the specified number of bits
 430  |
 431  | This routines appends the value of num to the bits already present in
 432  | the array buf. You need to give it the number of bits to use and you
 433  | better make sure that this number of bits is enough to hold the value
 434  | Also num must be positive.
 435  |
 436 */
 437
 438 static void sendbits(int buf[], int num_of_bits, int num) {
 439
 440     unsigned int cnt, lastbyte;
 441     int lastbits;
 442     unsigned char * cbuf;
 443
 444     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 445     cnt = (unsigned int) buf[0];
 446     lastbits = buf[1];
 447     lastbyte =(unsigned int) buf[2];
 448     while (num_of_bits >= 8) {
 449         lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
 450         cbuf[cnt++] = lastbyte >> lastbits;
 451         num_of_bits -= 8;
 452     }
 453     if (num_of_bits > 0) {
 454         lastbyte = (lastbyte << num_of_bits) | num;
 455         lastbits += num_of_bits;
 456         if (lastbits >= 8) {
 457             lastbits -= 8;
 458             cbuf[cnt++] = lastbyte >> lastbits;
 459         }
 460     }
 461     buf[0] = cnt;
 462     buf[1] = lastbits;
 463     buf[2] = lastbyte;
 464     if (lastbits>0) {
 465         cbuf[cnt] = lastbyte << (8 - lastbits);
 466     }
 467 }
 468
 469 /*_________________________________________________________________________
 470  |
 471  | sizeofint - calculate bitsize of an integer
 472  |
 473  | return the number of bits needed to store an integer with given max size
 474  |
 475 */
 476
 477 static int sizeofint(const int size) {
 478     int num = 1;
 479     int num_of_bits = 0;
 480
 481     while (size >= num && num_of_bits < 32) {
 482         num_of_bits++;
 483         num <<= 1;
 484     }
 485     return num_of_bits;
 486 }
 487
 488 /*___________________________________________________________________________
 489  |
 490  | sizeofints - calculate 'bitsize' of compressed ints
 491  |
 492  | given the number of small unsigned integers and the maximum value
 493  | return the number of bits needed to read or write them with the
 494  | routines receiveints and sendints. You need this parameter when
 495  | calling these routines. Note that for many calls I can use
 496  | the variable 'smallidx' which is exactly the number of bits, and
 497  | So I don't need to call 'sizeofints for those calls.
 498 */
 499
 500 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
 501     int i, num;
 502         int bytes[32];
 503     unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
 504     num_of_bytes = 1;
 505     bytes[0] = 1;
 506     num_of_bits = 0;
 507     for (i=0; i < num_of_ints; i++) {
 508         tmp = 0;
 509         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 510             tmp = bytes[bytecnt] * sizes[i] + tmp;
 511             bytes[bytecnt] = tmp & 0xff;
 512             tmp >>= 8;
 513         }
 514         while (tmp != 0) {
 515             bytes[bytecnt++] = tmp & 0xff;
 516             tmp >>= 8;
 517         }
 518         num_of_bytes = bytecnt;
 519     }
 520     num = 1;
 521     num_of_bytes--;
 522     while (bytes[num_of_bytes] >= num) {
 523         num_of_bits++;
 524         num *= 2;
 525     }
 526     return num_of_bits + num_of_bytes * 8;
 527
 528 }
 529
 530 /*____________________________________________________________________________
 531  |
 532  | sendints - send a small set of small integers in compressed format
 533  |
 534  | this routine is used internally by xdr3dfcoord, to send a set of
 535  | small integers to the buffer.
 536  | Multiplication with fixed (specified maximum ) sizes is used to get
 537  | to one big, multibyte integer. Allthough the routine could be
 538  | modified to handle sizes bigger than 16777216, or more than just
 539  | a few integers, this is not done, because the gain in compression
 540  | isn't worth the effort. Note that overflowing the multiplication
 541  | or the byte buffer (32 bytes) is unchecked and causes bad results.
 542  |
 543  */
 544
 545 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
 546         unsigned int sizes[], unsigned int nums[]) {
 547
 548     int i, num_of_bytes, bytecnt;
 549     unsigned int bytes[32], tmp;
 550
 551     tmp = nums[0];
 552     num_of_bytes = 0;
 553     do {
 554         bytes[num_of_bytes++] = tmp & 0xff;
 555         tmp >>= 8;
 556     } while (tmp != 0);
 557
 558     for (i = 1; i < num_of_ints; i++) {
 559         if (nums[i] >= sizes[i]) {
 560             fprintf(stderr,"major breakdown in sendints num %u doesn't "
 561                     "match size %u\n", nums[i], sizes[i]);
 562             exit(1);
 563         }
 564         /* use one step multiply */
 565         tmp = nums[i];
 566         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 567             tmp = bytes[bytecnt] * sizes[i] + tmp;
 568             bytes[bytecnt] = tmp & 0xff;
 569             tmp >>= 8;
 570         }
 571         while (tmp != 0) {
 572             bytes[bytecnt++] = tmp & 0xff;
 573             tmp >>= 8;
 574         }
 575         num_of_bytes = bytecnt;
 576     }
 577     if (num_of_bits >= num_of_bytes * 8) {
 578         for (i = 0; i < num_of_bytes; i++) {
 579             sendbits(buf, 8, bytes[i]);
 580         }
 581         sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
 582     } else {
 583         for (i = 0; i < num_of_bytes-1; i++) {
 584             sendbits(buf, 8, bytes[i]);
 585         }
 586         sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
 587     }
 588 }
 589
 590
 591 /*___________________________________________________________________________
 592  |
 593  | receivebits - decode number from buf using specified number of bits
 594  |
 595  | extract the number of bits from the array buf and construct an integer
 596  | from it. Return that value.
 597  |
 598 */
 599
 600 static int receivebits(int buf[], int num_of_bits) {
 601
 602     int cnt, num, lastbits;
 603     unsigned int lastbyte;
 604     unsigned char * cbuf;
 605     int mask = (1 << num_of_bits) -1;
 606
 607     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 608     cnt = buf[0];
 609     lastbits = (unsigned int) buf[1];
 610     lastbyte = (unsigned int) buf[2];
 611
 612     num = 0;
 613     while (num_of_bits >= 8) {
 614         lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
 615         num |=  (lastbyte >> lastbits) << (num_of_bits - 8);
 616         num_of_bits -=8;
 617     }
 618     if (num_of_bits > 0) {
 619         if (lastbits < num_of_bits) {
 620             lastbits += 8;
 621             lastbyte = (lastbyte << 8) | cbuf[cnt++];
 622         }
 623         lastbits -= num_of_bits;
 624         num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
 625     }
 626     num &= mask;
 627     buf[0] = cnt;
 628     buf[1] = lastbits;
 629     buf[2] = lastbyte;
 630     return num;
 631 }
 632
 633 /*____________________________________________________________________________
 634  |
 635  | receiveints - decode 'small' integers from the buf array
 636  |
 637  | this routine is the inverse from sendints() and decodes the small integers
 638  | written to buf by calculating the remainder and doing divisions with
 639  | the given sizes[]. You need to specify the total number of bits to be
 640  | used from buf in num_of_bits.
 641  |
 642 */
 643
 644 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
 645         unsigned int sizes[], int nums[]) {
 646     int bytes[32];
 647     int i, j, num_of_bytes, p, num;
 648
 649     bytes[1] = bytes[2] = bytes[3] = 0;
 650     num_of_bytes = 0;
 651     while (num_of_bits > 8) {
 652         bytes[num_of_bytes++] = receivebits(buf, 8);
 653         num_of_bits -= 8;
 654     }
 655     if (num_of_bits > 0) {
 656         bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
 657     }
 658     for (i = num_of_ints-1; i > 0; i--) {
 659         num = 0;
 660         for (j = num_of_bytes-1; j >=0; j--) {
 661             num = (num << 8) | bytes[j];
 662             p = num / sizes[i];
 663             bytes[j] = p;
 664             num = num - p * sizes[i];
 665         }
 666         nums[i] = num;
 667     }
 668     nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
 669 }
 670
 671 /*____________________________________________________________________________
 672  |
 673  | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
 674  |
 675  | this routine reads or writes (depending on how you opened the file with
 676  | xdropen() ) a large number of 3d coordinates (stored in *fp).
 677  | The number of coordinates triplets to write is given by *size. On
 678  | read this number may be zero, in which case it reads as many as were written
 679  | or it may specify the number if triplets to read (which should match the
 680  | number written).
 681  | Compression is achieved by first converting all floating numbers to integer
 682  | using multiplication by *precision and rounding to the nearest integer.
 683  | Then the minimum and maximum value are calculated to determine the range.
 684  | The limited range of integers so found, is used to compress the coordinates.
 685  | In addition the differences between succesive coordinates is calculated.
 686  | If the difference happens to be 'small' then only the difference is saved,
 687  | compressing the data even more. The notion of 'small' is changed dynamically
 688  | and is enlarged or reduced whenever needed or possible.
 689  | Extra compression is achieved in the case of GROMOS and coordinates of
 690  | water molecules. GROMOS first writes out the Oxygen position, followed by
 691  | the two hydrogens. In order to make the differences smaller (and thereby
 692  | compression the data better) the order is changed into first one hydrogen
 693  | then the oxygen, followed by the other hydrogen. This is rather special, but
 694  | it shouldn't harm in the general case.
 695  |
 696  */
 697
 698 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
 699
 700
 701     static int *ip = NULL;
 702     static int oldsize;
 703     static int *buf;
 704
 705     int minint[3], maxint[3], mindiff, *lip, diff;
 706     int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
 707     int minidx, maxidx;
 708     unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
 709     int flag, k;
 710     int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
 711     float *lfp, lf;
 712     int tmp, *thiscoord,  prevcoord[3];
 713     unsigned int tmpcoord[30];
 714
 715     int bufsize, xdrid, lsize;
 716     unsigned int bitsize;
 717     float inv_precision;
 718     int errval = 1;
 719     int rc;
 720
 721     bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
 722     prevcoord[0]  = prevcoord[1]  = prevcoord[2]  = 0;
 723
 724     /* find out if xdrs is opened for reading or for writing */
 725     xdrid = 0;
 726     while (xdridptr[xdrid] != xdrs) {
 727         xdrid++;
 728         if (xdrid >= MAXID) {
 729             fprintf(stderr, "xdr error. no open xdr stream\n");
 730             exit (1);
 731         }
 732     }
 733     if ((xdrmodes[xdrid] == 'w') || (xdrmodes[xdrid] == 'a')) {
 734
 735         /* xdrs is open for writing */
 736
 737         if (xdr_int(xdrs, size) == 0)
 738             return 0;
 739         size3 = *size * 3;
 740         /* when the number of coordinates is small, don't try to compress; just
 741          * write them as floats using xdr_vector
 742          */
 743         if (*size <= 9 ) {
 744             return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
 745                     (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
 746         }
 747
 748         if(xdr_float(xdrs, precision) == 0)
 749                 return 0;
 750
 751         if (ip == NULL) {
 752             ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
 753             if (ip == NULL) {
 754                 fprintf(stderr,"malloc failed\n");
 755                 exit(1);
 756             }
 757             bufsize = size3 * 1.2;
 758             buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
 759             if (buf == NULL) {
 760                 fprintf(stderr,"malloc failed\n");
 761                 exit(1);
 762             }
 763             oldsize = *size;
 764         } else if (*size > oldsize) {
 765             ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
 766             if (ip == NULL) {
 767                 fprintf(stderr,"malloc failed\n");
 768                 exit(1);
 769             }
 770             bufsize = size3 * 1.2;
 771             buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
 772             if (buf == NULL) {
 773                 fprintf(stderr,"realloc failed\n");
 774                 exit(1);
 775             }
 776             oldsize = *size;
 777         }
 778         /* buf[0-2] are special and do not contain actual data */
 779         buf[0] = buf[1] = buf[2] = 0;
 780         minint[0] = minint[1] = minint[2] = INT_MAX;
 781         maxint[0] = maxint[1] = maxint[2] = INT_MIN;
 782         prevrun = -1;
 783         lfp = fp;
 784         lip = ip;
 785         mindiff = INT_MAX;
 786         oldlint1 = oldlint2 = oldlint3 = 0;
 787         while(lfp < fp + size3 ) {
 788             /* find nearest integer */
 789             if (*lfp >= 0.0)
 790                 lf = *lfp * *precision + 0.5;
 791             else
 792                 lf = *lfp * *precision - 0.5;
 793             if (fabs(lf) > MAXABS) {
 794                 /* scaling would cause overflow */
 795                 errval = 0;
 796             }
 797             lint1 = lf;
 798             if (lint1 < minint[0]) minint[0] = lint1;
 799             if (lint1 > maxint[0]) maxint[0] = lint1;
 800             *lip++ = lint1;
 801             lfp++;
 802             if (*lfp >= 0.0)
 803                 lf = *lfp * *precision + 0.5;
 804             else
 805                 lf = *lfp * *precision - 0.5;
 806             if (fabs(lf) > MAXABS) {
 807                 /* scaling would cause overflow */
 808                 errval = 0;
 809             }
 810             lint2 = lf;
 811             if (lint2 < minint[1]) minint[1] = lint2;
 812             if (lint2 > maxint[1]) maxint[1] = lint2;
 813             *lip++ = lint2;
 814             lfp++;
 815             if (*lfp >= 0.0)
 816                 lf = *lfp * *precision + 0.5;
 817             else
 818                 lf = *lfp * *precision - 0.5;
 819             if (fabs(lf) > MAXABS) {
 820                 /* scaling would cause overflow */
 821                 errval = 0;
 822             }
 823             lint3 = lf;
 824             if (lint3 < minint[2]) minint[2] = lint3;
 825             if (lint3 > maxint[2]) maxint[2] = lint3;
 826             *lip++ = lint3;
 827             lfp++;
 828             diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
 829             if (diff < mindiff && lfp > fp + 3)
 830                 mindiff = diff;
 831             oldlint1 = lint1;
 832             oldlint2 = lint2;
 833             oldlint3 = lint3;
 834         }
 835         if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
 836                  (xdr_int(xdrs, &(minint[1])) == 0) ||
 837                  (xdr_int(xdrs, &(minint[2])) == 0) ||
 838              (xdr_int(xdrs, &(maxint[0])) == 0) ||
 839                  (xdr_int(xdrs, &(maxint[1])) == 0) ||
 840                  (xdr_int(xdrs, &(maxint[2])) == 0))
 841         {
 842                 return 0;
 843         }
 844
 845         if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
 846                 (float)maxint[1] - (float)minint[1] >= MAXABS ||
 847                 (float)maxint[2] - (float)minint[2] >= MAXABS) {
 848             /* turning value in unsigned by subtracting minint
 849              * would cause overflow
 850              */
 851             errval = 0;
 852         }
 853         sizeint[0] = maxint[0] - minint[0]+1;
 854         sizeint[1] = maxint[1] - minint[1]+1;
 855         sizeint[2] = maxint[2] - minint[2]+1;
 856
 857         /* check if one of the sizes is to big to be multiplied */
 858         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
 859             bitsizeint[0] = sizeofint(sizeint[0]);
 860             bitsizeint[1] = sizeofint(sizeint[1]);
 861             bitsizeint[2] = sizeofint(sizeint[2]);
 862             bitsize = 0; /* flag the use of large sizes */
 863         } else {
 864             bitsize = sizeofints(3, sizeint);
 865         }
 866         lip = ip;
 867         luip = (unsigned int *) ip;
 868         smallidx = FIRSTIDX;
 869         while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
 870             smallidx++;
 871         }
 872         if(xdr_int(xdrs, &smallidx) == 0)
 873                 return 0;
 874
 875         maxidx = MIN(LASTIDX, smallidx + 8) ;
 876         minidx = maxidx - 8; /* often this equal smallidx */
 877         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
 878         smallnum = magicints[smallidx] / 2;
 879         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 880         larger = magicints[maxidx] / 2;
 881         i = 0;
 882         while (i < *size) {
 883             is_small = 0;
 884             thiscoord = (int *)(luip) + i * 3;
 885             if (smallidx < maxidx && i >= 1 &&
 886                     abs(thiscoord[0] - prevcoord[0]) < larger &&
 887                     abs(thiscoord[1] - prevcoord[1]) < larger &&
 888                     abs(thiscoord[2] - prevcoord[2]) < larger) {
 889                 is_smaller = 1;
 890             } else if (smallidx > minidx) {
 891                 is_smaller = -1;
 892             } else {
 893                 is_smaller = 0;
 894             }
 895             if (i + 1 < *size) {
 896                 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
 897                         abs(thiscoord[1] - thiscoord[4]) < smallnum &&
 898                         abs(thiscoord[2] - thiscoord[5]) < smallnum) {
 899                     /* interchange first with second atom for better
 900                      * compression of water molecules
 901                      */
 902                     tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
 903                         thiscoord[3] = tmp;
 904                     tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
 905                         thiscoord[4] = tmp;
 906                     tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
 907                         thiscoord[5] = tmp;
 908                     is_small = 1;
 909                 }
 910
 911             }
 912             tmpcoord[0] = thiscoord[0] - minint[0];
 913             tmpcoord[1] = thiscoord[1] - minint[1];
 914             tmpcoord[2] = thiscoord[2] - minint[2];
 915             if (bitsize == 0) {
 916                 sendbits(buf, bitsizeint[0], tmpcoord[0]);
 917                 sendbits(buf, bitsizeint[1], tmpcoord[1]);
 918                 sendbits(buf, bitsizeint[2], tmpcoord[2]);
 919             } else {
 920                 sendints(buf, 3, bitsize, sizeint, tmpcoord);
 921             }
 922             prevcoord[0] = thiscoord[0];
 923             prevcoord[1] = thiscoord[1];
 924             prevcoord[2] = thiscoord[2];
 925             thiscoord = thiscoord + 3;
 926             i++;
 927
 928             run = 0;
 929             if (is_small == 0 && is_smaller == -1)
 930                 is_smaller = 0;
 931             while (is_small && run < 8*3) {
 932                 if (is_smaller == -1 && (
 933                         SQR(thiscoord[0] - prevcoord[0]) +
 934                         SQR(thiscoord[1] - prevcoord[1]) +
 935                         SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
 936                     is_smaller = 0;
 937                 }
 938
 939                 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
 940                 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
 941                 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
 942
 943                 prevcoord[0] = thiscoord[0];
 944                 prevcoord[1] = thiscoord[1];
 945                 prevcoord[2] = thiscoord[2];
 946
 947                 i++;
 948                 thiscoord = thiscoord + 3;
 949                 is_small = 0;
 950                 if (i < *size &&
 951                         abs(thiscoord[0] - prevcoord[0]) < smallnum &&
 952                         abs(thiscoord[1] - prevcoord[1]) < smallnum &&
 953                         abs(thiscoord[2] - prevcoord[2]) < smallnum) {
 954                     is_small = 1;
 955                 }
 956             }
 957             if (run != prevrun || is_smaller != 0) {
 958                 prevrun = run;
 959                 sendbits(buf, 1, 1); /* flag the change in run-length */
 960                 sendbits(buf, 5, run+is_smaller+1);
 961             } else {
 962                 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
 963             }
 964             for (k=0; k < run; k+=3) {
 965                 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
 966             }
 967             if (is_smaller != 0) {
 968                 smallidx += is_smaller;
 969                 if (is_smaller < 0) {
 970                     smallnum = smaller;
 971                     smaller = magicints[smallidx-1] / 2;
 972                 } else {
 973                     smaller = smallnum;
 974                     smallnum = magicints[smallidx] / 2;
 975                 }
 976                 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 977             }
 978         }
 979         if (buf[1] != 0) buf[0]++;
 980                 /* buf[0] holds the length in bytes */
 981         if(xdr_int(xdrs, &(buf[0])) == 0)
 982                 return 0;
 983
 984         return errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
 985     } else {
 986
 987         /* xdrs is open for reading */
 988
 989         if (xdr_int(xdrs, &lsize) == 0)
 990             return 0;
 991         if (*size != 0 && lsize != *size) {
 992             fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
 993                     "%d arg vs %d in file", *size, lsize);
 994         }
 995         *size = lsize;
 996         size3 = *size * 3;
 997         if (*size <= 9) {
 998             *precision = -1;
 999             return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
1000                     (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
1001         }
1002         if(xdr_float(xdrs, precision) == 0)
1003                 return 0;
1004
1005         if (ip == NULL) {
1006             ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
1007             if (ip == NULL) {
1008                 fprintf(stderr,"malloc failed\n");
1009                 exit(1);
1010             }
1011             bufsize = size3 * 1.2;
1012             buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
1013             if (buf == NULL) {
1014                 fprintf(stderr,"malloc failed\n");
1015                 exit(1);
1016             }
1017             oldsize = *size;
1018         } else if (*size > oldsize) {
1019             ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
1020             if (ip == NULL) {
1021                 fprintf(stderr,"malloc failed\n");
1022                 exit(1);
1023             }
1024             bufsize = size3 * 1.2;
1025             buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
1026             if (buf == NULL) {
1027                 fprintf(stderr,"malloc failed\n");
1028                 exit(1);
1029             }
1030             oldsize = *size;
1031         }
1032         buf[0] = buf[1] = buf[2] = 0;
1033
1034         if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
1035                  (xdr_int(xdrs, &(minint[1])) == 0) ||
1036                  (xdr_int(xdrs, &(minint[2])) == 0) ||
1037                  (xdr_int(xdrs, &(maxint[0])) == 0) ||
1038                  (xdr_int(xdrs, &(maxint[1])) == 0) ||
1039                  (xdr_int(xdrs, &(maxint[2])) == 0))
1040         {
1041                 return 0;
1042         }
1043
1044         sizeint[0] = maxint[0] - minint[0]+1;
1045         sizeint[1] = maxint[1] - minint[1]+1;
1046         sizeint[2] = maxint[2] - minint[2]+1;
1047
1048         /* check if one of the sizes is to big to be multiplied */
1049         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1050             bitsizeint[0] = sizeofint(sizeint[0]);
1051             bitsizeint[1] = sizeofint(sizeint[1]);
1052             bitsizeint[2] = sizeofint(sizeint[2]);
1053             bitsize = 0; /* flag the use of large sizes */
1054         } else {
1055             bitsize = sizeofints(3, sizeint);
1056         }
1057
1058         if (xdr_int(xdrs, &smallidx) == 0)
1059             return 0;
1060         maxidx = MIN(LASTIDX, smallidx + 8) ;
1061         minidx = maxidx - 8; /* often this equal smallidx */
1062         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1063         smallnum = magicints[smallidx] / 2;
1064         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1065         larger = magicints[maxidx];
1066
1067         /* buf[0] holds the length in bytes */
1068
1069         if (xdr_int(xdrs, &(buf[0])) == 0)
1070             return 0;
1071         if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
1072             return 0;
1073         buf[0] = buf[1] = buf[2] = 0;
1074
1075         lfp = fp;
1076         inv_precision = 1.0 / * precision;
1077         run = 0;
1078         i = 0;
1079         lip = ip;
1080         while ( i < lsize ) {
1081             thiscoord = (int *)(lip) + i * 3;
1082
1083             if (bitsize == 0) {
1084                 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1085                 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1086                 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1087             } else {
1088                 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1089             }
1090
1091             i++;
1092             thiscoord[0] += minint[0];
1093             thiscoord[1] += minint[1];
1094             thiscoord[2] += minint[2];
1095
1096             prevcoord[0] = thiscoord[0];
1097             prevcoord[1] = thiscoord[1];
1098             prevcoord[2] = thiscoord[2];
1099
1100
1101             flag = receivebits(buf, 1);
1102             is_smaller = 0;
1103             if (flag == 1) {
1104                 run = receivebits(buf, 5);
1105                 is_smaller = run % 3;
1106                 run -= is_smaller;
1107                 is_smaller--;
1108             }
1109             if (run > 0) {
1110                 thiscoord += 3;
1111                 for (k = 0; k < run; k+=3) {
1112                     receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1113                     i++;
1114                     thiscoord[0] += prevcoord[0] - smallnum;
1115                     thiscoord[1] += prevcoord[1] - smallnum;
1116                     thiscoord[2] += prevcoord[2] - smallnum;
1117                     if (k == 0) {
1118                         /* interchange first with second atom for better
1119                          * compression of water molecules
1120                          */
1121                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1122                                 prevcoord[0] = tmp;
1123                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1124                                 prevcoord[1] = tmp;
1125                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1126                                 prevcoord[2] = tmp;
1127                         *lfp++ = prevcoord[0] * inv_precision;
1128                         *lfp++ = prevcoord[1] * inv_precision;
1129                         *lfp++ = prevcoord[2] * inv_precision;
1130                     } else {
1131                         prevcoord[0] = thiscoord[0];
1132                         prevcoord[1] = thiscoord[1];
1133                         prevcoord[2] = thiscoord[2];
1134                     }
1135                     *lfp++ = thiscoord[0] * inv_precision;
1136                     *lfp++ = thiscoord[1] * inv_precision;
1137                     *lfp++ = thiscoord[2] * inv_precision;
1138                 }
1139             } else {
1140                 *lfp++ = thiscoord[0] * inv_precision;
1141                 *lfp++ = thiscoord[1] * inv_precision;
1142                 *lfp++ = thiscoord[2] * inv_precision;
1143             }
1144             smallidx += is_smaller;
1145             if (is_smaller < 0) {
1146                 smallnum = smaller;
1147                 if (smallidx > FIRSTIDX) {
1148                     smaller = magicints[smallidx - 1] /2;
1149                 } else {
1150                     smaller = 0;
1151                 }
1152             } else if (is_smaller > 0) {
1153                 smaller = smallnum;
1154                 smallnum = magicints[smallidx] / 2;
1155             }
1156             sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1157         }
1158     }
1159     return 1;
1160 }
1161
1162
1163
1164 /******************************************************************
1165
1166   XTC files have a relatively simple structure.
1167   They have a header of 16 bytes and the rest are
1168   the compressed coordinates of the files. Due to the
1169   compression 00 is not present in the coordinates.
1170   The first 4 bytes of the header are the magic number
1171   1995 (0x000007CB). If we find this number we are guaranteed
1172   to be in the header, due to the presence of so many zeros.
1173   The second 4 bytes are the number of atoms in the frame, and is
1174   assumed to be constant. The third 4 bytes are the frame number.
1175   The last 4 bytes are a floating point representation of the time.
1176
1177 ********************************************************************/
1178
1179 /* Must match definition in xtcio.c */
1180 #ifndef XTC_MAGIC
1181 #define XTC_MAGIC 1995
1182 #endif
1183
1184 static const int header_size = 16;
1185
1186 /* Check if we are at the header start.
1187    At the same time it will also read 1 int */
1188 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1189                                int natoms, int * timestep, float * time){
1190   int i_inp[3];
1191   float f_inp[10];
1192   int i;
1193   gmx_off_t off;
1194
1195
1196   if((off = gmx_ftell(fp)) < 0){
1197     return -1;
1198   }
1199   /* read magic natoms and timestep */
1200   for(i = 0;i<3;i++){
1201     if(!xdr_int(xdrs, &(i_inp[i]))){
1202       gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET);
1203       return -1;
1204     }
1205   }
1206   /* quick return */
1207   if(i_inp[0] != XTC_MAGIC){
1208     if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
1209       return -1;
1210     }
1211     return 0;
1212   }
1213   /* read time and box */
1214   for(i = 0;i<10;i++){
1215     if(!xdr_float(xdrs, &(f_inp[i]))){
1216       gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET);
1217       return -1;
1218     }
1219   }
1220   /* Make a rigourous check to see if we are in the beggining of a header
1221      Hopefully there are no ambiguous cases */
1222   /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1223      right triangle and that the first element must be nonzero unless the entire matrix is zero
1224   */
1225   if(i_inp[1] == natoms &&
1226      ((f_inp[1] != 0 && f_inp[6] == 0) ||
1227       (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0))){
1228     if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
1229       return -1;
1230     }
1231     *time = f_inp[0];
1232     *timestep = i_inp[2];
1233     return 1;
1234   }
1235   if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
1236     return -1;
1237   }
1238   return 0;
1239 }
1240
1241 static int
1242 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1243 {
1244     gmx_off_t off;
1245     int step;
1246     float time;
1247     int ret;
1248
1249     if((off = gmx_ftell(fp)) < 0){
1250       return -1;
1251     }
1252
1253     /* read one int just to make sure we dont read this frame but the next */
1254     xdr_int(xdrs,&step);
1255     while(1){
1256       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
1257       if(ret == 1){
1258         if(gmx_fseek(fp,off,SEEK_SET)){
1259           return -1;
1260         }
1261         return step;
1262       }else if(ret == -1){
1263         if(gmx_fseek(fp,off,SEEK_SET)){
1264           return -1;
1265         }
1266       }
1267     }
1268     return -1;
1269 }
1270
1271
1272 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1273                                      bool * bOK)
1274 {
1275     gmx_off_t off;
1276     float time;
1277     int step;
1278     int ret;
1279     *bOK = 0;
1280
1281     if ((off = gmx_ftell(fp)) < 0)
1282     {
1283         return -1;
1284     }
1285     /* read one int just to make sure we dont read this frame but the next */
1286     xdr_int(xdrs, &step);
1287     while (1)
1288     {
1289         ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1290         if (ret == 1)
1291         {
1292             *bOK = 1;
1293             if (gmx_fseek(fp,off,SEEK_SET))
1294             {
1295                 *bOK = 0;
1296                 return -1;
1297             }
1298             return time;
1299         }
1300         else if (ret == -1)
1301         {
1302             if (gmx_fseek(fp,off,SEEK_SET))
1303             {
1304                 return -1;
1305             }
1306             return -1;
1307         }
1308     }
1309     return -1;
1310 }
1311
1312
1313 static float
1314 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, bool * bOK)
1315 {
1316     gmx_off_t off;
1317     int step;
1318     float time;
1319     int ret;
1320     *bOK = 0;
1321
1322     if ((off = gmx_ftell(fp)) < 0)
1323     {
1324         return -1;
1325     }
1326
1327     while (1)
1328     {
1329         ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1330         if (ret == 1)
1331         {
1332             *bOK = 1;
1333             if (gmx_fseek(fp,off,SEEK_SET))
1334             {
1335                 *bOK = 0;
1336                 return -1;
1337             }
1338             return time;
1339         }
1340         else if (ret == -1)
1341         {
1342             if (gmx_fseek(fp,off,SEEK_SET))
1343             {
1344                 return -1;
1345             }
1346             return -1;
1347         }
1348         else if (ret == 0)
1349         {
1350             /*Go back.*/
1351             if (gmx_fseek(fp,-2*XDR_INT_SIZE,SEEK_CUR))
1352             {
1353                 return -1;
1354             }
1355         }
1356     }
1357     return -1;
1358 }
1359
1360 /* Currently not used, just for completeness */
1361 static int
1362 xtc_get_current_frame_number(FILE *fp,XDR *xdrs,int natoms, bool * bOK)
1363 {
1364     gmx_off_t off;
1365     int ret;
1366     int step;
1367     float time;
1368     *bOK = 0;
1369
1370     if((off = gmx_ftell(fp)) < 0){
1371       return -1;
1372     }
1373
1374
1375     while(1){
1376       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
1377       if(ret == 1){
1378         *bOK = 1;
1379         if(gmx_fseek(fp,off,SEEK_SET)){
1380                 *bOK = 0;
1381           return -1;
1382         }
1383         return step;
1384       }else if(ret == -1){
1385         if(gmx_fseek(fp,off,SEEK_SET)){
1386           return -1;
1387         }
1388         return -1;
1389
1390       }else if(ret == 0){
1391                   /*Go back.*/
1392                   if(gmx_fseek(fp,-2*XDR_INT_SIZE,SEEK_CUR)){
1393                           return -1;
1394                   }
1395       }
1396     }
1397     return -1;
1398 }
1399
1400
1401
1402 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1403 {
1404     int inp;
1405     gmx_off_t res;
1406     int ret;
1407     int step;
1408     float time;
1409     /* read one int just to make sure we dont read this frame but the next */
1410     xdr_int(xdrs,&step);
1411     while(1)
1412     {
1413       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
1414       if(ret == 1){
1415         if((res = gmx_ftell(fp)) >= 0){
1416           return res - XDR_INT_SIZE;
1417         }else{
1418           return res;
1419         }
1420       }else if(ret == -1){
1421         return -1;
1422       }
1423     }
1424     return -1;
1425 }
1426
1427
1428
1429 float
1430 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, bool * bOK)
1431 {
1432   float  res;
1433   float  tinit;
1434   gmx_off_t off;
1435
1436   if((off   = gmx_ftell(fp)) < 0){
1437     return -1;
1438   }
1439
1440     tinit = xtc_get_current_frame_time(fp,xdrs,natoms,bOK);
1441
1442     *bOK = 1;
1443
1444     if(!(*bOK))
1445     {
1446         return -1;
1447     }
1448
1449     res = xtc_get_next_frame_time(fp,xdrs,natoms,bOK);
1450
1451     if(!(*bOK))
1452     {
1453         return -1;
1454     }
1455
1456     res -= tinit;
1457     if(gmx_fseek(fp,off,SEEK_SET)){
1458       *bOK = 0;
1459       return -1;
1460     }
1461     return res;
1462 }
1463
1464
1465 int
1466 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1467 {
1468     gmx_off_t low = 0;
1469     gmx_off_t high,pos;
1470
1471
1472     /* round to 4 bytes */
1473     int fr;
1474     gmx_off_t  offset;
1475     if(gmx_fseek(fp,0,SEEK_END)){
1476       return -1;
1477     }
1478
1479     if((high = gmx_ftell(fp)) < 0){
1480       return -1;
1481     }
1482
1483     /* round to 4 bytes  */
1484     high /= XDR_INT_SIZE;
1485     high *= XDR_INT_SIZE;
1486     offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1487
1488     if(gmx_fseek(fp,offset,SEEK_SET)){
1489       return -1;
1490     }
1491
1492     while(1)
1493     {
1494         fr = xtc_get_next_frame_number(fp,xdrs,natoms);
1495         if(fr < 0)
1496         {
1497             return -1;
1498         }
1499         if(fr != frame && abs(low-high) > header_size)
1500         {
1501             if(fr < frame)
1502             {
1503                 low = offset;
1504             }
1505             else
1506             {
1507                 high = offset;
1508             }
1509             /* round to 4 bytes */
1510             offset = (((high+low)/2)/4)*4;
1511
1512             if(gmx_fseek(fp,offset,SEEK_SET)){
1513               return -1;
1514             }
1515         }
1516         else
1517         {
1518             break;
1519         }
1520     }
1521     if(offset <= header_size)
1522     {
1523         offset = low;
1524     }
1525
1526     if(gmx_fseek(fp,offset,SEEK_SET)){
1527       return -1;
1528     }
1529
1530     if((pos = xtc_get_next_frame_start(fp,xdrs,natoms))< 0){
1531     /* we probably hit an end of file */
1532       return -1;
1533     }
1534
1535     if(gmx_fseek(fp,pos,SEEK_SET)){
1536       return -1;
1537     }
1538
1539     return 0;
1540 }
1541
1542
1543
1544 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms)
1545 {
1546     float t;
1547     float dt;
1548     bool bOK;
1549     gmx_off_t low = 0;
1550     gmx_off_t high, offset, pos;
1551     int res;
1552     int dt_sign = 0;
1553
1554     if (gmx_fseek(fp,0,SEEK_END))
1555     {
1556         return -1;
1557     }
1558
1559     if ((high = gmx_ftell(fp)) < 0)
1560     {
1561         return -1;
1562     }
1563     /* round to int  */
1564     high /= XDR_INT_SIZE;
1565     high *= XDR_INT_SIZE;
1566     offset = ((high / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1567
1568     if (gmx_fseek(fp,offset,SEEK_SET))
1569     {
1570         return -1;
1571     }
1572
1573
1574     /*
1575      * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1576     dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1577
1578     if (!bOK)
1579     {
1580         return -1;
1581     }
1582     */
1583
1584     while (1)
1585     {
1586         dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1587         if (!bOK)
1588         {
1589             return -1;
1590         }
1591         else
1592         {
1593             if (dt > 0)
1594             {
1595                 if (dt_sign == -1)
1596                 {
1597                     /* Found a place in the trajectory that has positive time step while
1598                      other has negative time step */
1599                     return -2;
1600                 }
1601                 dt_sign = 1;
1602             }
1603             else if (dt < 0)
1604             {
1605                 if (dt_sign == 1)
1606                 {
1607                     /* Found a place in the trajectory that has positive time step while
1608                      other has negative time step */
1609                     return -2;
1610                 }
1611                 dt_sign = -1;
1612             }
1613         }
1614         t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1615         if (!bOK)
1616         {
1617             return -1;
1618         }
1619
1620         /* If we are before the target time and the time step is positive or 0, or we have
1621          after the target time and the time step is negative, or the difference between
1622          the current time and the target time is bigger than dt and above all the distance between high
1623          and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1624          if we reached the solution */
1625         if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
1626             - time) >= dt && dt_sign >= 0)
1627             || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
1628             > header_size))
1629         {
1630             if (dt >= 0 && dt_sign != -1)
1631             {
1632                 if (t < time)
1633                 {
1634                     low = offset;
1635                 }
1636                 else
1637                 {
1638                     high = offset;
1639                 }
1640             }
1641             else if (dt <= 0 && dt_sign == -1)
1642             {
1643                 if (t >= time)
1644                 {
1645                     low = offset;
1646                 }
1647                 else
1648                 {
1649                     high = offset;
1650                 }
1651             }
1652             else
1653             {
1654                 /* We should never reach here */
1655                 return -1;
1656             }
1657             /* round to 4 bytes and subtract header*/
1658             offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1659             if (gmx_fseek(fp,offset,SEEK_SET))
1660             {
1661                 return -1;
1662             }
1663         }
1664         else
1665         {
1666             if (abs(low - high) <= header_size)
1667             {
1668                 break;
1669             }
1670             /* re-estimate dt */
1671             if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1672             {
1673                 if (bOK)
1674                 {
1675                     dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1676                 }
1677             }
1678             if (t >= time && t - time < dt)
1679             {
1680                 break;
1681             }
1682         }
1683     }
1684
1685     if (offset <= header_size)
1686     {
1687         offset = low;
1688     }
1689
1690     gmx_fseek(fp,offset,SEEK_SET);
1691
1692     if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1693     {
1694         return -1;
1695     }
1696
1697     if (gmx_fseek(fp,pos,SEEK_SET))
1698     {
1699         return -1;
1700     }
1701     return 0;
1702 }
1703
1704 float
1705 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, bool * bOK)
1706 {
1707     float  time;
1708     gmx_off_t  off;
1709     int res;
1710     *bOK = 1;
1711     off = gmx_ftell(fp);
1712     if(off < 0){
1713       *bOK = 0;
1714       return -1;
1715     }
1716
1717     if( (res = gmx_fseek(fp,-3*XDR_INT_SIZE,SEEK_END)) != 0){
1718       *bOK = 0;
1719       return -1;
1720     }
1721
1722     time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1723     if(!(*bOK)){
1724       return -1;
1725     }
1726
1727     if( (res = gmx_fseek(fp,off,SEEK_SET)) != 0){
1728       *bOK = 0;
1729       return -1;
1730     }
1731     return time;
1732 }
1733
1734
1735 int
1736 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, bool * bOK)
1737 {
1738     int    frame;
1739     gmx_off_t  off;
1740     int res;
1741     *bOK = 1;
1742
1743     if((off = gmx_ftell(fp)) < 0){
1744       *bOK = 0;
1745       return -1;
1746     }
1747
1748
1749     if(gmx_fseek(fp,-3*XDR_INT_SIZE,SEEK_END)){
1750       *bOK = 0;
1751       return -1;
1752     }
1753
1754     frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1755     if(!bOK){
1756       return -1;
1757     }
1758
1759
1760     if(gmx_fseek(fp,off,SEEK_SET)){
1761       *bOK = 0;
1762       return -1;
1763     }
1764
1765     return frame;
1766 }