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