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, 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.
50 #include "gmx_fatal.h"
53 /* This is just for clarity - it can never be anything but 4! */
54 #define XDR_INT_SIZE 4
56 /* same order as the definition of xdr_datatype */
57 const char *xdr_datatype_names[] =
68 /*___________________________________________________________________________
70 | what follows are the C routine to read/write compressed coordinates together
71 | with some routines to assist in this task (those are marked
72 | static and cannot be called from user programs)
74 #define MAXABS INT_MAX-2
77 #define MIN(x, y) ((x) < (y) ? (x) : (y))
80 #define MAX(x, y) ((x) > (y) ? (x) : (y))
83 #define SQR(x) ((x)*(x))
85 static const int magicints[] = {
86 0, 0, 0, 0, 0, 0, 0, 0, 0,
87 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
88 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
89 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
90 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
91 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
92 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
93 8388607, 10568983, 13316085, 16777216
97 /* note that magicints[FIRSTIDX-1] == 0 */
98 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
101 /*____________________________________________________________________________
103 | sendbits - encode num into buf using the specified number of bits
105 | This routines appends the value of num to the bits already present in
106 | the array buf. You need to give it the number of bits to use and you
107 | better make sure that this number of bits is enough to hold the value
108 | Also num must be positive.
112 static void sendbits(int buf[], int num_of_bits, int num)
115 unsigned int cnt, lastbyte;
117 unsigned char * cbuf;
119 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
120 cnt = (unsigned int) buf[0];
122 lastbyte = (unsigned int) buf[2];
123 while (num_of_bits >= 8)
125 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
126 cbuf[cnt++] = lastbyte >> lastbits;
131 lastbyte = (lastbyte << num_of_bits) | num;
132 lastbits += num_of_bits;
136 cbuf[cnt++] = lastbyte >> lastbits;
144 cbuf[cnt] = lastbyte << (8 - lastbits);
148 /*_________________________________________________________________________
150 | sizeofint - calculate bitsize of an integer
152 | return the number of bits needed to store an integer with given max size
156 static int sizeofint(const int size)
161 while (size >= num && num_of_bits < 32)
169 /*___________________________________________________________________________
171 | sizeofints - calculate 'bitsize' of compressed ints
173 | given the number of small unsigned integers and the maximum value
174 | return the number of bits needed to read or write them with the
175 | routines receiveints and sendints. You need this parameter when
176 | calling these routines. Note that for many calls I can use
177 | the variable 'smallidx' which is exactly the number of bits, and
178 | So I don't need to call 'sizeofints for those calls.
181 static int sizeofints( const int num_of_ints, unsigned int sizes[])
185 unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
189 for (i = 0; i < num_of_ints; i++)
192 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
194 tmp = bytes[bytecnt] * sizes[i] + tmp;
195 bytes[bytecnt] = tmp & 0xff;
200 bytes[bytecnt++] = tmp & 0xff;
203 num_of_bytes = bytecnt;
207 while (bytes[num_of_bytes] >= num)
212 return num_of_bits + num_of_bytes * 8;
216 /*____________________________________________________________________________
218 | sendints - send a small set of small integers in compressed format
220 | this routine is used internally by xdr3dfcoord, to send a set of
221 | small integers to the buffer.
222 | Multiplication with fixed (specified maximum ) sizes is used to get
223 | to one big, multibyte integer. Allthough the routine could be
224 | modified to handle sizes bigger than 16777216, or more than just
225 | a few integers, this is not done, because the gain in compression
226 | isn't worth the effort. Note that overflowing the multiplication
227 | or the byte buffer (32 bytes) is unchecked and causes bad results.
231 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
232 unsigned int sizes[], unsigned int nums[])
235 int i, num_of_bytes, bytecnt;
236 unsigned int bytes[32], tmp;
242 bytes[num_of_bytes++] = tmp & 0xff;
247 for (i = 1; i < num_of_ints; i++)
249 if (nums[i] >= sizes[i])
251 fprintf(stderr, "major breakdown in sendints num %u doesn't "
252 "match size %u\n", nums[i], sizes[i]);
255 /* use one step multiply */
257 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
259 tmp = bytes[bytecnt] * sizes[i] + tmp;
260 bytes[bytecnt] = tmp & 0xff;
265 bytes[bytecnt++] = tmp & 0xff;
268 num_of_bytes = bytecnt;
270 if (num_of_bits >= num_of_bytes * 8)
272 for (i = 0; i < num_of_bytes; i++)
274 sendbits(buf, 8, bytes[i]);
276 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
280 for (i = 0; i < num_of_bytes-1; i++)
282 sendbits(buf, 8, bytes[i]);
284 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
289 /*___________________________________________________________________________
291 | receivebits - decode number from buf using specified number of bits
293 | extract the number of bits from the array buf and construct an integer
294 | from it. Return that value.
298 static int receivebits(int buf[], int num_of_bits)
301 int cnt, num, lastbits;
302 unsigned int lastbyte;
303 unsigned char * cbuf;
304 int mask = (1 << num_of_bits) -1;
306 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
308 lastbits = (unsigned int) buf[1];
309 lastbyte = (unsigned int) buf[2];
312 while (num_of_bits >= 8)
314 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
315 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
320 if (lastbits < num_of_bits)
323 lastbyte = (lastbyte << 8) | cbuf[cnt++];
325 lastbits -= num_of_bits;
326 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
335 /*____________________________________________________________________________
337 | receiveints - decode 'small' integers from the buf array
339 | this routine is the inverse from sendints() and decodes the small integers
340 | written to buf by calculating the remainder and doing divisions with
341 | the given sizes[]. You need to specify the total number of bits to be
342 | used from buf in num_of_bits.
346 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
347 unsigned int sizes[], int nums[])
350 int i, j, num_of_bytes, p, num;
352 bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
354 while (num_of_bits > 8)
356 bytes[num_of_bytes++] = receivebits(buf, 8);
361 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
363 for (i = num_of_ints-1; i > 0; i--)
366 for (j = num_of_bytes-1; j >= 0; j--)
368 num = (num << 8) | bytes[j];
371 num = num - p * sizes[i];
375 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
378 /*____________________________________________________________________________
380 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
382 | this routine reads or writes (depending on how you opened the file with
383 | xdropen() ) a large number of 3d coordinates (stored in *fp).
384 | The number of coordinates triplets to write is given by *size. On
385 | read this number may be zero, in which case it reads as many as were written
386 | or it may specify the number if triplets to read (which should match the
388 | Compression is achieved by first converting all floating numbers to integer
389 | using multiplication by *precision and rounding to the nearest integer.
390 | Then the minimum and maximum value are calculated to determine the range.
391 | The limited range of integers so found, is used to compress the coordinates.
392 | In addition the differences between succesive coordinates is calculated.
393 | If the difference happens to be 'small' then only the difference is saved,
394 | compressing the data even more. The notion of 'small' is changed dynamically
395 | and is enlarged or reduced whenever needed or possible.
396 | Extra compression is achieved in the case of GROMOS and coordinates of
397 | water molecules. GROMOS first writes out the Oxygen position, followed by
398 | the two hydrogens. In order to make the differences smaller (and thereby
399 | compression the data better) the order is changed into first one hydrogen
400 | then the oxygen, followed by the other hydrogen. This is rather special, but
401 | it shouldn't harm in the general case.
405 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
411 /* preallocate a small buffer and ip on the stack - if we need more
412 we can always malloc(). This is faster for small values of size: */
413 unsigned prealloc_size = 3*16;
414 int prealloc_ip[3*16], prealloc_buf[3*20];
415 int we_should_free = 0;
417 int minint[3], maxint[3], mindiff, *lip, diff;
418 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
420 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
422 int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
424 int tmp, *thiscoord, prevcoord[3];
425 unsigned int tmpcoord[30];
427 int bufsize, xdrid, lsize;
428 unsigned int bitsize;
433 bRead = (xdrs->x_op == XDR_DECODE);
434 bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
435 prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
439 /* xdrs is open for writing */
441 if (xdr_int(xdrs, size) == 0)
446 /* when the number of coordinates is small, don't try to compress; just
447 * write them as floats using xdr_vector
451 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
452 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
455 if (xdr_float(xdrs, precision) == 0)
460 if (size3 <= prealloc_size)
468 bufsize = size3 * 1.2;
469 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
470 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
471 if (ip == NULL || buf == NULL)
473 fprintf(stderr, "malloc failed\n");
477 /* buf[0-2] are special and do not contain actual data */
478 buf[0] = buf[1] = buf[2] = 0;
479 minint[0] = minint[1] = minint[2] = INT_MAX;
480 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
485 oldlint1 = oldlint2 = oldlint3 = 0;
486 while (lfp < fp + size3)
488 /* find nearest integer */
491 lf = *lfp * *precision + 0.5;
495 lf = *lfp * *precision - 0.5;
497 if (fabs(lf) > MAXABS)
499 /* scaling would cause overflow */
503 if (lint1 < minint[0])
507 if (lint1 > maxint[0])
515 lf = *lfp * *precision + 0.5;
519 lf = *lfp * *precision - 0.5;
521 if (fabs(lf) > MAXABS)
523 /* scaling would cause overflow */
527 if (lint2 < minint[1])
531 if (lint2 > maxint[1])
539 lf = *lfp * *precision + 0.5;
543 lf = *lfp * *precision - 0.5;
545 if (fabs(lf) > MAXABS)
547 /* scaling would cause overflow */
551 if (lint3 < minint[2])
555 if (lint3 > maxint[2])
561 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
562 if (diff < mindiff && lfp > fp + 3)
570 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
571 (xdr_int(xdrs, &(minint[1])) == 0) ||
572 (xdr_int(xdrs, &(minint[2])) == 0) ||
573 (xdr_int(xdrs, &(maxint[0])) == 0) ||
574 (xdr_int(xdrs, &(maxint[1])) == 0) ||
575 (xdr_int(xdrs, &(maxint[2])) == 0))
585 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
586 (float)maxint[1] - (float)minint[1] >= MAXABS ||
587 (float)maxint[2] - (float)minint[2] >= MAXABS)
589 /* turning value in unsigned by subtracting minint
590 * would cause overflow
594 sizeint[0] = maxint[0] - minint[0]+1;
595 sizeint[1] = maxint[1] - minint[1]+1;
596 sizeint[2] = maxint[2] - minint[2]+1;
598 /* check if one of the sizes is to big to be multiplied */
599 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
601 bitsizeint[0] = sizeofint(sizeint[0]);
602 bitsizeint[1] = sizeofint(sizeint[1]);
603 bitsizeint[2] = sizeofint(sizeint[2]);
604 bitsize = 0; /* flag the use of large sizes */
608 bitsize = sizeofints(3, sizeint);
611 luip = (unsigned int *) ip;
613 while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
617 if (xdr_int(xdrs, &smallidx) == 0)
627 maxidx = MIN(LASTIDX, smallidx + 8);
628 minidx = maxidx - 8; /* often this equal smallidx */
629 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
630 smallnum = magicints[smallidx] / 2;
631 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
632 larger = magicints[maxidx] / 2;
637 thiscoord = (int *)(luip) + i * 3;
638 if (smallidx < maxidx && i >= 1 &&
639 abs(thiscoord[0] - prevcoord[0]) < larger &&
640 abs(thiscoord[1] - prevcoord[1]) < larger &&
641 abs(thiscoord[2] - prevcoord[2]) < larger)
645 else if (smallidx > minidx)
655 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
656 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
657 abs(thiscoord[2] - thiscoord[5]) < smallnum)
659 /* interchange first with second atom for better
660 * compression of water molecules
662 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
664 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
666 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
672 tmpcoord[0] = thiscoord[0] - minint[0];
673 tmpcoord[1] = thiscoord[1] - minint[1];
674 tmpcoord[2] = thiscoord[2] - minint[2];
677 sendbits(buf, bitsizeint[0], tmpcoord[0]);
678 sendbits(buf, bitsizeint[1], tmpcoord[1]);
679 sendbits(buf, bitsizeint[2], tmpcoord[2]);
683 sendints(buf, 3, bitsize, sizeint, tmpcoord);
685 prevcoord[0] = thiscoord[0];
686 prevcoord[1] = thiscoord[1];
687 prevcoord[2] = thiscoord[2];
688 thiscoord = thiscoord + 3;
692 if (is_small == 0 && is_smaller == -1)
696 while (is_small && run < 8*3)
698 if (is_smaller == -1 && (
699 SQR(thiscoord[0] - prevcoord[0]) +
700 SQR(thiscoord[1] - prevcoord[1]) +
701 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller))
706 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
707 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
708 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
710 prevcoord[0] = thiscoord[0];
711 prevcoord[1] = thiscoord[1];
712 prevcoord[2] = thiscoord[2];
715 thiscoord = thiscoord + 3;
718 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
719 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
720 abs(thiscoord[2] - prevcoord[2]) < smallnum)
725 if (run != prevrun || is_smaller != 0)
728 sendbits(buf, 1, 1); /* flag the change in run-length */
729 sendbits(buf, 5, run+is_smaller+1);
733 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
735 for (k = 0; k < run; k += 3)
737 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
741 smallidx += is_smaller;
745 smaller = magicints[smallidx-1] / 2;
750 smallnum = magicints[smallidx] / 2;
752 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
759 /* buf[0] holds the length in bytes */
760 if (xdr_int(xdrs, &(buf[0])) == 0)
771 rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
783 /* xdrs is open for reading */
785 if (xdr_int(xdrs, &lsize) == 0)
789 if (*size != 0 && lsize != *size)
791 fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
792 "%d arg vs %d in file", *size, lsize);
799 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
800 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
802 if (xdr_float(xdrs, precision) == 0)
807 if (size3 <= prealloc_size)
815 bufsize = size3 * 1.2;
816 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
817 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
818 if (ip == NULL || buf == NULL)
820 fprintf(stderr, "malloc failed\n");
825 buf[0] = buf[1] = buf[2] = 0;
827 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
828 (xdr_int(xdrs, &(minint[1])) == 0) ||
829 (xdr_int(xdrs, &(minint[2])) == 0) ||
830 (xdr_int(xdrs, &(maxint[0])) == 0) ||
831 (xdr_int(xdrs, &(maxint[1])) == 0) ||
832 (xdr_int(xdrs, &(maxint[2])) == 0))
842 sizeint[0] = maxint[0] - minint[0]+1;
843 sizeint[1] = maxint[1] - minint[1]+1;
844 sizeint[2] = maxint[2] - minint[2]+1;
846 /* check if one of the sizes is to big to be multiplied */
847 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
849 bitsizeint[0] = sizeofint(sizeint[0]);
850 bitsizeint[1] = sizeofint(sizeint[1]);
851 bitsizeint[2] = sizeofint(sizeint[2]);
852 bitsize = 0; /* flag the use of large sizes */
856 bitsize = sizeofints(3, sizeint);
859 if (xdr_int(xdrs, &smallidx) == 0)
869 maxidx = MIN(LASTIDX, smallidx + 8);
870 minidx = maxidx - 8; /* often this equal smallidx */
871 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
872 smallnum = magicints[smallidx] / 2;
873 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
874 larger = magicints[maxidx];
876 /* buf[0] holds the length in bytes */
878 if (xdr_int(xdrs, &(buf[0])) == 0)
889 if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
901 buf[0] = buf[1] = buf[2] = 0;
904 inv_precision = 1.0 / *precision;
910 thiscoord = (int *)(lip) + i * 3;
914 thiscoord[0] = receivebits(buf, bitsizeint[0]);
915 thiscoord[1] = receivebits(buf, bitsizeint[1]);
916 thiscoord[2] = receivebits(buf, bitsizeint[2]);
920 receiveints(buf, 3, bitsize, sizeint, thiscoord);
924 thiscoord[0] += minint[0];
925 thiscoord[1] += minint[1];
926 thiscoord[2] += minint[2];
928 prevcoord[0] = thiscoord[0];
929 prevcoord[1] = thiscoord[1];
930 prevcoord[2] = thiscoord[2];
933 flag = receivebits(buf, 1);
937 run = receivebits(buf, 5);
938 is_smaller = run % 3;
945 for (k = 0; k < run; k += 3)
947 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
949 thiscoord[0] += prevcoord[0] - smallnum;
950 thiscoord[1] += prevcoord[1] - smallnum;
951 thiscoord[2] += prevcoord[2] - smallnum;
954 /* interchange first with second atom for better
955 * compression of water molecules
957 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
959 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
961 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
963 *lfp++ = prevcoord[0] * inv_precision;
964 *lfp++ = prevcoord[1] * inv_precision;
965 *lfp++ = prevcoord[2] * inv_precision;
969 prevcoord[0] = thiscoord[0];
970 prevcoord[1] = thiscoord[1];
971 prevcoord[2] = thiscoord[2];
973 *lfp++ = thiscoord[0] * inv_precision;
974 *lfp++ = thiscoord[1] * inv_precision;
975 *lfp++ = thiscoord[2] * inv_precision;
980 *lfp++ = thiscoord[0] * inv_precision;
981 *lfp++ = thiscoord[1] * inv_precision;
982 *lfp++ = thiscoord[2] * inv_precision;
984 smallidx += is_smaller;
988 if (smallidx > FIRSTIDX)
990 smaller = magicints[smallidx - 1] /2;
997 else if (is_smaller > 0)
1000 smallnum = magicints[smallidx] / 2;
1002 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1015 /******************************************************************
1017 XTC files have a relatively simple structure.
1018 They have a header of 16 bytes and the rest are
1019 the compressed coordinates of the files. Due to the
1020 compression 00 is not present in the coordinates.
1021 The first 4 bytes of the header are the magic number
1022 1995 (0x000007CB). If we find this number we are guaranteed
1023 to be in the header, due to the presence of so many zeros.
1024 The second 4 bytes are the number of atoms in the frame, and is
1025 assumed to be constant. The third 4 bytes are the frame number.
1026 The last 4 bytes are a floating point representation of the time.
1028 ********************************************************************/
1030 /* Must match definition in xtcio.c */
1032 #define XTC_MAGIC 1995
1035 static const int header_size = 16;
1037 /* Check if we are at the header start.
1038 At the same time it will also read 1 int */
1039 static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1040 int natoms, int * timestep, float * time)
1048 if ((off = gmx_ftell(fp)) < 0)
1052 /* read magic natoms and timestep */
1053 for (i = 0; i < 3; i++)
1055 if (!xdr_int(xdrs, &(i_inp[i])))
1057 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1062 if (i_inp[0] != XTC_MAGIC)
1064 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1070 /* read time and box */
1071 for (i = 0; i < 10; i++)
1073 if (!xdr_float(xdrs, &(f_inp[i])))
1075 gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
1079 /* Make a rigourous check to see if we are in the beggining of a header
1080 Hopefully there are no ambiguous cases */
1081 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1082 right triangle and that the first element must be nonzero unless the entire matrix is zero
1084 if (i_inp[1] == natoms &&
1085 ((f_inp[1] != 0 && f_inp[6] == 0) ||
1086 (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
1088 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1093 *timestep = i_inp[2];
1096 if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
1104 xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1111 if ((off = gmx_ftell(fp)) < 0)
1116 /* read one int just to make sure we dont read this frame but the next */
1117 xdr_int(xdrs, &step);
1120 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1123 if (gmx_fseek(fp, off, SEEK_SET))
1131 if (gmx_fseek(fp, off, SEEK_SET))
1141 static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1150 if ((off = gmx_ftell(fp)) < 0)
1154 /* read one int just to make sure we dont read this frame but the next */
1155 xdr_int(xdrs, &step);
1158 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1162 if (gmx_fseek(fp, off, SEEK_SET))
1171 if (gmx_fseek(fp, off, SEEK_SET))
1183 xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1191 if ((off = gmx_ftell(fp)) < 0)
1198 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1202 if (gmx_fseek(fp, off, SEEK_SET))
1211 if (gmx_fseek(fp, off, SEEK_SET))
1220 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1229 /* Currently not used, just for completeness */
1231 xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1239 if ((off = gmx_ftell(fp)) < 0)
1247 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1251 if (gmx_fseek(fp, off, SEEK_SET))
1260 if (gmx_fseek(fp, off, SEEK_SET))
1270 if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
1281 static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1288 /* read one int just to make sure we dont read this frame but the next */
1289 xdr_int(xdrs, &step);
1292 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1295 if ((res = gmx_ftell(fp)) >= 0)
1297 return res - XDR_INT_SIZE;
1315 xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1322 if ((off = gmx_ftell(fp)) < 0)
1327 tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1334 res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
1342 if (0 != gmx_fseek(fp, off, SEEK_SET))
1352 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1355 gmx_off_t high, pos;
1358 /* round to 4 bytes */
1361 if (gmx_fseek(fp, 0, SEEK_END))
1366 if ((high = gmx_ftell(fp)) < 0)
1371 /* round to 4 bytes */
1372 high /= XDR_INT_SIZE;
1373 high *= XDR_INT_SIZE;
1374 offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
1376 if (gmx_fseek(fp, offset, SEEK_SET))
1383 fr = xtc_get_next_frame_number(fp, xdrs, natoms);
1388 if (fr != frame && abs(low-high) > header_size)
1398 /* round to 4 bytes */
1399 offset = (((high+low)/2)/4)*4;
1401 if (gmx_fseek(fp, offset, SEEK_SET))
1411 if (offset <= header_size)
1416 if (gmx_fseek(fp, offset, SEEK_SET))
1421 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1423 /* we probably hit an end of file */
1427 if (gmx_fseek(fp, pos, SEEK_SET))
1437 int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
1441 gmx_bool bOK = FALSE;
1443 gmx_off_t high, offset, pos;
1447 if (bSeekForwardOnly)
1449 low = gmx_ftell(fp);
1451 if (gmx_fseek(fp, 0, SEEK_END))
1456 if ((high = gmx_ftell(fp)) < 0)
1461 high /= XDR_INT_SIZE;
1462 high *= XDR_INT_SIZE;
1463 offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1465 if (gmx_fseek(fp, offset, SEEK_SET))
1472 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1473 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1483 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1494 /* Found a place in the trajectory that has positive time step while
1495 other has negative time step */
1504 /* Found a place in the trajectory that has positive time step while
1505 other has negative time step */
1511 t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1517 /* If we are before the target time and the time step is positive or 0, or we have
1518 after the target time and the time step is negative, or the difference between
1519 the current time and the target time is bigger than dt and above all the distance between high
1520 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1521 if we reached the solution */
1522 if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
1523 - time) >= dt && dt_sign >= 0)
1524 || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
1527 if (dt >= 0 && dt_sign != -1)
1538 else if (dt <= 0 && dt_sign == -1)
1551 /* We should never reach here */
1554 /* round to 4 bytes and subtract header*/
1555 offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
1556 if (gmx_fseek(fp, offset, SEEK_SET))
1563 if (abs(low - high) <= header_size)
1567 /* re-estimate dt */
1568 if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1572 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1575 if (t >= time && t - time < dt)
1582 if (offset <= header_size)
1587 gmx_fseek(fp, offset, SEEK_SET);
1589 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1594 if (gmx_fseek(fp, pos, SEEK_SET))
1602 xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1608 off = gmx_ftell(fp);
1615 if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END)) != 0)
1621 time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1627 if ( (res = gmx_fseek(fp, off, SEEK_SET)) != 0)
1637 xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1644 if ((off = gmx_ftell(fp)) < 0)
1651 if (gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END))
1657 frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1664 if (gmx_fseek(fp, off, SEEK_SET))