src/gromacs/fileio/libxdrf.c

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