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.
45 #include "gromacs/fileio/xdr_datatype.h"
46 #include "gromacs/fileio/xdrf.h"
47 #include "gromacs/utility/futil.h"
49 /* This is just for clarity - it can never be anything but 4! */
50 #define XDR_INT_SIZE 4
52 /* same order as the definition of xdr_datatype */
53 const char *xdr_datatype_names[] =
64 /*___________________________________________________________________________
66 | what follows are the C routine to read/write compressed coordinates together
67 | with some routines to assist in this task (those are marked
68 | static and cannot be called from user programs)
70 #define MAXABS INT_MAX-2
73 #define MIN(x, y) ((x) < (y) ? (x) : (y))
76 #define MAX(x, y) ((x) > (y) ? (x) : (y))
79 #define SQR(x) ((x)*(x))
81 static const int magicints[] = {
82 0, 0, 0, 0, 0, 0, 0, 0, 0,
83 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
84 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
85 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
86 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
87 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
88 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
89 8388607, 10568983, 13316085, 16777216
93 /* note that magicints[FIRSTIDX-1] == 0 */
94 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
97 /*____________________________________________________________________________
99 | sendbits - encode num into buf using the specified number of bits
101 | This routines appends the value of num to the bits already present in
102 | the array buf. You need to give it the number of bits to use and you
103 | better make sure that this number of bits is enough to hold the value
104 | Also num must be positive.
108 static void sendbits(int buf[], int num_of_bits, int num)
111 unsigned int cnt, lastbyte;
113 unsigned char * cbuf;
115 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
116 cnt = (unsigned int) buf[0];
118 lastbyte = (unsigned int) buf[2];
119 while (num_of_bits >= 8)
121 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
122 cbuf[cnt++] = lastbyte >> lastbits;
127 lastbyte = (lastbyte << num_of_bits) | num;
128 lastbits += num_of_bits;
132 cbuf[cnt++] = lastbyte >> lastbits;
140 cbuf[cnt] = lastbyte << (8 - lastbits);
144 /*_________________________________________________________________________
146 | sizeofint - calculate bitsize of an integer
148 | return the number of bits needed to store an integer with given max size
152 static int sizeofint(const int size)
157 while (size >= num && num_of_bits < 32)
165 /*___________________________________________________________________________
167 | sizeofints - calculate 'bitsize' of compressed ints
169 | given the number of small unsigned integers and the maximum value
170 | return the number of bits needed to read or write them with the
171 | routines receiveints and sendints. You need this parameter when
172 | calling these routines. Note that for many calls I can use
173 | the variable 'smallidx' which is exactly the number of bits, and
174 | So I don't need to call 'sizeofints for those calls.
177 static int sizeofints( const int num_of_ints, unsigned int sizes[])
181 unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
185 for (i = 0; i < num_of_ints; i++)
188 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
190 tmp = bytes[bytecnt] * sizes[i] + tmp;
191 bytes[bytecnt] = tmp & 0xff;
196 bytes[bytecnt++] = tmp & 0xff;
199 num_of_bytes = bytecnt;
203 while (bytes[num_of_bytes] >= num)
208 return num_of_bits + num_of_bytes * 8;
212 /*____________________________________________________________________________
214 | sendints - send a small set of small integers in compressed format
216 | this routine is used internally by xdr3dfcoord, to send a set of
217 | small integers to the buffer.
218 | Multiplication with fixed (specified maximum ) sizes is used to get
219 | to one big, multibyte integer. Allthough the routine could be
220 | modified to handle sizes bigger than 16777216, or more than just
221 | a few integers, this is not done, because the gain in compression
222 | isn't worth the effort. Note that overflowing the multiplication
223 | or the byte buffer (32 bytes) is unchecked and causes bad results.
227 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
228 unsigned int sizes[], unsigned int nums[])
231 int i, num_of_bytes, bytecnt;
232 unsigned int bytes[32], tmp;
238 bytes[num_of_bytes++] = tmp & 0xff;
243 for (i = 1; i < num_of_ints; i++)
245 if (nums[i] >= sizes[i])
247 fprintf(stderr, "major breakdown in sendints num %u doesn't "
248 "match size %u\n", nums[i], sizes[i]);
251 /* use one step multiply */
253 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
255 tmp = bytes[bytecnt] * sizes[i] + tmp;
256 bytes[bytecnt] = tmp & 0xff;
261 bytes[bytecnt++] = tmp & 0xff;
264 num_of_bytes = bytecnt;
266 if (num_of_bits >= num_of_bytes * 8)
268 for (i = 0; i < num_of_bytes; i++)
270 sendbits(buf, 8, bytes[i]);
272 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
276 for (i = 0; i < num_of_bytes-1; i++)
278 sendbits(buf, 8, bytes[i]);
280 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
285 /*___________________________________________________________________________
287 | receivebits - decode number from buf using specified number of bits
289 | extract the number of bits from the array buf and construct an integer
290 | from it. Return that value.
294 static int receivebits(int buf[], int num_of_bits)
297 int cnt, num, lastbits;
298 unsigned int lastbyte;
299 unsigned char * cbuf;
300 int mask = (1 << num_of_bits) -1;
302 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
304 lastbits = (unsigned int) buf[1];
305 lastbyte = (unsigned int) buf[2];
308 while (num_of_bits >= 8)
310 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
311 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
316 if (lastbits < num_of_bits)
319 lastbyte = (lastbyte << 8) | cbuf[cnt++];
321 lastbits -= num_of_bits;
322 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
331 /*____________________________________________________________________________
333 | receiveints - decode 'small' integers from the buf array
335 | this routine is the inverse from sendints() and decodes the small integers
336 | written to buf by calculating the remainder and doing divisions with
337 | the given sizes[]. You need to specify the total number of bits to be
338 | used from buf in num_of_bits.
342 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
343 unsigned int sizes[], int nums[])
346 int i, j, num_of_bytes, p, num;
348 bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
350 while (num_of_bits > 8)
352 bytes[num_of_bytes++] = receivebits(buf, 8);
357 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
359 for (i = num_of_ints-1; i > 0; i--)
362 for (j = num_of_bytes-1; j >= 0; j--)
364 num = (num << 8) | bytes[j];
367 num = num - p * sizes[i];
371 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
374 /*____________________________________________________________________________
376 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
378 | this routine reads or writes (depending on how you opened the file with
379 | xdropen() ) a large number of 3d coordinates (stored in *fp).
380 | The number of coordinates triplets to write is given by *size. On
381 | read this number may be zero, in which case it reads as many as were written
382 | or it may specify the number if triplets to read (which should match the
384 | Compression is achieved by first converting all floating numbers to integer
385 | using multiplication by *precision and rounding to the nearest integer.
386 | Then the minimum and maximum value are calculated to determine the range.
387 | The limited range of integers so found, is used to compress the coordinates.
388 | In addition the differences between succesive coordinates is calculated.
389 | If the difference happens to be 'small' then only the difference is saved,
390 | compressing the data even more. The notion of 'small' is changed dynamically
391 | and is enlarged or reduced whenever needed or possible.
392 | Extra compression is achieved in the case of GROMOS and coordinates of
393 | water molecules. GROMOS first writes out the Oxygen position, followed by
394 | the two hydrogens. In order to make the differences smaller (and thereby
395 | compression the data better) the order is changed into first one hydrogen
396 | then the oxygen, followed by the other hydrogen. This is rather special, but
397 | it shouldn't harm in the general case.
401 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
407 /* preallocate a small buffer and ip on the stack - if we need more
408 we can always malloc(). This is faster for small values of size: */
409 unsigned prealloc_size = 3*16;
410 int prealloc_ip[3*16], prealloc_buf[3*20];
411 int we_should_free = 0;
413 int minint[3], maxint[3], mindiff, *lip, diff;
414 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
416 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
418 int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
420 int tmp, *thiscoord, prevcoord[3];
421 unsigned int tmpcoord[30];
423 int bufsize, xdrid, lsize;
424 unsigned int bitsize;
429 bRead = (xdrs->x_op == XDR_DECODE);
430 bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
431 prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
435 /* xdrs is open for writing */
437 if (xdr_int(xdrs, size) == 0)
442 /* when the number of coordinates is small, don't try to compress; just
443 * write them as floats using xdr_vector
447 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
448 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
451 if (xdr_float(xdrs, precision) == 0)
456 if (size3 <= prealloc_size)
464 bufsize = size3 * 1.2;
465 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
466 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
467 if (ip == NULL || buf == NULL)
469 fprintf(stderr, "malloc failed\n");
473 /* buf[0-2] are special and do not contain actual data */
474 buf[0] = buf[1] = buf[2] = 0;
475 minint[0] = minint[1] = minint[2] = INT_MAX;
476 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
481 oldlint1 = oldlint2 = oldlint3 = 0;
482 while (lfp < fp + size3)
484 /* find nearest integer */
487 lf = *lfp * *precision + 0.5;
491 lf = *lfp * *precision - 0.5;
493 if (fabs(lf) > MAXABS)
495 /* scaling would cause overflow */
499 if (lint1 < minint[0])
503 if (lint1 > maxint[0])
511 lf = *lfp * *precision + 0.5;
515 lf = *lfp * *precision - 0.5;
517 if (fabs(lf) > MAXABS)
519 /* scaling would cause overflow */
523 if (lint2 < minint[1])
527 if (lint2 > maxint[1])
535 lf = *lfp * *precision + 0.5;
539 lf = *lfp * *precision - 0.5;
541 if (fabs(lf) > MAXABS)
543 /* scaling would cause overflow */
547 if (lint3 < minint[2])
551 if (lint3 > maxint[2])
557 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
558 if (diff < mindiff && lfp > fp + 3)
566 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
567 (xdr_int(xdrs, &(minint[1])) == 0) ||
568 (xdr_int(xdrs, &(minint[2])) == 0) ||
569 (xdr_int(xdrs, &(maxint[0])) == 0) ||
570 (xdr_int(xdrs, &(maxint[1])) == 0) ||
571 (xdr_int(xdrs, &(maxint[2])) == 0))
581 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
582 (float)maxint[1] - (float)minint[1] >= MAXABS ||
583 (float)maxint[2] - (float)minint[2] >= MAXABS)
585 /* turning value in unsigned by subtracting minint
586 * would cause overflow
590 sizeint[0] = maxint[0] - minint[0]+1;
591 sizeint[1] = maxint[1] - minint[1]+1;
592 sizeint[2] = maxint[2] - minint[2]+1;
594 /* check if one of the sizes is to big to be multiplied */
595 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
597 bitsizeint[0] = sizeofint(sizeint[0]);
598 bitsizeint[1] = sizeofint(sizeint[1]);
599 bitsizeint[2] = sizeofint(sizeint[2]);
600 bitsize = 0; /* flag the use of large sizes */
604 bitsize = sizeofints(3, sizeint);
607 luip = (unsigned int *) ip;
609 while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
613 if (xdr_int(xdrs, &smallidx) == 0)
623 maxidx = MIN(LASTIDX, smallidx + 8);
624 minidx = maxidx - 8; /* often this equal smallidx */
625 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
626 smallnum = magicints[smallidx] / 2;
627 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
628 larger = magicints[maxidx] / 2;
633 thiscoord = (int *)(luip) + i * 3;
634 if (smallidx < maxidx && i >= 1 &&
635 abs(thiscoord[0] - prevcoord[0]) < larger &&
636 abs(thiscoord[1] - prevcoord[1]) < larger &&
637 abs(thiscoord[2] - prevcoord[2]) < larger)
641 else if (smallidx > minidx)
651 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
652 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
653 abs(thiscoord[2] - thiscoord[5]) < smallnum)
655 /* interchange first with second atom for better
656 * compression of water molecules
658 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
660 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
662 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
668 tmpcoord[0] = thiscoord[0] - minint[0];
669 tmpcoord[1] = thiscoord[1] - minint[1];
670 tmpcoord[2] = thiscoord[2] - minint[2];
673 sendbits(buf, bitsizeint[0], tmpcoord[0]);
674 sendbits(buf, bitsizeint[1], tmpcoord[1]);
675 sendbits(buf, bitsizeint[2], tmpcoord[2]);
679 sendints(buf, 3, bitsize, sizeint, tmpcoord);
681 prevcoord[0] = thiscoord[0];
682 prevcoord[1] = thiscoord[1];
683 prevcoord[2] = thiscoord[2];
684 thiscoord = thiscoord + 3;
688 if (is_small == 0 && is_smaller == -1)
692 while (is_small && run < 8*3)
694 if (is_smaller == -1 && (
695 SQR(thiscoord[0] - prevcoord[0]) +
696 SQR(thiscoord[1] - prevcoord[1]) +
697 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller))
702 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
703 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
704 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
706 prevcoord[0] = thiscoord[0];
707 prevcoord[1] = thiscoord[1];
708 prevcoord[2] = thiscoord[2];
711 thiscoord = thiscoord + 3;
714 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
715 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
716 abs(thiscoord[2] - prevcoord[2]) < smallnum)
721 if (run != prevrun || is_smaller != 0)
724 sendbits(buf, 1, 1); /* flag the change in run-length */
725 sendbits(buf, 5, run+is_smaller+1);
729 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
731 for (k = 0; k < run; k += 3)
733 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
737 smallidx += is_smaller;
741 smaller = magicints[smallidx-1] / 2;
746 smallnum = magicints[smallidx] / 2;
748 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
755 /* buf[0] holds the length in bytes */
756 if (xdr_int(xdrs, &(buf[0])) == 0)
767 rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
779 /* xdrs is open for reading */
781 if (xdr_int(xdrs, &lsize) == 0)
785 if (*size != 0 && lsize != *size)
787 fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
788 "%d arg vs %d in file", *size, lsize);
795 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
796 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
798 if (xdr_float(xdrs, precision) == 0)
803 if (size3 <= prealloc_size)
811 bufsize = size3 * 1.2;
812 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
813 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
814 if (ip == NULL || buf == NULL)
816 fprintf(stderr, "malloc failed\n");
821 buf[0] = buf[1] = buf[2] = 0;
823 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
824 (xdr_int(xdrs, &(minint[1])) == 0) ||
825 (xdr_int(xdrs, &(minint[2])) == 0) ||
826 (xdr_int(xdrs, &(maxint[0])) == 0) ||
827 (xdr_int(xdrs, &(maxint[1])) == 0) ||
828 (xdr_int(xdrs, &(maxint[2])) == 0))
838 sizeint[0] = maxint[0] - minint[0]+1;
839 sizeint[1] = maxint[1] - minint[1]+1;
840 sizeint[2] = maxint[2] - minint[2]+1;
842 /* check if one of the sizes is to big to be multiplied */
843 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
845 bitsizeint[0] = sizeofint(sizeint[0]);
846 bitsizeint[1] = sizeofint(sizeint[1]);
847 bitsizeint[2] = sizeofint(sizeint[2]);
848 bitsize = 0; /* flag the use of large sizes */
852 bitsize = sizeofints(3, sizeint);
855 if (xdr_int(xdrs, &smallidx) == 0)
865 maxidx = MIN(LASTIDX, smallidx + 8);
866 minidx = maxidx - 8; /* often this equal smallidx */
867 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
868 smallnum = magicints[smallidx] / 2;
869 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
870 larger = magicints[maxidx];
872 /* buf[0] holds the length in bytes */
874 if (xdr_int(xdrs, &(buf[0])) == 0)
885 if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
897 buf[0] = buf[1] = buf[2] = 0;
900 inv_precision = 1.0 / *precision;
906 thiscoord = (int *)(lip) + i * 3;
910 thiscoord[0] = receivebits(buf, bitsizeint[0]);
911 thiscoord[1] = receivebits(buf, bitsizeint[1]);
912 thiscoord[2] = receivebits(buf, bitsizeint[2]);
916 receiveints(buf, 3, bitsize, sizeint, thiscoord);
920 thiscoord[0] += minint[0];
921 thiscoord[1] += minint[1];
922 thiscoord[2] += minint[2];
924 prevcoord[0] = thiscoord[0];
925 prevcoord[1] = thiscoord[1];
926 prevcoord[2] = thiscoord[2];
929 flag = receivebits(buf, 1);
933 run = receivebits(buf, 5);
934 is_smaller = run % 3;
941 for (k = 0; k < run; k += 3)
943 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
945 thiscoord[0] += prevcoord[0] - smallnum;
946 thiscoord[1] += prevcoord[1] - smallnum;
947 thiscoord[2] += prevcoord[2] - smallnum;
950 /* interchange first with second atom for better
951 * compression of water molecules
953 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
955 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
957 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
959 *lfp++ = prevcoord[0] * inv_precision;
960 *lfp++ = prevcoord[1] * inv_precision;
961 *lfp++ = prevcoord[2] * inv_precision;
965 prevcoord[0] = thiscoord[0];
966 prevcoord[1] = thiscoord[1];
967 prevcoord[2] = thiscoord[2];
969 *lfp++ = thiscoord[0] * inv_precision;
970 *lfp++ = thiscoord[1] * inv_precision;
971 *lfp++ = thiscoord[2] * inv_precision;
976 *lfp++ = thiscoord[0] * inv_precision;
977 *lfp++ = thiscoord[1] * inv_precision;
978 *lfp++ = thiscoord[2] * inv_precision;
980 smallidx += is_smaller;
984 if (smallidx > FIRSTIDX)
986 smaller = magicints[smallidx - 1] /2;
993 else if (is_smaller > 0)
996 smallnum = magicints[smallidx] / 2;
998 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1011 /******************************************************************
1013 XTC files have a relatively simple structure.
1014 They have a header of 16 bytes and the rest are
1015 the compressed coordinates of the files. Due to the
1016 compression 00 is not present in the coordinates.
1017 The first 4 bytes of the header are the magic number
1018 1995 (0x000007CB). If we find this number we are guaranteed
1019 to be in the header, due to the presence of so many zeros.
1020 The second 4 bytes are the number of atoms in the frame, and is
1021 assumed to be constant. The third 4 bytes are the frame number.
1022 The last 4 bytes are a floating point representation of the time.
1024 ********************************************************************/
1026 /* Must match definition in xtcio.c */
1028 #define XTC_MAGIC 1995
1031 static const int header_size = 16;
1033 /* Check if we are at the header start.
1034 At the same time it will also read 1 int */
1035 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1036 int natoms, int * timestep, float * time)
1044 if ((off = gmx_ftell(fp)) < 0)
1048 /* read magic natoms and timestep */
1049 for (i = 0; i < 3; i++)
1051 if (!xdr_int(xdrs, &(i_inp[i])))
1053 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1058 if (i_inp[0] != XTC_MAGIC)
1060 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1066 /* read time and box */
1067 for (i = 0; i < 10; i++)
1069 if (!xdr_float(xdrs, &(f_inp[i])))
1071 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1075 /* Make a rigourous check to see if we are in the beggining of a header
1076 Hopefully there are no ambiguous cases */
1077 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1078 right triangle and that the first element must be nonzero unless the entire matrix is zero
1080 if (i_inp[1] == natoms &&
1081 ((f_inp[1] != 0 && f_inp[6] == 0) ||
1082 (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
1084 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1089 *timestep = i_inp[2];
1092 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1100 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1107 if ((off = gmx_ftell(fp)) < 0)
1112 /* read one int just to make sure we dont read this frame but the next */
1113 xdr_int(xdrs, &step);
1116 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1119 if (gmx_fseek(fp, off, SEEK_SET))
1127 if (gmx_fseek(fp, off, SEEK_SET))
1137 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1146 if ((off = gmx_ftell(fp)) < 0)
1150 /* read one int just to make sure we dont read this frame but the next */
1151 xdr_int(xdrs, &step);
1154 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1158 if (gmx_fseek(fp, off, SEEK_SET))
1167 if (gmx_fseek(fp, off, SEEK_SET))
1179 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1187 if ((off = gmx_ftell(fp)) < 0)
1194 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1198 if (gmx_fseek(fp, off, SEEK_SET))
1207 if (gmx_fseek(fp, off, SEEK_SET))
1216 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1225 /* Currently not used, just for completeness */
1227 xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1235 if ((off = gmx_ftell(fp)) < 0)
1243 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1247 if (gmx_fseek(fp, off, SEEK_SET))
1256 if (gmx_fseek(fp, off, SEEK_SET))
1266 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1277 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1284 /* read one int just to make sure we dont read this frame but the next */
1285 xdr_int(xdrs, &step);
1288 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1291 if ((res = gmx_ftell(fp)) >= 0)
1293 return res - XDR_INT_SIZE;
1311 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1318 if ((off = gmx_ftell(fp)) < 0)
1323 tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1330 res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
1338 if (0 != gmx_fseek(fp, off, SEEK_SET))
1347 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1350 gmx_off_t high, pos;
1353 /* round to 4 bytes */
1356 if (gmx_fseek(fp, 0, SEEK_END))
1361 if ((high = gmx_ftell(fp)) < 0)
1366 /* round to 4 bytes */
1367 high /= XDR_INT_SIZE;
1368 high *= XDR_INT_SIZE;
1369 offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1371 if (gmx_fseek(fp, offset, SEEK_SET))
1378 fr = xtc_get_next_frame_number(fp, xdrs, natoms);
1383 if (fr != frame && llabs(low-high) > header_size)
1393 /* round to 4 bytes */
1394 offset = (((high+low)/2)/4)*4;
1396 if (gmx_fseek(fp, offset, SEEK_SET))
1406 if (offset <= header_size)
1411 if (gmx_fseek(fp, offset, SEEK_SET))
1416 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1418 /* we probably hit an end of file */
1422 if (gmx_fseek(fp, pos, SEEK_SET))
1432 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
1436 gmx_bool bOK = FALSE;
1438 gmx_off_t high, offset, pos;
1442 if (bSeekForwardOnly)
1444 low = gmx_ftell(fp)-header_size;
1446 if (gmx_fseek(fp, 0, SEEK_END))
1451 if ((high = gmx_ftell(fp)) < 0)
1456 high /= XDR_INT_SIZE;
1457 high *= XDR_INT_SIZE;
1458 offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1460 if (gmx_fseek(fp, offset, SEEK_SET))
1467 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1468 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1478 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1489 /* Found a place in the trajectory that has positive time step while
1490 other has negative time step */
1499 /* Found a place in the trajectory that has positive time step while
1500 other has negative time step */
1506 t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1512 /* If we are before the target time and the time step is positive or 0, or we have
1513 after the target time and the time step is negative, or the difference between
1514 the current time and the target time is bigger than dt and above all the distance between high
1515 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1516 if we reached the solution */
1517 if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) ||
1518 ((t - time) >= dt && dt_sign >= 0) || ((time - t) >= -dt && dt_sign < 0)) &&
1519 (llabs(low - high) > header_size))
1521 if (dt >= 0 && dt_sign != -1)
1532 else if (dt <= 0 && dt_sign == -1)
1545 /* We should never reach here */
1548 /* round to 4 bytes and subtract header*/
1549 offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1550 if (gmx_fseek(fp, offset, SEEK_SET))
1557 if (llabs(low - high) <= header_size)
1561 /* re-estimate dt */
1562 if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1566 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1569 if (t >= time && t - time < dt)
1576 if (offset <= header_size)
1581 gmx_fseek(fp, offset, SEEK_SET);
1583 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1588 if (gmx_fseek(fp, pos, SEEK_SET))
1596 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1602 off = gmx_ftell(fp);
1609 if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END)) != 0)
1615 time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1621 if ( (res = gmx_fseek(fp, off, SEEK_SET)) != 0)
1631 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1638 if ((off = gmx_ftell(fp)) < 0)
1645 if (gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END))
1651 frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1658 if (gmx_fseek(fp, off, SEEK_SET))