← Back to branch summary

~ubuntu-branches/debian/stretch/openbabel/stretch

« back to all changes in this revision

Viewing changes to src/formats/opendxformat.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-07-22 23:54:58 UTC
  • mfrom: (3.1.10 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080722235458-3o606czluviz4akx
Tags: 2.2.0-2
* Upload to unstable.
* debian/control: Updated descriptions.
* debian/patches/gauss_cube_format.patch: New patch, makes the 
  gaussian cube format available again.
* debian/rules (DEB_DH_MAKESHLIBS_ARGS_libopenbabel3): Removed.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Likewise.
* debian/libopenbabel3.install: Adjust formats directory to single 
  version hierarchy.

Show diffs side-by-side

added added

removed removed

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