src/gromacs/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 #include "string2.h"
  47 #include "futil.h"
  48 #include "gmx_fatal.h"
  49
  50
  51 #if 0
  52 #ifdef HAVE_FSEEKO
  53 #  define gmx_fseek(A,B,C) fseeko(A,B,C)
  54 #  define gmx_ftell(A) ftello(A)
  55 #  define gmx_off_t off_t
  56 #else
  57 #  define gmx_fseek(A,B,C) fseek(A,B,C)
  58 #  define gmx_ftell(A) ftell(A)
  59 #  define gmx_off_t int
  60 #endif
  61 #endif
  62
  63
  64 /* This is just for clarity - it can never be anything but 4! */
  65 #define XDR_INT_SIZE 4
  66
  67 /* same order as the definition of xdr_datatype */
  68 const char *xdr_datatype_names[] =
  69 {
  70     "int",
  71     "float",
  72     "double",
  73     "large int",
  74     "char",
  75     "string"
  76 };
  77
  78
  79 /*___________________________________________________________________________
  80  |
  81  | what follows are the C routine to read/write compressed coordinates together
  82  | with some routines to assist in this task (those are marked
  83  | static and cannot be called from user programs)
  84 */
  85 #define MAXABS INT_MAX-2
  86
  87 #ifndef MIN
  88 #define MIN(x,y) ((x) < (y) ? (x):(y))
  89 #endif
  90 #ifndef MAX
  91 #define MAX(x,y) ((x) > (y) ? (x):(y))
  92 #endif
  93 #ifndef SQR
  94 #define SQR(x) ((x)*(x))
  95 #endif
  96 static const int magicints[] = {
  97     0, 0, 0, 0, 0, 0, 0, 0, 0,
  98     8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
  99     80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
 100     812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
 101     8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
 102     82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
 103     832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
 104     8388607, 10568983, 13316085, 16777216 };
 105
 106 #define FIRSTIDX 9
 107 /* note that magicints[FIRSTIDX-1] == 0 */
 108 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
 109
 110
 111 /*____________________________________________________________________________
 112  |
 113  | sendbits - encode num into buf using the specified number of bits
 114  |
 115  | This routines appends the value of num to the bits already present in
 116  | the array buf. You need to give it the number of bits to use and you
 117  | better make sure that this number of bits is enough to hold the value
 118  | Also num must be positive.
 119  |
 120 */
 121
 122 static void sendbits(int buf[], int num_of_bits, int num) {
 123
 124     unsigned int cnt, lastbyte;
 125     int lastbits;
 126     unsigned char * cbuf;
 127
 128     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 129     cnt = (unsigned int) buf[0];
 130     lastbits = buf[1];
 131     lastbyte =(unsigned int) buf[2];
 132     while (num_of_bits >= 8) {
 133         lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
 134         cbuf[cnt++] = lastbyte >> lastbits;
 135         num_of_bits -= 8;
 136     }
 137     if (num_of_bits > 0) {
 138         lastbyte = (lastbyte << num_of_bits) | num;
 139         lastbits += num_of_bits;
 140         if (lastbits >= 8) {
 141             lastbits -= 8;
 142             cbuf[cnt++] = lastbyte >> lastbits;
 143         }
 144     }
 145     buf[0] = cnt;
 146     buf[1] = lastbits;
 147     buf[2] = lastbyte;
 148     if (lastbits>0) {
 149         cbuf[cnt] = lastbyte << (8 - lastbits);
 150     }
 151 }
 152
 153 /*_________________________________________________________________________
 154  |
 155  | sizeofint - calculate bitsize of an integer
 156  |
 157  | return the number of bits needed to store an integer with given max size
 158  |
 159 */
 160
 161 static int sizeofint(const int size) {
 162     int num = 1;
 163     int num_of_bits = 0;
 164
 165     while (size >= num && num_of_bits < 32) {
 166         num_of_bits++;
 167         num <<= 1;
 168     }
 169     return num_of_bits;
 170 }
 171
 172 /*___________________________________________________________________________
 173  |
 174  | sizeofints - calculate 'bitsize' of compressed ints
 175  |
 176  | given the number of small unsigned integers and the maximum value
 177  | return the number of bits needed to read or write them with the
 178  | routines receiveints and sendints. You need this parameter when
 179  | calling these routines. Note that for many calls I can use
 180  | the variable 'smallidx' which is exactly the number of bits, and
 181  | So I don't need to call 'sizeofints for those calls.
 182 */
 183
 184 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
 185     int i, num;
 186         int bytes[32];
 187     unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
 188     num_of_bytes = 1;
 189     bytes[0] = 1;
 190     num_of_bits = 0;
 191     for (i=0; i < num_of_ints; i++) {
 192         tmp = 0;
 193         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 194             tmp = bytes[bytecnt] * sizes[i] + tmp;
 195             bytes[bytecnt] = tmp & 0xff;
 196             tmp >>= 8;
 197         }
 198         while (tmp != 0) {
 199             bytes[bytecnt++] = tmp & 0xff;
 200             tmp >>= 8;
 201         }
 202         num_of_bytes = bytecnt;
 203     }
 204     num = 1;
 205     num_of_bytes--;
 206     while (bytes[num_of_bytes] >= num) {
 207         num_of_bits++;
 208         num *= 2;
 209     }
 210     return num_of_bits + num_of_bytes * 8;
 211
 212 }
 213
 214 /*____________________________________________________________________________
 215  |
 216  | sendints - send a small set of small integers in compressed format
 217  |
 218  | this routine is used internally by xdr3dfcoord, to send a set of
 219  | small integers to the buffer.
 220  | Multiplication with fixed (specified maximum ) sizes is used to get
 221  | to one big, multibyte integer. Allthough the routine could be
 222  | modified to handle sizes bigger than 16777216, or more than just
 223  | a few integers, this is not done, because the gain in compression
 224  | isn't worth the effort. Note that overflowing the multiplication
 225  | or the byte buffer (32 bytes) is unchecked and causes bad results.
 226  |
 227  */
 228
 229 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
 230         unsigned int sizes[], unsigned int nums[]) {
 231
 232     int i, num_of_bytes, bytecnt;
 233     unsigned int bytes[32], tmp;
 234
 235     tmp = nums[0];
 236     num_of_bytes = 0;
 237     do {
 238         bytes[num_of_bytes++] = tmp & 0xff;
 239         tmp >>= 8;
 240     } while (tmp != 0);
 241
 242     for (i = 1; i < num_of_ints; i++) {
 243         if (nums[i] >= sizes[i]) {
 244             fprintf(stderr,"major breakdown in sendints num %u doesn't "
 245                     "match size %u\n", nums[i], sizes[i]);
 246             exit(1);
 247         }
 248         /* use one step multiply */
 249         tmp = nums[i];
 250         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
 251             tmp = bytes[bytecnt] * sizes[i] + tmp;
 252             bytes[bytecnt] = tmp & 0xff;
 253             tmp >>= 8;
 254         }
 255         while (tmp != 0) {
 256             bytes[bytecnt++] = tmp & 0xff;
 257             tmp >>= 8;
 258         }
 259         num_of_bytes = bytecnt;
 260     }
 261     if (num_of_bits >= num_of_bytes * 8) {
 262         for (i = 0; i < num_of_bytes; i++) {
 263             sendbits(buf, 8, bytes[i]);
 264         }
 265         sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
 266     } else {
 267         for (i = 0; i < num_of_bytes-1; i++) {
 268             sendbits(buf, 8, bytes[i]);
 269         }
 270         sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
 271     }
 272 }
 273
 274
 275 /*___________________________________________________________________________
 276  |
 277  | receivebits - decode number from buf using specified number of bits
 278  |
 279  | extract the number of bits from the array buf and construct an integer
 280  | from it. Return that value.
 281  |
 282 */
 283
 284 static int receivebits(int buf[], int num_of_bits) {
 285
 286     int cnt, num, lastbits;
 287     unsigned int lastbyte;
 288     unsigned char * cbuf;
 289     int mask = (1 << num_of_bits) -1;
 290
 291     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
 292     cnt = buf[0];
 293     lastbits = (unsigned int) buf[1];
 294     lastbyte = (unsigned int) buf[2];
 295
 296     num = 0;
 297     while (num_of_bits >= 8) {
 298         lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
 299         num |=  (lastbyte >> lastbits) << (num_of_bits - 8);
 300         num_of_bits -=8;
 301     }
 302     if (num_of_bits > 0) {
 303         if (lastbits < num_of_bits) {
 304             lastbits += 8;
 305             lastbyte = (lastbyte << 8) | cbuf[cnt++];
 306         }
 307         lastbits -= num_of_bits;
 308         num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
 309     }
 310     num &= mask;
 311     buf[0] = cnt;
 312     buf[1] = lastbits;
 313     buf[2] = lastbyte;
 314     return num;
 315 }
 316
 317 /*____________________________________________________________________________
 318  |
 319  | receiveints - decode 'small' integers from the buf array
 320  |
 321  | this routine is the inverse from sendints() and decodes the small integers
 322  | written to buf by calculating the remainder and doing divisions with
 323  | the given sizes[]. You need to specify the total number of bits to be
 324  | used from buf in num_of_bits.
 325  |
 326 */
 327
 328 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
 329         unsigned int sizes[], int nums[]) {
 330     int bytes[32];
 331     int i, j, num_of_bytes, p, num;
 332
 333     bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
 334     num_of_bytes = 0;
 335     while (num_of_bits > 8) {
 336         bytes[num_of_bytes++] = receivebits(buf, 8);
 337         num_of_bits -= 8;
 338     }
 339     if (num_of_bits > 0) {
 340         bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
 341     }
 342     for (i = num_of_ints-1; i > 0; i--) {
 343         num = 0;
 344         for (j = num_of_bytes-1; j >=0; j--) {
 345             num = (num << 8) | bytes[j];
 346             p = num / sizes[i];
 347             bytes[j] = p;
 348             num = num - p * sizes[i];
 349         }
 350         nums[i] = num;
 351     }
 352     nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
 353 }
 354
 355 /*____________________________________________________________________________
 356  |
 357  | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
 358  |
 359  | this routine reads or writes (depending on how you opened the file with
 360  | xdropen() ) a large number of 3d coordinates (stored in *fp).
 361  | The number of coordinates triplets to write is given by *size. On
 362  | read this number may be zero, in which case it reads as many as were written
 363  | or it may specify the number if triplets to read (which should match the
 364  | number written).
 365  | Compression is achieved by first converting all floating numbers to integer
 366  | using multiplication by *precision and rounding to the nearest integer.
 367  | Then the minimum and maximum value are calculated to determine the range.
 368  | The limited range of integers so found, is used to compress the coordinates.
 369  | In addition the differences between succesive coordinates is calculated.
 370  | If the difference happens to be 'small' then only the difference is saved,
 371  | compressing the data even more. The notion of 'small' is changed dynamically
 372  | and is enlarged or reduced whenever needed or possible.
 373  | Extra compression is achieved in the case of GROMOS and coordinates of
 374  | water molecules. GROMOS first writes out the Oxygen position, followed by
 375  | the two hydrogens. In order to make the differences smaller (and thereby
 376  | compression the data better) the order is changed into first one hydrogen
 377  | then the oxygen, followed by the other hydrogen. This is rather special, but
 378  | it shouldn't harm in the general case.
 379  |
 380  */
 381
 382 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
 383 {
 384     int *ip = NULL;
 385     int *buf = NULL;
 386     gmx_bool bRead;
 387
 388     /* preallocate a small buffer and ip on the stack - if we need more
 389        we can always malloc(). This is faster for small values of size: */
 390     unsigned prealloc_size=3*16;
 391     int prealloc_ip[3*16], prealloc_buf[3*20];
 392     int we_should_free=0;
 393
 394     int minint[3], maxint[3], mindiff, *lip, diff;
 395     int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
 396     int minidx, maxidx;
 397     unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
 398     int flag, k;
 399     int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
 400     float *lfp, lf;
 401     int tmp, *thiscoord,  prevcoord[3];
 402     unsigned int tmpcoord[30];
 403
 404     int bufsize, xdrid, lsize;
 405     unsigned int bitsize;
 406     float inv_precision;
 407     int errval = 1;
 408     int rc;
 409
 410     bRead = (xdrs->x_op == XDR_DECODE);
 411     bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
 412     prevcoord[0]  = prevcoord[1]  = prevcoord[2]  = 0;
 413
 414     if (!bRead)
 415     {
 416         /* xdrs is open for writing */
 417
 418         if (xdr_int(xdrs, size) == 0)
 419             return 0;
 420         size3 = *size * 3;
 421         /* when the number of coordinates is small, don't try to compress; just
 422          * write them as floats using xdr_vector
 423          */
 424         if (*size <= 9 ) {
 425             return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
 426                     (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
 427         }
 428
 429         if(xdr_float(xdrs, precision) == 0)
 430             return 0;
 431
 432         if (size3 <= prealloc_size)
 433         {
 434             ip=prealloc_ip;
 435             buf=prealloc_buf;
 436         }
 437         else
 438         {
 439             we_should_free=1;
 440             bufsize = size3 * 1.2;
 441             ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
 442             buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
 443             if (ip == NULL || buf==NULL)
 444             {
 445                 fprintf(stderr,"malloc failed\n");
 446                 exit(1);
 447             }
 448         }
 449         /* buf[0-2] are special and do not contain actual data */
 450         buf[0] = buf[1] = buf[2] = 0;
 451         minint[0] = minint[1] = minint[2] = INT_MAX;
 452         maxint[0] = maxint[1] = maxint[2] = INT_MIN;
 453         prevrun = -1;
 454         lfp = fp;
 455         lip = ip;
 456         mindiff = INT_MAX;
 457         oldlint1 = oldlint2 = oldlint3 = 0;
 458         while(lfp < fp + size3 ) {
 459             /* find nearest integer */
 460             if (*lfp >= 0.0)
 461                 lf = *lfp * *precision + 0.5;
 462             else
 463                 lf = *lfp * *precision - 0.5;
 464             if (fabs(lf) > MAXABS) {
 465                 /* scaling would cause overflow */
 466                 errval = 0;
 467             }
 468             lint1 = lf;
 469             if (lint1 < minint[0]) minint[0] = lint1;
 470             if (lint1 > maxint[0]) maxint[0] = lint1;
 471             *lip++ = lint1;
 472             lfp++;
 473             if (*lfp >= 0.0)
 474                 lf = *lfp * *precision + 0.5;
 475             else
 476                 lf = *lfp * *precision - 0.5;
 477             if (fabs(lf) > MAXABS) {
 478                 /* scaling would cause overflow */
 479                 errval = 0;
 480             }
 481             lint2 = lf;
 482             if (lint2 < minint[1]) minint[1] = lint2;
 483             if (lint2 > maxint[1]) maxint[1] = lint2;
 484             *lip++ = lint2;
 485             lfp++;
 486             if (*lfp >= 0.0)
 487                 lf = *lfp * *precision + 0.5;
 488             else
 489                 lf = *lfp * *precision - 0.5;
 490             if (fabs(lf) > MAXABS) {
 491                 /* scaling would cause overflow */
 492                 errval = 0;
 493             }
 494             lint3 = lf;
 495             if (lint3 < minint[2]) minint[2] = lint3;
 496             if (lint3 > maxint[2]) maxint[2] = lint3;
 497             *lip++ = lint3;
 498             lfp++;
 499             diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
 500             if (diff < mindiff && lfp > fp + 3)
 501                 mindiff = diff;
 502             oldlint1 = lint1;
 503             oldlint2 = lint2;
 504             oldlint3 = lint3;
 505         }
 506         if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
 507                  (xdr_int(xdrs, &(minint[1])) == 0) ||
 508                  (xdr_int(xdrs, &(minint[2])) == 0) ||
 509              (xdr_int(xdrs, &(maxint[0])) == 0) ||
 510                  (xdr_int(xdrs, &(maxint[1])) == 0) ||
 511                  (xdr_int(xdrs, &(maxint[2])) == 0))
 512         {
 513             if (we_should_free)
 514             {
 515                 free(ip);
 516                 free(buf);
 517             }
 518             return 0;
 519         }
 520
 521         if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
 522                 (float)maxint[1] - (float)minint[1] >= MAXABS ||
 523                 (float)maxint[2] - (float)minint[2] >= MAXABS) {
 524             /* turning value in unsigned by subtracting minint
 525              * would cause overflow
 526              */
 527             errval = 0;
 528         }
 529         sizeint[0] = maxint[0] - minint[0]+1;
 530         sizeint[1] = maxint[1] - minint[1]+1;
 531         sizeint[2] = maxint[2] - minint[2]+1;
 532
 533         /* check if one of the sizes is to big to be multiplied */
 534         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
 535             bitsizeint[0] = sizeofint(sizeint[0]);
 536             bitsizeint[1] = sizeofint(sizeint[1]);
 537             bitsizeint[2] = sizeofint(sizeint[2]);
 538             bitsize = 0; /* flag the use of large sizes */
 539         } else {
 540             bitsize = sizeofints(3, sizeint);
 541         }
 542         lip = ip;
 543         luip = (unsigned int *) ip;
 544         smallidx = FIRSTIDX;
 545         while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
 546             smallidx++;
 547         }
 548         if(xdr_int(xdrs, &smallidx) == 0)
 549         {
 550             if (we_should_free)
 551             {
 552                 free(ip);
 553                 free(buf);
 554             }
 555             return 0;
 556         }
 557
 558         maxidx = MIN(LASTIDX, smallidx + 8) ;
 559         minidx = maxidx - 8; /* often this equal smallidx */
 560         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
 561         smallnum = magicints[smallidx] / 2;
 562         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 563         larger = magicints[maxidx] / 2;
 564         i = 0;
 565         while (i < *size) {
 566             is_small = 0;
 567             thiscoord = (int *)(luip) + i * 3;
 568             if (smallidx < maxidx && i >= 1 &&
 569                     abs(thiscoord[0] - prevcoord[0]) < larger &&
 570                     abs(thiscoord[1] - prevcoord[1]) < larger &&
 571                     abs(thiscoord[2] - prevcoord[2]) < larger) {
 572                 is_smaller = 1;
 573             } else if (smallidx > minidx) {
 574                 is_smaller = -1;
 575             } else {
 576                 is_smaller = 0;
 577             }
 578             if (i + 1 < *size) {
 579                 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
 580                         abs(thiscoord[1] - thiscoord[4]) < smallnum &&
 581                         abs(thiscoord[2] - thiscoord[5]) < smallnum) {
 582                     /* interchange first with second atom for better
 583                      * compression of water molecules
 584                      */
 585                     tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
 586                         thiscoord[3] = tmp;
 587                     tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
 588                         thiscoord[4] = tmp;
 589                     tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
 590                         thiscoord[5] = tmp;
 591                     is_small = 1;
 592                 }
 593
 594             }
 595             tmpcoord[0] = thiscoord[0] - minint[0];
 596             tmpcoord[1] = thiscoord[1] - minint[1];
 597             tmpcoord[2] = thiscoord[2] - minint[2];
 598             if (bitsize == 0) {
 599                 sendbits(buf, bitsizeint[0], tmpcoord[0]);
 600                 sendbits(buf, bitsizeint[1], tmpcoord[1]);
 601                 sendbits(buf, bitsizeint[2], tmpcoord[2]);
 602             } else {
 603                 sendints(buf, 3, bitsize, sizeint, tmpcoord);
 604             }
 605             prevcoord[0] = thiscoord[0];
 606             prevcoord[1] = thiscoord[1];
 607             prevcoord[2] = thiscoord[2];
 608             thiscoord = thiscoord + 3;
 609             i++;
 610
 611             run = 0;
 612             if (is_small == 0 && is_smaller == -1)
 613                 is_smaller = 0;
 614             while (is_small && run < 8*3) {
 615                 if (is_smaller == -1 && (
 616                         SQR(thiscoord[0] - prevcoord[0]) +
 617                         SQR(thiscoord[1] - prevcoord[1]) +
 618                         SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
 619                     is_smaller = 0;
 620                 }
 621
 622                 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
 623                 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
 624                 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
 625
 626                 prevcoord[0] = thiscoord[0];
 627                 prevcoord[1] = thiscoord[1];
 628                 prevcoord[2] = thiscoord[2];
 629
 630                 i++;
 631                 thiscoord = thiscoord + 3;
 632                 is_small = 0;
 633                 if (i < *size &&
 634                         abs(thiscoord[0] - prevcoord[0]) < smallnum &&
 635                         abs(thiscoord[1] - prevcoord[1]) < smallnum &&
 636                         abs(thiscoord[2] - prevcoord[2]) < smallnum) {
 637                     is_small = 1;
 638                 }
 639             }
 640             if (run != prevrun || is_smaller != 0) {
 641                 prevrun = run;
 642                 sendbits(buf, 1, 1); /* flag the change in run-length */
 643                 sendbits(buf, 5, run+is_smaller+1);
 644             } else {
 645                 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
 646             }
 647             for (k=0; k < run; k+=3) {
 648                 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
 649             }
 650             if (is_smaller != 0) {
 651                 smallidx += is_smaller;
 652                 if (is_smaller < 0) {
 653                     smallnum = smaller;
 654                     smaller = magicints[smallidx-1] / 2;
 655                 } else {
 656                     smaller = smallnum;
 657                     smallnum = magicints[smallidx] / 2;
 658                 }
 659                 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
 660             }
 661         }
 662         if (buf[1] != 0) buf[0]++;
 663                 /* buf[0] holds the length in bytes */
 664         if(xdr_int(xdrs, &(buf[0])) == 0)
 665         {
 666             if (we_should_free)
 667             {
 668                 free(ip);
 669                 free(buf);
 670             }
 671             return 0;
 672         }
 673
 674
 675         rc=errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
 676         if (we_should_free)
 677         {
 678             free(ip);
 679             free(buf);
 680         }
 681         return rc;
 682
 683     } else {
 684
 685         /* xdrs is open for reading */
 686
 687         if (xdr_int(xdrs, &lsize) == 0)
 688             return 0;
 689         if (*size != 0 && lsize != *size) {
 690             fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
 691                     "%d arg vs %d in file", *size, lsize);
 692         }
 693         *size = lsize;
 694         size3 = *size * 3;
 695         if (*size <= 9) {
 696             *precision = -1;
 697             return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
 698                     (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
 699         }
 700         if(xdr_float(xdrs, precision) == 0)
 701                 return 0;
 702
 703         if (size3 <= prealloc_size)
 704         {
 705             ip=prealloc_ip;
 706             buf=prealloc_buf;
 707         }
 708         else
 709         {
 710             we_should_free=1;
 711             bufsize = size3 * 1.2;
 712             ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
 713             buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
 714             if (ip == NULL || buf==NULL)
 715             {
 716                 fprintf(stderr,"malloc failed\n");
 717                 exit(1);
 718             }
 719         }
 720
 721         buf[0] = buf[1] = buf[2] = 0;
 722
 723         if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
 724                  (xdr_int(xdrs, &(minint[1])) == 0) ||
 725                  (xdr_int(xdrs, &(minint[2])) == 0) ||
 726                  (xdr_int(xdrs, &(maxint[0])) == 0) ||
 727                  (xdr_int(xdrs, &(maxint[1])) == 0) ||
 728                  (xdr_int(xdrs, &(maxint[2])) == 0))
 729         {
 730             if (we_should_free)
 731             {
 732                 free(ip);
 733                 free(buf);
 734             }
 735             return 0;
 736         }
 737
 738         sizeint[0] = maxint[0] - minint[0]+1;
 739         sizeint[1] = maxint[1] - minint[1]+1;
 740         sizeint[2] = maxint[2] - minint[2]+1;
 741
 742         /* check if one of the sizes is to big to be multiplied */
 743         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
 744             bitsizeint[0] = sizeofint(sizeint[0]);
 745             bitsizeint[1] = sizeofint(sizeint[1]);
 746             bitsizeint[2] = sizeofint(sizeint[2]);
 747             bitsize = 0; /* flag the use of large sizes */
 748         } else {
 749             bitsize = sizeofints(3, sizeint);
 750         }
 751
 752         if (xdr_int(xdrs, &smallidx) == 0)
 753         {
 754             if (we_should_free)
 755             {
 756                 free(ip);
 757                 free(buf);
 758             }
 759             return 0;
 760         }
 761
 762         maxidx = MIN(LASTIDX, smallidx + 8) ;
 763         minidx = maxidx - 8; /* often this equal smallidx */
 764         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
 765         smallnum = magicints[smallidx] / 2;
 766         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
 767         larger = magicints[maxidx];
 768
 769         /* buf[0] holds the length in bytes */
 770
 771         if (xdr_int(xdrs, &(buf[0])) == 0)
 772         {
 773             if (we_should_free)
 774             {
 775                 free(ip);
 776                 free(buf);
 777             }
 778             return 0;
 779         }
 780
 781
 782         if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
 783         {
 784             if (we_should_free)
 785             {
 786                 free(ip);
 787                 free(buf);
 788             }
 789             return 0;
 790         }
 791
 792
 793
 794         buf[0] = buf[1] = buf[2] = 0;
 795
 796         lfp = fp;
 797         inv_precision = 1.0 / * precision;
 798         run = 0;
 799         i = 0;
 800         lip = ip;
 801         while ( i < lsize ) {
 802             thiscoord = (int *)(lip) + i * 3;
 803
 804             if (bitsize == 0) {
 805                 thiscoord[0] = receivebits(buf, bitsizeint[0]);
 806                 thiscoord[1] = receivebits(buf, bitsizeint[1]);
 807                 thiscoord[2] = receivebits(buf, bitsizeint[2]);
 808             } else {
 809                 receiveints(buf, 3, bitsize, sizeint, thiscoord);
 810             }
 811
 812             i++;
 813             thiscoord[0] += minint[0];
 814             thiscoord[1] += minint[1];
 815             thiscoord[2] += minint[2];
 816
 817             prevcoord[0] = thiscoord[0];
 818             prevcoord[1] = thiscoord[1];
 819             prevcoord[2] = thiscoord[2];
 820
 821
 822             flag = receivebits(buf, 1);
 823             is_smaller = 0;
 824             if (flag == 1) {
 825                 run = receivebits(buf, 5);
 826                 is_smaller = run % 3;
 827                 run -= is_smaller;
 828                 is_smaller--;
 829             }
 830             if (run > 0) {
 831                 thiscoord += 3;
 832                 for (k = 0; k < run; k+=3) {
 833                     receiveints(buf, 3, smallidx, sizesmall, thiscoord);
 834                     i++;
 835                     thiscoord[0] += prevcoord[0] - smallnum;
 836                     thiscoord[1] += prevcoord[1] - smallnum;
 837                     thiscoord[2] += prevcoord[2] - smallnum;
 838                     if (k == 0) {
 839                         /* interchange first with second atom for better
 840                          * compression of water molecules
 841                          */
 842                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
 843                                 prevcoord[0] = tmp;
 844                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
 845                                 prevcoord[1] = tmp;
 846                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
 847                                 prevcoord[2] = tmp;
 848                         *lfp++ = prevcoord[0] * inv_precision;
 849                         *lfp++ = prevcoord[1] * inv_precision;
 850                         *lfp++ = prevcoord[2] * inv_precision;
 851                     } else {
 852                         prevcoord[0] = thiscoord[0];
 853                         prevcoord[1] = thiscoord[1];
 854                         prevcoord[2] = thiscoord[2];
 855                     }
 856                     *lfp++ = thiscoord[0] * inv_precision;
 857                     *lfp++ = thiscoord[1] * inv_precision;
 858                     *lfp++ = thiscoord[2] * inv_precision;
 859                 }
 860             } else {
 861                 *lfp++ = thiscoord[0] * inv_precision;
 862                 *lfp++ = thiscoord[1] * inv_precision;
 863                 *lfp++ = thiscoord[2] * inv_precision;
 864             }
 865             smallidx += is_smaller;
 866             if (is_smaller < 0) {
 867                 smallnum = smaller;
 868                 if (smallidx > FIRSTIDX) {
 869                     smaller = magicints[smallidx - 1] /2;
 870                 } else {
 871                     smaller = 0;
 872                 }
 873             } else if (is_smaller > 0) {
 874                 smaller = smallnum;
 875                 smallnum = magicints[smallidx] / 2;
 876             }
 877             sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
 878         }
 879     }
 880     if (we_should_free)
 881     {
 882         free(ip);
 883         free(buf);
 884     }
 885     return 1;
 886 }
 887
 888
 889
 890 /******************************************************************
 891
 892   XTC files have a relatively simple structure.
 893   They have a header of 16 bytes and the rest are
 894   the compressed coordinates of the files. Due to the
 895   compression 00 is not present in the coordinates.
 896   The first 4 bytes of the header are the magic number
 897   1995 (0x000007CB). If we find this number we are guaranteed
 898   to be in the header, due to the presence of so many zeros.
 899   The second 4 bytes are the number of atoms in the frame, and is
 900   assumed to be constant. The third 4 bytes are the frame number.
 901   The last 4 bytes are a floating point representation of the time.
 902
 903 ********************************************************************/
 904
 905 /* Must match definition in xtcio.c */
 906 #ifndef XTC_MAGIC
 907 #define XTC_MAGIC 1995
 908 #endif
 909
 910 static const int header_size = 16;
 911
 912 /* Check if we are at the header start.
 913    At the same time it will also read 1 int */
 914 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
 915                                int natoms, int * timestep, float * time){
 916   int i_inp[3];
 917   float f_inp[10];
 918   int i;
 919   gmx_off_t off;
 920
 921
 922   if((off = gmx_ftell(fp)) < 0){
 923     return -1;
 924   }
 925   /* read magic natoms and timestep */
 926   for(i = 0;i<3;i++){
 927     if(!xdr_int(xdrs, &(i_inp[i]))){
 928       gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET);
 929       return -1;
 930     }
 931   }
 932   /* quick return */
 933   if(i_inp[0] != XTC_MAGIC){
 934     if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
 935       return -1;
 936     }
 937     return 0;
 938   }
 939   /* read time and box */
 940   for(i = 0;i<10;i++){
 941     if(!xdr_float(xdrs, &(f_inp[i]))){
 942       gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET);
 943       return -1;
 944     }
 945   }
 946   /* Make a rigourous check to see if we are in the beggining of a header
 947      Hopefully there are no ambiguous cases */
 948   /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
 949      right triangle and that the first element must be nonzero unless the entire matrix is zero
 950   */
 951   if(i_inp[1] == natoms &&
 952      ((f_inp[1] != 0 && f_inp[6] == 0) ||
 953       (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0))){
 954     if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
 955       return -1;
 956     }
 957     *time = f_inp[0];
 958     *timestep = i_inp[2];
 959     return 1;
 960   }
 961   if(gmx_fseek(fp,off+XDR_INT_SIZE,SEEK_SET)){
 962     return -1;
 963   }
 964   return 0;
 965 }
 966
 967 static int
 968 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
 969 {
 970     gmx_off_t off;
 971     int step;
 972     float time;
 973     int ret;
 974
 975     if((off = gmx_ftell(fp)) < 0){
 976       return -1;
 977     }
 978
 979     /* read one int just to make sure we dont read this frame but the next */
 980     xdr_int(xdrs,&step);
 981     while(1){
 982       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
 983       if(ret == 1){
 984         if(gmx_fseek(fp,off,SEEK_SET)){
 985           return -1;
 986         }
 987         return step;
 988       }else if(ret == -1){
 989         if(gmx_fseek(fp,off,SEEK_SET)){
 990           return -1;
 991         }
 992       }
 993     }
 994     return -1;
 995 }
 996
 997
 998 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
 999                                      gmx_bool * bOK)
