src/gmxlib/libxdrf.c

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