LOOS  v2.3.2
trr.hpp
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 #if !defined(LOOS_TRR_HPP)
32 #define LOOS_TRR_HPP
33 #include <ios>
34 #include <iostream>
35 #include <string>
36 #include <stdexcept>
37 #include <vector>
38 
39 #include <loos_defs.hpp>
40 
41 #include <Coord.hpp>
42 #include <xdr.hpp>
43 #include <AtomicGroup.hpp>
44 #include <Trajectory.hpp>
45 
46 #include <boost/format.hpp>
47 
48 
49 namespace loos {
50 
52 
78  class TRR : public Trajectory {
79  static const int magic; // various internal constants
80  // GROMACS uses...
81  static const int DIM;
82 
83  // This is the internal header for a TRR frame
84  struct Header {
85  bool bDouble;
86  int ir_size;
87  int e_size;
88  int box_size;
89  int vir_size;
90  int pres_size;
91  int top_size;
92  int sym_size;
93  int x_size;
94  int v_size;
95  int f_size;
96 
97  int natoms;
98  int step;
99  int nre;
100  float tf;
101  float lambdaf;
102  double td;
103  double lambdad;
104  };
105 
106  public:
107  explicit TRR(const std::string& s) : Trajectory(s), xdr_file(ifs()) {
108  init();
109  }
110 
111  explicit TRR(const char* p) : Trajectory(p), xdr_file(ifs()) {
112  init();
113  }
114  explicit TRR(std::istream& is) : Trajectory(is), xdr_file(ifs()) {
115  init();
116  }
117 
118  std::string description() const { return("Gromacs TRR"); }
119  static pTraj create(const std::string& fname, const AtomicGroup& model) {
120  return(pTraj(new TRR(fname)));
121  }
122 
123 
124  uint natoms(void) const { return(hdr_.natoms); }
125 
129  float timestep(void) const { return(0); }
130 
131  uint nframes(void) const { return(frame_indices.size()); }
132  bool hasPeriodicBox(void) const { return(hdr_.box_size != 0); }
133  GCoord periodicBox(void) const { return(box); }
134 
135 
136  std::vector<GCoord> coords(void) { return(coords_); }
137 
138  // TRR specific attributes...
139  std::vector<double> virial(void) const { return(vir_); }
140  std::vector<double> pressure(void) const { return(pres_); }
141  std::vector<GCoord> velocities(void) const { return(velo_); }
142  std::vector<GCoord> forces(void) const { return(forc_); }
143 
144  bool isDouble(void) const { return(hdr_.bDouble); }
145  bool hasVirial(void) const { return(hdr_.vir_size != 0); }
146  bool hasPressure(void) const { return(hdr_.pres_size != 0); }
147  bool hasCoords(void) const { return(hdr_.x_size != 0); }
148  bool hasVelocities(void) const { return(hdr_.v_size != 0); }
149  bool hasForces(void) const { return(hdr_.f_size != 0); }
150 
151  double time(void) const { return( hdr_.bDouble ? hdr_.td : hdr_.tf ); }
152  double lambda(void) const { return( hdr_.bDouble ? hdr_.lambdad : hdr_.lambdaf); }
153 
154  int step(void) const { return(hdr_.step); }
155 
156  bool parseFrame(void) {
157  if (!readHeader(hdr_)) // Catch EOF
158  return(false);
159 
160  if (hdr_.bDouble)
161  return(readRawFrame<double>());
162 
163  return(readRawFrame<float>());
164  }
165 
166 
167  private:
168  void init(void);
169  int floatSize(Header& h);
170  bool readHeader(Header& h);
171 
172  // These simplify reading of blocks of data from the TRR file.
173  // Templatized since GROMACS can store data in either single or
174  // double-precision...
175 
176  template<typename T>
177  void readBlock(std::vector<double>& v, const uint n, const std::string& msg) {
178  T buf[n];
179  uint i = xdr_file.read(buf, n);
180  if (i != n)
181  throw(FileReadError(_filename, "Unable to read " + msg));
182 
183  for (i=0; i<n; ++i)
184  v.push_back(buf[i]);
185  }
186 
187 
188  // This assumes the block of data are triplets and converts them
189  // into GCoords, scaling from nm to Angstroms along the way...
190  template<typename T>
191  void readBlock(std::vector<GCoord>& v, const uint n, const std::string& msg) {
192  T* buf = new T[n];
193  if (buf == 0)
194  throw(std::runtime_error("Out of memory"));
195 
196  if (xdr_file.read(buf, n) != n) {
197  delete[] buf;
198  throw(FileReadError(_filename, "Unable to read " + msg));
199  }
200  for (uint i=0; i<n; i += DIM)
201  v.push_back(GCoord(buf[i], buf[i+1], buf[i+2]) * 10.0);
202 
203  delete[] buf;
204  }
205 
206 
207  // Note: Assumes that the object Header has already been read...
208  template<typename T>
209  bool readRawFrame() {
210 
211  // Clear data first...
212  box_.clear();
213  vir_.clear();
214  pres_.clear();
215  velo_.clear();
216  forc_.clear();
217  coords_.clear();
218 
219  if (hdr_.box_size) {
220  readBlock<T>(box_, DIM*DIM, "box");
221  box = GCoord(box_[0], box_[4], box_[8]) * 10.0; // Convert
222  // to angstroms
223  }
224 
225  if (hdr_.vir_size)
226  readBlock<T>(vir_, DIM*DIM, "virial");
227  if (hdr_.pres_size)
228  readBlock<T>(pres_, DIM*DIM, "pressure");
229 
230  if (hdr_.x_size)
231  readBlock<T>(coords_, hdr_.natoms * DIM, "Coordinates");
232 
233  if (hdr_.v_size)
234  readBlock<T>(velo_, hdr_.natoms * DIM, "Velocities");
235 
236  if (hdr_.f_size)
237  readBlock<T>(forc_, hdr_.natoms * DIM, "Forces");
238 
239 
240  return(! ((xdr_file.get())->fail() || (xdr_file.get())->eof()) );
241  }
242 
243  void rewindImpl(void) { ifs()->clear(); ifs()->seekg(0, std::ios_base::beg); }
244  void seekNextFrameImpl(void) { }
245  void seekFrameImpl(uint);
246  void updateGroupCoordsImpl(AtomicGroup& g);
247 
248 
249 
250  private:
251  internal::XDRReader xdr_file;
252  std::vector<GCoord> coords_;
253  GCoord box;
254  std::vector<size_t> frame_indices; // Index into file for start
255  // of frame header
256 
257  std::vector<double> box_;
258  std::vector<double> vir_;
259  std::vector<double> pres_;
260  std::vector<GCoord> velo_;
261  std::vector<GCoord> forc_;
262 
263  Header hdr_;
264  };
265 
266 
267 
268 };
269 
270 
271 #endif
272 
Errors that occur while reading a file.
Definition: exceptions.hpp:180
uint read(T *p)
Read a single datum.
Definition: xdr.hpp:67
bool parseFrame(void)
Parse an actual frame.
Definition: trr.hpp:156
GCoord periodicBox(void) const
Returns the periodic box for the current frame/trajectory.
Definition: trr.hpp:133
Class representing the GROMACS TRR trajectory files.
Definition: trr.hpp:78
uint natoms(void) const
of atoms per frame
Definition: trr.hpp:124
bool hasPeriodicBox(void) const
Definition: trr.hpp:132
std::istream * get(void)
Returns the stored istream pointer.
Definition: xdr.hpp:64
uint nframes(void) const
Number of frames in the trajectory.
Definition: trr.hpp:131
float timestep(void) const
Definition: trr.hpp:129
std::string description() const
Return a string describing trajectory format.
Definition: trr.hpp:118
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
Base-class for polymorphic trajectories.
Definition: Trajectory.hpp:64
std::vector< GCoord > coords(void)
Returns the current frames coordinates as a vector of GCoords.
Definition: trr.hpp:136