1000 {
1001     gmx_off_t off;
1002     float time;
1003     int step;
1004     int ret;
1005     *bOK = 0;
1006
1007     if ((off = gmx_ftell(fp)) < 0)
1008     {
1009         return -1;
1010     }
1011     /* read one int just to make sure we dont read this frame but the next */
1012     xdr_int(xdrs, &step);
1013     while (1)
1014     {
1015         ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1016         if (ret == 1)
1017         {
1018             *bOK = 1;
1019             if (gmx_fseek(fp,off,SEEK_SET))
1020             {
1021                 *bOK = 0;
1022                 return -1;
1023             }
1024             return time;
1025         }
1026         else if (ret == -1)
1027         {
1028             if (gmx_fseek(fp,off,SEEK_SET))
1029             {
1030                 return -1;
1031             }
1032             return -1;
1033         }
1034     }
1035     return -1;
1036 }
1037
1038
1039 static float
1040 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1041 {
1042     gmx_off_t off;
1043     int step;
1044     float time;
1045     int ret;
1046     *bOK = 0;
1047
1048     if ((off = gmx_ftell(fp)) < 0)
1049     {
1050         return -1;
1051     }
1052
1053     while (1)
1054     {
1055         ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1056         if (ret == 1)
1057         {
1058             *bOK = 1;
1059             if (gmx_fseek(fp,off,SEEK_SET))
1060             {
1061                 *bOK = 0;
1062                 return -1;
1063             }
1064             return time;
1065         }
1066         else if (ret == -1)
1067         {
1068             if (gmx_fseek(fp,off,SEEK_SET))
1069             {
1070                 return -1;
1071             }
1072             return -1;
1073         }
1074         else if (ret == 0)
1075         {
1076             /*Go back.*/
1077             if (gmx_fseek(fp,-2*XDR_INT_SIZE,SEEK_CUR))
1078             {
1079                 return -1;
1080             }
1081         }
1082     }
1083     return -1;
1084 }
1085
1086 /* Currently not used, just for completeness */
1087 static int
1088 xtc_get_current_frame_number(FILE *fp,XDR *xdrs,int natoms, gmx_bool * bOK)
1089 {
1090     gmx_off_t off;
1091     int ret;
1092     int step;
1093     float time;
1094     *bOK = 0;
1095
1096     if((off = gmx_ftell(fp)) < 0){
1097       return -1;
1098     }
1099
1100
1101     while(1){
1102       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
1103       if(ret == 1){
1104         *bOK = 1;
1105         if(gmx_fseek(fp,off,SEEK_SET)){
1106                 *bOK = 0;
1107           return -1;
1108         }
1109         return step;
1110       }else if(ret == -1){
1111         if(gmx_fseek(fp,off,SEEK_SET)){
1112           return -1;
1113         }
1114         return -1;
1115
1116       }else if(ret == 0){
1117                   /*Go back.*/
1118                   if(gmx_fseek(fp,-2*XDR_INT_SIZE,SEEK_CUR)){
1119                           return -1;
1120                   }
1121       }
1122     }
1123     return -1;
1124 }
1125
1126
1127
1128 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1129 {
1130     int inp;
1131     gmx_off_t res;
1132     int ret;
1133     int step;
1134     float time;
1135     /* read one int just to make sure we dont read this frame but the next */
1136     xdr_int(xdrs,&step);
1137     while(1)
1138     {
1139       ret = xtc_at_header_start(fp,xdrs,natoms,&step,&time);
1140       if(ret == 1){
1141         if((res = gmx_ftell(fp)) >= 0){
1142           return res - XDR_INT_SIZE;
1143         }else{
1144           return res;
1145         }
1146       }else if(ret == -1){
1147         return -1;
1148       }
1149     }
1150     return -1;
1151 }
1152
1153
1154 static
1155 float
1156 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1157 {
1158   float  res;
1159   float  tinit;
1160   gmx_off_t off;
1161
1162   *bOK = 0;
1163   if((off   = gmx_ftell(fp)) < 0){
1164     return -1;
1165   }
1166
1167     tinit = xtc_get_current_frame_time(fp,xdrs,natoms,bOK);
1168
1169     if(!(*bOK))
1170     {
1171         return -1;
1172     }
1173
1174     res = xtc_get_next_frame_time(fp,xdrs,natoms,bOK);
1175
1176     if(!(*bOK))
1177     {
1178         return -1;
1179     }
1180
1181     res -= tinit;
1182     if (0 != gmx_fseek(fp,off,SEEK_SET)) {
1183       *bOK = 0;
1184       return -1;
1185     }
1186     return res;
1187 }
1188
1189
1190 int
1191 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1192 {
1193     gmx_off_t low = 0;
1194     gmx_off_t high,pos;
1195
1196
1197     /* round to 4 bytes */
1198     int fr;
1199     gmx_off_t  offset;
1200     if(gmx_fseek(fp,0,SEEK_END)){
1201       return -1;
1202     }
1203
1204     if((high = gmx_ftell(fp)) < 0){
1205       return -1;
1206     }
1207
1208     /* round to 4 bytes  */
1209     high /= XDR_INT_SIZE;
1210     high *= XDR_INT_SIZE;
1211     offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1212
1213     if(gmx_fseek(fp,offset,SEEK_SET)){
1214       return -1;
1215     }
1216
1217     while(1)
1218     {
1219         fr = xtc_get_next_frame_number(fp,xdrs,natoms);
1220         if(fr < 0)
1221         {
1222             return -1;
1223         }
1224         if(fr != frame && abs(low-high) > header_size)
1225         {
1226             if(fr < frame)
1227             {
1228                 low = offset;
1229             }
1230             else
1231             {
1232                 high = offset;
1233             }
1234             /* round to 4 bytes */
1235             offset = (((high+low)/2)/4)*4;
1236
1237             if(gmx_fseek(fp,offset,SEEK_SET)){
1238               return -1;
1239             }
1240         }
1241         else
1242         {
1243             break;
1244         }
1245     }
1246     if(offset <= header_size)
1247     {
1248         offset = low;
1249     }
1250
1251     if(gmx_fseek(fp,offset,SEEK_SET)){
1252       return -1;
1253     }
1254
1255     if((pos = xtc_get_next_frame_start(fp,xdrs,natoms))< 0){
1256     /* we probably hit an end of file */
1257       return -1;
1258     }
1259
1260     if(gmx_fseek(fp,pos,SEEK_SET)){
1261       return -1;
1262     }
1263
1264     return 0;
1265 }
1266
1267
1268
1269 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms,gmx_bool bSeekForwardOnly)
1270 {
1271     float t;
1272     float dt;
1273     gmx_bool bOK = FALSE;
1274     gmx_off_t low = 0;
1275     gmx_off_t high, offset, pos;
1276     int res;
1277     int dt_sign = 0;
1278
1279     if (bSeekForwardOnly)
1280     {
1281         low = gmx_ftell(fp);
1282     }
1283     if (gmx_fseek(fp,0,SEEK_END))
1284     {
1285         return -1;
1286     }
1287
1288     if ((high = gmx_ftell(fp)) < 0)
1289     {
1290         return -1;
1291     }
1292     /* round to int  */
1293     high /= XDR_INT_SIZE;
1294     high *= XDR_INT_SIZE;
1295     offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1296
1297     if (gmx_fseek(fp,offset,SEEK_SET))
1298     {
1299         return -1;
1300     }
1301
1302
1303     /*
1304      * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1305     dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1306
1307     if (!bOK)
1308     {
1309         return -1;
1310     }
1311     */
1312
1313     while (1)
1314     {
1315         dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1316         if (!bOK)
1317         {
1318             return -1;
1319         }
1320         else
1321         {
1322             if (dt > 0)
1323             {
1324                 if (dt_sign == -1)
1325                 {
1326                     /* Found a place in the trajectory that has positive time step while
1327                      other has negative time step */
1328                     return -2;
1329                 }
1330                 dt_sign = 1;
1331             }
1332             else if (dt < 0)
1333             {
1334                 if (dt_sign == 1)
1335                 {
1336                     /* Found a place in the trajectory that has positive time step while
1337                      other has negative time step */
1338                     return -2;
1339                 }
1340                 dt_sign = -1;
1341             }
1342         }
1343         t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1344         if (!bOK)
1345         {
1346             return -1;
1347         }
1348
1349         /* If we are before the target time and the time step is positive or 0, or we have
1350          after the target time and the time step is negative, or the difference between
1351          the current time and the target time is bigger than dt and above all the distance between high
1352          and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1353          if we reached the solution */
1354         if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
1355             - time) >= dt && dt_sign >= 0)
1356             || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
1357             > header_size))
1358         {
1359             if (dt >= 0 && dt_sign != -1)
1360             {
1361                 if (t < time)
1362                 {
1363                     low = offset;
1364                 }
1365                 else
1366                 {
1367                     high = offset;
1368                 }
1369             }
1370             else if (dt <= 0 && dt_sign == -1)
1371             {
1372                 if (t >= time)
1373                 {
1374                     low = offset;
1375                 }
1376                 else
1377                 {
1378                     high = offset;
1379                 }
1380             }
1381             else
1382             {
1383                 /* We should never reach here */
1384                 return -1;
1385             }
1386             /* round to 4 bytes and subtract header*/
1387             offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1388             if (gmx_fseek(fp,offset,SEEK_SET))
1389             {
1390                 return -1;
1391             }
1392         }
1393         else
1394         {
1395             if (abs(low - high) <= header_size)
1396             {
1397                 break;
1398             }
1399             /* re-estimate dt */
1400             if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1401             {
1402                 if (bOK)
1403                 {
1404                     dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1405                 }
1406             }
1407             if (t >= time && t - time < dt)
1408             {
1409                 break;
1410             }
1411         }
1412     }
1413
1414     if (offset <= header_size)
1415     {
1416         offset = low;
1417     }
1418
1419     gmx_fseek(fp,offset,SEEK_SET);
1420
1421     if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1422     {
1423         return -1;
1424     }
1425
1426     if (gmx_fseek(fp,pos,SEEK_SET))
1427     {
1428         return -1;
1429     }
1430     return 0;
1431 }
1432
1433 float
1434 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1435 {
1436     float  time;
1437     gmx_off_t  off;
1438     int res;
1439     *bOK = 1;
1440     off = gmx_ftell(fp);
1441     if(off < 0){
1442       *bOK = 0;
1443       return -1;
1444     }
1445
1446     if( (res = gmx_fseek(fp,-3*XDR_INT_SIZE,SEEK_END)) != 0){
1447       *bOK = 0;
1448       return -1;
1449     }
1450
1451     time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1452     if(!(*bOK)){
1453       return -1;
1454     }
1455
1456     if( (res = gmx_fseek(fp,off,SEEK_SET)) != 0){
1457       *bOK = 0;
1458       return -1;
1459     }
1460     return time;
1461 }
1462
1463
1464 int
1465 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1466 {
1467     int    frame;
1468     gmx_off_t  off;
1469     int res;
1470     *bOK = 1;
1471
1472     if((off = gmx_ftell(fp)) < 0){
1473       *bOK = 0;
1474       return -1;
1475     }
1476
1477
1478     if(gmx_fseek(fp,-3*XDR_INT_SIZE,SEEK_END)){
1479       *bOK = 0;
1480       return -1;
1481     }
1482
1483     frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1484     if(!bOK){
1485       return -1;
1486     }
1487
1488
1489     if(gmx_fseek(fp,off,SEEK_SET)){
1490       *bOK = 0;
1491       return -1;
1492     }
1493
1494     return frame;
1495 }