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 3446 by cli2, Wed Sep 10 19:51:45 2008 UTC

# Line 19 | Line 19 | GNU General Public License for more details.
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"
# Line 53 | Line 54 | namespace OpenBabel
54        return NOTREADABLE | WRITEONEONLY;
55      }
56      
56 –    //*** This section identical for most OBMol conversions ***
57 –    ////////////////////////////////////////////////////
58 –    /// 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);
63 –    //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 <  /////////////////////////////////////////////////////////////////
72 <
73 <
74 <  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
75 <  {
73 >
74 >  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
75      OBMol* pmol = dynamic_cast<OBMol*>(pOb);
76      if(pmol==NULL)
77        return false;
# Line 84 | 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;
90 <        }
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 <            {
105 <              sameMolTypes.push_back(j);
106 <              indices.insert(indices.end(), fragmentLists[j].begin(),
107 <                             fragmentLists[j].end());
108 <              used[j]=true;
109 <            }
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      
115 –    //
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      size_t pos = OutputFileName.rfind(".");
# Line 129 | Line 121 | namespace OpenBabel
121        OutputFileName = OutputFileName.substr(0, pos) + ".md";      
122      else
123        OutputFileName += ".md";
124 <
124 >    
125      ofstream ofs(OutputFileName.c_str());
126 <    if(!ofs)
135 <      {
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;
145 <      }
133 >    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) {
134 >      delete *i;
135 >    }
136      
147 –    //    
148 –    
137      return(true);
138    }
139    
140    bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
141 <                                     vector<int>& frag2)
154 <  {
141 >                                     vector<int>& frag2) {
142      if (frag1.size() != frag2.size())
143 <      {
157 <        return false;
158 <      }
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 <          }
171 <      }
154 >        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
155 >          return false;
156 >        
157 >    }
158      return true;
159    }
160    
161 <  struct SameAngle
176 <  {
161 >  struct SameAngle {
162      bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
163 <                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
179 <    {
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  
184 –  /*  
185 –    void OOPSEFormat::findAngles(OBMol& mol) {
168  
169 <    //if already has data return
170 <    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) {
169 >  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol,
170 >                                            vector<int>& fragment) {
171      
172      OBMol* newMol = new OBMol();
230 –    OBChainsParser* chainParser = new OBChainsParser();
173      newMol->ReserveAtoms(fragment.size());
174      newMol->BeginModify();
175      for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
176        OBAtom* newAtom = newMol->NewAtom();
177        *newAtom = *mol.GetAtom(*i);
178      }
179 +
180      newMol->EndModify();
181      newMol->ConnectTheDots();
182      newMol->PerceiveBondOrders();
240 –    newMol->FindAngles();
241 –    newMol->FindTorsions();
242 –    //chainParser->PerceiveChains(*newMol, false);
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 <    OBChainsParser* chainParser = new OBChainsParser();
198 >    bool molIsWater = false;
199      OBResidue *r;
200 <    int resKey;
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      
264 –    std::cerr << "yo\n";
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];
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 –      
212        map<OBAtom*, int> atomMap;
279 –      os << "molecule {\n";
280 –      sprintf(buffer, "%d", i);
281 –      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
282 –      
213  
214 <      //atom
215 <      
216 <      int ai = 0;
217 <      FOR_RESIDUES_OF_MOL(res, *pmol) {
218 <
219 <        std::cerr << "found residue" << res->GetName() << "\n";
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 < //       FOR_RESIDUES_OF_MOL(res, *pmol) {
224 <        
225 < //      resName = res->GetName();
226 < //      resKey = res->GetResKey();
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 <        
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 >          if (r == NULL)
239 >            resName = "NULL";
240 >          else
241 >            resName = r->GetName();
242 >          
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 < //      std::cerr << "found residue " << resName << "\n";
258 < //      std::cerr << "  num = " << res->GetNum() << "\n";
259 < //      std::cerr << "  numAtoms = " << res->GetNumAtoms() << "\n";
260 < //      std::cerr << "  num = " << res->GetNum() << "\n";
261 < //      std::cerr << "  chain = " << res->GetChain() << "\n";
262 < //      std::cerr << "  chainnum = " << res->GetChainNum() << "\n";
263 < //      std::cerr << "  idx = " << res->GetIdx() << "\n";
264 < //      std::cerr << "  key = " << res->GetResKey() << "\n";
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 < //      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 < //       }      
345 >            }
346  
347 <      ai = 0;
348 <      FOR_ATOMS_OF_MOL(atom, *pmol ) {
349 <        str = atom->GetType();
350 <        r = atom->GetResidue();
351 <        
352 <        if (r == NULL)
353 <          resName = "NULL";
354 <        else
355 <          resName = r->GetName();
356 <                
357 <        if (resName.compare("NULL") ==0 || resName.compare("LIG") == 0 || resName.compare("UNK") == 0) {
358 <          // Either couldn't find a residue at all or couldn't find a reasonable
359 <          // residue name to use.  We'll punt and use OpenBabel's internal atom typing:
360 <          ttab.SetFromType("INT");
361 <          ttab.SetToType("INT");
362 <          ttab.Translate(str1,str);
363 <        } else {
364 <                
365 <
366 <          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 <          }
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 <          //torsion
440 <          pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
441 <          OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
442 <          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 <          }
368 >          os << "  bond { ";
369            
370 <          os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
371 <          int torsionIndex = 0;
372 <          for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
373 <          {
374 <          os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
375 <          os << indentLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
376 <          os << indentLevel1 << "}\n";          
377 <          }
378 <      */
379 <      os << "}" << endl;
380 <      os << endl;
381 <      
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 479 | 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      
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 –    
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