LOOS  v2.3.2
xtcwriter.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2014, Tod D. Romo, Alan Grossfield
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 #include <xtcwriter.hpp>
24 #include <xtc.hpp>
25 
26 namespace loos
27 {
28 
29 
30  // -- The following is adapted from xdrfile-1.1b --
31  // Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
32  // Coordinate compression (c) by Frans van Hoesel.
33 
34  // http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library
35 
36  /******************************************************************
37  GNU LESSER GENERAL PUBLIC LICENSE
38  Version 3, 29 June 2007
39 
40  Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
41  Everyone is permitted to copy and distribute verbatim copies
42  of this license document, but changing it is not allowed.
43 
44 
45  This version of the GNU Lesser General Public License incorporates
46 the terms and conditions of version 3 of the GNU General Public
47 License, supplemented by the additional permissions listed below.
48 
49  0. Additional Definitions.
50 
51  As used herein, "this License" refers to version 3 of the GNU Lesser
52 General Public License, and the "GNU GPL" refers to version 3 of the GNU
53 General Public License.
54 
55  "The Library" refers to a covered work governed by this License,
56 other than an Application or a Combined Work as defined below.
57 
58  An "Application" is any work that makes use of an interface provided
59 by the Library, but which is not otherwise based on the Library.
60 Defining a subclass of a class defined by the Library is deemed a mode
61 of using an interface provided by the Library.
62 
63  A "Combined Work" is a work produced by combining or linking an
64 Application with the Library. The particular version of the Library
65 with which the Combined Work was made is also called the "Linked
66 Version".
67 
68  The "Minimal Corresponding Source" for a Combined Work means the
69 Corresponding Source for the Combined Work, excluding any source code
70 for portions of the Combined Work that, considered in isolation, are
71 based on the Application, and not on the Linked Version.
72 
73  The "Corresponding Application Code" for a Combined Work means the
74 object code and/or source code for the Application, including any data
75 and utility programs needed for reproducing the Combined Work from the
76 Application, but excluding the System Libraries of the Combined Work.
77 
78  1. Exception to Section 3 of the GNU GPL.
79 
80  You may convey a covered work under sections 3 and 4 of this License
81 without being bound by section 3 of the GNU GPL.
82 
83  2. Conveying Modified Versions.
84 
85  If you modify a copy of the Library, and, in your modifications, a
86 facility refers to a function or data to be supplied by an Application
87 that uses the facility (other than as an argument passed when the
88 facility is invoked), then you may convey a copy of the modified
89 version:
90 
91  a) under this License, provided that you make a good faith effort to
92  ensure that, in the event an Application does not supply the
93  function or data, the facility still operates, and performs
94  whatever part of its purpose remains meaningful, or
95 
96  b) under the GNU GPL, with none of the additional permissions of
97  this License applicable to that copy.
98 
99  3. Object Code Incorporating Material from Library Header Files.
100 
101  The object code form of an Application may incorporate material from
102 a header file that is part of the Library. You may convey such object
103 code under terms of your choice, provided that, if the incorporated
104 material is not limited to numerical parameters, data structure
105 layouts and accessors, or small macros, inline functions and templates
106 (ten or fewer lines in length), you do both of the following:
107 
108  a) Give prominent notice with each copy of the object code that the
109  Library is used in it and that the Library and its use are
110  covered by this License.
111 
112  b) Accompany the object code with a copy of the GNU GPL and this license
113  document.
114 
115  4. Combined Works.
116 
117  You may convey a Combined Work under terms of your choice that,
118 taken together, effectively do not restrict modification of the
119 portions of the Library contained in the Combined Work and reverse
120 engineering for debugging such modifications, if you also do each of
121 the following:
122 
123  a) Give prominent notice with each copy of the Combined Work that
124  the Library is used in it and that the Library and its use are
125  covered by this License.
126 
127  b) Accompany the Combined Work with a copy of the GNU GPL and this license
128  document.
129 
130  c) For a Combined Work that displays copyright notices during
131  execution, include the copyright notice for the Library among
132  these notices, as well as a reference directing the user to the
133  copies of the GNU GPL and this license document.
134 
135  d) Do one of the following:
136 
137  0) Convey the Minimal Corresponding Source under the terms of this
138  License, and the Corresponding Application Code in a form
139  suitable for, and under terms that permit, the user to
140  recombine or relink the Application with a modified version of
141  the Linked Version to produce a modified Combined Work, in the
142  manner specified by section 6 of the GNU GPL for conveying
143  Corresponding Source.
144 
145  1) Use a suitable shared library mechanism for linking with the
146  Library. A suitable mechanism is one that (a) uses at run time
147  a copy of the Library already present on the user's computer
148  system, and (b) will operate properly with a modified version
149  of the Library that is interface-compatible with the Linked
150  Version.
151 
152  e) Provide Installation Information, but only if you would otherwise
153  be required to provide such information under section 6 of the
154  GNU GPL, and only to the extent that such information is
155  necessary to install and execute a modified version of the
156  Combined Work produced by recombining or relinking the
157  Application with a modified version of the Linked Version. (If
158  you use option 4d0, the Installation Information must accompany
159  the Minimal Corresponding Source and Corresponding Application
160  Code. If you use option 4d1, you must provide the Installation
161  Information in the manner specified by section 6 of the GNU GPL
162  for conveying Corresponding Source.)
163 
164  5. Combined Libraries.
165 
166  You may place library facilities that are a work based on the
167 Library side by side in a single library together with other library
168 facilities that are not Applications and are not covered by this
169 License, and convey such a combined library under terms of your
170 choice, if you do both of the following:
171 
172  a) Accompany the combined library with a copy of the same work based
173  on the Library, uncombined with any other library facilities,
174  conveyed under the terms of this License.
175 
176  b) Give prominent notice with the combined library that part of it
177  is a work based on the Library, and explaining where to find the
178  accompanying uncombined form of the same work.
179 
180  6. Revised Versions of the GNU Lesser General Public License.
181 
182  The Free Software Foundation may publish revised and/or new versions
183 of the GNU Lesser General Public License from time to time. Such new
184 versions will be similar in spirit to the present version, but may
185 differ in detail to address new problems or concerns.
186 
187  Each version is given a distinguishing version number. If the
188 Library as you received it specifies that a certain numbered version
189 of the GNU Lesser General Public License "or any later version"
190 applies to it, you have the option of following the terms and
191 conditions either of that published version or of any later version
192 published by the Free Software Foundation. If the Library as you
193 received it does not specify a version number of the GNU Lesser
194 General Public License, you may choose any version of the GNU Lesser
195 General Public License ever published by the Free Software Foundation.
196 
197  If the Library as you received it specifies that a proxy can decide
198 whether future versions of the GNU Lesser General Public License shall
199 apply, that proxy's public statement of acceptance of any version is
200 permanent authorization for you to choose that version for the
201 Library.
202 
203 ******************************************************************/
204 
205  // TDR - This needs to move?
206  const int XTCWriter::magicints[] = {
207  0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
208  80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
209  1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003,
210  16384, 20642, 26007, 32768, 41285, 52015, 65536,82570, 104031,
211  131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
212  832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
213  4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216
214  };
215 
216 
217  const int XTCWriter::firstidx = 9;
218  const int XTCWriter::lastidx = sizeof(magicints) / sizeof(*magicints);
219 
220  const int XTCWriter::DIM = 3;
221 
222  /* Internal support routines for reading/writing compressed coordinates
223  * sizeofint - calculate smallest number of bits necessary
224  * to represent a certain integer.
225  */
226  int XTCWriter::sizeofint(const int size) const {
227  int num = 1;
228  int num_of_bits = 0;
229 
230  while (size >= num && num_of_bits < 32) {
231  num_of_bits++;
232  num <<= 1;
233  }
234  return num_of_bits;
235  }
236 
237 
238  /*
239  * sizeofints - calculate 'bitsize' of compressed ints
240  *
241  * given a number of small unsigned integers and the maximum value
242  * return the number of bits needed to read or write them with the
243  * routines encodeints/decodeints. You need this parameter when
244  * calling those routines.
245  * (However, in some cases we can just use the variable 'smallidx'
246  * which is the exact number of bits, and them we dont need to call
247  * this routine).
248  */
249  int XTCWriter::sizeofints(const int num_of_ints, const unsigned int sizes[]) const
250  {
251  int i, num;
252  unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
253  num_of_bytes = 1;
254  bytes[0] = 1;
255  num_of_bits = 0;
256  for (i=0; i < num_of_ints; i++) {
257  tmp = 0;
258  for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
259  tmp = bytes[bytecnt] * sizes[i] + tmp;
260  bytes[bytecnt] = tmp & 0xff;
261  tmp >>= 8;
262  }
263  while (tmp != 0) {
264  bytes[bytecnt++] = tmp & 0xff;
265  tmp >>= 8;
266  }
267  num_of_bytes = bytecnt;
268  }
269  num = 1;
270  num_of_bytes--;
271  while (bytes[num_of_bytes] >= num) {
272  num_of_bits++;
273  num *= 2;
274  }
275  return(num_of_bits + num_of_bytes * 8);
276  }
277 
278 
279  /*
280  * encodebits - encode num into buf using the specified number of bits
281  *
282  * This routines appends the value of num to the bits already present in
283  * the array buf. You need to give it the number of bits to use and you had
284  * better make sure that this number of bits is enough to hold the value.
285  * Num must also be positive.
286  */
287  void XTCWriter::encodebits(int buf[], int num_of_bits, const int num) const
288  {
289 
290  unsigned int cnt, lastbyte;
291  int lastbits;
292  unsigned char * cbuf;
293 
294  cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
295  cnt = (unsigned int) buf[0];
296  lastbits = buf[1];
297  lastbyte =(unsigned int) buf[2];
298  while (num_of_bits >= 8) {
299  lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
300  cbuf[cnt++] = lastbyte >> lastbits;
301  num_of_bits -= 8;
302  }
303  if (num_of_bits > 0) {
304  lastbyte = (lastbyte << num_of_bits) | num;
305  lastbits += num_of_bits;
306  if (lastbits >= 8) {
307  lastbits -= 8;
308  cbuf[cnt++] = lastbyte >> lastbits;
309  }
310  }
311  buf[0] = cnt;
312  buf[1] = lastbits;
313  buf[2] = lastbyte;
314  if (lastbits>0) {
315  cbuf[cnt] = lastbyte << (8 - lastbits);
316  }
317  }
318 
319  /*
320  * encodeints - encode a small set of small integers in compressed format
321  *
322  * this routine is used internally by xdr3dfcoord, to encode a set of
323  * small integers to the buffer for writing to a file.
324  * Multiplication with fixed (specified maximum) sizes is used to get
325  * to one big, multibyte integer. Allthough the routine could be
326  * modified to handle sizes bigger than 16777216, or more than just
327  * a few integers, this is not done because the gain in compression
328  * isn't worth the effort. Note that overflowing the multiplication
329  * or the byte buffer (32 bytes) is unchecked and whould cause bad results.
330  * THese things are checked in the calling routines, so make sure not
331  * to remove those checks...
332  */
333 
334  void XTCWriter::encodeints(int buf[], const int num_of_ints, const int num_of_bits,
335  const unsigned int sizes[], const unsigned int nums[]) const
336  {
337 
338  int i;
339  unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
340 
341  tmp = nums[0];
342  num_of_bytes = 0;
343  do {
344  bytes[num_of_bytes++] = tmp & 0xff;
345  tmp >>= 8;
346  } while (tmp != 0);
347 
348  for (i = 1; i < num_of_ints; i++) {
349  if (nums[i] >= sizes[i]) {
350  std::ostringstream oss;
351  oss << boost::format("Major breakdown in XTCWriter::encodeints() - num %u doesn't match size %u")
352  % nums[i]
353  % sizes[i];
354  throw(LOOSError(oss.str()));
355  }
356  /* use one step multiply */
357  tmp = nums[i];
358  for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
359  tmp = bytes[bytecnt] * sizes[i] + tmp;
360  bytes[bytecnt] = tmp & 0xff;
361  tmp >>= 8;
362  }
363  while (tmp != 0) {
364  bytes[bytecnt++] = tmp & 0xff;
365  tmp >>= 8;
366  }
367  num_of_bytes = bytecnt;
368  }
369  if (num_of_bits >= num_of_bytes * 8) {
370  for (i = 0; i < num_of_bytes; i++) {
371  encodebits(buf, 8, bytes[i]);
372  }
373  encodebits(buf, num_of_bits - num_of_bytes * 8, 0);
374  }
375  else {
376  for (i = 0; i < num_of_bytes-1; i++){
377  encodebits(buf, 8, bytes[i]);
378  }
379  encodebits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
380  }
381  }
382 
383 
384 
385 
386 
387  void XTCWriter::writeCompressedCoordsFloat(float* ptr, int size, float precision)
388  {
389  int minint[3], maxint[3], mindiff, *lip, diff;
390  int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
391  int minidx, maxidx;
392  unsigned sizeint[3], sizesmall[3], bitsizeint[3], *luip;
393  int k;
394  int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
395  float *lfp, lf;
396  int tmp, tmpsum, *thiscoord, prevcoord[3];
397  unsigned int tmpcoord[30];
398  unsigned int bitsize;
399 
400  uint size3 = size * 3;
401 
402  bitsizeint[0] = 0;
403  bitsizeint[1] = 0;
404  bitsizeint[2] = 0;
405 
406  allocateBuffers(size);
407  if (!xdr.write(size))
408  throw(FileWriteError(_filename, "Could not write size to XTC file"));
409 
410  /* Dont bother with compression for three atoms or less */
411  if(size<=9)
412  {
413  xdr.write(ptr, size3);
414  return;
415  }
416  /* Compression-time if we got here. Write precision first */
417  if (precision <= 0)
418  precision = 1000;
419 
420  xdr.write(precision);
421  /* buf2[0-2] are special and do not contain actual data */
422  buf2[0] = buf2[1] = buf2[2] = 0;
423  minint[0] = minint[1] = minint[2] = INT_MAX;
424  maxint[0] = maxint[1] = maxint[2] = INT_MIN;
425  prevrun = -1;
426  lfp = ptr;
427  lip = buf1;
428  mindiff = INT_MAX;
429  oldlint1 = oldlint2 = oldlint3 = 0;
430  while(lfp < ptr + size3 )
431  {
432  /* find nearest integer */
433  if (*lfp >= 0.0)
434  lf = *lfp * precision + 0.5;
435  else
436  lf = *lfp * precision - 0.5;
437  if (fabs(lf) > INT_MAX-2)
438  {
439  /* scaling would cause overflow */
440  throw(LOOSError("Internal overflow compressing coordinates"));
441  }
442  lint1 = lf;
443  if (lint1 < minint[0]) minint[0] = lint1;
444  if (lint1 > maxint[0]) maxint[0] = lint1;
445  *lip++ = lint1;
446  lfp++;
447  if (*lfp >= 0.0)
448  lf = *lfp * precision + 0.5;
449  else
450  lf = *lfp * precision - 0.5;
451  if (fabs(lf) > INT_MAX-2)
452  {
453  /* scaling would cause overflow */
454  throw(LOOSError("Internal overflow compressing coordinates"));
455  }
456  lint2 = lf;
457  if (lint2 < minint[1]) minint[1] = lint2;
458  if (lint2 > maxint[1]) maxint[1] = lint2;
459  *lip++ = lint2;
460  lfp++;
461  if (*lfp >= 0.0)
462  lf = *lfp * precision + 0.5;
463  else
464  lf = *lfp * precision - 0.5;
465 
466  // *** TDR - This is not actually used
467  // if (fabs(lf) > INT_MAX-2)
468  // {
469  // errval=0;
470  // }
471  lint3 = lf;
472  if (lint3 < minint[2]) minint[2] = lint3;
473  if (lint3 > maxint[2]) maxint[2] = lint3;
474  *lip++ = lint3;
475  lfp++;
476  diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
477  if (diff < mindiff && lfp > ptr + 3)
478  mindiff = diff;
479  oldlint1 = lint1;
480  oldlint2 = lint2;
481  oldlint3 = lint3;
482  }
483  xdr.write(minint, 3);
484  xdr.write(maxint, 3);
485 
486  if ((float)maxint[0] - (float)minint[0] >= INT_MAX-2 ||
487  (float)maxint[1] - (float)minint[1] >= INT_MAX-2 ||
488  (float)maxint[2] - (float)minint[2] >= INT_MAX-2) {
489  /* turning value in unsigned by subtracting minint
490  * would cause overflow
491  */
492  throw(LOOSError("Internal overflow compressing internal coordinates"));
493  }
494  sizeint[0] = maxint[0] - minint[0]+1;
495  sizeint[1] = maxint[1] - minint[1]+1;
496  sizeint[2] = maxint[2] - minint[2]+1;
497 
498  /* check if one of the sizes is to big to be multiplied */
499  if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
500  {
501  bitsizeint[0] = sizeofint(sizeint[0]);
502  bitsizeint[1] = sizeofint(sizeint[1]);
503  bitsizeint[2] = sizeofint(sizeint[2]);
504  bitsize = 0; /* flag the use of large sizes */
505  }
506  else
507  {
508  bitsize = sizeofints(3, sizeint);
509  }
510  lip = buf1;
511  luip = (unsigned int *) buf1;
512  smallidx = firstidx;
513  while (smallidx < lastidx && magicints[smallidx] < mindiff)
514  {
515  smallidx++;
516  }
517  xdr.write(smallidx);
518  tmp=smallidx+8;
519  maxidx = (lastidx<tmp) ? lastidx : tmp;
520  minidx = maxidx - 8; /* often this equal smallidx */
521  tmp=smallidx-1;
522  tmp= (firstidx>tmp) ? firstidx : tmp;
523  smaller = magicints[tmp] / 2;
524  smallnum = magicints[smallidx] / 2;
525  sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
526  larger = magicints[maxidx] / 2;
527  i = 0;
528  while (i < size)
529  {
530  is_small = 0;
531  thiscoord = (int *)(luip) + i * 3;
532  if (smallidx < maxidx && i >= 1 &&
533  abs(thiscoord[0] - prevcoord[0]) < larger &&
534  abs(thiscoord[1] - prevcoord[1]) < larger &&
535  abs(thiscoord[2] - prevcoord[2]) < larger) {
536  is_smaller = 1;
537  }
538  else if (smallidx > minidx)
539  {
540  is_smaller = -1;
541  }
542  else
543  {
544  is_smaller = 0;
545  }
546  if (i + 1 < size)
547  {
548  if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
549  abs(thiscoord[1] - thiscoord[4]) < smallnum &&
550  abs(thiscoord[2] - thiscoord[5]) < smallnum)
551  {
552  /* interchange first with second atom for better
553  * compression of water molecules
554  */
555  tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
556  thiscoord[3] = tmp;
557  tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
558  thiscoord[4] = tmp;
559  tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
560  thiscoord[5] = tmp;
561  is_small = 1;
562  }
563  }
564  tmpcoord[0] = thiscoord[0] - minint[0];
565  tmpcoord[1] = thiscoord[1] - minint[1];
566  tmpcoord[2] = thiscoord[2] - minint[2];
567  if (bitsize == 0)
568  {
569  encodebits(buf2, bitsizeint[0], tmpcoord[0]);
570  encodebits(buf2, bitsizeint[1], tmpcoord[1]);
571  encodebits(buf2, bitsizeint[2], tmpcoord[2]);
572  }
573  else
574  {
575  encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
576  }
577  prevcoord[0] = thiscoord[0];
578  prevcoord[1] = thiscoord[1];
579  prevcoord[2] = thiscoord[2];
580  thiscoord = thiscoord + 3;
581  i++;
582 
583  run = 0;
584  if (is_small == 0 && is_smaller == -1)
585  is_smaller = 0;
586  while (is_small && run < 8*3)
587  {
588  tmpsum=0;
589  for(j=0;j<3;j++)
590  {
591  tmp=thiscoord[j] - prevcoord[j];
592  tmpsum+=tmp*tmp;
593  }
594  if (is_smaller == -1 && tmpsum >= smaller * smaller)
595  {
596  is_smaller = 0;
597  }
598 
599  tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
600  tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
601  tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
602 
603  prevcoord[0] = thiscoord[0];
604  prevcoord[1] = thiscoord[1];
605  prevcoord[2] = thiscoord[2];
606 
607  i++;
608  thiscoord = thiscoord + 3;
609  is_small = 0;
610  if (i < size &&
611  abs(thiscoord[0] - prevcoord[0]) < smallnum &&
612  abs(thiscoord[1] - prevcoord[1]) < smallnum &&
613  abs(thiscoord[2] - prevcoord[2]) < smallnum)
614  {
615  is_small = 1;
616  }
617  }
618  if (run != prevrun || is_smaller != 0)
619  {
620  prevrun = run;
621  encodebits(buf2, 1, 1); /* flag the change in run-length */
622  encodebits(buf2, 5, run+is_smaller+1);
623  }
624  else
625  {
626  encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
627  }
628  for (k=0; k < run; k+=3)
629  {
630  encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
631  }
632  if (is_smaller != 0)
633  {
634  smallidx += is_smaller;
635  if (is_smaller < 0)
636  {
637  smallnum = smaller;
638  smaller = magicints[smallidx-1] / 2;
639  }
640  else
641  {
642  smaller = smallnum;
643  smallnum = magicints[smallidx] / 2;
644  }
645  sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
646  }
647  }
648  if (buf2[1] != 0) buf2[0]++;
649  xdr.write(buf2[0]);
650  tmp=xdr.write((char *)&(buf2[3]),(unsigned int)buf2[0]);
651  if(tmp!=(unsigned int)buf2[0])
652  throw(FileWriteError(_filename, "Error while writing compressed coordinates to XTC file"));
653  }
654 
655 
656 
657  // -- End of code from xdrfile library --
658 
659 
660  // Handle allocation of buffers (would be handle by system xdr lib)
661  void XTCWriter::allocateBuffers(const size_t size) {
662  size_t size3 = size * 3;
663  if (size3 > buf1size) {
664  if (buf1)
665  delete[] buf1;
666  if (buf2)
667  delete[] buf2;
668 
669  buf1 = new int[size3];
670  size_t size3plus = size3 * 1.2;
671  buf2 = new int[size3plus];
672 
673  buf1size = size3;
674  buf2size = size3 * 1.2;
675  }
676  }
677 
678 
679  // Write a frame header
680  void XTCWriter::writeHeader(const int natoms, const int step, const float time) {
681  int magic = 1995;
682 
683  xdr.write(magic);
684  xdr.write(natoms);
685  xdr.write(step);
686  xdr.write(time);
687  }
688 
689 
690  // Write a periodic box, translating from A to nm
691  void XTCWriter::writeBox(const GCoord& box) {
692  float outbox[DIM*DIM];
693  for (uint i=0; i < DIM*DIM; ++i)
694  outbox[i] = 0.0;
695  outbox[0] = box[0] / 10.0;
696  outbox[4] = box[1] / 10.0;
697  outbox[8] = box[2] / 10.0;
698 
699  xdr.write(outbox, DIM*DIM);
700  }
701 
702 
703 
704  // Write a frame, converting units from A to nm. Will allocate a temp array to hold coords...
705  void XTCWriter::writeFrame(const AtomicGroup& model, const uint step, const double time) {
706 
707  writeHeader(model.size(), step, time);
708  writeBox(model.periodicBox());
709  uint n = model.size();
710 
711  if (n > crds_size_) {
712  delete[] crds_;
713  crds_ = new float[n * 3];
714  crds_size_ = n;
715  }
716 
717  for (uint i=0,k=0; i<n; ++i) {
718  GCoord c = model[i]->coords();
719  crds_[k++] = c.x() / 10.0; // Convert to nm
720  crds_[k++] = c.y() / 10.0;
721  crds_[k++] = c.z() / 10.0;
722  }
723  writeCompressedCoordsFloat(crds_, n, precision_);
724 
725  ++current_;
726  }
727 
728 
729 
730  void XTCWriter::writeFrame(const AtomicGroup& model) {
731  writeFrame(model, step_, dt_ * step_);
732  step_ += steps_per_frame_;
733  }
734 
735 
736  // Read existing XTC to get frame count...
737  void XTCWriter::prepareToAppend() {
738  stream_->seekg(0);
739  XTC xtc(*stream_);
740  current_ = xtc.nframes();
741  stream_->seekp(0, std::ios_base::end);
742  }
743 
744 };
void writeFrame(const AtomicGroup &model)
Write a frame to the trajectory.
Definition: xtcwriter.cpp:730
uint write(const T &p)
Writes a single datum.
Definition: xdr.hpp:192
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
Class representing GROMACS reduced precision, compressed trajectories.
Definition: xtc.hpp:61
GCoord periodicBox(void) const
Fetch the periodic boundary conditions.