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