src/gromacs/fileio/libxdrf.c

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