src/gromacs/fileio/libxdrf.cpp

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