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