ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/atom2md/oopseformat.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/atom2md/oopseformat.cpp (file contents):
Revision 3319 by gezelter, Wed Jan 23 03:45:33 2008 UTC vs.
Revision 3446 by cli2, Wed Sep 10 19:51:45 2008 UTC

# Line 16 | Line 16 | GNU General Public License for more details.
16  
17   #include <openbabel/babelconfig.h>
18   #include <openbabel/obmolecformat.h>
19 + #include <openbabel/obiter.h>
20 + #include <openbabel/mol.h>
21 + #include <openbabel/chains.h>
22 + #include <openbabel/data.h>
23   #include <fstream>
24  
25 + #include "utils/StringUtils.hpp"
26 +
27   using namespace std;
28   namespace OpenBabel
29   {
# Line 48 | Line 54 | namespace OpenBabel
54        return NOTREADABLE | WRITEONEONLY;
55      }
56      
51 –    //*** This section identical for most OBMol conversions ***
52 –    ////////////////////////////////////////////////////
53 –    /// The "API" interface functions
57      virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
58  
59    private:
60      bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
58 –    void findAngles(OBMol& mol);
61      OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
62 <    void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices);
62 >    void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os,
63 >                     OBMol& mol, vector<int>& indices);
64 >    void CalcBoundingBox(OBMol &mol,
65 >                         double &min_x, double &max_x,
66 >                         double &min_y, double &max_y,
67 >                         double &min_z, double &max_z);
68 >      
69    };
70 <
70 >  
71    //Make an instance of the format class
72    OOPSEFormat theOOPSEFormat;
73 <  
74 <  /////////////////////////////////////////////////////////////////
67 <
68 <
69 <  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
70 <  {
73 >
74 >  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
75      OBMol* pmol = dynamic_cast<OBMol*>(pOb);
76      if(pmol==NULL)
77        return false;
# Line 79 | Line 83 | namespace OpenBabel
83      vector<vector<int> > molecules;
84      vector<int> indices;
85      for(int i =0; i < used.size(); ++i) {
86 <      if (used[i])
87 <        {
88 <          continue;
85 <        }
86 >      if (used[i])
87 >        continue;
88 >
89        used[i] = true;
90        vector<int> sameMolTypes;
91        sameMolTypes.push_back(i);
92        indices.insert(indices.end(), fragmentLists[i].begin(),
93                       fragmentLists[i].end());
94 <      for (int j = i + 1;j < used.size(); ++j)
95 <        {
96 <          if (used[j])
97 <            {
98 <              continue;
99 <            }
100 <          
101 <          if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
102 <            {
100 <              sameMolTypes.push_back(j);
101 <              indices.insert(indices.end(), fragmentLists[j].begin(),
102 <                             fragmentLists[j].end());
103 <              used[j]=true;
104 <            }
94 >      for (int j = i + 1;j < used.size(); ++j) {
95 >        if (used[j])
96 >          continue;
97 >                
98 >        if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j])) {
99 >          sameMolTypes.push_back(j);
100 >          indices.insert(indices.end(), fragmentLists[j].begin(),
101 >                         fragmentLists[j].end());
102 >          used[j]=true;
103          }
104 <      molecules.push_back(sameMolTypes);
105 <      
104 >      }
105 >      molecules.push_back(sameMolTypes);      
106      }
107      
110 –    //
108      vector<OBMol*> mdMols;    
109      vector<int> numMols;
110      for(vector<vector<int> >::iterator  i = molecules.begin();
111 <        i != molecules.end(); ++i)
112 <      {
113 <        mdMols.push_back(createMolFromFragment(*pmol,
114 <                                               fragmentLists[i->front()]));
115 <        numMols.push_back((*i).size());
116 <      }
111 >        i != molecules.end(); ++i) {
112 >      
113 >      mdMols.push_back(createMolFromFragment(*pmol,
114 >                                             fragmentLists[i->front()]));
115 >      numMols.push_back((*i).size());
116 >    }
117      
118      string OutputFileName = pConv->GetInFilename();
119 <    unsigned int pos = OutputFileName.rfind(".");
120 <    if(pos==string::npos)
124 <      OutputFileName += ".md";
125 <    else
119 >    size_t pos = OutputFileName.rfind(".");
120 >    if(pos!=string::npos)
121        OutputFileName = OutputFileName.substr(0, pos) + ".md";      
122 +    else
123 +      OutputFileName += ".md";
124 +    
125      ofstream ofs(OutputFileName.c_str());
126 <    if(!ofs)
129 <      {
126 >    if(!ofs) {
127          cerr << "Cannot write to " << OutputFileName <<endl;
128          return false;
129 <      }
129 >    }
130      
131      WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
132      
133 <    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i)
134 <      {
135 <        delete *i;
139 <      }
133 >    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) {
134 >      delete *i;
135 >    }
136      
141 –    //    
142 –    
137      return(true);
138    }
139    
140    bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
141 <                                     vector<int>& frag2)
148 <  {
141 >                                     vector<int>& frag2) {
142      if (frag1.size() != frag2.size())
143 <      {
151 <        return false;
152 <      }
143 >      return false;
144      
145 <    //exact graph matching is a NP complete problem
146 <    /** @todo using sparse matrix to store the connectivities*/
147 <    for (unsigned int i =0 ; i < frag1.size(); ++i)
148 <      {
145 >    // Exact graph matching is an NP complete problem.
146 >    // This just matches all of the atom atomic numbers and may falsely
147 >    // detect identical fragments which aren't really identical.
148 >    // @todo using sparse matrix to store the connectivities
149 >
150 >    for (unsigned int i =0 ; i < frag1.size(); ++i) {
151          OBAtom* atom1 = mol.GetAtom(frag1[i]);
152          OBAtom* atom2 = mol.GetAtom(frag2[i]);
153          
154 <        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
155 <          {
156 <            return false;
157 <          }
165 <      }
154 >        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
155 >          return false;
156 >        
157 >    }
158      return true;
159    }
160    
161 <  struct SameAngle
170 <  {
161 >  struct SameAngle {
162      bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
163 <                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
173 <    {
163 >                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const {
164        return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) || (t1.first == t2.third && t1.third == t2.first));
165      }
166    };
167 <  
178 <  void OOPSEFormat::findAngles(OBMol& mol) {
179 <    /*
180 <    //if already has data return
181 <    if(mol.HasData(OBGenericDataType::AngleData))
182 <    return;
183 <    
184 <    vector<OBEdgeBase*>::iterator bi1,bi2;
185 <    OBBond* bond;
186 <    OBAtom *a,*b,*c;
187 <    
188 <    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle> uniqueAngles;
189 <    //loop through all bonds generating torsions
190 <    for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
191 <    {
192 <        b = bond->GetBeginAtom();
193 <        c = bond->GetEndAtom();
194 <        if(b->IsHydrogen())
195 <            continue;
167 >
168  
169 <        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
170 <        {
199 <            if(a == c)
200 <                continue;          
201 <            
202 <            uniqueAngles.insert(triple<OBAtom*,OBAtom*,OBAtom*>(a, b, c));
203 <        }
204 <    }
205 <
206 <    //get new data and attach it to molecule
207 <    OBAngleData *angles = new OBAngleData;
208 <    mol.SetData(angles);
209 <    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle>::iterator i;
210 <
211 <    for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
212 <        OBAngle angle;
213 <        angle.SetAtoms(i->first, i->second, i->second);
214 <        angles->SetData(angle);
215 <    }
216 <    */
217 <  }
218 <  
219 <  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment) {
169 >  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol,
170 >                                            vector<int>& fragment) {
171      
172      OBMol* newMol = new OBMol();
173      newMol->ReserveAtoms(fragment.size());
# Line 225 | Line 176 | namespace OpenBabel
176        OBAtom* newAtom = newMol->NewAtom();
177        *newAtom = *mol.GetAtom(*i);
178      }
179 +
180      newMol->EndModify();
181      newMol->ConnectTheDots();
182 <    findAngles(*newMol);
183 <    newMol->FindTorsions();
182 >    newMol->PerceiveBondOrders();
183 >
184      return newMol;
185    }
186    
187 <  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices) {
188 <    std::string indentLevel1("  ");
189 <    std::string indentLevel2("    ");
187 >  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols,
188 >                                ostream& os, OBMol& mol,
189 >                                vector<int>& indices) {
190 >    
191      std::string molPrefix("MolName");
192 +    std::string resName;
193      unsigned int i;
194      const int BUFFLEN = 1024;
195      char buffer[BUFFLEN];
196 <    string str, str1;
196 >    string str, str1, str2, str3;
197 >    OBAtom *a, *b, *c, *d;    
198 >    bool molIsWater = false;
199 >    OBResidue *r;
200 >    int resKey, myserial;
201 >    char type_name[10];
202 >    char *element_name;
203 >    int res_num;
204 >    OBChainsParser* chainParser = new OBChainsParser();  
205 >    double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
206      
244 –    
207      os << "<OOPSE version=4>" << endl;
208      os << "  <MetaData>" << endl << endl;
209      
210      for(i = 0; i < mols.size(); ++i) {
211        OBMol* pmol = mols[i];
250 –      
251 –      pmol->ConnectTheDots();
252 –      pmol->PerceiveBondOrders();
253 –      //pmol->FindSSSR();
254 –      //pmol->SetAromaticPerceived();
255 –      //pmol->Kekulize();
256 –      
212        map<OBAtom*, int> atomMap;
213 <      os << "molecule {\n";
214 <      sprintf(buffer, "%d", i);
215 <      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
213 >
214 >      chainParser->PerceiveChains(*pmol, false);
215 >      molIsWater = false;
216 >      FOR_RESIDUES_OF_MOL(residue, *pmol) {
217 >        std::cerr << "residue = " << residue->GetName() << "\n";
218 >        if (residue->GetName().compare("HOH") == 0) {
219 >          molIsWater = true;
220 >        }
221 >      }
222        
223 <      //atom
224 <      int ai = 0;
225 <      FOR_ATOMS_OF_MOL(atom, *pmol ) {
226 <        str = atom->GetType();
227 <        ttab.SetFromType("INT");
228 <        ttab.SetToType("INT");
229 <        ttab.Translate(str1,str);
230 <        os << indentLevel1 << "atom[" << ai << "] {\n";
231 <        // os << indentLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
232 <        os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
233 <        os << indentLevel1 << "}\n";
234 <        atomMap[&(*atom)] = ai++;
235 <      }        
236 <      os << "\n";
276 <      
277 <      //bond
278 <      FOR_BONDS_OF_MOL(bond, *pmol ) {
279 <        os << indentLevel1 << "bond {\n";
280 <        os << indentLevel2 << "members(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << ");\n";
281 <        os << indentLevel1 << "}\n";
282 <      }  
283 <      /*
284 <      //bend
285 <      OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
286 <      OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
287 <      vector<OBAngle> angles = pAngleData->GetData();
288 <      
289 <      os << indentLevel1 << "nBends = " << angles.size() << ";\n";        
290 <      int bendIndex = 0;
291 <      for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
292 <      {
293 <      triple<OBAtom*, OBAtom*, OBAtom*> bendAtoms = ti->getAtoms();
294 <      os << indentLevel1 << "bend[" << bendIndex++ << "] {\n";
295 <          os << indentLevel2 << "member(" << atomMap[bendAtoms.first] <<  ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
296 <          os << indentLevel1 << "}\n";            
297 <          }
223 >      if (molIsWater) {
224 >        // water include files define all of the known water types
225 >        os << "#include \"water.md\";\n";
226 >        pmol->SetTitle("HOH");
227 >      } else {
228 >
229 >        os << "molecule {\n";
230 >        sprintf(buffer, "%d", i);
231 >        os << "  name = \"" << molPrefix << buffer << "\";\n";
232 >        
233 >        int ai = 0;
234 >        FOR_ATOMS_OF_MOL(atom, *pmol ) {
235 >          str = atom->GetType();
236 >          r = atom->GetResidue();
237            
238 <          //torsion
239 <          pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
240 <          OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
241 <          vector<OBTorsion> torsions = pTorsionData->GetData();
303 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > torsionArray;
304 <          for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
305 <          {
306 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > tmpTorsions = ti->getTorsions();
307 <          torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());            
308 <          }
238 >          if (r == NULL)
239 >            resName = "NULL";
240 >          else
241 >            resName = r->GetName();
242            
243 <          os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
244 <          int torsionIndex = 0;
245 <          for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
246 <          {
247 <          os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
248 <          os << indentLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
249 <          os << indentLevel1 << "}\n";          
250 <          }
251 <      */
252 <      os << "}" << endl;
253 <      os << endl;
254 <      
243 >          if (resName.compare("NULL") ==0 ||
244 >              resName.compare("LIG") == 0 ||
245 >              resName.compare("UNK") == 0) {
246 >            // Either couldn't find a residue at all or couldn't find a
247 >            // reasonable residue name to use.  We'll punt and use
248 >            // OpenBabel's internal atom typing:
249 >            ttab.SetFromType("INT");
250 >            ttab.SetToType("INT");
251 >            ttab.Translate(str1, str);
252 >          } else {                        
253 >            
254 >            // If we know what residue we've got, the specific atom name can
255 >            // be used to help specify partial charges.
256 >
257 >            //resdat.SetResName(resName);
258 >            
259 >            // atom type from residue:
260 >            str = r->GetAtomID(&*atom);
261 >          
262 >            // arginine has separate indices for chemically-identical
263 >            // nitrogen atoms:
264 >            if (resName.compare("ARG") == 0) {
265 >              if (str.compare("NH1") == 0 || str.compare("NH2") == 0) {
266 >                str = "NH";
267 >              }
268 >            }
269 >            if (resName.compare("VAL") == 0) {
270 >              if (str.compare("CG1") == 0 || str.compare("CG2") == 0) {
271 >                str = "CG";
272 >              }
273 >            }
274 >            if (resName.compare("LEU") == 0) {
275 >              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
276 >                str = "CD";
277 >              }
278 >            }
279 >            if (resName.compare("ASP") == 0) {
280 >              if (str.compare("OD1") == 0 || str.compare("OD2") == 0) {
281 >                str = "OD";
282 >              }
283 >            }
284 >            if (resName.compare("GLU") == 0) {
285 >              if (str.compare("OE1") == 0 || str.compare("OE2") == 0) {
286 >                str = "OE";
287 >              }
288 >            }
289 >            if (resName.compare("TYR") == 0) {
290 >              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
291 >                str = "CD";
292 >              }
293 >              if (str.compare("CE1") == 0 || str.compare("CE2") == 0) {
294 >                str = "CE";
295 >              }
296 >            }
297 >            
298 >
299 >            if ((&*atom)->IsHydrogen()) {
300 >               FOR_NBORS_OF_ATOM(nbr, *atom) {
301 >                 str2 = r->GetAtomID(&*nbr);
302 >                 size_t startpos = str2.find_first_not_of(" ");
303 >                 size_t endpos = str2.find_last_not_of(" ");
304 >                 if ((endpos - startpos) < 1) {
305 >                   // if the bonded atom type has only one character (i.e. N)
306 >                   // then the hydrogen will be labeled "HN" to show what
307 >                   // kind of proton it is:
308 >                   str3 = str2;
309 >                 } else {
310 >                   if (str2.compare("OH") == 0) {
311 >                      str3 = "O";
312 >                   } else {
313 >                     // When the bonded atom type is more specific, we drop
314 >                     // the first character:  i.e. H bonded to OG1 is HG1 type:
315 >                     str3 = str2.substr(startpos+1, endpos-startpos);
316 >                   }
317 >                 }
318 >                str = "H" + str3;
319 >               }
320 >               // same problem with arginine NH atoms, but now for connected hydrogens
321 >               if (resName.compare("ARG") == 0) {
322 >                 if (str.compare("HH1") == 0 || str.compare("HH2") == 0) {
323 >                   str = "HH";
324 >                 }
325 >               }
326 >               if (resName.compare("VAL") == 0) {
327 >                 if (str.compare("HG1") == 0 || str.compare("HG2") == 0) {
328 >                   str = "HG";
329 >                 }
330 >               }
331 >               if (resName.compare("LEU") == 0) {
332 >                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
333 >                   str = "HD";
334 >                 }
335 >               }
336 >               if (resName.compare("TYR") == 0) {
337 >                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
338 >                   str = "HD";
339 >                 }
340 >                 if (str.compare("HE1") == 0 || str.compare("HE2") == 0) {
341 >                   str = "HE";
342 >                 }
343 >               }
344 >
345 >            }
346 >
347 >            // atom type from residue table:
348 >            //resdat.LookupType(str, str2, hyb);
349 >            size_t startpos = str.find_first_not_of(" ");
350 >            size_t endpos = str.find_last_not_of(" ");
351 >            str = str.substr( startpos, endpos-startpos+1 );
352 >            str1 = resName + "-" + str;
353 >          }      
354 >          os << "  atom[" << ai << "] { ";
355 >          os << "type = " << "\"" << str1 << "\"" << "; ";
356 >          os << "}\n";
357 >          atomMap[&(*atom)] = ai++;
358 >        }        
359 >        os << "\n";
360 >        
361 >        //bond
362 >        
363 >        int b1, b2;
364 >        FOR_BONDS_OF_MOL(bond, *pmol ) {
365 >          b1 = atomMap[bond->GetBeginAtom()];
366 >          b2 = atomMap[bond->GetEndAtom()];
367 >          
368 >          os << "  bond { ";
369 >          
370 >          if (b1 < b2)
371 >            os << "members(" << b1 <<  ", " << b2 << "); ";
372 >          else
373 >            os << "members(" << b2 <<  ", " << b1 << "); ";
374 >          
375 >          os << "}" << endl;
376 >        }
377 >        
378 >        os << endl;
379 >      
380 >        os << "}" << endl;
381 >        os << endl;
382 >      }
383      }
384      
385      os << endl;
386 <    
387 <    
388 <    for(i=0; i < mols.size(); ++i) {      
386 >        
387 >    for(i=0; i < mols.size(); ++i) {
388 >      OBMol* pmol = mols[i];      
389        os << "component{" << endl;
390 <      sprintf(buffer, "%d", i);
391 <      os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
392 <      os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
390 >      if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
391 >        os << "  type = " << "HOH" << ";" << endl;
392 >      } else {
393 >        sprintf(buffer, "%d", i);
394 >        os << "  type = " << molPrefix << buffer << ";" << endl;
395 >      }
396 >      os << "  nMol = " << numMols[i]<< ";" << endl;
397        os << "}" << endl;
398      }
399      
# Line 339 | Line 404 | namespace OpenBabel
404      sprintf(buffer, "        Time: %.10g", 0.0);
405      
406      os << buffer << endl;
407 +
408 +    CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);
409 +
410 +    // still to do: should compute a bounding box here
411 +    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}",
412 +            max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0, max_z - min_z);
413      
343 –    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}", 100.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 100.0);
344 –    
414      os << buffer << endl;
415      os << "    </FrameData>" << endl;
416      os << "    <StuntDoubles>" << endl;
417      
418      OBAtom *atom;
419      
420 +    // still to do: intercept waters and recompute pvqj lines
421 +
422      for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {    
423        atom = mol.GetAtom(*i);
424 <      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1, "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
424 >      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1,
425 >              "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
426        os << buffer << endl;
427      }
428      os << "    </StuntDoubles>" << endl;
429      os << "  </Snapshot>" << endl;
430      os << "</OOPSE>" << endl;
431    }
432 <  
432 >
433 >  void OOPSEFormat::CalcBoundingBox(OBMol &mol,
434 >                                    double &min_x, double &max_x,
435 >                                    double &min_y, double &max_y,
436 >                                    double &min_z, double &max_z
437 >                                    )
438 >  {
439 >    /* ---- Init bounding-box variables ---- */
440 >    min_x = (double) 0.0;
441 >    max_x = (double) 0.0;
442 >    min_y = (double) 0.0;
443 >    max_y = (double) 0.0;
444 >    min_z = (double) 0.0;
445 >    max_z = (double) 0.0;
446 >    
447 >    /* ---- Check all atoms ---- */
448 >    for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
449 >      {
450 >        
451 >        /* ---- Get a pointer to ith atom ---- */
452 >        OBAtom *atom = mol.GetAtom(i);
453 >        
454 >        /* ---- Check for minimal/maximal x-position ---- */
455 >        if (atom -> GetX() < min_x)
456 >          min_x = atom -> GetX();
457 >        if (atom -> GetX() > max_x)
458 >          max_x = atom -> GetX();
459 >        
460 >        /* ---- Check for minimal/maximal y-position ---- */
461 >        if (atom -> GetY() < min_y)
462 >          min_y = atom -> GetY();
463 >        if (atom -> GetY() > max_y)
464 >          max_y = atom -> GetY();
465 >        
466 >        /* ---- Check for minimal/maximal z-position ---- */
467 >        if (atom -> GetZ() < min_z)
468 >          min_z = atom -> GetZ();
469 >        if (atom -> GetZ() > max_z)
470 >          max_z = atom -> GetZ();
471 >        
472 >      }    
473 >  }
474   } //namespace OpenBabel
475  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines