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,2015,2016,2017, 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.
47 #include "gromacs/fileio/xdr_datatype.h"
48 #include "gromacs/fileio/xdrf.h"
49 #include "gromacs/utility/futil.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 SQR(x) ((x)*(x))
77 static const int magicints
[] = {
78 0, 0, 0, 0, 0, 0, 0, 0, 0,
79 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
80 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
81 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
82 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
83 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
84 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
85 8388607, 10568983, 13316085, 16777216
89 /* note that magicints[FIRSTIDX-1] == 0 */
90 #define LASTIDX static_cast<int>((sizeof(magicints) / sizeof(*magicints)))
93 /*____________________________________________________________________________
95 | sendbits - encode num into buf using the specified number of bits
97 | This routines appends the value of num to the bits already present in
98 | the array buf. You need to give it the number of bits to use and you
99 | better make sure that this number of bits is enough to hold the value
100 | Also num must be positive.
104 static void sendbits(int buf
[], int num_of_bits
, int num
)
107 unsigned int cnt
, lastbyte
;
109 unsigned char * cbuf
;
111 cbuf
= (reinterpret_cast<unsigned char *>(buf
)) + 3 * sizeof(*buf
);
112 cnt
= static_cast<unsigned int>(buf
[0]);
114 lastbyte
= static_cast<unsigned int>(buf
[2]);
115 while (num_of_bits
>= 8)
117 lastbyte
= (lastbyte
<< 8) | ((num
>> (num_of_bits
-8)) /* & 0xff*/);
118 cbuf
[cnt
++] = lastbyte
>> lastbits
;
123 lastbyte
= (lastbyte
<< num_of_bits
) | num
;
124 lastbits
+= num_of_bits
;
128 cbuf
[cnt
++] = lastbyte
>> lastbits
;
136 cbuf
[cnt
] = lastbyte
<< (8 - lastbits
);
140 /*_________________________________________________________________________
142 | sizeofint - calculate bitsize of an integer
144 | return the number of bits needed to store an integer with given max size
148 static int sizeofint(const int size
)
153 while (size
>= num
&& num_of_bits
< 32)
161 /*___________________________________________________________________________
163 | sizeofints - calculate 'bitsize' of compressed ints
165 | given the number of small unsigned integers and the maximum value
166 | return the number of bits needed to read or write them with the
167 | routines receiveints and sendints. You need this parameter when
168 | calling these routines. Note that for many calls I can use
169 | the variable 'smallidx' which is exactly the number of bits, and
170 | So I don't need to call 'sizeofints for those calls.
173 static int sizeofints( const int num_of_ints
, unsigned int sizes
[])
177 unsigned int num_of_bytes
, num_of_bits
, bytecnt
, tmp
;
181 for (i
= 0; i
< num_of_ints
; i
++)
184 for (bytecnt
= 0; bytecnt
< num_of_bytes
; bytecnt
++)
186 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
187 bytes
[bytecnt
] = tmp
& 0xff;
192 bytes
[bytecnt
++] = tmp
& 0xff;
195 num_of_bytes
= bytecnt
;
199 while (bytes
[num_of_bytes
] >= num
)
204 return num_of_bits
+ num_of_bytes
* 8;
208 /*____________________________________________________________________________
210 | sendints - send a small set of small integers in compressed format
212 | this routine is used internally by xdr3dfcoord, to send a set of
213 | small integers to the buffer.
214 | Multiplication with fixed (specified maximum ) sizes is used to get
215 | to one big, multibyte integer. Allthough the routine could be
216 | modified to handle sizes bigger than 16777216, or more than just
217 | a few integers, this is not done, because the gain in compression
218 | isn't worth the effort. Note that overflowing the multiplication
219 | or the byte buffer (32 bytes) is unchecked and causes bad results.
223 static void sendints(int buf
[], const int num_of_ints
, const int num_of_bits
,
224 unsigned int sizes
[], unsigned int nums
[])
227 int i
, num_of_bytes
, bytecnt
;
228 unsigned int bytes
[32], tmp
;
234 bytes
[num_of_bytes
++] = tmp
& 0xff;
239 for (i
= 1; i
< num_of_ints
; i
++)
241 if (nums
[i
] >= sizes
[i
])
243 fprintf(stderr
, "major breakdown in sendints num %u doesn't "
244 "match size %u\n", nums
[i
], sizes
[i
]);
247 /* use one step multiply */
249 for (bytecnt
= 0; bytecnt
< num_of_bytes
; bytecnt
++)
251 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
252 bytes
[bytecnt
] = tmp
& 0xff;
257 bytes
[bytecnt
++] = tmp
& 0xff;
260 num_of_bytes
= bytecnt
;
262 if (num_of_bits
>= num_of_bytes
* 8)
264 for (i
= 0; i
< num_of_bytes
; i
++)
266 sendbits(buf
, 8, bytes
[i
]);
268 sendbits(buf
, num_of_bits
- num_of_bytes
* 8, 0);
272 for (i
= 0; i
< num_of_bytes
-1; i
++)
274 sendbits(buf
, 8, bytes
[i
]);
276 sendbits(buf
, num_of_bits
- (num_of_bytes
-1) * 8, bytes
[i
]);
281 /*___________________________________________________________________________
283 | receivebits - decode number from buf using specified number of bits
285 | extract the number of bits from the array buf and construct an integer
286 | from it. Return that value.
290 static int receivebits(int buf
[], int num_of_bits
)
293 int cnt
, num
, lastbits
;
294 unsigned int lastbyte
;
295 unsigned char * cbuf
;
296 int mask
= (1 << num_of_bits
) -1;
298 cbuf
= reinterpret_cast<unsigned char *>(buf
) + 3 * sizeof(*buf
);
300 lastbits
= static_cast<unsigned int>(buf
[1]);
301 lastbyte
= static_cast<unsigned int>(buf
[2]);
304 while (num_of_bits
>= 8)
306 lastbyte
= ( lastbyte
<< 8 ) | cbuf
[cnt
++];
307 num
|= (lastbyte
>> lastbits
) << (num_of_bits
- 8);
312 if (lastbits
< num_of_bits
)
315 lastbyte
= (lastbyte
<< 8) | cbuf
[cnt
++];
317 lastbits
-= num_of_bits
;
318 num
|= (lastbyte
>> lastbits
) & ((1 << num_of_bits
) -1);
327 /*____________________________________________________________________________
329 | receiveints - decode 'small' integers from the buf array
331 | this routine is the inverse from sendints() and decodes the small integers
332 | written to buf by calculating the remainder and doing divisions with
333 | the given sizes[]. You need to specify the total number of bits to be
334 | used from buf in num_of_bits.
338 static void receiveints(int buf
[], const int num_of_ints
, int num_of_bits
,
339 unsigned int sizes
[], int nums
[])
342 int i
, j
, num_of_bytes
, p
, num
;
344 bytes
[0] = bytes
[1] = bytes
[2] = bytes
[3] = 0;
346 while (num_of_bits
> 8)
348 bytes
[num_of_bytes
++] = receivebits(buf
, 8);
353 bytes
[num_of_bytes
++] = receivebits(buf
, num_of_bits
);
355 for (i
= num_of_ints
-1; i
> 0; i
--)
358 for (j
= num_of_bytes
-1; j
>= 0; j
--)
360 num
= (num
<< 8) | bytes
[j
];
363 num
= num
- p
* sizes
[i
];
367 nums
[0] = bytes
[0] | (bytes
[1] << 8) | (bytes
[2] << 16) | (bytes
[3] << 24);
370 /*____________________________________________________________________________
372 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
374 | this routine reads or writes (depending on how you opened the file with
375 | xdropen() ) a large number of 3d coordinates (stored in *fp).
376 | The number of coordinates triplets to write is given by *size. On
377 | read this number may be zero, in which case it reads as many as were written
378 | or it may specify the number if triplets to read (which should match the
380 | Compression is achieved by first converting all floating numbers to integer
381 | using multiplication by *precision and rounding to the nearest integer.
382 | Then the minimum and maximum value are calculated to determine the range.
383 | The limited range of integers so found, is used to compress the coordinates.
384 | In addition the differences between succesive coordinates is calculated.
385 | If the difference happens to be 'small' then only the difference is saved,
386 | compressing the data even more. The notion of 'small' is changed dynamically
387 | and is enlarged or reduced whenever needed or possible.
388 | Extra compression is achieved in the case of GROMOS and coordinates of
389 | water molecules. GROMOS first writes out the Oxygen position, followed by
390 | the two hydrogens. In order to make the differences smaller (and thereby
391 | compression the data better) the order is changed into first one hydrogen
392 | then the oxygen, followed by the other hydrogen. This is rather special, but
393 | it shouldn't harm in the general case.
397 int xdr3dfcoord(XDR
*xdrs
, float *fp
, int *size
, float *precision
)
403 /* preallocate a small buffer and ip on the stack - if we need more
404 we can always malloc(). This is faster for small values of size: */
405 unsigned prealloc_size
= 3*16;
406 int prealloc_ip
[3*16], prealloc_buf
[3*20];
407 int we_should_free
= 0;
409 int minint
[3], maxint
[3], mindiff
, *lip
, diff
;
410 int lint1
, lint2
, lint3
, oldlint1
, oldlint2
, oldlint3
, smallidx
;
412 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
, *luip
;
414 int smallnum
, smaller
, larger
, i
, is_small
, is_smaller
, run
, prevrun
;
416 int tmp
, *thiscoord
, prevcoord
[3];
417 unsigned int tmpcoord
[30];
420 unsigned int bitsize
;
425 bRead
= (xdrs
->x_op
== XDR_DECODE
);
426 bitsizeint
[0] = bitsizeint
[1] = bitsizeint
[2] = 0;
427 prevcoord
[0] = prevcoord
[1] = prevcoord
[2] = 0;
429 // The static analyzer warns about garbage values for thiscoord[] further
430 // down. It might be thrown off by all the reinterpret_casts, but we might
431 // as well make sure the small preallocated buffer is zero-initialized.
432 for (i
= 0; i
< static_cast<int>(prealloc_size
); i
++)
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
, reinterpret_cast<char *>(fp
), static_cast<unsigned int>(size3
),
452 static_cast<unsigned int>(sizeof(*fp
)), (xdrproc_t
)xdr_float
));
455 if (xdr_float(xdrs
, precision
) == 0)
460 if (size3
<= prealloc_size
)
468 bufsize
= static_cast<int>(size3
* 1.2);
469 ip
= reinterpret_cast<int *>(malloc(size3
* sizeof(*ip
)));
470 buf
= reinterpret_cast<int *>(malloc(bufsize
* sizeof(*buf
)));
471 if (ip
== nullptr || buf
== nullptr)
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 */
502 lint1
= static_cast<int>(lf
);
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 */
526 lint2
= static_cast<int>(lf
);
527 if (lint2
< minint
[1])
531 if (lint2
> maxint
[1])
539 lf
= *lfp
* *precision
+ 0.5;
543 lf
= *lfp
* *precision
- 0.5;
545 if (std::abs(lf
) > MAXABS
)
547 /* scaling would cause overflow */
550 lint3
= static_cast<int>(lf
);
551 if (lint3
< minint
[2])
555 if (lint3
> maxint
[2])
561 diff
= std::abs(oldlint1
-lint1
)+std::abs(oldlint2
-lint2
)+std::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
);
610 luip
= reinterpret_cast<unsigned int *>(ip
);
612 while (smallidx
< LASTIDX
&& magicints
[smallidx
] < mindiff
)
616 if (xdr_int(xdrs
, &smallidx
) == 0)
626 maxidx
= std::min(LASTIDX
, smallidx
+ 8);
627 minidx
= maxidx
- 8; /* often this equal smallidx */
628 smaller
= magicints
[std::max(FIRSTIDX
, smallidx
-1)] / 2;
629 smallnum
= magicints
[smallidx
] / 2;
630 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
631 larger
= magicints
[maxidx
] / 2;
636 thiscoord
= reinterpret_cast<int *>(luip
) + i
* 3;
637 if (smallidx
< maxidx
&& i
>= 1 &&
638 std::abs(thiscoord
[0] - prevcoord
[0]) < larger
&&
639 std::abs(thiscoord
[1] - prevcoord
[1]) < larger
&&
640 std::abs(thiscoord
[2] - prevcoord
[2]) < larger
)
644 else if (smallidx
> minidx
)
654 if (std::abs(thiscoord
[0] - thiscoord
[3]) < smallnum
&&
655 std::abs(thiscoord
[1] - thiscoord
[4]) < smallnum
&&
656 std::abs(thiscoord
[2] - thiscoord
[5]) < smallnum
)
658 /* interchange first with second atom for better
659 * compression of water molecules
661 tmp
= thiscoord
[0]; thiscoord
[0] = thiscoord
[3];
663 tmp
= thiscoord
[1]; thiscoord
[1] = thiscoord
[4];
665 tmp
= thiscoord
[2]; thiscoord
[2] = thiscoord
[5];
671 tmpcoord
[0] = thiscoord
[0] - minint
[0];
672 tmpcoord
[1] = thiscoord
[1] - minint
[1];
673 tmpcoord
[2] = thiscoord
[2] - minint
[2];
676 sendbits(buf
, bitsizeint
[0], tmpcoord
[0]);
677 sendbits(buf
, bitsizeint
[1], tmpcoord
[1]);
678 sendbits(buf
, bitsizeint
[2], tmpcoord
[2]);
682 sendints(buf
, 3, bitsize
, sizeint
, tmpcoord
);
684 prevcoord
[0] = thiscoord
[0];
685 prevcoord
[1] = thiscoord
[1];
686 prevcoord
[2] = thiscoord
[2];
687 thiscoord
= thiscoord
+ 3;
691 if (is_small
== 0 && is_smaller
== -1)
695 while (is_small
&& run
< 8*3)
697 if (is_smaller
== -1 && (
698 SQR(thiscoord
[0] - prevcoord
[0]) +
699 SQR(thiscoord
[1] - prevcoord
[1]) +
700 SQR(thiscoord
[2] - prevcoord
[2]) >= smaller
* smaller
))
705 tmpcoord
[run
++] = thiscoord
[0] - prevcoord
[0] + smallnum
;
706 tmpcoord
[run
++] = thiscoord
[1] - prevcoord
[1] + smallnum
;
707 tmpcoord
[run
++] = thiscoord
[2] - prevcoord
[2] + smallnum
;
709 prevcoord
[0] = thiscoord
[0];
710 prevcoord
[1] = thiscoord
[1];
711 prevcoord
[2] = thiscoord
[2];
714 thiscoord
= thiscoord
+ 3;
717 abs(thiscoord
[0] - prevcoord
[0]) < smallnum
&&
718 abs(thiscoord
[1] - prevcoord
[1]) < smallnum
&&
719 abs(thiscoord
[2] - prevcoord
[2]) < smallnum
)
724 if (run
!= prevrun
|| is_smaller
!= 0)
727 sendbits(buf
, 1, 1); /* flag the change in run-length */
728 sendbits(buf
, 5, run
+is_smaller
+1);
732 sendbits(buf
, 1, 0); /* flag the fact that runlength did not change */
734 for (k
= 0; k
< run
; k
+= 3)
736 sendints(buf
, 3, smallidx
, sizesmall
, &tmpcoord
[k
]);
740 smallidx
+= is_smaller
;
744 smaller
= magicints
[smallidx
-1] / 2;
749 smallnum
= magicints
[smallidx
] / 2;
751 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
758 /* buf[0] holds the length in bytes */
759 if (xdr_int(xdrs
, &(buf
[0])) == 0)
770 rc
= errval
* (xdr_opaque(xdrs
, reinterpret_cast<char *>(&(buf
[3])), static_cast<unsigned int>(buf
[0])));
782 /* xdrs is open for reading */
784 if (xdr_int(xdrs
, &lsize
) == 0)
788 if (*size
!= 0 && lsize
!= *size
)
790 fprintf(stderr
, "wrong number of coordinates in xdr3dfcoord; "
791 "%d arg vs %d in file", *size
, lsize
);
798 return (xdr_vector(xdrs
, reinterpret_cast<char *>(fp
), static_cast<unsigned int>(size3
),
799 static_cast<unsigned int>(sizeof(*fp
)), (xdrproc_t
)xdr_float
));
801 if (xdr_float(xdrs
, precision
) == 0)
806 if (size3
<= prealloc_size
)
814 bufsize
= static_cast<int>(size3
* 1.2);
815 ip
= reinterpret_cast<int *>(malloc(size3
* sizeof(*ip
)));
816 buf
= reinterpret_cast<int *>(malloc(bufsize
* sizeof(*buf
)));
817 if (ip
== nullptr || buf
== nullptr)
819 fprintf(stderr
, "malloc failed\n");
824 buf
[0] = buf
[1] = buf
[2] = 0;
826 if ( (xdr_int(xdrs
, &(minint
[0])) == 0) ||
827 (xdr_int(xdrs
, &(minint
[1])) == 0) ||
828 (xdr_int(xdrs
, &(minint
[2])) == 0) ||
829 (xdr_int(xdrs
, &(maxint
[0])) == 0) ||
830 (xdr_int(xdrs
, &(maxint
[1])) == 0) ||
831 (xdr_int(xdrs
, &(maxint
[2])) == 0))
841 sizeint
[0] = maxint
[0] - minint
[0]+1;
842 sizeint
[1] = maxint
[1] - minint
[1]+1;
843 sizeint
[2] = maxint
[2] - minint
[2]+1;
845 /* check if one of the sizes is to big to be multiplied */
846 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff)
848 bitsizeint
[0] = sizeofint(sizeint
[0]);
849 bitsizeint
[1] = sizeofint(sizeint
[1]);
850 bitsizeint
[2] = sizeofint(sizeint
[2]);
851 bitsize
= 0; /* flag the use of large sizes */
855 bitsize
= sizeofints(3, sizeint
);
858 if (xdr_int(xdrs
, &smallidx
) == 0)
868 smaller
= magicints
[std::max(FIRSTIDX
, smallidx
-1)] / 2;
869 smallnum
= magicints
[smallidx
] / 2;
870 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
872 /* buf[0] holds the length in bytes */
874 if (xdr_int(xdrs
, &(buf
[0])) == 0)
885 if (xdr_opaque(xdrs
, reinterpret_cast<char *>(&(buf
[3])), static_cast<unsigned int>(buf
[0])) == 0)
897 buf
[0] = buf
[1] = buf
[2] = 0;
900 inv_precision
= 1.0 / *precision
;
906 thiscoord
= reinterpret_cast<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
))
1136 static float xtc_get_next_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
,
1145 if ((off
= gmx_ftell(fp
)) < 0)
1149 /* read one int just to make sure we dont read this frame but the next */
1150 xdr_int(xdrs
, &step
);
1153 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1157 if (gmx_fseek(fp
, off
, SEEK_SET
))
1166 if (gmx_fseek(fp
, off
, SEEK_SET
))
1177 xtc_get_current_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1185 if ((off
= gmx_ftell(fp
)) < 0)
1192 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1196 if (gmx_fseek(fp
, off
, SEEK_SET
))
1205 if (gmx_fseek(fp
, off
, SEEK_SET
))
1214 if (gmx_fseek(fp
, -2*XDR_INT_SIZE
, SEEK_CUR
))
1222 /* Currently not used, just for completeness */
1224 xtc_get_current_frame_number(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1232 if ((off
= gmx_ftell(fp
)) < 0)
1240 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1244 if (gmx_fseek(fp
, off
, SEEK_SET
))
1253 if (gmx_fseek(fp
, off
, SEEK_SET
))
1263 if (gmx_fseek(fp
, -2*XDR_INT_SIZE
, SEEK_CUR
))
1273 static gmx_off_t
xtc_get_next_frame_start(FILE *fp
, XDR
*xdrs
, int natoms
)
1279 /* read one int just to make sure we dont read this frame but the next */
1280 xdr_int(xdrs
, &step
);
1283 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1286 if ((res
= gmx_ftell(fp
)) >= 0)
1288 return res
- XDR_INT_SIZE
;
1305 xdr_xtc_estimate_dt(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1312 if ((off
= gmx_ftell(fp
)) < 0)
1317 tinit
= xtc_get_current_frame_time(fp
, xdrs
, natoms
, bOK
);
1324 res
= xtc_get_next_frame_time(fp
, xdrs
, natoms
, bOK
);
1332 if (0 != gmx_fseek(fp
, off
, SEEK_SET
))
1341 xdr_xtc_seek_frame(int frame
, FILE *fp
, XDR
*xdrs
, int natoms
)
1344 gmx_off_t high
, pos
;
1347 /* round to 4 bytes */
1350 if (gmx_fseek(fp
, 0, SEEK_END
))
1355 if ((high
= gmx_ftell(fp
)) < 0)
1360 /* round to 4 bytes */
1361 high
/= XDR_INT_SIZE
;
1362 high
*= XDR_INT_SIZE
;
1363 offset
= ((high
/2)/XDR_INT_SIZE
)*XDR_INT_SIZE
;
1365 if (gmx_fseek(fp
, offset
, SEEK_SET
))
1372 fr
= xtc_get_next_frame_number(fp
, xdrs
, natoms
);
1377 if (fr
!= frame
&& llabs(low
-high
) > header_size
)
1387 /* round to 4 bytes */
1388 offset
= (((high
+low
)/2)/4)*4;
1390 if (gmx_fseek(fp
, offset
, SEEK_SET
))
1400 if (offset
<= header_size
)
1405 if (gmx_fseek(fp
, offset
, SEEK_SET
))
1410 if ((pos
= xtc_get_next_frame_start(fp
, xdrs
, natoms
)) < 0)
1412 /* we probably hit an end of file */
1416 if (gmx_fseek(fp
, pos
, SEEK_SET
))
1426 int xdr_xtc_seek_time(real time
, FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool bSeekForwardOnly
)
1430 gmx_bool bOK
= FALSE
;
1432 gmx_off_t high
, offset
, pos
;
1435 if (bSeekForwardOnly
)
1437 low
= gmx_ftell(fp
)-header_size
;
1439 if (gmx_fseek(fp
, 0, SEEK_END
))
1444 if ((high
= gmx_ftell(fp
)) < 0)
1449 high
/= XDR_INT_SIZE
;
1450 high
*= XDR_INT_SIZE
;
1451 offset
= (((high
-low
) / 2) / XDR_INT_SIZE
) * XDR_INT_SIZE
;
1453 if (gmx_fseek(fp
, offset
, SEEK_SET
))
1460 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1461 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1471 dt
= xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
);
1482 /* Found a place in the trajectory that has positive time step while
1483 other has negative time step */
1492 /* Found a place in the trajectory that has positive time step while
1493 other has negative time step */
1499 t
= xtc_get_next_frame_time(fp
, xdrs
, natoms
, &bOK
);
1505 /* If we are before the target time and the time step is positive or 0, or we have
1506 after the target time and the time step is negative, or the difference between
1507 the current time and the target time is bigger than dt and above all the distance between high
1508 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1509 if we reached the solution */
1510 if ((((t
< time
&& dt_sign
>= 0) || (t
> time
&& dt_sign
== -1)) ||
1511 ((t
- time
) >= dt
&& dt_sign
>= 0) || ((time
- t
) >= -dt
&& dt_sign
< 0)) &&
1512 (llabs(low
- high
) > header_size
))
1514 if (dt
>= 0 && dt_sign
!= -1)
1525 else if (dt
<= 0 && dt_sign
== -1)
1538 /* We should never reach here */
1541 /* round to 4 bytes and subtract header*/
1542 offset
= (((high
+ low
) / 2) / XDR_INT_SIZE
) * XDR_INT_SIZE
;
1543 if (gmx_fseek(fp
, offset
, SEEK_SET
))
1550 if (llabs(low
- high
) <= header_size
)
1554 /* re-estimate dt */
1555 if (xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
) != dt
)
1559 dt
= xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
);
1562 if (t
>= time
&& t
- time
< dt
)
1569 if (offset
<= header_size
)
1574 gmx_fseek(fp
, offset
, SEEK_SET
);
1576 if ((pos
= xtc_get_next_frame_start(fp
, xdrs
, natoms
)) < 0)
1581 if (gmx_fseek(fp
, pos
, SEEK_SET
))
1589 xdr_xtc_get_last_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1594 off
= gmx_ftell(fp
);
1601 if (gmx_fseek(fp
, -3*XDR_INT_SIZE
, SEEK_END
) != 0)
1607 time
= xtc_get_current_frame_time(fp
, xdrs
, natoms
, bOK
);
1613 if (gmx_fseek(fp
, off
, SEEK_SET
) != 0)
1623 xdr_xtc_get_last_frame_number(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1629 if ((off
= gmx_ftell(fp
)) < 0)
1636 if (gmx_fseek(fp
, -3*XDR_INT_SIZE
, SEEK_END
))
1642 frame
= xtc_get_current_frame_number(fp
, xdrs
, natoms
, bOK
);
1649 if (gmx_fseek(fp
, off
, SEEK_SET
))