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