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