LOOS  v2.3.2
xtc.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2009, 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 /*
24  NOTE: This code is based on the xdrfile library authored by:
25  David van der Spoel <spoel@gromacs.org>
26  Erik Lindahl <lindahl@gromacs.org>
27  and covered by the GLPL-v3 license
28 */
29 
30 
31 #include <xtc.hpp>
32 
33 
34 namespace loos {
35 
36 
37  // Globals...
38  const int XTC::magicints[] = {
39  0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
40  80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290, // 31
41  1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, // 41
42  16384, 20642, 26007, 32768, 41285, 52015, 65536,82570, 104031, // 50
43  131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561, // 58
44  832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, // 65
45  4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216 // 72
46  };
47 
48  const int XTC::firstidx = 9;
49 
50  const int XTC::lastidx = 73;
51 
52  const int XTC::magic = 1995;
53 
54  const uint XTC::min_compressed_system_size = 9;
55 
56 
57 
58  // The following are largely from the xdrlib...
59  int XTC::sizeofint(int size) {
60  int n = 0;
61  while ( size > 0 ) {
62  size >>= 1;
63  ++n;
64  }
65  return(n);
66  }
67 
68 
69  int XTC::sizeofints(uint* data, const uint n) {
70  uint nbytes = 1;
71  uint bytes[32];
72  uint nbits = 0;
73 
74  bytes[0] = 1;
75  for (uint i = 0; i < n; ++i) {
76  uint tmp = 0;
77  uint bytecnt;
78  for (bytecnt = 0; bytecnt < nbytes; ++bytecnt) {
79  tmp += bytes[bytecnt] * data[i];
80  bytes[bytecnt] = tmp & 0xff;
81  tmp >>= 8;
82  }
83  while (tmp != 0) {
84  bytes[bytecnt++] = tmp & 0xff;
85  tmp >>= 8;
86  }
87  nbytes = bytecnt;
88  }
89 
90  uint num = 1;
91  --nbytes;
92  while (bytes[nbytes] >= num) {
93  ++nbits;
94  num *= 2;
95  }
96 
97  return(nbits + nbytes*8);
98  }
99 
100 
101  int XTC::decodebits(int* buf, uint nbits) {
102 
103  int mask = (1 << nbits) -1;
104 
105  unsigned char *cbuf = reinterpret_cast<unsigned char*>(buf) + 3*sizeof(*buf);
106  int cnt = buf[0];
107  uint lastbits = static_cast<uint>(buf[1]);
108  uint lastbyte = static_cast<uint>(buf[2]);
109 
110  int num = 0;
111  while (nbits >= 8) {
112  lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
113  num |= (lastbyte >> lastbits) << (nbits - 8);
114  nbits -=8;
115  }
116  if (nbits > 0) {
117  if (lastbits < nbits) {
118  lastbits += 8;
119  lastbyte = (lastbyte << 8) | cbuf[cnt++];
120  }
121  lastbits -= nbits;
122  num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
123  }
124  num &= mask;
125  buf[0] = cnt;
126  buf[1] = static_cast<int>(lastbits);
127  buf[2] = static_cast<int>(lastbyte);
128 
129  return(num);
130  }
131 
132 
133  void XTC::decodeints(int* buf, const int nints, int nbits,
134  uint* sizes, int* nums) {
135  int bytes[32];
136  int i, j, num_of_bytes, p, num;
137 
138  bytes[1] = bytes[2] = bytes[3] = 0;
139  num_of_bytes = 0;
140  while (nbits > 8) {
141  bytes[num_of_bytes++] = decodebits(buf, 8);
142  nbits -= 8;
143  }
144  if (nbits > 0) {
145  bytes[num_of_bytes++] = decodebits(buf, nbits);
146  }
147  for (i = nints-1; i > 0; i--) {
148  num = 0;
149  for (j = num_of_bytes-1; j >=0; j--) {
150  num = (num << 8) | bytes[j];
151  p = num / sizes[i];
152  bytes[j] = p;
153  num = num - p * sizes[i];
154  }
155  nums[i] = num;
156  }
157  nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
158  }
159 
160 
161 
162  // Coordinates are converted into GCoords and stored in the object's
163  // coords_ vector
164 
165  bool XTC::readCompressedCoords(void)
166  {
167  int minint[3], maxint[3], *lip;
168  int smallidx;
169  uint sizeint[3], sizesmall[3], bitsizeint[3] = {0,0,0}, size3;
170  int k, *buf1, *buf2, lsize, flag;
171  int smallnum, smaller, i, is_smaller, run;
172  xtc_t precision, inv_precision;
173  int tmp, *thiscoord, prevcoord[3];
174  unsigned int bitsize;
175 
176 
177  if (!xdr_file.read(lsize))
178  return(false);
179 
180  size3 = lsize * 3;
181 
182  /* Dont bother with compression for three atoms or less */
183  if(lsize<=9) {
184  float* tmp = new xtc_t[size3];
185  xdr_file.read(tmp, size3);
186  for (uint i=0; i<size3; i += 3)
187  coords_.push_back(GCoord(tmp[i], tmp[i+1], tmp[i+2]) * 10.0);
188  delete[] tmp;
189  return(true);
190  }
191 
192  /* Compression-time if we got here. Read precision first */
193  xdr_file.read(precision);
194  precision_ = precision;
195 
196  int size3padded = static_cast<int>(size3 * 1.2);
197  buf1 = new int[size3padded];
198  buf2 = new int[size3padded];
199  /* buf2[0-2] are special and do not contain actual data */
200  buf2[0] = buf2[1] = buf2[2] = 0;
201  xdr_file.read(minint, 3);
202  xdr_file.read(maxint, 3);
203 
204  sizeint[0] = maxint[0] - minint[0]+1;
205  sizeint[1] = maxint[1] - minint[1]+1;
206  sizeint[2] = maxint[2] - minint[2]+1;
207 
208  /* check if one of the sizes is to big to be multiplied */
209  if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
210  bitsizeint[0] = sizeofint(sizeint[0]);
211  bitsizeint[1] = sizeofint(sizeint[1]);
212  bitsizeint[2] = sizeofint(sizeint[2]);
213  bitsize = 0; /* flag the use of large sizes */
214  } else {
215  bitsize = sizeofints(sizeint, 3);
216  }
217 
218  if (!xdr_file.read(smallidx)) {
219  delete[] buf1;
220  delete[] buf2;
221  return(false);
222  }
223 
224  tmp=smallidx+8;
225  tmp = smallidx-1;
226  tmp = (firstidx>tmp) ? firstidx : tmp;
227  smaller = magicints[tmp] / 2;
228  smallnum = magicints[smallidx] / 2;
229  sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
230 
231  /* buf2[0] holds the length in bytes */
232 
233  if (!xdr_file.read(buf2, 1)) {
234  delete[] buf1;
235  delete[] buf2;
236  return(false);
237  }
238 
239  if (!xdr_file.read(reinterpret_cast<char*>(&(buf2[3])), static_cast<uint>(buf2[0]))) {
240  delete[] buf1;
241  delete[] buf2;
242  return(false);
243  }
244  buf2[0] = buf2[1] = buf2[2] = 0;
245 
246  inv_precision = 1.0 / precision;
247  run = 0;
248  i = 0;
249  lip = buf1;
250  while ( i < lsize ) {
251  thiscoord = (int *)(lip) + i * 3;
252 
253  if (bitsize == 0) {
254  thiscoord[0] = decodebits(buf2, bitsizeint[0]);
255  thiscoord[1] = decodebits(buf2, bitsizeint[1]);
256  thiscoord[2] = decodebits(buf2, bitsizeint[2]);
257  } else {
258  decodeints(buf2, 3, bitsize, sizeint, thiscoord);
259  }
260 
261  i++;
262  thiscoord[0] += minint[0];
263  thiscoord[1] += minint[1];
264  thiscoord[2] += minint[2];
265 
266  prevcoord[0] = thiscoord[0];
267  prevcoord[1] = thiscoord[1];
268  prevcoord[2] = thiscoord[2];
269 
270  flag = decodebits(buf2, 1);
271  is_smaller = 0;
272  if (flag == 1) {
273  run = decodebits(buf2, 5);
274  is_smaller = run % 3;
275  run -= is_smaller;
276  is_smaller--;
277  }
278  if (run > 0) {
279  thiscoord += 3;
280  for (k = 0; k < run; k+=3) {
281  decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
282  i++;
283  thiscoord[0] += prevcoord[0] - smallnum;
284  thiscoord[1] += prevcoord[1] - smallnum;
285  thiscoord[2] += prevcoord[2] - smallnum;
286  if (k == 0) {
287  /* interchange first with second atom for better
288  * compression of water molecules
289  */
290  tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
291  prevcoord[0] = tmp;
292  tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
293  prevcoord[1] = tmp;
294  tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
295  prevcoord[2] = tmp;
296 
297  coords_.push_back(GCoord(prevcoord[0] * inv_precision,
298  prevcoord[1] * inv_precision,
299  prevcoord[2] * inv_precision) * 10.0);
300  } else {
301  prevcoord[0] = thiscoord[0];
302  prevcoord[1] = thiscoord[1];
303  prevcoord[2] = thiscoord[2];
304  }
305  coords_.push_back(GCoord(thiscoord[0] * inv_precision,
306  thiscoord[1] * inv_precision,
307  thiscoord[2] * inv_precision) * 10.0);
308  }
309  } else {
310  coords_.push_back(GCoord(thiscoord[0] * inv_precision,
311  thiscoord[1] * inv_precision,
312  thiscoord[2] * inv_precision) * 10.0);
313  }
314  smallidx += is_smaller;
315  if (is_smaller < 0) {
316  smallnum = smaller;
317  if (smallidx > firstidx) {
318  smaller = magicints[smallidx - 1] /2;
319  } else {
320  smaller = 0;
321  }
322  } else if (is_smaller > 0) {
323  smaller = smallnum;
324  smallnum = magicints[smallidx] / 2;
325  }
326  sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
327  }
328 
329  delete[] buf1;
330  delete[] buf2;
331  return(true);
332  }
333 
334 
335 
336  bool XTC::readUncompressedCoords(void)
337  {
338  uint lsize;
339 
340  if (!xdr_file.read(lsize))
341  return(false);
342 
343  uint size3 = lsize * 3;
344  coords_ = std::vector<GCoord>(lsize);
345  float* tmp_coords = new float[size3];
346  uint n = xdr_file.read(tmp_coords, size3);
347  if (n != size3)
348  throw(FileReadError(_filename, "XTC Error: number of uncompressed coords read did not match number expected"));
349 
350  uint i = 0;
351  for (uint j=0; j<lsize; ++j, i += 3)
352  coords_[j] = GCoord(tmp_coords[i], tmp_coords[i+1], tmp_coords[i+2]) * 10.0;
353 
354  delete[] tmp_coords;
355  return(true);
356  }
357 
358 
359 
360  void XTC::updateGroupCoordsImpl(AtomicGroup& g) {
361 
362  for (AtomicGroup::iterator i = g.begin(); i != g.end(); ++i) {
363  uint idx = (*i)->index();
364  if (idx > natoms_)
365  throw(LOOSError(**i, "atom index into trajectory frame is out of range"));
366  (*i)->coords(coords_[idx]);
367  }
368 
369  // XTC files *always* have a periodic box...
370  g.periodicBox(box);
371  }
372 
373 
374  bool XTC::parseFrame(void) {
375  if (ifs()->eof())
376  return(false);
377 
378  // First, clear out existing coords... A read error after this
379  // point will invalidate the current object's coord state
380 
381  coords_.clear();
382  if (!readFrameHeader(current_header_))
383  return(false);
384 
385  box = GCoord(current_header_.box[0],
386  current_header_.box[4],
387  current_header_.box[8]) * 10.0; // Convert to Angstroms
388  if (natoms_ <= min_compressed_system_size)
389  return(readUncompressedCoords());
390  else
391  return(readCompressedCoords());
392  }
393 
394 
395  bool XTC::readFrameHeader(XTC::Header& hdr) {
396  int magic_no;
397  int ok = xdr_file.read(magic_no);
398  if (!ok)
399  return(false);
400  if (magic_no != magic) {
401  std::ostringstream oss;
402  oss << "Invalid XTC magic number (got " << magic_no << " but expected " << magic << ")";
403  throw(FileReadError(_filename, oss.str()));
404  }
405 
406  // Defer error-checks until the end...
407  xdr_file.read(hdr.natoms);
408 
409  xdr_file.read(hdr.step);
410  xdr_file.read(hdr.time);
411  ok = xdr_file.read(hdr.box, 9);
412  if (!ok)
413  throw(FileReadError(_filename, "Problem reading XTC header"));
414 
415  return(true);
416  }
417 
418 
419  // Scan the trajectory file, skipping each compressed frame. In the
420  // process, we build up an index relating file-pos to frame index.
421  // This permits fast seeking of indivual frames.
422  void XTC::scanFrames(void) {
423  frame_indices.clear();
424 
425  rewindImpl();
426 
427  Header h;
428 
429  while (! ifs()->eof()) {
430  size_t pos = ifs()->tellg();
431 
432  bool ok = readFrameHeader(h);
433  if (!ok) {
434  rewindImpl();
435  return;
436  }
437 
438  frame_indices.push_back(pos);
439  if (natoms_ == 0)
440  natoms_ = h.natoms;
441  else if (natoms_ != h.natoms)
442  throw(FileOpenError(_filename, "XTC frames have differing numbers of atoms"));
443 
444  uint block_size = sizeof(internal::XDRReader::block_type);
445 
446  // Always update estimated timestep...
447  if (h.step != 0)
448  timestep_ = h.time / h.step;
449 
450  size_t offset = 0;
451  uint nbytes = 0;
452 
453  if (natoms_ <= min_compressed_system_size) {
454  nbytes = natoms_ * 3 * sizeof(float);
455  uint dummy;
456  xdr_file.read(dummy);
457  if (dummy != natoms_)
458  throw(FileOpenError(_filename, "XTC small system vector size is not what was expected"));
459  } else {
460  offset = 9 * block_size;
461  ifs()->seekg(offset, std::ios_base::cur);
462  xdr_file.read(nbytes);
463  }
464 
465  uint nblocks = nbytes / block_size;
466 
467  if (nbytes % block_size != 0)
468  ++nblocks; // round up
469  offset = nblocks * block_size;
470  ifs()->seekg(offset, std::ios_base::cur);
471  }
472 
473  // Catch-all for I/O errors
474  if (ifs()->fail())
475  throw(FileOpenError(_filename, "Problem scanning XTC trajectory to build frame indices"));
476 
477  }
478 
479 
480  void XTC::seekFrameImpl(const uint i) {
481  if (i >= frame_indices.size())
482  throw(FileError(_filename, "Requested XTC frame is out of range"));
483 
484  ifs()->clear();
485  ifs()->seekg(frame_indices[i], std::ios_base::beg);
486  }
487 
488 }
489 
uint read(T *p)
Read a single datum.
Definition: xdr.hpp:67
double precision(void) const
Return the stored file's precision.
Definition: xtc.hpp:111
unsigned int block_type
Type (and hence size) of the external block.
Definition: xdr.hpp:51
Namespace for most things not already encapsulated within a class.