1 /**********************************************************************
2 opendxformat.cpp - Read in OpenDX cube format files from APBS
3
4 Copyright (c) 2007-2008 by Geoffrey R. Hutchison
5 Some Portions Copyright (C) 2008 by Marcus D. Hanwell
6
7 This file is part of the Open Babel project.
8 For more information, see <http://openbabel.org/>
9
10 This program is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation: either version 2 of the License, or (at
13 your option) any later version.
14
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19 ***********************************************************************/
20
21 #include <sstream>
22
23 #include <openbabel/babelconfig.h>
24 #include <openbabel/obmolecformat.h>
25 #include <openbabel/mol.h>
26 #include <openbabel/atom.h>
27 #include <openbabel/elements.h>
28
29 #include <openbabel/griddata.h>
30 #include <cstdlib>
31
32 using namespace std;
33
34 namespace OpenBabel
35 {
36
37 // Ultimately, this should just read in an OBFloatGrid, rather than a "molecule"
38 // since there's no actual molecular data. For now, we'll read in an OBMol with no atoms
39 // and attach an OBGridData
40 class OBOpenDXCubeFormat : public OpenBabel::OBMoleculeFormat
41 {
42 public:
43 // Constructor: register "dx" extension code
OBOpenDXCubeFormat()44 OBOpenDXCubeFormat()
45 {
46 OpenBabel::OBConversion::RegisterFormat( "dx", this );
47 }
48
49 // Return description.
Description()50 virtual const char* Description() //required
51 {
52 return
53 "OpenDX cube format for APBS\n"
54 "A volume data format for IBM's Open Source visualization software\n"
55 "The OpenDX support is currently designed to read the OpenDX cube\n"
56 "files from APBS. \n\n";
57 }
58
59 // Return a specification url, not really a specification since
60 // I couldn't find it but close enough.
SpecificationURL()61 virtual const char* SpecificationURL()
62 {
63 return "http://apbs.sourceforge.net/doc/user-guide/index.html#id504516";
64 }
65
66 // Return MIME type, NULL in this case.
GetMIMEType()67 virtual const char* GetMIMEType() { return nullptr; }
68
69 // Skip to object: used for multi-object file formats.
SkipObjects(int n,OpenBabel::OBConversion * pConv)70 virtual int SkipObjects( int n, OpenBabel::OBConversion* pConv ) { return 0; }
71
Flags()72 virtual unsigned int Flags()
73 {
74 return READONEONLY | WRITEONEONLY | ZEROATOMSOK;
75 };
76
77 /// The "API" interface functions
78 virtual bool ReadMolecule( OpenBabel::OBBase* pOb, OpenBabel::OBConversion* pConv );
79 virtual bool WriteMolecule( OpenBabel::OBBase* pOb, OpenBabel::OBConversion* pConv );
80
81 };
82
83 //------------------------------------------------------------------------------
84
85 // Global variable used to register OpenDX cube format.
86 OBOpenDXCubeFormat theOpenDXCubeFormat;
87
88 //------------------------------------------------------------------------------
ReadMolecule(OBBase * pOb,OBConversion * pConv)89 bool OBOpenDXCubeFormat::ReadMolecule( OBBase* pOb, OBConversion* pConv )
90 {
91 OBMol* pmol = pOb->CastAndClear<OBMol>();
92 if (pmol == nullptr)
93 return false;
94
95 istream& ifs = *pConv->GetInStream();
96
97 const char* title = pConv->GetTitle();
98 char buffer[BUFF_SIZE];
99
100 stringstream errorMsg;
101
102 if (!ifs)
103 return false; // We are attempting to read past the end of the file
104
105 pmol->SetTitle(title);
106
107 while (ifs.good() && ifs.getline(buffer,BUFF_SIZE)) {
108 if (buffer[0] == '#')
109 continue; // comment line
110 if (EQn(buffer, "object", 6))
111 break;
112 }
113 if (!ifs)
114 return false; // ran out of lines
115
116 vector<string> vs;
117 tokenize(vs, buffer);
118
119 // Number of grid points (voxels)
120 vector<int> voxels(3);
121 if (!EQn(buffer, "object", 6) || vs.size() != 8)
122 return false;
123 else {
124 voxels[0] = atoi(vs[5].c_str());
125 voxels[1] = atoi(vs[6].c_str());
126 voxels[2] = atoi(vs[7].c_str());
127 }
128
129 double x, y, z;
130 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "origin", 6))
131 return false;
132 else {
133 tokenize(vs, buffer);
134 if (vs.size() != 4)
135 return false;
136 x = atof(vs[1].c_str());
137 y = atof(vs[2].c_str());
138 z = atof(vs[3].c_str());
139 }
140 vector3 origin(x, y, z);
141
142 // now three lines with the x, y, and z axes
143 vector<vector3> axes;
144 for (unsigned int i = 0; i < 3; ++i) {
145 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "delta", 5))
146 return false;
147 else {
148 tokenize(vs, buffer);
149 if (vs.size() != 4)
150 return false;
151 x = atof(vs[1].c_str());
152 y = atof(vs[2].c_str());
153 z = atof(vs[3].c_str());
154 axes.push_back(vector3(x, y, z));
155 }
156 }
157
158 // Two remaining header lines before the data:
159 /*
160 object 2 class gridconnections counts nx ny nz
161 object 3 class array type double rank 0 times n data follows
162 */
163 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
164 return false;
165 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
166 return false;
167
168 pmol->BeginModify();
169 pmol->SetDimension(3);
170
171 OBGridData *gd = new OBGridData;
172 gd->SetAttribute("OpenDX");
173
174 // get all values as one vector<double>
175 char *endptr;
176 vector<double> values;
177 int n = voxels[0]*voxels[1]*voxels[2];
178 int line = 0;
179 values.reserve(n);
180 while (ifs.getline(buffer, BUFF_SIZE))
181 {
182 ++line;
183 if (EQn(buffer, "attribute", 9))
184 break; // we're finished with reading data -- although we should probably have a voxel check in here too
185
186 tokenize(vs, buffer);
187 if (vs.size() == 0)
188 {
189 errorMsg << "Problem reading the OpenDX grid file: cannot"
190 << " read line " << line
191 << ", there does not appear to be any data in it.\n"
192 << buffer << "\n";
193 obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obError);
194 return false;
195 }
196
197 for (unsigned int l = 0; l < vs.size(); ++l)
198 {
199 values.push_back(strtod(static_cast<const char*>(vs[l].c_str()), &endptr));
200 }
201 }
202
203 gd->SetNumberOfPoints(voxels[0], voxels[1], voxels[2]);
204 gd->SetLimits(origin, axes[0], axes[1], axes[2]);
205 gd->SetUnit(OBGridData::ANGSTROM);
206 gd->SetOrigin(fileformatInput); // i.e., is this data from a file or determined by Open Babel
207 gd->SetValues(values); // set the values
208 pmol->SetData(gd); // store the grids in the OBMol
209 pmol->EndModify();
210
211 // Trailing control lines
212 /*
213 attribute "dep" string "positions"
214 object "regular positions regular connections" class field
215 component "positions" value 1
216 component "connections" value 2
217 component "data" value 3
218 */
219 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "object", 6))
220 return false;
221 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
222 return false;
223 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
224 return false;
225 if (!ifs.getline(buffer, BUFF_SIZE) || !EQn(buffer, "component", 9))
226 return false;
227
228 // clean out any remaining blank lines
229 std::streampos ipos;
230 do
231 {
232 ipos = ifs.tellg();
233 ifs.getline(buffer,BUFF_SIZE);
234 }
235 while(strlen(buffer) == 0 && !ifs.eof() );
236 ifs.seekg(ipos);
237
238 return true;
239 }
240
241 //------------------------------------------------------------------------------
WriteMolecule(OBBase * pOb,OBConversion * pConv)242 bool OBOpenDXCubeFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
243 {
244 OBMol* pmol = dynamic_cast<OBMol*>(pOb);
245 if (pmol == nullptr)
246 return false;
247
248 ostream &ofs = *pConv->GetOutStream();
249 OBMol &mol = *pmol;
250
251 char buffer[BUFF_SIZE];
252 string str;
253 stringstream errorMsg;
254
255 OBGridData *gd = (OBGridData*)mol.GetData(OBGenericDataType::GridData);
256 if (gd == nullptr) {
257 errorMsg << "The molecule has no grid.";
258 obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);
259 return false;
260 }
261
262 // APBS-style OpenDX Multigrid
263 // First some comments
264 ofs << "# Data from Open Babel " << BABEL_VERSION << "\n";
265 str = mol.GetTitle();
266 if (str.empty())
267 ofs << "# Molecule Title: *****" << "\n";
268 else
269 ofs << "# Molecule Title: " << str << "\n";
270
271 int nx, ny, nz;
272 double origin[3], xAxis[3], yAxis[3], zAxis[3];
273 gd->GetAxes(xAxis, yAxis, zAxis);
274 gd->GetNumberOfPoints(nx, ny, nz);
275 gd->GetOriginVector(origin);
276
277 // data line 1: # of points in x, y, z (nx, ny, nz)
278 snprintf(buffer, BUFF_SIZE, "object 1 class gridpositions counts %5d %5d %5d", nx, ny, nz);
279 ofs << buffer << "\n";
280
281 // data line 2: origin (x, y, z)
282 snprintf(buffer, BUFF_SIZE,"origin %12.6f %12.6f %12.6f",
283 origin[0], origin[1], origin[2]);
284 ofs << buffer << "\n";
285
286 // data line 3: x-displacement
287 snprintf(buffer, BUFF_SIZE,"delta %12.6f %12.6f %12.6f",
288 xAxis[0], xAxis[1], xAxis[2]);
289 ofs << buffer << "\n";
290
291 // data line 4: y-displacement
292 snprintf(buffer, BUFF_SIZE,"delta %12.6f %12.6f %12.6f",
293 yAxis[0], yAxis[1], yAxis[2]);
294 ofs << buffer << "\n";
295
296 // data line 5: z-displacement
297 snprintf(buffer, BUFF_SIZE,"delta %12.6f %12.6f %12.6f",
298 zAxis[0], zAxis[1], zAxis[2]);
299 ofs << buffer << "\n";
300
301 // data line 6: # of points in x, y, z (nx, ny, nz)
302 snprintf(buffer, BUFF_SIZE, "object 2 class gridconnections counts %5d %5d %5d", nx, ny, nz);
303 ofs << buffer << "\n";
304
305 // data line 7: total # of points
306 snprintf(buffer, BUFF_SIZE, "object 3 class array type double rank 0 items %5d data follows", nx*ny*nz);
307 ofs << buffer << "\n";
308
309 // The cube(s)
310 double value;
311 unsigned int count = 1;
312 for (int i = 0; i < nx; ++i)
313 {
314 for (int j = 0; j < ny; ++j)
315 {
316 for (int k = 0; k < nz; ++k)
317 {
318 value = gd->GetValue(i, j, k);
319 snprintf(buffer, BUFF_SIZE," %12.5E", value);
320 if (count % 3 == 0)
321 ofs << buffer << "\n";
322 else
323 ofs << buffer;
324 count++;
325 } // z-axis
326 } // y-axis
327 } // x-axis
328
329 if (count % 3 != 0)
330 ofs << "\n";
331 ofs << "attribute \"dep\" string \"positions\"\n";
332 ofs << "object \"regular positions regular connections\" class field\n";
333 ofs << "component \"positions\" value 1\n";
334 ofs << "component \"connections\" value 2\n";
335 ofs << "component \"data\" value 3\n";
336
337 return true;
338 }
339 }
340