LOOS  v2.3.2
amber_netcdf.cpp
1 // (c) 2012 Tod D. Romo, Grossfield Lab, URMC
2 
3 #include <amber_netcdf.hpp>
4 #include <AtomicGroup.hpp>
5 
6 namespace loos {
7 
8 
9  bool isFileNetCDF(const std::string& fname) {
10  std::ifstream ifs(fname.c_str());
11 
12  char buf[4];
13  ifs.read(buf, 4);
14  return (buf[0] == 'C' && buf[1] == 'D' && buf[2] == 'F' && (buf[3] = 0x01 || buf[3] == 0x02));
15  }
16 
17 
18 
19 
20  void AmberNetcdf::init(const char* name, const uint natoms) {
21  int retval;
22 
23 
24  retval = nc_open(name, NC_NOWRITE, &_ncid);
25  if (retval)
26  throw(FileOpenError(name, "", retval)); // May want to preserve code here in the future...
27 
28  // Read and validate global attributes...
29  readGlobalAttributes();
30  if (_conventions.empty() || _conventionVersion.empty())
31  throw(FileOpenError(name, "Unable to find convention global attributes. Is this really an Amber NetCDF trajectory?"));
32 
33  if (_conventions.find("AMBER") == std::string::npos)
34  throw(FileOpenError(name, "Cannot find AMBER tag in global attribues. Is this really an Amber NetCDF trajectory?"));
35 
36  if (_conventionVersion != std::string("1.0"))
37  throw(FileOpenError(name, "Convention version is '" + _conventionVersion + "', but only 1.0 is supported for Amber NetCDF trajectories."));
38 
39  // Verify # of atoms match...
40  int atom_id;
41  retval = nc_inq_dimid(_ncid, "atom", &atom_id);
42  if (retval)
43  throw(FileOpenError(name, "Cannot read atom id"), retval);
44  retval = nc_inq_dimlen(_ncid, atom_id, &_natoms);
45  if (retval)
46  throw(FileOpenError(name, "Cannot read atom length", retval));
47  if (_natoms != natoms)
48  throw(FileOpenError(name, "AmberNetcdf has different number of atoms than expected"));
49 
50 
51  // Get nframes
52  int frame_id;
53  retval = nc_inq_dimid(_ncid, "frame", &frame_id);
54  if (retval)
55  throw(FileOpenError(name, "Cannot read frame information", retval));
56  retval = nc_inq_dimlen(_ncid, frame_id, &_nframes);
57 
58  // Check for periodic cells...
59  retval = nc_inq_varid(_ncid, "cell_lengths", &_cell_lengths_id);
60  _periodic = !retval;
61 
62  // Any additional validation should go here...
63 
64  // Get coord-id for later use...
65  retval = nc_inq_varid(_ncid, "coordinates", &_coord_id);
66  if (retval)
67  throw(FileOpenError(name, "Cannot get id for coordinates", retval));
68 
69  // Attempt to determine timestep by looking at dT between frames 1 & 2
70  if (_nframes >= 2) {
71  int time_id;
72  retval = nc_inq_varid(_ncid, "time", &time_id);
73  if (!retval) {
74  float t0, t1;
75  size_t idx[1];
76 
77  idx[0] = 0;
78  retval = nc_get_var1_float(_ncid, time_id, idx, &t0);
79  if (retval)
80  throw(FileOpenError(name, "Cannot get first time point", retval));
81 
82  idx[0] = 1;
83  retval = nc_get_var1_float(_ncid, time_id, idx, &t1);
84  if (retval)
85  throw(FileOpenError(name, "Cannot get second time point", retval));
86 
87  // Assume units are in picoseconds
88  _timestep = (t1-t0)*1e-12;
89  }
90  }
91 
92 
93  // Now cache the first frame...
94  readRawFrame(0);
95  cached_first = true;
96 
97  }
98 
99 
100  // Given a frame number, read the coord data into the internal array
101  // and retrieve the corresponding periodic box (if present)
102  void AmberNetcdf::readRawFrame(const uint frameno) {
103  size_t start[3] = {0, 0, 0};
104  size_t count[3] = {1, 1, 3};
105 
106 
107  // Read coordinates first...
108  start[0] = frameno;
109  count[1] = _natoms;
110 
111 
112  int retval = VarTypeDecider<GCoord::element_type>::read(_ncid, _coord_id, start, count, _coord_data);
113  if (retval)
114  throw(FileReadError(_filename, "Cannot read Amber netcdf frame", retval));
115 
116  // Now get box if present...
117  if (_periodic) {
118  start[1] = 0;
119  count[1] = 3;
120 
121  retval = VarTypeDecider<GCoord::element_type>::read(_ncid, _cell_lengths_id, start, count, _box_data);
122  if (retval)
123  throw(FileReadError(_filename, "Cannot read Amber netcdf periodic box", retval));
124  }
125 
126  }
127 
128 
129  void AmberNetcdf::seekNextFrameImpl() {
130  ++_current_frame;
131  }
132 
133  void AmberNetcdf::seekFrameImpl(const uint i) {
134  if (i >= _nframes)
135  throw(FileError(_filename, "Attempting seek frame beyond end of trajectory"));
136  _current_frame = i;
137  }
138 
139  void AmberNetcdf::rewindImpl() {
140  _current_frame = 0;
141  }
142 
143  bool AmberNetcdf::parseFrame() {
144  if (_current_frame >= _nframes)
145  return(false);
146 
147  readRawFrame(_current_frame);
148  return(true);
149  }
150 
151 void AmberNetcdf::updateGroupCoordsImpl(AtomicGroup& g) {
152 
153  for (AtomicGroup::iterator i = g.begin(); i != g.end(); ++i) {
154  uint idx = (*i)->index();
155  if (idx >= _natoms)
156  throw(LOOSError(**i, "Atom index into trajectory frame is out of bounds"));
157  idx *= 3;
158  (*i)->coords(GCoord(_coord_data[idx], _coord_data[idx+1], _coord_data[idx+2]));
159  }
160 
161  if (_periodic)
162  g.periodicBox(GCoord(_box_data[0], _box_data[1], _box_data[2]));
163  }
164 
165 
166  void AmberNetcdf::readGlobalAttributes() {
167 
168  _title = readGlobalAttribute("title");
169  _application = readGlobalAttribute("application");
170  _program = readGlobalAttribute("program");
171  _programVersion = readGlobalAttribute("programVersion");
172  _conventions = readGlobalAttribute("Conventions");
173  _conventionVersion = readGlobalAttribute("ConventionVersion");
174  }
175 
176 
177  // Will return an emptry string if the attribute is not found
178  // May want to consider special exception for type errors...
179  std::string AmberNetcdf::readGlobalAttribute(const std::string& name) {
180  size_t len;
181 
182  int retval = nc_inq_attlen(_ncid, NC_GLOBAL, name.c_str(), &len);
183  if (retval)
184  return(std::string());
185 
186  nc_type type;
187  retval = nc_inq_atttype(_ncid, NC_GLOBAL, name.c_str(), &type);
188  if (type != NC_CHAR)
189  throw(FileOpenError(_filename, "Only character data is supported for global attributes", retval));
190 
191  char* buf = new char[len+1];
192  retval = nc_get_att_text(_ncid, NC_GLOBAL, name.c_str(), buf);
193  if (retval) {
194  delete[] buf;
195  throw(FileOpenError(_filename, "Cannot read attribute " + name, retval));
196  }
197  buf[len]='\0';
198 
199  return(std::string(buf));
200  }
201 
202 
203 
204 
205 
206 };
bool isFileNetCDF(const std::string &fname)
Returns true if the file is a NetCDF file.
Definition: amber_netcdf.cpp:9
Namespace for most things not already encapsulated within a class.