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