LOOS  v2.3.2
dcd.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, 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 
25 #include <iostream>
26 #include <fstream>
27 #include <sstream>
28 #include <exception>
29 #include <stdexcept>
30 #include <vector>
31 
32 #include <stdio.h>
33 #include <string.h>
34 #include <assert.h>
35 
36 #include <dcd.hpp>
37 #include <AtomicGroup.hpp>
38 
39 
40 namespace loos {
41 
42  uint DCD::natoms(void) const { return(_natoms); }
43  bool DCD::hasPeriodicBox(void) const { return(_icntrl[10] == 1); }
44  GCoord DCD::periodicBox(void) const { return(GCoord(qcrys[0], qcrys[1], qcrys[2])); }
45 
46 
47  bool DCD::suppress_warnings = false;
48 
49 
50  std::vector<std::string> DCD::titles(void) const { return(_titles); }
51 
52  int DCD::icntrl(const int i) const { assert(i>=0 && i<20); return(_icntrl[i]); }
53  void DCD::icntrl(const int i, const int val) { assert(i>=0 && i<20); _icntrl[i] = val; }
54 
55  // * legacy *
56  std::vector<double> DCD::crystalParams(void) const { return(qcrys); }
57  bool DCD::hasCrystalParams(void) const { return(_icntrl[10] == 1); }
58 
59  float DCD::timestep(void) const { return(_delta); }
60  uint DCD::nframes(void) const { return(_icntrl[0]); }
61 
62  std::vector<dcd_real> DCD::xcoords(void) const { return(xcrds); }
63  std::vector<dcd_real> DCD::ycoords(void) const { return(ycrds); }
64  std::vector<dcd_real> DCD::zcoords(void) const { return(zcrds); }
65 
66  // The following track CHARMm names (more or less...)
67  unsigned int DCD::nsteps(void) const { return(_icntrl[3]); }
68  float DCD::delta(void) const { return(_delta); }
69  int DCD::nsavc(void) const { return(_icntrl[2]); }
70  int DCD::nfile(void) const { return(_icntrl[0]); }
71  int DCD::nfixed(void) const { return(_icntrl[8]); }
72 
73  bool DCD::nativeFormat(void) const { return(!swabbing); }
74 
75 
76  // Allocate the input buffer and pre-size the coordinate vectors.
77  // n = number of atoms
78 
79  void DCD::allocateSpace(const int n) {
80  xcrds = std::vector<dcd_real>(n);
81  ycrds = std::vector<dcd_real>(n);
82  zcrds = std::vector<dcd_real>(n);
83  }
84 
85 
86 
87  // Read the F77 record length from the file stream
88  // Returns 0 at EOF
89 
90  unsigned int DCD::readRecordLen(void) {
91  DataOverlay o;
92 
93  ifs()->read(o.c, 4);
94 
95  if (ifs()->eof())
96  return(0);
97  if (ifs()->fail())
98  throw(FileReadError(_filename, "Unable to read DCD record length"));
99 
100  uint data = o.ui;
101  if (swabbing)
102  data = swab(data);
103 
104  return(data);
105  }
106 
107 
108  // Check for endian-ness
109  void DCD::endianMatch(StreamWrapper& fsw) {
110  unsigned long curpos = fsw()->tellg();
111  unsigned int datum;
112  ifs()->read((char *)(&datum), sizeof(datum));
113  fsw()->seekg(curpos);
114 
115  if (ifs()->eof() || ifs()->fail())
116  throw(FileReadError(_filename, "Unable to read first datum from DCD file"));
117 
118  if (datum != 0x54) {
119  datum = swab(datum);
120  if (datum != 0x54)
121  throw(FileReadError(_filename, "Unable to determine endian-ness of DCD file"));
122  swabbing = true;
123  } else
124  swabbing = false;
125  }
126 
127 
128  // Read a full line of F77-formatted data.
129  // Returns a pointer to the read data and puts the # of bytes read into *len
130  // Returns a null pointer and 0 length at EOF
131  // Note: It is up to the caller to swab individual elements...
132 
133  DCD::DataOverlay* DCD::readF77Line(unsigned int *len) {
134  DataOverlay* ptr;
135  unsigned int n, n2;
136 
137  *len = 0;
138  n = readRecordLen();
139  if (n == 0)
140  return(0);
141 
142  ptr = new DataOverlay[n];
143 
144  ifs()->read((char *)ptr, n);
145  if (ifs()->fail())
146  throw(FileReadError(_filename, "Error reading data record from DCD"));
147 
148  n2 = readRecordLen();
149  if (n != n2)
150  throw(FileReadError(_filename, "Mismatch in record length while reading from DCD"));
151 
152  *len = n;
153  return(ptr);
154  }
155 
156 
157 
158  // Read in the DCD header...
159 
160  void DCD::readHeader(void) {
161  DataOverlay *ptr;
162  unsigned int len;
163  char *cp;
164  int i;
165 
166  endianMatch(ifs);
167  ptr = readF77Line(&len);
168  if (len != 84)
169  throw(FileReadError(_filename, "Error while reading DCD header"));
170 
171  // Check for the magic name. Ignore swabbing for now...
172  DataOverlay o = ptr[0];
173  if (!(o.c[0] == 'C' && o.c[1] == 'O' && o.c[2] == 'R' && o.c[3] == 'D'))
174  throw(FileReadError(_filename, "DCD is missing CORD magic marker"));
175 
176  // Copy in the ICNTRL data...
177  for (i=0; i<20; i++) {
178  if (swabbing)
179  _icntrl[i] = swab(ptr[i+1].i);
180  else
181  _icntrl[i] = ptr[i+1].i;
182  }
183 
184  // Extract the delta value...
185  if (swabbing)
186  _delta = swab(ptr[10].f);
187  else
188  _delta = ptr[10].f;
189 
190  if (nfixed() != 0)
191  throw(LOOSError("Fixed atoms not yet supported by LOOS DCD reader"));
192 
193  delete[] ptr;
194 
195  // Now read in the TITLE info...
196 
197  ptr = readF77Line(&len);
198  if (!ptr)
199  throw(FileReadError(_filename, "Unexpected EOF reading DCD TITLE"));
200  char sbuff[81];
201  int ntitle = ptr[0].i;
202  if (swabbing)
203  ntitle = swab(ntitle);
204 
205  cp = (char *)(ptr + 1);
206 
207  for (i=0; i<ntitle; i++) {
208  memcpy(sbuff, cp + 80*i, 80);
209  sbuff[80] = '\0';
210  std::string s(sbuff);
211  _titles.push_back(s);
212  }
213  delete[] ptr;
214 
215  // get the NATOMS...
216  ptr = readF77Line(&len);
217  if (len != 4)
218  throw(FileReadError(_filename, "Error reading number of atoms from DCD"));
219  if (swabbing)
220  _natoms = swab(ptr->i);
221  else
222  _natoms = ptr->i;
223  delete[] ptr;
224 
225 
226  // Finally, set internal variables and allocate space for a frame...
227  first_frame_pos = ifs()->tellg();
228 
229  frame_size = 12 * (2 + _natoms);
230  if (hasCrystalParams())
231  frame_size += 56;
232 
233  allocateSpace(_natoms);
234 
235  // Issue warnings
236  if (nframes() == 0 && !suppress_warnings)
237  std::cerr << "Warning- DCD '" << _filename << "' appears empty; verify with dcdinfo and fix with fixdcd" << std::endl;
238  }
239 
240 
241 
242  // Read in and reorder the crystal parameters...
243  // NOTE: This is already double!
244 
245 
246  bool DCD::readCrystalParams(void) {
247  unsigned int len;
248  DataOverlay* o;
249 
250  o = readF77Line(&len);
251  if (!o)
252  return(false);
253 
254  if (len != 48)
255  throw(FileReadError(_filename, "Cannot read crystal parameters"));
256 
257  double* dp = reinterpret_cast<double*>(o);
258 
259  qcrys[0] = dp[0];
260  qcrys[1] = dp[2];
261  qcrys[2] = dp[5];
262  qcrys[3] = dp[1];
263  qcrys[4] = dp[3];
264  qcrys[5] = dp[4];
265 
266  if (swabbing)
267  for (int i=0; i<6; ++i)
268  qcrys[i] = swab(qcrys[i]);
269 
270  delete[] o;
271 
272  return(true);
273  }
274 
275 
276 
277  // Read a line of coordinates into the specified vector.
278 
279  bool DCD::readCoordLine(std::vector<dcd_real>& v) {
280  DataOverlay *op;
281  int n = _natoms * sizeof(dcd_real);
282  unsigned int len;
283 
284 
285  op = readF77Line(&len);
286  if (!op)
287  return(false);
288 
289  if (len != (unsigned int)n)
290  throw(FileReadError(_filename, "Size of coords stored in frame does not match model size"));
291 
292  // Recast the values as floats and store them...
293  uint i;
294  for (i=0; i<_natoms; i++)
295  if (swabbing)
296  v[i] = swab(op[i].f);
297  else
298  v[i] = op[i].f;
299 
300  delete[] op;
301 
302  return(true);
303  }
304 
305 
306  void DCD::seekFrameImpl(const uint i) {
307 
308  if (first_frame_pos == 0)
309  throw(FileError(_filename, "Trying to seek to a DCD frame without having first read the header"));
310 
311  if (i >= nframes())
312  throw(FileError(_filename, "Requested DCD frame is out of range"));
313 
314  ifs()->clear();
315  ifs()->seekg(first_frame_pos + i * frame_size);
316  if (ifs()->fail() || ifs()->bad())
317  throw(FileError(_filename, "Cannot seek to requested frame"));
318  }
319 
320 
321  // Read in a frame of data.
322  // Returns TRUE if success or FALSE if at EOF.
323  // Throws an exception if there was an error... (Should we throw EOF
324  // instead?)
325 
326  bool DCD::parseFrame(void) {
327 
328  if (first_frame_pos == 0)
329  throw(FileReadError(_filename, "Trying to read a DCD frame without first having read the header."));
330 
331  // This will not catch most cases of reading to the end of the file...
332  if (ifs()->eof())
333  return(false);
334 
335  if (hasCrystalParams())
336  if (!readCrystalParams())
337  return(false);
338 
339  if (!readCoordLine(xcrds))
340  return(false);
341 
342  if (!readCoordLine(ycrds))
343  throw(FileReadError(_filename, "Unexpected EOF reading Y-coordinates from DCD"));
344  if (!readCoordLine(zcrds))
345  throw(FileReadError(_filename, "Unexepcted EOF reading Z-coordinates from DCD"));
346 
347  return(true);
348  }
349 
350 
351  void DCD::rewindImpl(void) {
352  ifs()->clear();
353  ifs()->seekg(first_frame_pos);
354  if (ifs()->fail() || ifs()->bad())
355  throw(FileError(_filename, "Cannot rewind DCD trajectory"));
356  }
357 
358 
359  // ----------------------------------------------------------
360 
361 
362  std::vector<GCoord> DCD::coords(void) {
363  std::vector<GCoord> crds(_natoms);
364 
365  for (uint i=0; i<_natoms; i++) {
366  crds[i].x(xcrds[i]);
367  crds[i].y(ycrds[i]);
368  crds[i].z(zcrds[i]);
369  }
370 
371  return(crds);
372  }
373 
374  std::vector<GCoord> DCD::mappedCoords(const std::vector<int>& indices) {
375  std::vector<int>::const_iterator iter;
376  std::vector<GCoord> crds(indices.size());
377 
378  int j = 0;
379  for (iter = indices.begin(); iter != indices.end(); iter++, j++) {
380  int index = *iter;
381  crds[j].x(xcrds[index]);
382  crds[j].y(ycrds[index]);
383  crds[j].z(zcrds[index]);
384  }
385 
386  return(crds);
387  }
388 
389 
390  void DCD::updateGroupCoordsImpl(AtomicGroup& g) {
391  for (AtomicGroup::iterator i = g.begin(); i != g.end(); ++i) {
392  uint idx = (*i)->index();
393  if (idx >= _natoms)
394  throw(LOOSError(**i, "Atom index into the trajectory frame is out of bounds"));
395  (*i)->coords(GCoord(xcrds[idx], ycrds[idx], zcrds[idx]));
396  }
397 
398  // Handle periodic boundary conditions (if present)
399  if (hasPeriodicBox()) {
401  }
402  }
403 
404 
405 
406  void DCD::initTrajectory() {
407  readHeader();
408  bool b = parseFrame();
409  if (!b)
410  throw(LOOSError("Cannot read first frame of DCD during initialization"));
411  cached_first = true;
412  }
413 
414 
415 }
std::vector< dcd_real > xcoords(void) const
Return the raw coords...
Definition: dcd.cpp:62
Errors related to File I/O.
Definition: exceptions.hpp:110
Errors that occur while reading a file.
Definition: exceptions.hpp:180
std::vector< dcd_real > zcoords(void) const
Return the raw coords...
Definition: dcd.cpp:64
virtual float timestep(void) const
Timestep per frame.
Definition: dcd.cpp:59
bool nativeFormat(void) const
Returns true if the DCD file being read is in the native endian format.
Definition: dcd.cpp:73
virtual std::vector< GCoord > coords(void)
Auto-interleave the coords into a vector of GCoord()'s.
Definition: dcd.cpp:362
virtual bool hasPeriodicBox(void) const
Definition: dcd.cpp:43
T swab(const T &datum)
Returns a byte-swapped copy of an arbitrary type.
Definition: utils.hpp:287
std::vector< GCoord > mappedCoords(const std::vector< int > &map)
Interleave coords, selecting entries indexed by map.
Definition: dcd.cpp:374
std::vector< dcd_real > ycoords(void) const
Return the raw coords...
Definition: dcd.cpp:63
Generic LOOS exception.
Definition: exceptions.hpp:40
virtual GCoord periodicBox(void) const
Returns the periodic box for the current frame/trajectory.
Definition: dcd.cpp:44
virtual uint nframes(void) const
Number of frames in the trajectory.
Definition: dcd.cpp:60
virtual uint natoms(void) const
of atoms per frame
Definition: dcd.cpp:42
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
virtual bool parseFrame(void)
Parse a frame of the DCD.
Definition: dcd.cpp:326
GCoord periodicBox(void) const
Fetch the periodic boundary conditions.