2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "xdr_datatype.h"
51 /* This is just for clarity - it can never be anything but 4! */
52 #define XDR_INT_SIZE 4
54 /* same order as the definition of xdr_datatype */
55 const char *xdr_datatype_names[] =
66 /*___________________________________________________________________________
68 | what follows are the C routine to read/write compressed coordinates together
69 | with some routines to assist in this task (those are marked
70 | static and cannot be called from user programs)
72 #define MAXABS INT_MAX-2
75 #define MIN(x, y) ((x) < (y) ? (x) : (y))
78 #define MAX(x, y) ((x) > (y) ? (x) : (y))
81 #define SQR(x) ((x)*(x))
83 static const int magicints[] = {
84 0, 0, 0, 0, 0, 0, 0, 0, 0,
85 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
86 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
87 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
88 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
89 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
90 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
91 8388607, 10568983, 13316085, 16777216
95 /* note that magicints[FIRSTIDX-1] == 0 */
96 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
99 /*____________________________________________________________________________
101 | sendbits - encode num into buf using the specified number of bits
103 | This routines appends the value of num to the bits already present in
104 | the array buf. You need to give it the number of bits to use and you
105 | better make sure that this number of bits is enough to hold the value
106 | Also num must be positive.
110 static void sendbits(int buf[], int num_of_bits, int num)
113 unsigned int cnt, lastbyte;
115 unsigned char * cbuf;
117 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
118 cnt = (unsigned int) buf[0];
120 lastbyte = (unsigned int) buf[2];
121 while (num_of_bits >= 8)
123 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
124 cbuf[cnt++] = lastbyte >> lastbits;
129 lastbyte = (lastbyte << num_of_bits) | num;
130 lastbits += num_of_bits;
134 cbuf[cnt++] = lastbyte >> lastbits;
142 cbuf[cnt] = lastbyte << (8 - lastbits);
146 /*_________________________________________________________________________
148 | sizeofint - calculate bitsize of an integer
150 | return the number of bits needed to store an integer with given max size
154 static int sizeofint(const int size)
159 while (size >= num && num_of_bits < 32)
167 /*___________________________________________________________________________
169 | sizeofints - calculate 'bitsize' of compressed ints
171 | given the number of small unsigned integers and the maximum value
172 | return the number of bits needed to read or write them with the
173 | routines receiveints and sendints. You need this parameter when
174 | calling these routines. Note that for many calls I can use
175 | the variable 'smallidx' which is exactly the number of bits, and
176 | So I don't need to call 'sizeofints for those calls.
179 static int sizeofints( const int num_of_ints, unsigned int sizes[])
183 unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
187 for (i = 0; i < num_of_ints; i++)
190 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
192 tmp = bytes[bytecnt] * sizes[i] + tmp;
193 bytes[bytecnt] = tmp & 0xff;
198 bytes[bytecnt++] = tmp & 0xff;
201 num_of_bytes = bytecnt;
205 while (bytes[num_of_bytes] >= num)
210 return num_of_bits + num_of_bytes * 8;
214 /*____________________________________________________________________________
216 | sendints - send a small set of small integers in compressed format
218 | this routine is used internally by xdr3dfcoord, to send a set of
219 | small integers to the buffer.
220 | Multiplication with fixed (specified maximum ) sizes is used to get
221 | to one big, multibyte integer. Allthough the routine could be
222 | modified to handle sizes bigger than 16777216, or more than just
223 | a few integers, this is not done, because the gain in compression
224 | isn't worth the effort. Note that overflowing the multiplication
225 | or the byte buffer (32 bytes) is unchecked and causes bad results.
229 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
230 unsigned int sizes[], unsigned int nums[])
233 int i, num_of_bytes, bytecnt;
234 unsigned int bytes[32], tmp;
240 bytes[num_of_bytes++] = tmp & 0xff;
245 for (i = 1; i < num_of_ints; i++)
247 if (nums[i] >= sizes[i])
249 fprintf(stderr, "major breakdown in sendints num %u doesn't "
250 "match size %u\n", nums[i], sizes[i]);
253 /* use one step multiply */
255 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
257 tmp = bytes[bytecnt] * sizes[i] + tmp;
258 bytes[bytecnt] = tmp & 0xff;
263 bytes[bytecnt++] = tmp & 0xff;
266 num_of_bytes = bytecnt;
268 if (num_of_bits >= num_of_bytes * 8)
270 for (i = 0; i < num_of_bytes; i++)
272 sendbits(buf, 8, bytes[i]);
274 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
278 for (i = 0; i < num_of_bytes-1; i++)
280 sendbits(buf, 8, bytes[i]);
282 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
287 /*___________________________________________________________________________
289 | receivebits - decode number from buf using specified number of bits
291 | extract the number of bits from the array buf and construct an integer
292 | from it. Return that value.
296 static int receivebits(int buf[], int num_of_bits)
299 int cnt, num, lastbits;
300 unsigned int lastbyte;
301 unsigned char * cbuf;
302 int mask = (1 << num_of_bits) -1;
304 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
306 lastbits = (unsigned int) buf[1];
307 lastbyte = (unsigned int) buf[2];
310 while (num_of_bits >= 8)
312 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
313 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
318 if (lastbits < num_of_bits)
321 lastbyte = (lastbyte << 8) | cbuf[cnt++];
323 lastbits -= num_of_bits;
324 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
333 /*____________________________________________________________________________
335 | receiveints - decode 'small' integers from the buf array
337 | this routine is the inverse from sendints() and decodes the small integers
338 | written to buf by calculating the remainder and doing divisions with
339 | the given sizes[]. You need to specify the total number of bits to be
340 | used from buf in num_of_bits.
344 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
345 unsigned int sizes[], int nums[])
348 int i, j, num_of_bytes, p, num;
350 bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
352 while (num_of_bits > 8)
354 bytes[num_of_bytes++] = receivebits(buf, 8);
359 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
361 for (i = num_of_ints-1; i > 0; i--)
364 for (j = num_of_bytes-1; j >= 0; j--)
366 num = (num << 8) | bytes[j];
369 num = num - p * sizes[i];
373 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
376 /*____________________________________________________________________________
378 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
380 | this routine reads or writes (depending on how you opened the file with
381 | xdropen() ) a large number of 3d coordinates (stored in *fp).
382 | The number of coordinates triplets to write is given by *size. On
383 | read this number may be zero, in which case it reads as many as were written
384 | or it may specify the number if triplets to read (which should match the
386 | Compression is achieved by first converting all floating numbers to integer
387 | using multiplication by *precision and rounding to the nearest integer.
388 | Then the minimum and maximum value are calculated to determine the range.
389 | The limited range of integers so found, is used to compress the coordinates.
390 | In addition the differences between succesive coordinates is calculated.
391 | If the difference happens to be 'small' then only the difference is saved,
392 | compressing the data even more. The notion of 'small' is changed dynamically
393 | and is enlarged or reduced whenever needed or possible.
394 | Extra compression is achieved in the case of GROMOS and coordinates of
395 | water molecules. GROMOS first writes out the Oxygen position, followed by
396 | the two hydrogens. In order to make the differences smaller (and thereby
397 | compression the data better) the order is changed into first one hydrogen
398 | then the oxygen, followed by the other hydrogen. This is rather special, but
399 | it shouldn't harm in the general case.
403 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
409 /* preallocate a small buffer and ip on the stack - if we need more
410 we can always malloc(). This is faster for small values of size: */
411 unsigned prealloc_size = 3*16;
412 int prealloc_ip[3*16], prealloc_buf[3*20];
413 int we_should_free = 0;
415 int minint[3], maxint[3], mindiff, *lip, diff;
416 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
418 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
420 int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
422 int tmp, *thiscoord, prevcoord[3];
423 unsigned int tmpcoord[30];
425 int bufsize, xdrid, lsize;
426 unsigned int bitsize;
431 bRead = (xdrs->x_op == XDR_DECODE);
432 bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
433 prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
437 /* xdrs is open for writing */
439 if (xdr_int(xdrs, size) == 0)
444 /* when the number of coordinates is small, don't try to compress; just
445 * write them as floats using xdr_vector
449 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
450 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
453 if (xdr_float(xdrs, precision) == 0)
458 if (size3 <= prealloc_size)
466 bufsize = size3 * 1.2;
467 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
468 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
469 if (ip == NULL || buf == NULL)
471 fprintf(stderr, "malloc failed\n");
475 /* buf[0-2] are special and do not contain actual data */
476 buf[0] = buf[1] = buf[2] = 0;
477 minint[0] = minint[1] = minint[2] = INT_MAX;
478 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
483 oldlint1 = oldlint2 = oldlint3 = 0;
484 while (lfp < fp + size3)
486 /* find nearest integer */
489 lf = *lfp * *precision + 0.5;
493 lf = *lfp * *precision - 0.5;
495 if (fabs(lf) > MAXABS)
497 /* scaling would cause overflow */
501 if (lint1 < minint[0])
505 if (lint1 > maxint[0])
513 lf = *lfp * *precision + 0.5;
517 lf = *lfp * *precision - 0.5;
519 if (fabs(lf) > MAXABS)
521 /* scaling would cause overflow */
525 if (lint2 < minint[1])
529 if (lint2 > maxint[1])
537 lf = *lfp * *precision + 0.5;
541 lf = *lfp * *precision - 0.5;
543 if (fabs(lf) > MAXABS)
545 /* scaling would cause overflow */
549 if (lint3 < minint[2])
553 if (lint3 > maxint[2])
559 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
560 if (diff < mindiff && lfp > fp + 3)
568 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
569 (xdr_int(xdrs, &(minint[1])) == 0) ||
570 (xdr_int(xdrs, &(minint[2])) == 0) ||
571 (xdr_int(xdrs, &(maxint[0])) == 0) ||
572 (xdr_int(xdrs, &(maxint[1])) == 0) ||
573 (xdr_int(xdrs, &(maxint[2])) == 0))
583 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
584 (float)maxint[1] - (float)minint[1] >= MAXABS ||
585 (float)maxint[2] - (float)minint[2] >= MAXABS)
587 /* turning value in unsigned by subtracting minint
588 * would cause overflow
592 sizeint[0] = maxint[0] - minint[0]+1;
593 sizeint[1] = maxint[1] - minint[1]+1;
594 sizeint[2] = maxint[2] - minint[2]+1;
596 /* check if one of the sizes is to big to be multiplied */
597 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
599 bitsizeint[0] = sizeofint(sizeint[0]);
600 bitsizeint[1] = sizeofint(sizeint[1]);
601 bitsizeint[2] = sizeofint(sizeint[2]);
602 bitsize = 0; /* flag the use of large sizes */
606 bitsize = sizeofints(3, sizeint);
609 luip = (unsigned int *) ip;
611 while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
615 if (xdr_int(xdrs, &smallidx) == 0)
625 maxidx = MIN(LASTIDX, smallidx + 8);
626 minidx = maxidx - 8; /* often this equal smallidx */
627 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
628 smallnum = magicints[smallidx] / 2;
629 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
630 larger = magicints[maxidx] / 2;
635 thiscoord = (int *)(luip) + i * 3;
636 if (smallidx < maxidx && i >= 1 &&
637 abs(thiscoord[0] - prevcoord[0]) < larger &&
638 abs(thiscoord[1] - prevcoord[1]) < larger &&
639 abs(thiscoord[2] - prevcoord[2]) < larger)
643 else if (smallidx > minidx)
653 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
654 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
655 abs(thiscoord[2] - thiscoord[5]) < smallnum)
657 /* interchange first with second atom for better
658 * compression of water molecules
660 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
662 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
664 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
670 tmpcoord[0] = thiscoord[0] - minint[0];
671 tmpcoord[1] = thiscoord[1] - minint[1];
672 tmpcoord[2] = thiscoord[2] - minint[2];
675 sendbits(buf, bitsizeint[0], tmpcoord[0]);
676 sendbits(buf, bitsizeint[1], tmpcoord[1]);
677 sendbits(buf, bitsizeint[2], tmpcoord[2]);
681 sendints(buf, 3, bitsize, sizeint, tmpcoord);
683 prevcoord[0] = thiscoord[0];
684 prevcoord[1] = thiscoord[1];
685 prevcoord[2] = thiscoord[2];
686 thiscoord = thiscoord + 3;
690 if (is_small == 0 && is_smaller == -1)
694 while (is_small && run < 8*3)
696 if (is_smaller == -1 && (
697 SQR(thiscoord[0] - prevcoord[0]) +
698 SQR(thiscoord[1] - prevcoord[1]) +
699 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller))
704 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
705 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
706 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
708 prevcoord[0] = thiscoord[0];
709 prevcoord[1] = thiscoord[1];
710 prevcoord[2] = thiscoord[2];
713 thiscoord = thiscoord + 3;
716 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
717 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
718 abs(thiscoord[2] - prevcoord[2]) < smallnum)
723 if (run != prevrun || is_smaller != 0)
726 sendbits(buf, 1, 1); /* flag the change in run-length */
727 sendbits(buf, 5, run+is_smaller+1);
731 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
733 for (k = 0; k < run; k += 3)
735 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
739 smallidx += is_smaller;
743 smaller = magicints[smallidx-1] / 2;
748 smallnum = magicints[smallidx] / 2;
750 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
757 /* buf[0] holds the length in bytes */
758 if (xdr_int(xdrs, &(buf[0])) == 0)
769 rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
781 /* xdrs is open for reading */
783 if (xdr_int(xdrs, &lsize) == 0)
787 if (*size != 0 && lsize != *size)
789 fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
790 "%d arg vs %d in file", *size, lsize);
797 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
798 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
800 if (xdr_float(xdrs, precision) == 0)
805 if (size3 <= prealloc_size)
813 bufsize = size3 * 1.2;
814 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
815 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
816 if (ip == NULL || buf == NULL)
818 fprintf(stderr, "malloc failed\n");
823 buf[0] = buf[1] = buf[2] = 0;
825 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
826 (xdr_int(xdrs, &(minint[1])) == 0) ||
827 (xdr_int(xdrs, &(minint[2])) == 0) ||
828 (xdr_int(xdrs, &(maxint[0])) == 0) ||
829 (xdr_int(xdrs, &(maxint[1])) == 0) ||
830 (xdr_int(xdrs, &(maxint[2])) == 0))
840 sizeint[0] = maxint[0] - minint[0]+1;
841 sizeint[1] = maxint[1] - minint[1]+1;
842 sizeint[2] = maxint[2] - minint[2]+1;
844 /* check if one of the sizes is to big to be multiplied */
845 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
847 bitsizeint[0] = sizeofint(sizeint[0]);
848 bitsizeint[1] = sizeofint(sizeint[1]);
849 bitsizeint[2] = sizeofint(sizeint[2]);
850 bitsize = 0; /* flag the use of large sizes */
854 bitsize = sizeofints(3, sizeint);
857 if (xdr_int(xdrs, &smallidx) == 0)
867 maxidx = MIN(LASTIDX, smallidx + 8);
868 minidx = maxidx - 8; /* often this equal smallidx */
869 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
870 smallnum = magicints[smallidx] / 2;
871 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
872 larger = magicints[maxidx];
874 /* buf[0] holds the length in bytes */
876 if (xdr_int(xdrs, &(buf[0])) == 0)
887 if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
899 buf[0] = buf[1] = buf[2] = 0;
902 inv_precision = 1.0 / *precision;
908 thiscoord = (int *)(lip) + i * 3;
912 thiscoord[0] = receivebits(buf, bitsizeint[0]);
913 thiscoord[1] = receivebits(buf, bitsizeint[1]);
914 thiscoord[2] = receivebits(buf, bitsizeint[2]);
918 receiveints(buf, 3, bitsize, sizeint, thiscoord);
922 thiscoord[0] += minint[0];
923 thiscoord[1] += minint[1];
924 thiscoord[2] += minint[2];
926 prevcoord[0] = thiscoord[0];
927 prevcoord[1] = thiscoord[1];
928 prevcoord[2] = thiscoord[2];
931 flag = receivebits(buf, 1);
935 run = receivebits(buf, 5);
936 is_smaller = run % 3;
943 for (k = 0; k < run; k += 3)
945 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
947 thiscoord[0] += prevcoord[0] - smallnum;
948 thiscoord[1] += prevcoord[1] - smallnum;
949 thiscoord[2] += prevcoord[2] - smallnum;
952 /* interchange first with second atom for better
953 * compression of water molecules
955 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
957 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
959 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
961 *lfp++ = prevcoord[0] * inv_precision;
962 *lfp++ = prevcoord[1] * inv_precision;
963 *lfp++ = prevcoord[2] * inv_precision;
967 prevcoord[0] = thiscoord[0];
968 prevcoord[1] = thiscoord[1];
969 prevcoord[2] = thiscoord[2];
971 *lfp++ = thiscoord[0] * inv_precision;
972 *lfp++ = thiscoord[1] * inv_precision;
973 *lfp++ = thiscoord[2] * inv_precision;
978 *lfp++ = thiscoord[0] * inv_precision;
979 *lfp++ = thiscoord[1] * inv_precision;
980 *lfp++ = thiscoord[2] * inv_precision;
982 smallidx += is_smaller;
986 if (smallidx > FIRSTIDX)
988 smaller = magicints[smallidx - 1] /2;
995 else if (is_smaller > 0)
998 smallnum = magicints[smallidx] / 2;
1000 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1013 /******************************************************************
1015 XTC files have a relatively simple structure.
1016 They have a header of 16 bytes and the rest are
1017 the compressed coordinates of the files. Due to the
1018 compression 00 is not present in the coordinates.
1019 The first 4 bytes of the header are the magic number
1020 1995 (0x000007CB). If we find this number we are guaranteed
1021 to be in the header, due to the presence of so many zeros.
1022 The second 4 bytes are the number of atoms in the frame, and is
1023 assumed to be constant. The third 4 bytes are the frame number.
1024 The last 4 bytes are a floating point representation of the time.
1026 ********************************************************************/
1028 /* Must match definition in xtcio.c */
1030 #define XTC_MAGIC 1995
1033 static const int header_size = 16;
1035 /* Check if we are at the header start.
1036 At the same time it will also read 1 int */
1037 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1038 int natoms, int * timestep, float * time)
1046 if ((off = gmx_ftell(fp)) < 0)
1050 /* read magic natoms and timestep */
1051 for (i = 0; i < 3; i++)
1053 if (!xdr_int(xdrs, &(i_inp[i])))
1055 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1060 if (i_inp[0] != XTC_MAGIC)
1062 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1068 /* read time and box */
1069 for (i = 0; i < 10; i++)
1071 if (!xdr_float(xdrs, &(f_inp[i])))
1073 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1077 /* Make a rigourous check to see if we are in the beggining of a header
1078 Hopefully there are no ambiguous cases */
1079 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1080 right triangle and that the first element must be nonzero unless the entire matrix is zero
1082 if (i_inp[1] == natoms &&
1083 ((f_inp[1] != 0 && f_inp[6] == 0) ||
1084 (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
1086 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1091 *timestep = i_inp[2];
1094 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1102 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1109 if ((off = gmx_ftell(fp)) < 0)
1114 /* read one int just to make sure we dont read this frame but the next */
1115 xdr_int(xdrs, &step);
1118 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1121 if (gmx_fseek(fp, off, SEEK_SET))
1129 if (gmx_fseek(fp, off, SEEK_SET))
1139 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1148 if ((off = gmx_ftell(fp)) < 0)
1152 /* read one int just to make sure we dont read this frame but the next */
1153 xdr_int(xdrs, &step);
1156 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1160 if (gmx_fseek(fp, off, SEEK_SET))
1169 if (gmx_fseek(fp, off, SEEK_SET))
1181 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1189 if ((off = gmx_ftell(fp)) < 0)
1196 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1200 if (gmx_fseek(fp, off, SEEK_SET))
1209 if (gmx_fseek(fp, off, SEEK_SET))
1218 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1227 /* Currently not used, just for completeness */
1229 xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1237 if ((off = gmx_ftell(fp)) < 0)
1245 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1249 if (gmx_fseek(fp, off, SEEK_SET))
1258 if (gmx_fseek(fp, off, SEEK_SET))
1268 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1279 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1286 /* read one int just to make sure we dont read this frame but the next */
1287 xdr_int(xdrs, &step);
1290 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1293 if ((res = gmx_ftell(fp)) >= 0)
1295 return res - XDR_INT_SIZE;
1313 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1320 if ((off = gmx_ftell(fp)) < 0)
1325 tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1332 res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
1340 if (0 != gmx_fseek(fp, off, SEEK_SET))
1349 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1352 gmx_off_t high, pos;
1355 /* round to 4 bytes */
1358 if (gmx_fseek(fp, 0, SEEK_END))
1363 if ((high = gmx_ftell(fp)) < 0)
1368 /* round to 4 bytes */
1369 high /= XDR_INT_SIZE;
1370 high *= XDR_INT_SIZE;
1371 offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1373 if (gmx_fseek(fp, offset, SEEK_SET))
1380 fr = xtc_get_next_frame_number(fp, xdrs, natoms);
1385 if (fr != frame && llabs(low-high) > header_size)
1395 /* round to 4 bytes */
1396 offset = (((high+low)/2)/4)*4;
1398 if (gmx_fseek(fp, offset, SEEK_SET))
1408 if (offset <= header_size)
1413 if (gmx_fseek(fp, offset, SEEK_SET))
1418 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1420 /* we probably hit an end of file */
1424 if (gmx_fseek(fp, pos, SEEK_SET))
1434 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
1438 gmx_bool bOK = FALSE;
1440 gmx_off_t high, offset, pos;
1444 if (bSeekForwardOnly)
1446 low = gmx_ftell(fp);
1448 if (gmx_fseek(fp, 0, SEEK_END))
1453 if ((high = gmx_ftell(fp)) < 0)
1458 high /= XDR_INT_SIZE;
1459 high *= XDR_INT_SIZE;
1460 offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1462 if (gmx_fseek(fp, offset, SEEK_SET))
1469 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1470 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1480 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1491 /* Found a place in the trajectory that has positive time step while
1492 other has negative time step */
1501 /* Found a place in the trajectory that has positive time step while
1502 other has negative time step */
1508 t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1514 /* If we are before the target time and the time step is positive or 0, or we have
1515 after the target time and the time step is negative, or the difference between
1516 the current time and the target time is bigger than dt and above all the distance between high
1517 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1518 if we reached the solution */
1519 if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) ||
1520 ((t - time) >= dt && dt_sign >= 0) || ((time - t) >= -dt && dt_sign < 0)) &&
1521 (llabs(low - high) > header_size))
1523 if (dt >= 0 && dt_sign != -1)
1534 else if (dt <= 0 && dt_sign == -1)
1547 /* We should never reach here */
1550 /* round to 4 bytes and subtract header*/
1551 offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1552 if (gmx_fseek(fp, offset, SEEK_SET))
1559 if (llabs(low - high) <= header_size)
1563 /* re-estimate dt */
1564 if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1568 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1571 if (t >= time && t - time < dt)
1578 if (offset <= header_size)
1583 gmx_fseek(fp, offset, SEEK_SET);
1585 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1590 if (gmx_fseek(fp, pos, SEEK_SET))
1598 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1604 off = gmx_ftell(fp);
1611 if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END)) != 0)
1617 time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1623 if ( (res = gmx_fseek(fp, off, SEEK_SET)) != 0)
1633 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1640 if ((off = gmx_ftell(fp)) < 0)
1647 if (gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END))
1653 frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1660 if (gmx_fseek(fp, off, SEEK_SET))