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