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