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