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