3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
48 #include "gmx_fatal.h"
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
57 # define gmx_fseek(A, B, C) fseek(A, B, C)
58 # define gmx_ftell(A) ftell(A)
59 # define gmx_off_t int
64 /* This is just for clarity - it can never be anything but 4! */
65 #define XDR_INT_SIZE 4
67 /* same order as the definition of xdr_datatype */
68 const char *xdr_datatype_names[] =
79 /*___________________________________________________________________________
81 | what follows are the C routine to read/write compressed coordinates together
82 | with some routines to assist in this task (those are marked
83 | static and cannot be called from user programs)
85 #define MAXABS INT_MAX-2
88 #define MIN(x, y) ((x) < (y) ? (x) : (y))
91 #define MAX(x, y) ((x) > (y) ? (x) : (y))
94 #define SQR(x) ((x)*(x))
96 static const int magicints[] = {
97 0, 0, 0, 0, 0, 0, 0, 0, 0,
98 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
99 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
100 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
101 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
102 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
103 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
104 8388607, 10568983, 13316085, 16777216
108 /* note that magicints[FIRSTIDX-1] == 0 */
109 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
112 /*____________________________________________________________________________
114 | sendbits - encode num into buf using the specified number of bits
116 | This routines appends the value of num to the bits already present in
117 | the array buf. You need to give it the number of bits to use and you
118 | better make sure that this number of bits is enough to hold the value
119 | Also num must be positive.
123 static void sendbits(int buf[], int num_of_bits, int num)
126 unsigned int cnt, lastbyte;
128 unsigned char * cbuf;
130 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
131 cnt = (unsigned int) buf[0];
133 lastbyte = (unsigned int) buf[2];
134 while (num_of_bits >= 8)
136 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
137 cbuf[cnt++] = lastbyte >> lastbits;
142 lastbyte = (lastbyte << num_of_bits) | num;
143 lastbits += num_of_bits;
147 cbuf[cnt++] = lastbyte >> lastbits;
155 cbuf[cnt] = lastbyte << (8 - lastbits);
159 /*_________________________________________________________________________
161 | sizeofint - calculate bitsize of an integer
163 | return the number of bits needed to store an integer with given max size
167 static int sizeofint(const int size)
172 while (size >= num && num_of_bits < 32)
180 /*___________________________________________________________________________
182 | sizeofints - calculate 'bitsize' of compressed ints
184 | given the number of small unsigned integers and the maximum value
185 | return the number of bits needed to read or write them with the
186 | routines receiveints and sendints. You need this parameter when
187 | calling these routines. Note that for many calls I can use
188 | the variable 'smallidx' which is exactly the number of bits, and
189 | So I don't need to call 'sizeofints for those calls.
192 static int sizeofints( const int num_of_ints, unsigned int sizes[])
196 unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
200 for (i = 0; i < num_of_ints; i++)
203 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
205 tmp = bytes[bytecnt] * sizes[i] + tmp;
206 bytes[bytecnt] = tmp & 0xff;
211 bytes[bytecnt++] = tmp & 0xff;
214 num_of_bytes = bytecnt;
218 while (bytes[num_of_bytes] >= num)
223 return num_of_bits + num_of_bytes * 8;
227 /*____________________________________________________________________________
229 | sendints - send a small set of small integers in compressed format
231 | this routine is used internally by xdr3dfcoord, to send a set of
232 | small integers to the buffer.
233 | Multiplication with fixed (specified maximum ) sizes is used to get
234 | to one big, multibyte integer. Allthough the routine could be
235 | modified to handle sizes bigger than 16777216, or more than just
236 | a few integers, this is not done, because the gain in compression
237 | isn't worth the effort. Note that overflowing the multiplication
238 | or the byte buffer (32 bytes) is unchecked and causes bad results.
242 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
243 unsigned int sizes[], unsigned int nums[])
246 int i, num_of_bytes, bytecnt;
247 unsigned int bytes[32], tmp;
253 bytes[num_of_bytes++] = tmp & 0xff;
258 for (i = 1; i < num_of_ints; i++)
260 if (nums[i] >= sizes[i])
262 fprintf(stderr, "major breakdown in sendints num %u doesn't "
263 "match size %u\n", nums[i], sizes[i]);
266 /* use one step multiply */
268 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
270 tmp = bytes[bytecnt] * sizes[i] + tmp;
271 bytes[bytecnt] = tmp & 0xff;
276 bytes[bytecnt++] = tmp & 0xff;
279 num_of_bytes = bytecnt;
281 if (num_of_bits >= num_of_bytes * 8)
283 for (i = 0; i < num_of_bytes; i++)
285 sendbits(buf, 8, bytes[i]);
287 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
291 for (i = 0; i < num_of_bytes-1; i++)
293 sendbits(buf, 8, bytes[i]);
295 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
300 /*___________________________________________________________________________
302 | receivebits - decode number from buf using specified number of bits
304 | extract the number of bits from the array buf and construct an integer
305 | from it. Return that value.
309 static int receivebits(int buf[], int num_of_bits)
312 int cnt, num, lastbits;
313 unsigned int lastbyte;
314 unsigned char * cbuf;
315 int mask = (1 << num_of_bits) -1;
317 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
319 lastbits = (unsigned int) buf[1];
320 lastbyte = (unsigned int) buf[2];
323 while (num_of_bits >= 8)
325 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
326 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
331 if (lastbits < num_of_bits)
334 lastbyte = (lastbyte << 8) | cbuf[cnt++];
336 lastbits -= num_of_bits;
337 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
346 /*____________________________________________________________________________
348 | receiveints - decode 'small' integers from the buf array
350 | this routine is the inverse from sendints() and decodes the small integers
351 | written to buf by calculating the remainder and doing divisions with
352 | the given sizes[]. You need to specify the total number of bits to be
353 | used from buf in num_of_bits.
357 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
358 unsigned int sizes[], int nums[])
361 int i, j, num_of_bytes, p, num;
363 bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
365 while (num_of_bits > 8)
367 bytes[num_of_bytes++] = receivebits(buf, 8);
372 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
374 for (i = num_of_ints-1; i > 0; i--)
377 for (j = num_of_bytes-1; j >= 0; j--)
379 num = (num << 8) | bytes[j];
382 num = num - p * sizes[i];
386 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
389 /*____________________________________________________________________________
391 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
393 | this routine reads or writes (depending on how you opened the file with
394 | xdropen() ) a large number of 3d coordinates (stored in *fp).
395 | The number of coordinates triplets to write is given by *size. On
396 | read this number may be zero, in which case it reads as many as were written
397 | or it may specify the number if triplets to read (which should match the
399 | Compression is achieved by first converting all floating numbers to integer
400 | using multiplication by *precision and rounding to the nearest integer.
401 | Then the minimum and maximum value are calculated to determine the range.
402 | The limited range of integers so found, is used to compress the coordinates.
403 | In addition the differences between succesive coordinates is calculated.
404 | If the difference happens to be 'small' then only the difference is saved,
405 | compressing the data even more. The notion of 'small' is changed dynamically
406 | and is enlarged or reduced whenever needed or possible.
407 | Extra compression is achieved in the case of GROMOS and coordinates of
408 | water molecules. GROMOS first writes out the Oxygen position, followed by
409 | the two hydrogens. In order to make the differences smaller (and thereby
410 | compression the data better) the order is changed into first one hydrogen
411 | then the oxygen, followed by the other hydrogen. This is rather special, but
412 | it shouldn't harm in the general case.
416 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
422 /* preallocate a small buffer and ip on the stack - if we need more
423 we can always malloc(). This is faster for small values of size: */
424 unsigned prealloc_size = 3*16;
425 int prealloc_ip[3*16], prealloc_buf[3*20];
426 int we_should_free = 0;
428 int minint[3], maxint[3], mindiff, *lip, diff;
429 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
431 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
433 int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
435 int tmp, *thiscoord, prevcoord[3];
436 unsigned int tmpcoord[30];
438 int bufsize, xdrid, lsize;
439 unsigned int bitsize;
444 bRead = (xdrs->x_op == XDR_DECODE);
445 bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
446 prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
450 /* xdrs is open for writing */
452 if (xdr_int(xdrs, size) == 0)
457 /* when the number of coordinates is small, don't try to compress; just
458 * write them as floats using xdr_vector
462 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
463 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
466 if (xdr_float(xdrs, precision) == 0)
471 if (size3 <= prealloc_size)
479 bufsize = size3 * 1.2;
480 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
481 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
482 if (ip == NULL || buf == NULL)
484 fprintf(stderr, "malloc failed\n");
488 /* buf[0-2] are special and do not contain actual data */
489 buf[0] = buf[1] = buf[2] = 0;
490 minint[0] = minint[1] = minint[2] = INT_MAX;
491 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
496 oldlint1 = oldlint2 = oldlint3 = 0;
497 while (lfp < fp + size3)
499 /* find nearest integer */
502 lf = *lfp * *precision + 0.5;
506 lf = *lfp * *precision - 0.5;
508 if (fabs(lf) > MAXABS)
510 /* scaling would cause overflow */
514 if (lint1 < minint[0])
518 if (lint1 > maxint[0])
526 lf = *lfp * *precision + 0.5;
530 lf = *lfp * *precision - 0.5;
532 if (fabs(lf) > MAXABS)
534 /* scaling would cause overflow */
538 if (lint2 < minint[1])
542 if (lint2 > maxint[1])
550 lf = *lfp * *precision + 0.5;
554 lf = *lfp * *precision - 0.5;
556 if (fabs(lf) > MAXABS)
558 /* scaling would cause overflow */
562 if (lint3 < minint[2])
566 if (lint3 > maxint[2])
572 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
573 if (diff < mindiff && lfp > fp + 3)
581 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
582 (xdr_int(xdrs, &(minint[1])) == 0) ||
583 (xdr_int(xdrs, &(minint[2])) == 0) ||
584 (xdr_int(xdrs, &(maxint[0])) == 0) ||
585 (xdr_int(xdrs, &(maxint[1])) == 0) ||
586 (xdr_int(xdrs, &(maxint[2])) == 0))
596 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
597 (float)maxint[1] - (float)minint[1] >= MAXABS ||
598 (float)maxint[2] - (float)minint[2] >= MAXABS)
600 /* turning value in unsigned by subtracting minint
601 * would cause overflow
605 sizeint[0] = maxint[0] - minint[0]+1;
606 sizeint[1] = maxint[1] - minint[1]+1;
607 sizeint[2] = maxint[2] - minint[2]+1;
609 /* check if one of the sizes is to big to be multiplied */
610 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
612 bitsizeint[0] = sizeofint(sizeint[0]);
613 bitsizeint[1] = sizeofint(sizeint[1]);
614 bitsizeint[2] = sizeofint(sizeint[2]);
615 bitsize = 0; /* flag the use of large sizes */
619 bitsize = sizeofints(3, sizeint);
622 luip = (unsigned int *) ip;
624 while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
628 if (xdr_int(xdrs, &smallidx) == 0)
638 maxidx = MIN(LASTIDX, smallidx + 8);
639 minidx = maxidx - 8; /* often this equal smallidx */
640 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
641 smallnum = magicints[smallidx] / 2;
642 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
643 larger = magicints[maxidx] / 2;
648 thiscoord = (int *)(luip) + i * 3;
649 if (smallidx < maxidx && i >= 1 &&
650 abs(thiscoord[0] - prevcoord[0]) < larger &&
651 abs(thiscoord[1] - prevcoord[1]) < larger &&
652 abs(thiscoord[2] - prevcoord[2]) < larger)
656 else if (smallidx > minidx)
666 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
667 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
668 abs(thiscoord[2] - thiscoord[5]) < smallnum)
670 /* interchange first with second atom for better
671 * compression of water molecules
673 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
675 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
677 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
683 tmpcoord[0] = thiscoord[0] - minint[0];
684 tmpcoord[1] = thiscoord[1] - minint[1];
685 tmpcoord[2] = thiscoord[2] - minint[2];
688 sendbits(buf, bitsizeint[0], tmpcoord[0]);
689 sendbits(buf, bitsizeint[1], tmpcoord[1]);
690 sendbits(buf, bitsizeint[2], tmpcoord[2]);
694 sendints(buf, 3, bitsize, sizeint, tmpcoord);
696 prevcoord[0] = thiscoord[0];
697 prevcoord[1] = thiscoord[1];
698 prevcoord[2] = thiscoord[2];
699 thiscoord = thiscoord + 3;
703 if (is_small == 0 && is_smaller == -1)
707 while (is_small && run < 8*3)
709 if (is_smaller == -1 && (
710 SQR(thiscoord[0] - prevcoord[0]) +
711 SQR(thiscoord[1] - prevcoord[1]) +
712 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller))
717 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
718 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
719 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
721 prevcoord[0] = thiscoord[0];
722 prevcoord[1] = thiscoord[1];
723 prevcoord[2] = thiscoord[2];
726 thiscoord = thiscoord + 3;
729 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
730 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
731 abs(thiscoord[2] - prevcoord[2]) < smallnum)
736 if (run != prevrun || is_smaller != 0)
739 sendbits(buf, 1, 1); /* flag the change in run-length */
740 sendbits(buf, 5, run+is_smaller+1);
744 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
746 for (k = 0; k < run; k += 3)
748 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
752 smallidx += is_smaller;
756 smaller = magicints[smallidx-1] / 2;
761 smallnum = magicints[smallidx] / 2;
763 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
770 /* buf[0] holds the length in bytes */
771 if (xdr_int(xdrs, &(buf[0])) == 0)
782 rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
794 /* xdrs is open for reading */
796 if (xdr_int(xdrs, &lsize) == 0)
800 if (*size != 0 && lsize != *size)
802 fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
803 "%d arg vs %d in file", *size, lsize);
810 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
811 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
813 if (xdr_float(xdrs, precision) == 0)
818 if (size3 <= prealloc_size)
826 bufsize = size3 * 1.2;
827 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
828 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
829 if (ip == NULL || buf == NULL)
831 fprintf(stderr, "malloc failed\n");
836 buf[0] = buf[1] = buf[2] = 0;
838 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
839 (xdr_int(xdrs, &(minint[1])) == 0) ||
840 (xdr_int(xdrs, &(minint[2])) == 0) ||
841 (xdr_int(xdrs, &(maxint[0])) == 0) ||
842 (xdr_int(xdrs, &(maxint[1])) == 0) ||
843 (xdr_int(xdrs, &(maxint[2])) == 0))
853 sizeint[0] = maxint[0] - minint[0]+1;
854 sizeint[1] = maxint[1] - minint[1]+1;
855 sizeint[2] = maxint[2] - minint[2]+1;
857 /* check if one of the sizes is to big to be multiplied */
858 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 */
867 bitsize = sizeofints(3, sizeint);
870 if (xdr_int(xdrs, &smallidx) == 0)
880 maxidx = MIN(LASTIDX, smallidx + 8);
881 minidx = maxidx - 8; /* often this equal smallidx */
882 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
883 smallnum = magicints[smallidx] / 2;
884 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
885 larger = magicints[maxidx];
887 /* buf[0] holds the length in bytes */
889 if (xdr_int(xdrs, &(buf[0])) == 0)
900 if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
912 buf[0] = buf[1] = buf[2] = 0;
915 inv_precision = 1.0 / *precision;
921 thiscoord = (int *)(lip) + i * 3;
925 thiscoord[0] = receivebits(buf, bitsizeint[0]);
926 thiscoord[1] = receivebits(buf, bitsizeint[1]);
927 thiscoord[2] = receivebits(buf, bitsizeint[2]);
931 receiveints(buf, 3, bitsize, sizeint, thiscoord);
935 thiscoord[0] += minint[0];
936 thiscoord[1] += minint[1];
937 thiscoord[2] += minint[2];
939 prevcoord[0] = thiscoord[0];
940 prevcoord[1] = thiscoord[1];
941 prevcoord[2] = thiscoord[2];
944 flag = receivebits(buf, 1);
948 run = receivebits(buf, 5);
949 is_smaller = run % 3;
956 for (k = 0; k < run; k += 3)
958 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
960 thiscoord[0] += prevcoord[0] - smallnum;
961 thiscoord[1] += prevcoord[1] - smallnum;
962 thiscoord[2] += prevcoord[2] - smallnum;
965 /* interchange first with second atom for better
966 * compression of water molecules
968 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
970 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
972 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
974 *lfp++ = prevcoord[0] * inv_precision;
975 *lfp++ = prevcoord[1] * inv_precision;
976 *lfp++ = prevcoord[2] * inv_precision;
980 prevcoord[0] = thiscoord[0];
981 prevcoord[1] = thiscoord[1];
982 prevcoord[2] = thiscoord[2];
984 *lfp++ = thiscoord[0] * inv_precision;
985 *lfp++ = thiscoord[1] * inv_precision;
986 *lfp++ = thiscoord[2] * inv_precision;
991 *lfp++ = thiscoord[0] * inv_precision;
992 *lfp++ = thiscoord[1] * inv_precision;
993 *lfp++ = thiscoord[2] * inv_precision;
995 smallidx += is_smaller;
999 if (smallidx > FIRSTIDX)
1001 smaller = magicints[smallidx - 1] /2;
1008 else if (is_smaller > 0)
1011 smallnum = magicints[smallidx] / 2;
1013 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1026 /******************************************************************
1028 XTC files have a relatively simple structure.
1029 They have a header of 16 bytes and the rest are
1030 the compressed coordinates of the files. Due to the
1031 compression 00 is not present in the coordinates.
1032 The first 4 bytes of the header are the magic number
1033 1995 (0x000007CB). If we find this number we are guaranteed
1034 to be in the header, due to the presence of so many zeros.
1035 The second 4 bytes are the number of atoms in the frame, and is
1036 assumed to be constant. The third 4 bytes are the frame number.
1037 The last 4 bytes are a floating point representation of the time.
1039 ********************************************************************/
1041 /* Must match definition in xtcio.c */
1043 #define XTC_MAGIC 1995
1046 static const int header_size = 16;
1048 /* Check if we are at the header start.
1049 At the same time it will also read 1 int */
1050 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1051 int natoms, int * timestep, float * time)
1059 if ((off = gmx_ftell(fp)) < 0)
1063 /* read magic natoms and timestep */
1064 for (i = 0; i < 3; i++)
1066 if (!xdr_int(xdrs, &(i_inp[i])))
1068 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1073 if (i_inp[0] != XTC_MAGIC)
1075 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1081 /* read time and box */
1082 for (i = 0; i < 10; i++)
1084 if (!xdr_float(xdrs, &(f_inp[i])))
1086 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1090 /* Make a rigourous check to see if we are in the beggining of a header
1091 Hopefully there are no ambiguous cases */
1092 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1093 right triangle and that the first element must be nonzero unless the entire matrix is zero
1095 if (i_inp[1] == natoms &&
1096 ((f_inp[1] != 0 && f_inp[6] == 0) ||
1097 (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
1099 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1104 *timestep = i_inp[2];
1107 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1115 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1122 if ((off = gmx_ftell(fp)) < 0)
1127 /* read one int just to make sure we dont read this frame but the next */
1128 xdr_int(xdrs, &step);
1131 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1134 if (gmx_fseek(fp, off, SEEK_SET))
1142 if (gmx_fseek(fp, off, SEEK_SET))
1152 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1161 if ((off = gmx_ftell(fp)) < 0)
1165 /* read one int just to make sure we dont read this frame but the next */
1166 xdr_int(xdrs, &step);
1169 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1173 if (gmx_fseek(fp, off, SEEK_SET))
1182 if (gmx_fseek(fp, off, SEEK_SET))
1194 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1202 if ((off = gmx_ftell(fp)) < 0)
1209 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1213 if (gmx_fseek(fp, off, SEEK_SET))
1222 if (gmx_fseek(fp, off, SEEK_SET))
1231 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1240 /* Currently not used, just for completeness */
1242 xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1250 if ((off = gmx_ftell(fp)) < 0)
1258 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1262 if (gmx_fseek(fp, off, SEEK_SET))
1271 if (gmx_fseek(fp, off, SEEK_SET))
1281 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1292 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1299 /* read one int just to make sure we dont read this frame but the next */
1300 xdr_int(xdrs, &step);
1303 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1306 if ((res = gmx_ftell(fp)) >= 0)
1308 return res - XDR_INT_SIZE;
1326 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1333 if ((off = gmx_ftell(fp)) < 0)
1338 tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1345 res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
1353 if (0 != gmx_fseek(fp, off, SEEK_SET))
1363 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1366 gmx_off_t high, pos;
1369 /* round to 4 bytes */
1372 if (gmx_fseek(fp, 0, SEEK_END))
1377 if ((high = gmx_ftell(fp)) < 0)
1382 /* round to 4 bytes */
1383 high /= XDR_INT_SIZE;
1384 high *= XDR_INT_SIZE;
1385 offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1387 if (gmx_fseek(fp, offset, SEEK_SET))
1394 fr = xtc_get_next_frame_number(fp, xdrs, natoms);
1399 if (fr != frame && abs(low-high) > header_size)
1409 /* round to 4 bytes */
1410 offset = (((high+low)/2)/4)*4;
1412 if (gmx_fseek(fp, offset, SEEK_SET))
1422 if (offset <= header_size)
1427 if (gmx_fseek(fp, offset, SEEK_SET))
1432 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1434 /* we probably hit an end of file */
1438 if (gmx_fseek(fp, pos, SEEK_SET))
1448 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
1452 gmx_bool bOK = FALSE;
1454 gmx_off_t high, offset, pos;
1458 if (bSeekForwardOnly)
1460 low = gmx_ftell(fp);
1462 if (gmx_fseek(fp, 0, SEEK_END))
1467 if ((high = gmx_ftell(fp)) < 0)
1472 high /= XDR_INT_SIZE;
1473 high *= XDR_INT_SIZE;
1474 offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1476 if (gmx_fseek(fp, offset, SEEK_SET))
1483 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1484 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1494 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1505 /* Found a place in the trajectory that has positive time step while
1506 other has negative time step */
1515 /* Found a place in the trajectory that has positive time step while
1516 other has negative time step */
1522 t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1528 /* If we are before the target time and the time step is positive or 0, or we have
1529 after the target time and the time step is negative, or the difference between
1530 the current time and the target time is bigger than dt and above all the distance between high
1531 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1532 if we reached the solution */
1533 if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
1534 - time) >= dt && dt_sign >= 0)
1535 || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
1538 if (dt >= 0 && dt_sign != -1)
1549 else if (dt <= 0 && dt_sign == -1)
1562 /* We should never reach here */
1565 /* round to 4 bytes and subtract header*/
1566 offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1567 if (gmx_fseek(fp, offset, SEEK_SET))
1574 if (abs(low - high) <= header_size)
1578 /* re-estimate dt */
1579 if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1583 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1586 if (t >= time && t - time < dt)
1593 if (offset <= header_size)
1598 gmx_fseek(fp, offset, SEEK_SET);
1600 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1605 if (gmx_fseek(fp, pos, SEEK_SET))
1613 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1619 off = gmx_ftell(fp);
1626 if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END)) != 0)
1632 time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1638 if ( (res = gmx_fseek(fp, off, SEEK_SET)) != 0)
1648 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1655 if ((off = gmx_ftell(fp)) < 0)
1662 if (gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END))
1668 frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1675 if (gmx_fseek(fp, off, SEEK_SET))