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 3424 by gezelter, Tue Jul 1 13:28:23 2008 UTC vs.
Revision 3432 by gezelter, Mon Jul 14 12:35:58 2008 UTC

# Line 53 | Line 53 | namespace OpenBabel
53        return NOTREADABLE | WRITEONEONLY;
54      }
55      
56 –    //*** This section identical for most OBMol conversions ***
57 –    ////////////////////////////////////////////////////
58 –    /// The "API" interface functions
56      virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
57  
58    private:
59      bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
63 –    //void findAngles(OBMol& mol);
60      OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
61 <    void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices);
61 >    void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os,
62 >                     OBMol& mol, vector<int>& indices);
63 >    void CalcBoundingBox(OBMol &mol,
64 >                         double &min_x, double &max_x,
65 >                         double &min_y, double &max_y,
66 >                         double &min_z, double &max_z);
67 >      
68    };
69 <
69 >  
70    //Make an instance of the format class
71    OOPSEFormat theOOPSEFormat;
72 <  
73 <  /////////////////////////////////////////////////////////////////
72 <
73 <
74 <  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
75 <  {
72 >
73 >  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
74      OBMol* pmol = dynamic_cast<OBMol*>(pOb);
75      if(pmol==NULL)
76        return false;
# Line 84 | Line 82 | namespace OpenBabel
82      vector<vector<int> > molecules;
83      vector<int> indices;
84      for(int i =0; i < used.size(); ++i) {
85 <      if (used[i])
86 <        {
87 <          continue;
90 <        }
85 >      if (used[i])
86 >        continue;
87 >
88        used[i] = true;
89        vector<int> sameMolTypes;
90        sameMolTypes.push_back(i);
91        indices.insert(indices.end(), fragmentLists[i].begin(),
92                       fragmentLists[i].end());
93 <      for (int j = i + 1;j < used.size(); ++j)
94 <        {
95 <          if (used[j])
96 <            {
97 <              continue;
98 <            }
99 <          
100 <          if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
101 <            {
105 <              sameMolTypes.push_back(j);
106 <              indices.insert(indices.end(), fragmentLists[j].begin(),
107 <                             fragmentLists[j].end());
108 <              used[j]=true;
109 <            }
93 >      for (int j = i + 1;j < used.size(); ++j) {
94 >        if (used[j])
95 >          continue;
96 >                
97 >        if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j])) {
98 >          sameMolTypes.push_back(j);
99 >          indices.insert(indices.end(), fragmentLists[j].begin(),
100 >                         fragmentLists[j].end());
101 >          used[j]=true;
102          }
103 <      molecules.push_back(sameMolTypes);
104 <      
103 >      }
104 >      molecules.push_back(sameMolTypes);      
105      }
106      
115 –    //
107      vector<OBMol*> mdMols;    
108      vector<int> numMols;
109      for(vector<vector<int> >::iterator  i = molecules.begin();
110 <        i != molecules.end(); ++i)
111 <      {
112 <        mdMols.push_back(createMolFromFragment(*pmol,
113 <                                               fragmentLists[i->front()]));
114 <        numMols.push_back((*i).size());
115 <      }
110 >        i != molecules.end(); ++i) {
111 >      
112 >      mdMols.push_back(createMolFromFragment(*pmol,
113 >                                             fragmentLists[i->front()]));
114 >      numMols.push_back((*i).size());
115 >    }
116      
117      string OutputFileName = pConv->GetInFilename();
118      size_t pos = OutputFileName.rfind(".");
# Line 129 | Line 120 | namespace OpenBabel
120        OutputFileName = OutputFileName.substr(0, pos) + ".md";      
121      else
122        OutputFileName += ".md";
123 <
123 >    
124      ofstream ofs(OutputFileName.c_str());
125 <    if(!ofs)
135 <      {
125 >    if(!ofs) {
126          cerr << "Cannot write to " << OutputFileName <<endl;
127          return false;
128 <      }
128 >    }
129      
130      WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
131      
132 <    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i)
133 <      {
134 <        delete *i;
145 <      }
132 >    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) {
133 >      delete *i;
134 >    }
135      
147 –    //    
148 –    
136      return(true);
137    }
138    
139    bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
140 <                                     vector<int>& frag2)
154 <  {
140 >                                     vector<int>& frag2) {
141      if (frag1.size() != frag2.size())
142 <      {
157 <        return false;
158 <      }
142 >      return false;
143      
144 <    //exact graph matching is a NP complete problem
145 <    /** @todo using sparse matrix to store the connectivities*/
146 <    for (unsigned int i =0 ; i < frag1.size(); ++i)
147 <      {
144 >    // Exact graph matching is an NP complete problem.
145 >    // This just matches all of the atom atomic numbers and may falsely
146 >    // detect identical fragments which aren't really identical.
147 >    // @todo using sparse matrix to store the connectivities
148 >
149 >    for (unsigned int i =0 ; i < frag1.size(); ++i) {
150          OBAtom* atom1 = mol.GetAtom(frag1[i]);
151          OBAtom* atom2 = mol.GetAtom(frag2[i]);
152          
153 <        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
154 <          {
155 <            return false;
156 <          }
171 <      }
153 >        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
154 >          return false;
155 >        
156 >    }
157      return true;
158    }
159    
160 <  struct SameAngle
176 <  {
160 >  struct SameAngle {
161      bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
162 <                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
179 <    {
162 >                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const {
163        return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) || (t1.first == t2.third && t1.third == t2.first));
164      }
165    };
166  
184 –  /*  
185 –    void OOPSEFormat::findAngles(OBMol& mol) {
167  
168 <    //if already has data return
169 <    if(mol.HasData(OBGenericDataType::AngleData))
189 <    return;
190 <    
191 <    vector<OBEdgeBase*>::iterator bi1,bi2;
192 <    OBBond* bond;
193 <    OBAtom *a,*b,*c;
194 <    
195 <    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle> uniqueAngles;
196 <    //loop through all bonds generating torsions
197 <    for(bond = mol.BeginBond(bi1);bond;bond = mol.NextBond(bi1))
198 <    {
199 <        b = bond->GetBeginAtom();
200 <        c = bond->GetEndAtom();
201 <        if(b->IsHydrogen())
202 <            continue;
203 <
204 <        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
205 <        {
206 <            if(a == c)
207 <                continue;          
208 <            
209 <            uniqueAngles.insert(triple<OBAtom*,OBAtom*,OBAtom*>(a, b, c));
210 <        }
211 <    }
212 <
213 <    //get new data and attach it to molecule
214 <    OBAngleData *angles = new OBAngleData;
215 <    mol.SetData(angles);
216 <    set<triple<OBAtom*,OBAtom*,OBAtom*>, SameAngle>::iterator i;
217 <
218 <    for (i = uniqueAngles.begin(); i != uniqueAngles.end(); ++i) {
219 <        OBAngle angle;
220 <        angle.SetAtoms(i->first, i->second, i->second);
221 <        angles->SetData(angle);
222 <    }
223 <
224 <  }
225 <  */
226 <  
227 <  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment) {
168 >  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol,
169 >                                            vector<int>& fragment) {
170      
171      OBMol* newMol = new OBMol();
230 –    OBChainsParser* chainParser = new OBChainsParser();
172      newMol->ReserveAtoms(fragment.size());
173      newMol->BeginModify();
174      for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
175        OBAtom* newAtom = newMol->NewAtom();
176        *newAtom = *mol.GetAtom(*i);
177      }
178 +
179      newMol->EndModify();
180      newMol->ConnectTheDots();
181      newMol->PerceiveBondOrders();
240 –    newMol->FindAngles();
241 –    newMol->FindTorsions();
242 –    //chainParser->PerceiveChains(*newMol, false);
182  
183      return newMol;
184    }
185    
186 <  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices) {
187 <    std::string indentLevel1("  ");
188 <    std::string indentLevel2("    ");
186 >  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols,
187 >                                ostream& os, OBMol& mol,
188 >                                vector<int>& indices) {
189 >    
190      std::string molPrefix("MolName");
191      std::string resName;
192      unsigned int i;
# Line 254 | Line 194 | namespace OpenBabel
194      char buffer[BUFFLEN];
195      string str, str1;
196      OBAtom *a, *b, *c, *d;    
197 <    OBChainsParser* chainParser = new OBChainsParser();
197 >    bool molIsWater = false;
198      OBResidue *r;
199      int resKey;
200      char type_name[10];
201      char *element_name;
202      int res_num;
203 +    OBChainsParser* chainParser = new OBChainsParser();  
204 +    double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
205      
264 –    std::cerr << "yo\n";
206      os << "<OOPSE version=4>" << endl;
207      os << "  <MetaData>" << endl << endl;
208      
209      for(i = 0; i < mols.size(); ++i) {
210        OBMol* pmol = mols[i];
270 –      std::cerr << "yo1\n";
271 –      pmol->ConnectTheDots();
272 –      pmol->PerceiveBondOrders();
273 –      chainParser->PerceiveChains(*pmol, false);
274 –      //pmol->FindSSSR();
275 –      //pmol->SetAromaticPerceived();
276 –      //pmol->Kekulize();
277 –      
211        map<OBAtom*, int> atomMap;
279 –      os << "molecule {\n";
280 –      sprintf(buffer, "%d", i);
281 –      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
282 –      
212  
213 <      //atom
214 <      
215 <      int ai = 0;
216 <      FOR_RESIDUES_OF_MOL(res, *pmol) {
217 <
218 <        std::cerr << "found residue" << res->GetName() << "\n";
213 >      chainParser->PerceiveChains(*pmol, false);
214 >      molIsWater = false;
215 >      FOR_RESIDUES_OF_MOL(residue, *pmol) {
216 >        std::cerr << "residue = " << residue->GetName() << "\n";
217 >        if (residue->GetName().compare("HOH") == 0) {
218 >          molIsWater = true;
219 >        }
220        }
221        
222 < //       FOR_RESIDUES_OF_MOL(res, *pmol) {
223 <        
224 < //      resName = res->GetName();
225 < //      resKey = res->GetResKey();
222 >      if (molIsWater) {
223 >        // water include files define all of the known water types
224 >        os << "#include \"water.md\";\n";
225 >        pmol->SetTitle("HOH");
226 >      } else {
227  
228 <        
229 <
230 < //      std::cerr << "found residue " << resName << "\n";
231 < //      std::cerr << "  num = " << res->GetNum() << "\n";
232 < //      std::cerr << "  numAtoms = " << res->GetNumAtoms() << "\n";
233 < //      std::cerr << "  num = " << res->GetNum() << "\n";
234 < //      std::cerr << "  chain = " << res->GetChain() << "\n";
235 < //      std::cerr << "  chainnum = " << res->GetChainNum() << "\n";
305 < //      std::cerr << "  idx = " << res->GetIdx() << "\n";
306 < //      std::cerr << "  key = " << res->GetResKey() << "\n";
307 <
308 <
309 < //      FOR_ATOMS_OF_RESIDUE(atom, &*res) {
310 < //        str = atom->GetType();
311 < //        ttab.SetFromType("INT");
312 < //        ttab.SetToType("INT");
313 < //        ttab.Translate(str1,str);
314 < //        os << indentLevel1 << "atom[" << ai << "] {\n";
315 < //        os << indentLevel2 << "type = " << "\"" << resName << "-" << str1 << "\"" << ";\n";
316 < //        os << indentLevel1 << "}\n";
317 < //        atomMap[&(*atom)] = ai++;
318 < //      }        
319 < //      os << "\n";
320 <        
321 < //       }      
322 <
323 <      ai = 0;
324 <      FOR_ATOMS_OF_MOL(atom, *pmol ) {
325 <        str = atom->GetType();
326 <        r = atom->GetResidue();
327 <        
328 <        if (r == NULL)
329 <          resName = "NULL";
330 <        else
331 <          resName = r->GetName();
332 <                
333 <        if (resName.compare("NULL") ==0 || resName.compare("LIG") == 0 || resName.compare("UNK") == 0) {
334 <          // Either couldn't find a residue at all or couldn't find a reasonable
335 <          // residue name to use.  We'll punt and use OpenBabel's internal atom typing:
336 <          ttab.SetFromType("INT");
337 <          ttab.SetToType("INT");
338 <          ttab.Translate(str1,str);
339 <        } else {
340 <                
341 <
342 <          if (resName.compare("HOH") == 0) {
343 <            // HOH is a special residue name for water, and the standard atom types
344 <            // are OW and HW, so just append W to the string for the atom type:
345 <            ttab.SetFromType("INT");
346 <            ttab.SetToType("XYZ");
347 <            ttab.Translate(str1,str);
348 <            str1 += "W";
349 <          } else {
350 <
351 <            std::cerr << "found residue " << resName << "\n";
352 <
353 <            // If we know what residue we've got, the specific atom name can
354 <            // be used to help specify partial charges.
355 <
356 <            str = r->GetAtomID(&*atom);
357 <            size_t startpos = str.find_first_not_of(" ");
358 <            size_t endpos = str.find_last_not_of(" ");
359 <            str = str.substr( startpos, endpos-startpos+1 );
360 <            str1 = resName + "-" + str;
361 <          }
362 <
363 <        }
364 <
365 < //         os << indentLevel1 << "atom[" << ai << "] {\n";
366 < //         os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
367 < //         os << indentLevel1 << "}\n";
368 <        os << "atom[" << ai << "] { ";
369 <        os << "type = " << "\"" << str1 << "\"" << "; ";
370 <        os << "}\n";
371 <        atomMap[&(*atom)] = ai++;
372 <      }        
373 <      os << "\n";
374 <      
375 <      //bond
376 <      FOR_BONDS_OF_MOL(bond, *pmol ) {
377 < //         os << indentLevel1 << "bond {\n";
378 < //         os << indentLevel2 << "members(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << ");\n";
379 < //         os << indentLevel1 << "}\n";
380 <        os << "bond { ";
381 <        os << "members(" << atomMap[bond->GetBeginAtom()] <<  ", " << atomMap[bond->GetEndAtom()] << "); ";
382 <        os << "}\n";
383 <      }
384 <
385 <      FOR_ANGLES_OF_MOL(angle, *pmol) {
386 <
387 <        // OpenBabel's getAtoms returns the 3 atom pointer for the
388 <        // angle with the vertex first.  These need to be reordered
389 <        // for vertex-second ordering for OOPSE.
390 <
391 <         b = pmol->GetAtom((*angle)[0] + 1);
392 <         a = pmol->GetAtom((*angle)[1] + 1);
393 <         c = pmol->GetAtom((*angle)[2] + 1);
394 <
395 < //       os << indentLevel1 << "bend {\n";
396 < //       os << indentLevel2 << "members(" << atomMap[a] <<  ", " << atomMap[b] << ", " << atomMap[c] << ");\n";
397 < //       os << indentLevel1 << "}\n";
398 <         os << "bend { ";
399 <         os << "members(" << atomMap[a] <<  ", " << atomMap[b] << ", " << atomMap[c] << "); ";
400 <         os << "}\n";
401 <      }
402 <
403 <      FOR_TORSIONS_OF_MOL(torsion, *pmol) {
404 <
405 <        // OpenBabel's getAtoms returns the 3 atom pointer for the
406 <        // angle with the vertex first.  These need to be reordered
407 <        // for vertex-second ordering for OOPSE.
408 <
409 <         a = pmol->GetAtom((*torsion)[0] + 1);
410 <         b = pmol->GetAtom((*torsion)[1] + 1);
411 <         c = pmol->GetAtom((*torsion)[2] + 1);
412 <         d = pmol->GetAtom((*torsion)[3] + 1);
413 <
414 < //         os << indentLevel1 << "torsion {\n";
415 < //         os << indentLevel2 << "members(" << atomMap[a] <<  ", " << atomMap[b] << ", " << atomMap[c] << ", " << atomMap[d] << ");\n";
416 < //         os << indentLevel1 << "}\n";
417 <
418 <        os << "torsion { ";
419 <        os << "members(" << atomMap[a] <<  ", " << atomMap[b] << ", " << atomMap[c] << ", " << atomMap[d] << "); ";
420 <        os << "}\n";
421 <      }
422 <
423 <      /*
424 <      //bend
425 <      OBGenericData* pGenericData = pmol->GetData(OBGenericDataType::AngleData);
426 <      OBAngleData* pAngleData = dynamic_cast<OBAngleData*>(pGenericData);
427 <      vector<OBAngle> angles = pAngleData->GetData();
428 <      
429 <      os << indentLevel1 << "nBends = " << angles.size() << ";\n";        
430 <      int bendIndex = 0;
431 <      for (vector<OBAngle>::iterator ti = angles.begin(); ti != angles.end(); ++ti)
432 <      {
433 <      triple<OBAtom*, OBAtom*, OBAtom*> bendAtoms = ti->getAtoms();
434 <      os << indentLevel1 << "bend[" << bendIndex++ << "] {\n";
435 <          os << indentLevel2 << "member(" << atomMap[bendAtoms.first] <<  ", " << atomMap[bendAtoms.second] << atomMap[bendAtoms.third] <<");\n";
436 <          os << indentLevel1 << "}\n";            
437 <          }
228 >        os << "molecule {\n";
229 >        sprintf(buffer, "%d", i);
230 >        os << "  name = \"" << molPrefix << buffer << "\";\n";
231 >        
232 >        int ai = 0;
233 >        FOR_ATOMS_OF_MOL(atom, *pmol ) {
234 >          str = atom->GetType();
235 >          r = atom->GetResidue();
236            
237 <          //torsion
238 <          pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
239 <          OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
240 <          vector<OBTorsion> torsions = pTorsionData->GetData();
443 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > torsionArray;
444 <          for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
445 <          {
446 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > tmpTorsions = ti->getTorsions();
447 <          torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());            
448 <          }
237 >          if (r == NULL)
238 >            resName = "NULL";
239 >          else
240 >            resName = r->GetName();
241            
242 <          os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
243 <          int torsionIndex = 0;
244 <          for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
245 <          {
246 <          os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
247 <          os << indentLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
248 <          os << indentLevel1 << "}\n";          
249 <          }
250 <      */
251 <      os << "}" << endl;
252 <      os << endl;
253 <      
242 >          if (resName.compare("NULL") ==0 ||
243 >              resName.compare("LIG") == 0 ||
244 >              resName.compare("UNK") == 0) {
245 >            // Either couldn't find a residue at all or couldn't find a
246 >            // reasonable residue name to use.  We'll punt and use
247 >            // OpenBabel's internal atom typing:
248 >            ttab.SetFromType("INT");
249 >            ttab.SetToType("INT");
250 >            ttab.Translate(str1, str);
251 >          } else {                        
252 >            
253 >            // If we know what residue we've got, the specific atom name can
254 >            // be used to help specify partial charges.
255 >            
256 >            str = r->GetAtomID(&*atom);
257 >            size_t startpos = str.find_first_not_of(" ");
258 >            size_t endpos = str.find_last_not_of(" ");
259 >            str = str.substr( startpos, endpos-startpos+1 );
260 >            str1 = resName + "-" + str;
261 >          }      
262 >          os << "  atom[" << ai << "] { ";
263 >          os << "type = " << "\"" << str1 << "\"" << "; ";
264 >          os << "}\n";
265 >          atomMap[&(*atom)] = ai++;
266 >        }        
267 >        os << "\n";
268 >        
269 >        //bond
270 >        
271 >        int b1, b2;
272 >        FOR_BONDS_OF_MOL(bond, *pmol ) {
273 >          b1 = atomMap[bond->GetBeginAtom()];
274 >          b2 = atomMap[bond->GetEndAtom()];
275 >          
276 >          os << "  bond { ";
277 >          
278 >          if (b1 < b2)
279 >            os << "members(" << b1 <<  ", " << b2 << "); ";
280 >          else
281 >            os << "members(" << b2 <<  ", " << b1 << "); ";
282 >          
283 >          os << "}" << endl;
284 >        }
285 >        
286 >        os << endl;
287 >      
288 >        os << "}" << endl;
289 >        os << endl;
290 >      }
291      }
292      
293      os << endl;
294 <    
295 <    
296 <    for(i=0; i < mols.size(); ++i) {      
294 >        
295 >    for(i=0; i < mols.size(); ++i) {
296 >      OBMol* pmol = mols[i];      
297        os << "component{" << endl;
298 <      sprintf(buffer, "%d", i);
299 <      os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
300 <      os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
298 >      if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
299 >        os << "  type = " << "HOH" << ";" << endl;
300 >      } else {
301 >        sprintf(buffer, "%d", i);
302 >        os << "  type = " << molPrefix << buffer << ";" << endl;
303 >      }
304 >      os << "  nMol = " << numMols[i]<< ";" << endl;
305        os << "}" << endl;
306      }
307      
# Line 479 | Line 312 | namespace OpenBabel
312      sprintf(buffer, "        Time: %.10g", 0.0);
313      
314      os << buffer << endl;
315 +
316 +    CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);
317 +
318 +    // still to do: should compute a bounding box here
319 +    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}",
320 +            max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0, max_z - min_z);
321      
483 –    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);
484 –    
322      os << buffer << endl;
323      os << "    </FrameData>" << endl;
324      os << "    <StuntDoubles>" << endl;
325      
326      OBAtom *atom;
327      
328 +    // still to do: intercept waters and recompute pvqj lines
329 +
330      for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {    
331        atom = mol.GetAtom(*i);
332 <      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);
332 >      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1,
333 >              "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
334        os << buffer << endl;
335      }
336      os << "    </StuntDoubles>" << endl;
337      os << "  </Snapshot>" << endl;
338      os << "</OOPSE>" << endl;
339    }
340 <  
340 >
341 >  void OOPSEFormat::CalcBoundingBox(OBMol &mol,
342 >                                    double &min_x, double &max_x,
343 >                                    double &min_y, double &max_y,
344 >                                    double &min_z, double &max_z
345 >                                    )
346 >  {
347 >    /* ---- Init bounding-box variables ---- */
348 >    min_x = (double) 0.0;
349 >    max_x = (double) 0.0;
350 >    min_y = (double) 0.0;
351 >    max_y = (double) 0.0;
352 >    min_z = (double) 0.0;
353 >    max_z = (double) 0.0;
354 >    
355 >    /* ---- Check all atoms ---- */
356 >    for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
357 >      {
358 >        
359 >        /* ---- Get a pointer to ith atom ---- */
360 >        OBAtom *atom = mol.GetAtom(i);
361 >        
362 >        /* ---- Check for minimal/maximal x-position ---- */
363 >        if (atom -> GetX() < min_x)
364 >          min_x = atom -> GetX();
365 >        if (atom -> GetX() > max_x)
366 >          max_x = atom -> GetX();
367 >        
368 >        /* ---- Check for minimal/maximal y-position ---- */
369 >        if (atom -> GetY() < min_y)
370 >          min_y = atom -> GetY();
371 >        if (atom -> GetY() > max_y)
372 >          max_y = atom -> GetY();
373 >        
374 >        /* ---- Check for minimal/maximal z-position ---- */
375 >        if (atom -> GetZ() < min_z)
376 >          min_z = atom -> GetZ();
377 >        if (atom -> GetZ() > max_z)
378 >          max_z = atom -> GetZ();
379 >        
380 >      }    
381 >  }
382   } //namespace OpenBabel
383  

Diff Legend

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