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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gmx_fatal.h"
56 # define gmx_fseek(A,B,C) fseeko(A,B,C)
57 # define gmx_ftell(A) ftello(A)
58 # define gmx_off_t off_t
60 # define gmx_fseek(A,B,C) fseek(A,B,C)
61 # define gmx_ftell(A) ftell(A)
62 # define gmx_off_t int
67 /* This is just for clarity - it can never be anything but 4! */
68 #define XDR_INT_SIZE 4
70 /* same order as the definition of xdr_datatype */
71 const char *xdr_datatype_names
[] =
82 /*___________________________________________________________________________
84 | what follows are the C routine to read/write compressed coordinates together
85 | with some routines to assist in this task (those are marked
86 | static and cannot be called from user programs)
88 #define MAXABS INT_MAX-2
91 #define MIN(x,y) ((x) < (y) ? (x):(y))
94 #define MAX(x,y) ((x) > (y) ? (x):(y))
97 #define SQR(x) ((x)*(x))
99 static const int magicints
[] = {
100 0, 0, 0, 0, 0, 0, 0, 0, 0,
101 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
102 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
103 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
104 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
105 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
106 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
107 8388607, 10568983, 13316085, 16777216 };
110 /* note that magicints[FIRSTIDX-1] == 0 */
111 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
114 /*____________________________________________________________________________
116 | sendbits - encode num into buf using the specified number of bits
118 | This routines appends the value of num to the bits already present in
119 | the array buf. You need to give it the number of bits to use and you
120 | better make sure that this number of bits is enough to hold the value
121 | Also num must be positive.
125 static void sendbits(int buf
[], int num_of_bits
, int num
) {
127 unsigned int cnt
, lastbyte
;
129 unsigned char * cbuf
;
131 cbuf
= ((unsigned char *)buf
) + 3 * sizeof(*buf
);
132 cnt
= (unsigned int) buf
[0];
134 lastbyte
=(unsigned int) buf
[2];
135 while (num_of_bits
>= 8) {
136 lastbyte
= (lastbyte
<< 8) | ((num
>> (num_of_bits
-8)) /* & 0xff*/);
137 cbuf
[cnt
++] = lastbyte
>> lastbits
;
140 if (num_of_bits
> 0) {
141 lastbyte
= (lastbyte
<< num_of_bits
) | num
;
142 lastbits
+= num_of_bits
;
145 cbuf
[cnt
++] = lastbyte
>> lastbits
;
152 cbuf
[cnt
] = lastbyte
<< (8 - lastbits
);
156 /*_________________________________________________________________________
158 | sizeofint - calculate bitsize of an integer
160 | return the number of bits needed to store an integer with given max size
164 static int sizeofint(const int size
) {
168 while (size
>= num
&& num_of_bits
< 32) {
175 /*___________________________________________________________________________
177 | sizeofints - calculate 'bitsize' of compressed ints
179 | given the number of small unsigned integers and the maximum value
180 | return the number of bits needed to read or write them with the
181 | routines receiveints and sendints. You need this parameter when
182 | calling these routines. Note that for many calls I can use
183 | the variable 'smallidx' which is exactly the number of bits, and
184 | So I don't need to call 'sizeofints for those calls.
187 static int sizeofints( const int num_of_ints
, unsigned int sizes
[]) {
190 unsigned int num_of_bytes
, num_of_bits
, bytecnt
, tmp
;
194 for (i
=0; i
< num_of_ints
; i
++) {
196 for (bytecnt
= 0; bytecnt
< num_of_bytes
; bytecnt
++) {
197 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
198 bytes
[bytecnt
] = tmp
& 0xff;
202 bytes
[bytecnt
++] = tmp
& 0xff;
205 num_of_bytes
= bytecnt
;
209 while (bytes
[num_of_bytes
] >= num
) {
213 return num_of_bits
+ num_of_bytes
* 8;
217 /*____________________________________________________________________________
219 | sendints - send a small set of small integers in compressed format
221 | this routine is used internally by xdr3dfcoord, to send a set of
222 | small integers to the buffer.
223 | Multiplication with fixed (specified maximum ) sizes is used to get
224 | to one big, multibyte integer. Allthough the routine could be
225 | modified to handle sizes bigger than 16777216, or more than just
226 | a few integers, this is not done, because the gain in compression
227 | isn't worth the effort. Note that overflowing the multiplication
228 | or the byte buffer (32 bytes) is unchecked and causes bad results.
232 static void sendints(int buf
[], const int num_of_ints
, const int num_of_bits
,
233 unsigned int sizes
[], unsigned int nums
[]) {
235 int i
, num_of_bytes
, bytecnt
;
236 unsigned int bytes
[32], tmp
;
241 bytes
[num_of_bytes
++] = tmp
& 0xff;
245 for (i
= 1; i
< num_of_ints
; i
++) {
246 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
++) {
254 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
255 bytes
[bytecnt
] = tmp
& 0xff;
259 bytes
[bytecnt
++] = tmp
& 0xff;
262 num_of_bytes
= bytecnt
;
264 if (num_of_bits
>= num_of_bytes
* 8) {
265 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);
270 for (i
= 0; i
< num_of_bytes
-1; i
++) {
271 sendbits(buf
, 8, bytes
[i
]);
273 sendbits(buf
, num_of_bits
- (num_of_bytes
-1) * 8, bytes
[i
]);
278 /*___________________________________________________________________________
280 | receivebits - decode number from buf using specified number of bits
282 | extract the number of bits from the array buf and construct an integer
283 | from it. Return that value.
287 static int receivebits(int buf
[], int num_of_bits
) {
289 int cnt
, num
, lastbits
;
290 unsigned int lastbyte
;
291 unsigned char * cbuf
;
292 int mask
= (1 << num_of_bits
) -1;
294 cbuf
= ((unsigned char *)buf
) + 3 * sizeof(*buf
);
296 lastbits
= (unsigned int) buf
[1];
297 lastbyte
= (unsigned int) buf
[2];
300 while (num_of_bits
>= 8) {
301 lastbyte
= ( lastbyte
<< 8 ) | cbuf
[cnt
++];
302 num
|= (lastbyte
>> lastbits
) << (num_of_bits
- 8);
305 if (num_of_bits
> 0) {
306 if (lastbits
< num_of_bits
) {
308 lastbyte
= (lastbyte
<< 8) | cbuf
[cnt
++];
310 lastbits
-= num_of_bits
;
311 num
|= (lastbyte
>> lastbits
) & ((1 << num_of_bits
) -1);
320 /*____________________________________________________________________________
322 | receiveints - decode 'small' integers from the buf array
324 | this routine is the inverse from sendints() and decodes the small integers
325 | written to buf by calculating the remainder and doing divisions with
326 | the given sizes[]. You need to specify the total number of bits to be
327 | used from buf in num_of_bits.
331 static void receiveints(int buf
[], const int num_of_ints
, int num_of_bits
,
332 unsigned int sizes
[], int nums
[]) {
334 int i
, j
, num_of_bytes
, p
, num
;
336 bytes
[0] = bytes
[1] = bytes
[2] = bytes
[3] = 0;
338 while (num_of_bits
> 8) {
339 bytes
[num_of_bytes
++] = receivebits(buf
, 8);
342 if (num_of_bits
> 0) {
343 bytes
[num_of_bytes
++] = receivebits(buf
, num_of_bits
);
345 for (i
= num_of_ints
-1; i
> 0; i
--) {
347 for (j
= num_of_bytes
-1; j
>=0; j
--) {
348 num
= (num
<< 8) | bytes
[j
];
351 num
= num
- p
* sizes
[i
];
355 nums
[0] = bytes
[0] | (bytes
[1] << 8) | (bytes
[2] << 16) | (bytes
[3] << 24);
358 /*____________________________________________________________________________
360 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
362 | this routine reads or writes (depending on how you opened the file with
363 | xdropen() ) a large number of 3d coordinates (stored in *fp).
364 | The number of coordinates triplets to write is given by *size. On
365 | read this number may be zero, in which case it reads as many as were written
366 | or it may specify the number if triplets to read (which should match the
368 | Compression is achieved by first converting all floating numbers to integer
369 | using multiplication by *precision and rounding to the nearest integer.
370 | Then the minimum and maximum value are calculated to determine the range.
371 | The limited range of integers so found, is used to compress the coordinates.
372 | In addition the differences between succesive coordinates is calculated.
373 | If the difference happens to be 'small' then only the difference is saved,
374 | compressing the data even more. The notion of 'small' is changed dynamically
375 | and is enlarged or reduced whenever needed or possible.
376 | Extra compression is achieved in the case of GROMOS and coordinates of
377 | water molecules. GROMOS first writes out the Oxygen position, followed by
378 | the two hydrogens. In order to make the differences smaller (and thereby
379 | compression the data better) the order is changed into first one hydrogen
380 | then the oxygen, followed by the other hydrogen. This is rather special, but
381 | it shouldn't harm in the general case.
385 int xdr3dfcoord(XDR
*xdrs
, float *fp
, int *size
, float *precision
)
391 /* preallocate a small buffer and ip on the stack - if we need more
392 we can always malloc(). This is faster for small values of size: */
393 unsigned prealloc_size
=3*16;
394 int prealloc_ip
[3*16], prealloc_buf
[3*20];
395 int we_should_free
=0;
397 int minint
[3], maxint
[3], mindiff
, *lip
, diff
;
398 int lint1
, lint2
, lint3
, oldlint1
, oldlint2
, oldlint3
, smallidx
;
400 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
, *luip
;
402 int smallnum
, smaller
, larger
, i
, is_small
, is_smaller
, run
, prevrun
;
404 int tmp
, *thiscoord
, prevcoord
[3];
405 unsigned int tmpcoord
[30];
407 int bufsize
, xdrid
, lsize
;
408 unsigned int bitsize
;
413 bRead
= (xdrs
->x_op
== XDR_DECODE
);
414 bitsizeint
[0] = bitsizeint
[1] = bitsizeint
[2] = 0;
415 prevcoord
[0] = prevcoord
[1] = prevcoord
[2] = 0;
419 /* xdrs is open for writing */
421 if (xdr_int(xdrs
, size
) == 0)
424 /* when the number of coordinates is small, don't try to compress; just
425 * write them as floats using xdr_vector
428 return (xdr_vector(xdrs
, (char *) fp
, (unsigned int)size3
,
429 (unsigned int)sizeof(*fp
), (xdrproc_t
)xdr_float
));
432 if(xdr_float(xdrs
, precision
) == 0)
435 if (size3
<= prealloc_size
)
443 bufsize
= size3
* 1.2;
444 ip
= (int *)malloc((size_t)(size3
* sizeof(*ip
)));
445 buf
= (int *)malloc((size_t)(bufsize
* sizeof(*buf
)));
446 if (ip
== NULL
|| buf
==NULL
)
448 fprintf(stderr
,"malloc failed\n");
452 /* buf[0-2] are special and do not contain actual data */
453 buf
[0] = buf
[1] = buf
[2] = 0;
454 minint
[0] = minint
[1] = minint
[2] = INT_MAX
;
455 maxint
[0] = maxint
[1] = maxint
[2] = INT_MIN
;
460 oldlint1
= oldlint2
= oldlint3
= 0;
461 while(lfp
< fp
+ size3
) {
462 /* find nearest integer */
464 lf
= *lfp
* *precision
+ 0.5;
466 lf
= *lfp
* *precision
- 0.5;
467 if (fabs(lf
) > MAXABS
) {
468 /* scaling would cause overflow */
472 if (lint1
< minint
[0]) minint
[0] = lint1
;
473 if (lint1
> maxint
[0]) maxint
[0] = lint1
;
477 lf
= *lfp
* *precision
+ 0.5;
479 lf
= *lfp
* *precision
- 0.5;
480 if (fabs(lf
) > MAXABS
) {
481 /* scaling would cause overflow */
485 if (lint2
< minint
[1]) minint
[1] = lint2
;
486 if (lint2
> maxint
[1]) maxint
[1] = lint2
;
490 lf
= *lfp
* *precision
+ 0.5;
492 lf
= *lfp
* *precision
- 0.5;
493 if (fabs(lf
) > MAXABS
) {
494 /* scaling would cause overflow */
498 if (lint3
< minint
[2]) minint
[2] = lint3
;
499 if (lint3
> maxint
[2]) maxint
[2] = lint3
;
502 diff
= abs(oldlint1
-lint1
)+abs(oldlint2
-lint2
)+abs(oldlint3
-lint3
);
503 if (diff
< mindiff
&& lfp
> fp
+ 3)
509 if ( (xdr_int(xdrs
, &(minint
[0])) == 0) ||
510 (xdr_int(xdrs
, &(minint
[1])) == 0) ||
511 (xdr_int(xdrs
, &(minint
[2])) == 0) ||
512 (xdr_int(xdrs
, &(maxint
[0])) == 0) ||
513 (xdr_int(xdrs
, &(maxint
[1])) == 0) ||
514 (xdr_int(xdrs
, &(maxint
[2])) == 0))
524 if ((float)maxint
[0] - (float)minint
[0] >= MAXABS
||
525 (float)maxint
[1] - (float)minint
[1] >= MAXABS
||
526 (float)maxint
[2] - (float)minint
[2] >= MAXABS
) {
527 /* turning value in unsigned by subtracting minint
528 * would cause overflow
532 sizeint
[0] = maxint
[0] - minint
[0]+1;
533 sizeint
[1] = maxint
[1] - minint
[1]+1;
534 sizeint
[2] = maxint
[2] - minint
[2]+1;
536 /* check if one of the sizes is to big to be multiplied */
537 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff) {
538 bitsizeint
[0] = sizeofint(sizeint
[0]);
539 bitsizeint
[1] = sizeofint(sizeint
[1]);
540 bitsizeint
[2] = sizeofint(sizeint
[2]);
541 bitsize
= 0; /* flag the use of large sizes */
543 bitsize
= sizeofints(3, sizeint
);
546 luip
= (unsigned int *) ip
;
548 while (smallidx
< LASTIDX
&& magicints
[smallidx
] < mindiff
) {
551 if(xdr_int(xdrs
, &smallidx
) == 0)
561 maxidx
= MIN(LASTIDX
, smallidx
+ 8) ;
562 minidx
= maxidx
- 8; /* often this equal smallidx */
563 smaller
= magicints
[MAX(FIRSTIDX
, smallidx
-1)] / 2;
564 smallnum
= magicints
[smallidx
] / 2;
565 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
566 larger
= magicints
[maxidx
] / 2;
570 thiscoord
= (int *)(luip
) + i
* 3;
571 if (smallidx
< maxidx
&& i
>= 1 &&
572 abs(thiscoord
[0] - prevcoord
[0]) < larger
&&
573 abs(thiscoord
[1] - prevcoord
[1]) < larger
&&
574 abs(thiscoord
[2] - prevcoord
[2]) < larger
) {
576 } else if (smallidx
> minidx
) {
582 if (abs(thiscoord
[0] - thiscoord
[3]) < smallnum
&&
583 abs(thiscoord
[1] - thiscoord
[4]) < smallnum
&&
584 abs(thiscoord
[2] - thiscoord
[5]) < smallnum
) {
585 /* interchange first with second atom for better
586 * compression of water molecules
588 tmp
= thiscoord
[0]; thiscoord
[0] = thiscoord
[3];
590 tmp
= thiscoord
[1]; thiscoord
[1] = thiscoord
[4];
592 tmp
= thiscoord
[2]; thiscoord
[2] = thiscoord
[5];
598 tmpcoord
[0] = thiscoord
[0] - minint
[0];
599 tmpcoord
[1] = thiscoord
[1] - minint
[1];
600 tmpcoord
[2] = thiscoord
[2] - minint
[2];
602 sendbits(buf
, bitsizeint
[0], tmpcoord
[0]);
603 sendbits(buf
, bitsizeint
[1], tmpcoord
[1]);
604 sendbits(buf
, bitsizeint
[2], tmpcoord
[2]);
606 sendints(buf
, 3, bitsize
, sizeint
, tmpcoord
);
608 prevcoord
[0] = thiscoord
[0];
609 prevcoord
[1] = thiscoord
[1];
610 prevcoord
[2] = thiscoord
[2];
611 thiscoord
= thiscoord
+ 3;
615 if (is_small
== 0 && is_smaller
== -1)
617 while (is_small
&& run
< 8*3) {
618 if (is_smaller
== -1 && (
619 SQR(thiscoord
[0] - prevcoord
[0]) +
620 SQR(thiscoord
[1] - prevcoord
[1]) +
621 SQR(thiscoord
[2] - prevcoord
[2]) >= smaller
* smaller
)) {
625 tmpcoord
[run
++] = thiscoord
[0] - prevcoord
[0] + smallnum
;
626 tmpcoord
[run
++] = thiscoord
[1] - prevcoord
[1] + smallnum
;
627 tmpcoord
[run
++] = thiscoord
[2] - prevcoord
[2] + smallnum
;
629 prevcoord
[0] = thiscoord
[0];
630 prevcoord
[1] = thiscoord
[1];
631 prevcoord
[2] = thiscoord
[2];
634 thiscoord
= thiscoord
+ 3;
637 abs(thiscoord
[0] - prevcoord
[0]) < smallnum
&&
638 abs(thiscoord
[1] - prevcoord
[1]) < smallnum
&&
639 abs(thiscoord
[2] - prevcoord
[2]) < smallnum
) {
643 if (run
!= prevrun
|| is_smaller
!= 0) {
645 sendbits(buf
, 1, 1); /* flag the change in run-length */
646 sendbits(buf
, 5, run
+is_smaller
+1);
648 sendbits(buf
, 1, 0); /* flag the fact that runlength did not change */
650 for (k
=0; k
< run
; k
+=3) {
651 sendints(buf
, 3, smallidx
, sizesmall
, &tmpcoord
[k
]);
653 if (is_smaller
!= 0) {
654 smallidx
+= is_smaller
;
655 if (is_smaller
< 0) {
657 smaller
= magicints
[smallidx
-1] / 2;
660 smallnum
= magicints
[smallidx
] / 2;
662 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
665 if (buf
[1] != 0) buf
[0]++;
666 /* buf[0] holds the length in bytes */
667 if(xdr_int(xdrs
, &(buf
[0])) == 0)
678 rc
=errval
* (xdr_opaque(xdrs
, (char *)&(buf
[3]), (unsigned int)buf
[0]));
688 /* xdrs is open for reading */
690 if (xdr_int(xdrs
, &lsize
) == 0)
692 if (*size
!= 0 && lsize
!= *size
) {
693 fprintf(stderr
, "wrong number of coordinates in xdr3dfcoord; "
694 "%d arg vs %d in file", *size
, lsize
);
700 return (xdr_vector(xdrs
, (char *) fp
, (unsigned int)size3
,
701 (unsigned int)sizeof(*fp
), (xdrproc_t
)xdr_float
));
703 if(xdr_float(xdrs
, precision
) == 0)
706 if (size3
<= prealloc_size
)
714 bufsize
= size3
* 1.2;
715 ip
= (int *)malloc((size_t)(size3
* sizeof(*ip
)));
716 buf
= (int *)malloc((size_t)(bufsize
* sizeof(*buf
)));
717 if (ip
== NULL
|| buf
==NULL
)
719 fprintf(stderr
,"malloc failed\n");
724 buf
[0] = buf
[1] = buf
[2] = 0;
726 if ( (xdr_int(xdrs
, &(minint
[0])) == 0) ||
727 (xdr_int(xdrs
, &(minint
[1])) == 0) ||
728 (xdr_int(xdrs
, &(minint
[2])) == 0) ||
729 (xdr_int(xdrs
, &(maxint
[0])) == 0) ||
730 (xdr_int(xdrs
, &(maxint
[1])) == 0) ||
731 (xdr_int(xdrs
, &(maxint
[2])) == 0))
741 sizeint
[0] = maxint
[0] - minint
[0]+1;
742 sizeint
[1] = maxint
[1] - minint
[1]+1;
743 sizeint
[2] = maxint
[2] - minint
[2]+1;
745 /* check if one of the sizes is to big to be multiplied */
746 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff) {
747 bitsizeint
[0] = sizeofint(sizeint
[0]);
748 bitsizeint
[1] = sizeofint(sizeint
[1]);
749 bitsizeint
[2] = sizeofint(sizeint
[2]);
750 bitsize
= 0; /* flag the use of large sizes */
752 bitsize
= sizeofints(3, sizeint
);
755 if (xdr_int(xdrs
, &smallidx
) == 0)
765 maxidx
= MIN(LASTIDX
, smallidx
+ 8) ;
766 minidx
= maxidx
- 8; /* often this equal smallidx */
767 smaller
= magicints
[MAX(FIRSTIDX
, smallidx
-1)] / 2;
768 smallnum
= magicints
[smallidx
] / 2;
769 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
770 larger
= magicints
[maxidx
];
772 /* buf[0] holds the length in bytes */
774 if (xdr_int(xdrs
, &(buf
[0])) == 0)
785 if (xdr_opaque(xdrs
, (char *)&(buf
[3]), (unsigned int)buf
[0]) == 0)
797 buf
[0] = buf
[1] = buf
[2] = 0;
800 inv_precision
= 1.0 / * precision
;
804 while ( i
< lsize
) {
805 thiscoord
= (int *)(lip
) + i
* 3;
808 thiscoord
[0] = receivebits(buf
, bitsizeint
[0]);
809 thiscoord
[1] = receivebits(buf
, bitsizeint
[1]);
810 thiscoord
[2] = receivebits(buf
, bitsizeint
[2]);
812 receiveints(buf
, 3, bitsize
, sizeint
, thiscoord
);
816 thiscoord
[0] += minint
[0];
817 thiscoord
[1] += minint
[1];
818 thiscoord
[2] += minint
[2];
820 prevcoord
[0] = thiscoord
[0];
821 prevcoord
[1] = thiscoord
[1];
822 prevcoord
[2] = thiscoord
[2];
825 flag
= receivebits(buf
, 1);
828 run
= receivebits(buf
, 5);
829 is_smaller
= run
% 3;
835 for (k
= 0; k
< run
; k
+=3) {
836 receiveints(buf
, 3, smallidx
, sizesmall
, thiscoord
);
838 thiscoord
[0] += prevcoord
[0] - smallnum
;
839 thiscoord
[1] += prevcoord
[1] - smallnum
;
840 thiscoord
[2] += prevcoord
[2] - smallnum
;
842 /* interchange first with second atom for better
843 * compression of water molecules
845 tmp
= thiscoord
[0]; thiscoord
[0] = prevcoord
[0];
847 tmp
= thiscoord
[1]; thiscoord
[1] = prevcoord
[1];
849 tmp
= thiscoord
[2]; thiscoord
[2] = prevcoord
[2];
851 *lfp
++ = prevcoord
[0] * inv_precision
;
852 *lfp
++ = prevcoord
[1] * inv_precision
;
853 *lfp
++ = prevcoord
[2] * inv_precision
;
855 prevcoord
[0] = thiscoord
[0];
856 prevcoord
[1] = thiscoord
[1];
857 prevcoord
[2] = thiscoord
[2];
859 *lfp
++ = thiscoord
[0] * inv_precision
;
860 *lfp
++ = thiscoord
[1] * inv_precision
;
861 *lfp
++ = thiscoord
[2] * inv_precision
;
864 *lfp
++ = thiscoord
[0] * inv_precision
;
865 *lfp
++ = thiscoord
[1] * inv_precision
;
866 *lfp
++ = thiscoord
[2] * inv_precision
;
868 smallidx
+= is_smaller
;
869 if (is_smaller
< 0) {
871 if (smallidx
> FIRSTIDX
) {
872 smaller
= magicints
[smallidx
- 1] /2;
876 } else if (is_smaller
> 0) {
878 smallnum
= magicints
[smallidx
] / 2;
880 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
893 /******************************************************************
895 XTC files have a relatively simple structure.
896 They have a header of 16 bytes and the rest are
897 the compressed coordinates of the files. Due to the
898 compression 00 is not present in the coordinates.
899 The first 4 bytes of the header are the magic number
900 1995 (0x000007CB). If we find this number we are guaranteed
901 to be in the header, due to the presence of so many zeros.
902 The second 4 bytes are the number of atoms in the frame, and is
903 assumed to be constant. The third 4 bytes are the frame number.
904 The last 4 bytes are a floating point representation of the time.
906 ********************************************************************/
908 /* Must match definition in xtcio.c */
910 #define XTC_MAGIC 1995
913 static const int header_size
= 16;
915 /* Check if we are at the header start.
916 At the same time it will also read 1 int */
917 static int xtc_at_header_start(FILE *fp
, XDR
*xdrs
,
918 int natoms
, int * timestep
, float * time
){
925 if((off
= gmx_ftell(fp
)) < 0){
928 /* read magic natoms and timestep */
930 if(!xdr_int(xdrs
, &(i_inp
[i
]))){
931 gmx_fseek(fp
,off
+XDR_INT_SIZE
,SEEK_SET
);
936 if(i_inp
[0] != XTC_MAGIC
){
937 if(gmx_fseek(fp
,off
+XDR_INT_SIZE
,SEEK_SET
)){
942 /* read time and box */
944 if(!xdr_float(xdrs
, &(f_inp
[i
]))){
945 gmx_fseek(fp
,off
+XDR_INT_SIZE
,SEEK_SET
);
949 /* Make a rigourous check to see if we are in the beggining of a header
950 Hopefully there are no ambiguous cases */
951 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
952 right triangle and that the first element must be nonzero unless the entire matrix is zero
954 if(i_inp
[1] == natoms
&&
955 ((f_inp
[1] != 0 && f_inp
[6] == 0) ||
956 (f_inp
[1] == 0 && f_inp
[5] == 0 && f_inp
[9] == 0))){
957 if(gmx_fseek(fp
,off
+XDR_INT_SIZE
,SEEK_SET
)){
961 *timestep
= i_inp
[2];
964 if(gmx_fseek(fp
,off
+XDR_INT_SIZE
,SEEK_SET
)){
971 xtc_get_next_frame_number(FILE *fp
, XDR
*xdrs
, int natoms
)
978 if((off
= gmx_ftell(fp
)) < 0){
982 /* read one int just to make sure we dont read this frame but the next */
985 ret
= xtc_at_header_start(fp
,xdrs
,natoms
,&step
,&time
);
987 if(gmx_fseek(fp
,off
,SEEK_SET
)){
992 if(gmx_fseek(fp
,off
,SEEK_SET
)){
1001 static float xtc_get_next_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
,
1010 if ((off
= gmx_ftell(fp
)) < 0)
1014 /* read one int just to make sure we dont read this frame but the next */
1015 xdr_int(xdrs
, &step
);
1018 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1022 if (gmx_fseek(fp
,off
,SEEK_SET
))
1031 if (gmx_fseek(fp
,off
,SEEK_SET
))
1043 xtc_get_current_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1051 if ((off
= gmx_ftell(fp
)) < 0)
1058 ret
= xtc_at_header_start(fp
, xdrs
, natoms
, &step
, &time
);
1062 if (gmx_fseek(fp
,off
,SEEK_SET
))
1071 if (gmx_fseek(fp
,off
,SEEK_SET
))
1080 if (gmx_fseek(fp
,-2*XDR_INT_SIZE
,SEEK_CUR
))
1089 /* Currently not used, just for completeness */
1091 xtc_get_current_frame_number(FILE *fp
,XDR
*xdrs
,int natoms
, gmx_bool
* bOK
)
1099 if((off
= gmx_ftell(fp
)) < 0){
1105 ret
= xtc_at_header_start(fp
,xdrs
,natoms
,&step
,&time
);
1108 if(gmx_fseek(fp
,off
,SEEK_SET
)){
1113 }else if(ret
== -1){
1114 if(gmx_fseek(fp
,off
,SEEK_SET
)){
1121 if(gmx_fseek(fp
,-2*XDR_INT_SIZE
,SEEK_CUR
)){
1131 static gmx_off_t
xtc_get_next_frame_start(FILE *fp
, XDR
*xdrs
, int natoms
)
1138 /* read one int just to make sure we dont read this frame but the next */
1139 xdr_int(xdrs
,&step
);
1142 ret
= xtc_at_header_start(fp
,xdrs
,natoms
,&step
,&time
);
1144 if((res
= gmx_ftell(fp
)) >= 0){
1145 return res
- XDR_INT_SIZE
;
1149 }else if(ret
== -1){
1159 xdr_xtc_estimate_dt(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1166 if((off
= gmx_ftell(fp
)) < 0){
1170 tinit
= xtc_get_current_frame_time(fp
,xdrs
,natoms
,bOK
);
1177 res
= xtc_get_next_frame_time(fp
,xdrs
,natoms
,bOK
);
1185 if (0 != gmx_fseek(fp
,off
,SEEK_SET
)) {
1194 xdr_xtc_seek_frame(int frame
, FILE *fp
, XDR
*xdrs
, int natoms
)
1200 /* round to 4 bytes */
1203 if(gmx_fseek(fp
,0,SEEK_END
)){
1207 if((high
= gmx_ftell(fp
)) < 0){
1211 /* round to 4 bytes */
1212 high
/= XDR_INT_SIZE
;
1213 high
*= XDR_INT_SIZE
;
1214 offset
= ((high
/2)/XDR_INT_SIZE
)*XDR_INT_SIZE
;
1216 if(gmx_fseek(fp
,offset
,SEEK_SET
)){
1222 fr
= xtc_get_next_frame_number(fp
,xdrs
,natoms
);
1227 if(fr
!= frame
&& abs(low
-high
) > header_size
)
1237 /* round to 4 bytes */
1238 offset
= (((high
+low
)/2)/4)*4;
1240 if(gmx_fseek(fp
,offset
,SEEK_SET
)){
1249 if(offset
<= header_size
)
1254 if(gmx_fseek(fp
,offset
,SEEK_SET
)){
1258 if((pos
= xtc_get_next_frame_start(fp
,xdrs
,natoms
))< 0){
1259 /* we probably hit an end of file */
1263 if(gmx_fseek(fp
,pos
,SEEK_SET
)){
1272 int xdr_xtc_seek_time(real time
, FILE *fp
, XDR
*xdrs
, int natoms
,gmx_bool bSeekForwardOnly
)
1276 gmx_bool bOK
= FALSE
;
1278 gmx_off_t high
, offset
, pos
;
1282 if (bSeekForwardOnly
)
1284 low
= gmx_ftell(fp
);
1286 if (gmx_fseek(fp
,0,SEEK_END
))
1291 if ((high
= gmx_ftell(fp
)) < 0)
1296 high
/= XDR_INT_SIZE
;
1297 high
*= XDR_INT_SIZE
;
1298 offset
= (((high
-low
) / 2) / XDR_INT_SIZE
) * XDR_INT_SIZE
;
1300 if (gmx_fseek(fp
,offset
,SEEK_SET
))
1307 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1308 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1318 dt
= xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
);
1329 /* Found a place in the trajectory that has positive time step while
1330 other has negative time step */
1339 /* Found a place in the trajectory that has positive time step while
1340 other has negative time step */
1346 t
= xtc_get_next_frame_time(fp
, xdrs
, natoms
, &bOK
);
1352 /* If we are before the target time and the time step is positive or 0, or we have
1353 after the target time and the time step is negative, or the difference between
1354 the current time and the target time is bigger than dt and above all the distance between high
1355 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1356 if we reached the solution */
1357 if ((((t
< time
&& dt_sign
>= 0) || (t
> time
&& dt_sign
== -1)) || ((t
1358 - time
) >= dt
&& dt_sign
>= 0)
1359 || ((time
- t
) >= -dt
&& dt_sign
< 0)) && (abs(low
- high
)
1362 if (dt
>= 0 && dt_sign
!= -1)
1373 else if (dt
<= 0 && dt_sign
== -1)
1386 /* We should never reach here */
1389 /* round to 4 bytes and subtract header*/
1390 offset
= (((high
+ low
) / 2) / XDR_INT_SIZE
) * XDR_INT_SIZE
;
1391 if (gmx_fseek(fp
,offset
,SEEK_SET
))
1398 if (abs(low
- high
) <= header_size
)
1402 /* re-estimate dt */
1403 if (xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
) != dt
)
1407 dt
= xdr_xtc_estimate_dt(fp
, xdrs
, natoms
, &bOK
);
1410 if (t
>= time
&& t
- time
< dt
)
1417 if (offset
<= header_size
)
1422 gmx_fseek(fp
,offset
,SEEK_SET
);
1424 if ((pos
= xtc_get_next_frame_start(fp
, xdrs
, natoms
)) < 0)
1429 if (gmx_fseek(fp
,pos
,SEEK_SET
))
1437 xdr_xtc_get_last_frame_time(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1443 off
= gmx_ftell(fp
);
1449 if( (res
= gmx_fseek(fp
,-3*XDR_INT_SIZE
,SEEK_END
)) != 0){
1454 time
= xtc_get_current_frame_time(fp
, xdrs
, natoms
, bOK
);
1459 if( (res
= gmx_fseek(fp
,off
,SEEK_SET
)) != 0){
1468 xdr_xtc_get_last_frame_number(FILE *fp
, XDR
*xdrs
, int natoms
, gmx_bool
* bOK
)
1475 if((off
= gmx_ftell(fp
)) < 0){
1481 if(gmx_fseek(fp
,-3*XDR_INT_SIZE
,SEEK_END
)){
1486 frame
= xtc_get_current_frame_number(fp
, xdrs
, natoms
, bOK
);
1492 if(gmx_fseek(fp
,off
,SEEK_SET
)){