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