LOOS  v2.3.2
trr.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 <trr.hpp>
32 
33 
34 namespace loos {
35 
36  // GROMACS magic constants...
37  const int TRR::magic = 1993;
38  const int TRR::DIM = 3;
39 
40 
41  // GROMACS voodoo to determine precision of stored data...
42  int TRR::floatSize(Header& hdr) {
43  int i;
44 
45  if (hdr.box_size)
46  i = hdr.box_size/(DIM*DIM);
47  else if (hdr.x_size)
48  i = hdr.x_size/(hdr.natoms * DIM);
49  else if (hdr.v_size)
50  i = hdr.v_size/(hdr.natoms * DIM);
51  else if (hdr.f_size)
52  i = hdr.f_size/(hdr.natoms * DIM);
53  else
54  throw(FileReadError(_filename, "Cannot determine float size"));
55 
56  if (i != sizeof(float) && i != sizeof(double))
57  throw(FileReadError(_filename, "Float size does not match native sizes"));
58 
59  return(i);
60  }
61 
62 
63  bool TRR::readHeader(Header& hdr) {
64  uint magic_no;
65 
66  xdr_file.read(magic_no);
67  if ((xdr_file.get())->eof())
68  return(false);
69 
70  if (magic_no != magic) {
71  std::ostringstream oss;
72  oss << "Invalid magic number in TRR file...expected " << magic << ", but found " << magic_no;
73 
74  throw(FileReadError(_filename, oss.str()));
75  }
76 
77 
78  std::string version;
79  xdr_file.read(version);
80  xdr_file.read(hdr.ir_size);
81  xdr_file.read(hdr.e_size);
82  xdr_file.read(hdr.box_size);
83  xdr_file.read(hdr.vir_size);
84  xdr_file.read(hdr.pres_size);
85  xdr_file.read(hdr.top_size);
86  xdr_file.read(hdr.sym_size);
87  xdr_file.read(hdr.x_size);
88  xdr_file.read(hdr.v_size);
89  xdr_file.read(hdr.f_size);
90  xdr_file.read(hdr.natoms);
91 
92  int flsz = floatSize(hdr);
93  hdr.bDouble = (flsz == sizeof(double));
94 
95  xdr_file.read(hdr.step);
96  xdr_file.read(hdr.nre);
97 
98  if (hdr.bDouble) {
99  xdr_file.read(hdr.td);
100  hdr.tf = hdr.td;
101  xdr_file.read(hdr.lambdad);
102  hdr.lambdaf = hdr.lambdad;
103  } else {
104  xdr_file.read(hdr.tf);
105  hdr.td = hdr.tf;
106  xdr_file.read(hdr.lambdaf);
107  hdr.lambdad = hdr.lambdaf;
108  }
109 
110  if ((xdr_file.get())->fail())
111  throw(FileReadError(_filename, "Cannot read TRR header"));
112 
113  return(true);
114  }
115 
116 
117  // Initialize the object, along with scanning file for frames to
118  // build the frame index, and finally caches the first frame.
119  void TRR::init(void) {
120  Header h;
121  h.natoms = 0;
122 
123  // Since the trajectory can have different number of atoms per
124  // frame, we track the max and then reserve that space...
125  int maxatoms = 0;
126 
127  rewindImpl();
128  frame_indices.clear();
129 
130  size_t frame_start = (xdr_file.get())->tellg();
131  while (readHeader(h)) {
132  frame_indices.push_back(frame_start);
133  if (h.natoms > maxatoms)
134  maxatoms = h.natoms;
135 
136  uint b = sizeof(internal::XDRReader::block_type);
137  // Correct if double-precision TRR file...
138  if (h.bDouble)
139  b = sizeof(double);
140 
141  size_t offset = (h.box_size ? DIM*DIM*b : 0) + (h.vir_size ? DIM*DIM*b : 0)
142  + (h.pres_size ? DIM*DIM*b : 0) + (h.x_size ? h.natoms*DIM*b : 0)
143  + (h.v_size ? h.natoms*DIM*b : 0) + (h.f_size ? h.natoms*DIM*b : 0);
144 
145  (xdr_file.get())->seekg(offset, std::ios_base::cur);
146  frame_start = (xdr_file.get())->tellg();
147  }
148 
149  coords_.reserve(maxatoms);
150  velo_.reserve(maxatoms);
151  forc_.reserve(maxatoms);
152 
153  rewindImpl();
154 
155  parseFrame();
156  cached_first = true;
157  hdr_ = h;
158 
159  }
160 
161  void TRR::updateGroupCoordsImpl(AtomicGroup& g) {
162 
163  for (AtomicGroup::iterator i = g.begin(); i != g.end(); ++i) {
164  uint idx = (*i)->index();
165  if (static_cast<uint>(idx) >= natoms())
166  throw(LOOSError(**i, "atom index into trajectory frame is out of range"));
167  (*i)->coords(coords_[idx]);
168  }
169 
170  if (hdr_.box_size)
171  g.periodicBox(box);
172  }
173 
174  void TRR::seekFrameImpl(uint i) {
175  if (i >= frame_indices.size())
176  throw(FileError(_filename, "Requested TRR frame is out of range"));
177 
178  ifs()->clear();
179  ifs()->seekg(frame_indices[i], std::ios_base::beg);
180  }
181 
182 
183 }
uint read(T *p)
Read a single datum.
Definition: xdr.hpp:67
bool parseFrame(void)
Parse an actual frame.
Definition: trr.hpp:156
uint natoms(void) const
of atoms per frame
Definition: trr.hpp:124
unsigned int block_type
Type (and hence size) of the external block.
Definition: xdr.hpp:51
std::istream * get(void)
Returns the stored istream pointer.
Definition: xdr.hpp:64
Namespace for most things not already encapsulated within a class.