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