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 3432 by gezelter, Mon Jul 14 12:35:58 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 <fstream>
23  
24 + #include "utils/StringUtils.hpp"
25 +
26   using namespace std;
27   namespace OpenBabel
28   {
# Line 48 | Line 53 | namespace OpenBabel
53        return NOTREADABLE | WRITEONEONLY;
54      }
55      
51 –    //*** This section identical for most OBMol conversions ***
52 –    ////////////////////////////////////////////////////
53 –    /// 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);
58 –    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 <  /////////////////////////////////////////////////////////////////
67 <
68 <
69 <  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
70 <  {
72 >
73 >  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
74      OBMol* pmol = dynamic_cast<OBMol*>(pOb);
75      if(pmol==NULL)
76        return false;
# Line 79 | 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;
85 <        }
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 <            {
100 <              sameMolTypes.push_back(j);
101 <              indices.insert(indices.end(), fragmentLists[j].begin(),
102 <                             fragmentLists[j].end());
103 <              used[j]=true;
104 <            }
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      
110 –    //
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 <    unsigned int pos = OutputFileName.rfind(".");
119 <    if(pos==string::npos)
124 <      OutputFileName += ".md";
125 <    else
118 >    size_t pos = OutputFileName.rfind(".");
119 >    if(pos!=string::npos)
120        OutputFileName = OutputFileName.substr(0, pos) + ".md";      
121 +    else
122 +      OutputFileName += ".md";
123 +    
124      ofstream ofs(OutputFileName.c_str());
125 <    if(!ofs)
129 <      {
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;
139 <      }
132 >    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) {
133 >      delete *i;
134 >    }
135      
141 –    //    
142 –    
136      return(true);
137    }
138    
139    bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
140 <                                     vector<int>& frag2)
148 <  {
140 >                                     vector<int>& frag2) {
141      if (frag1.size() != frag2.size())
142 <      {
151 <        return false;
152 <      }
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 <          }
165 <      }
153 >        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
154 >          return false;
155 >        
156 >    }
157      return true;
158    }
159    
160 <  struct SameAngle
170 <  {
160 >  struct SameAngle {
161      bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
162 <                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
173 <    {
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    };
177 –  
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;
166  
197 –        for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
198 –        {
199 –            if(a == c)
200 –                continue;          
201 –            
202 –            uniqueAngles.insert(triple<OBAtom*,OBAtom*,OBAtom*>(a, b, c));
203 –        }
204 –    }
167  
168 <    //get new data and attach it to molecule
169 <    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) {
168 >  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol,
169 >                                            vector<int>& fragment) {
170      
171      OBMol* newMol = new OBMol();
172      newMol->ReserveAtoms(fragment.size());
# Line 225 | Line 175 | namespace OpenBabel
175        OBAtom* newAtom = newMol->NewAtom();
176        *newAtom = *mol.GetAtom(*i);
177      }
178 +
179      newMol->EndModify();
180      newMol->ConnectTheDots();
181 <    findAngles(*newMol);
182 <    newMol->FindTorsions();
181 >    newMol->PerceiveBondOrders();
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;
193      const int BUFFLEN = 1024;
194      char buffer[BUFFLEN];
195      string str, str1;
196 +    OBAtom *a, *b, *c, *d;    
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      
244 –    
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];
250 –      
251 –      pmol->ConnectTheDots();
252 –      pmol->PerceiveBondOrders();
253 –      //pmol->FindSSSR();
254 –      //pmol->SetAromaticPerceived();
255 –      //pmol->Kekulize();
256 –      
211        map<OBAtom*, int> atomMap;
212 <      os << "molecule {\n";
213 <      sprintf(buffer, "%d", i);
214 <      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
212 >
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 <      //atom
223 <      int ai = 0;
224 <      FOR_ATOMS_OF_MOL(atom, *pmol ) {
225 <        str = atom->GetType();
226 <        ttab.SetFromType("INT");
227 <        ttab.SetToType("INT");
228 <        ttab.Translate(str1,str);
229 <        os << indentLevel1 << "atom[" << ai << "] {\n";
230 <        // os << indentLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
231 <        os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
232 <        os << indentLevel1 << "}\n";
233 <        atomMap[&(*atom)] = ai++;
234 <      }        
235 <      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 <          }
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 >        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();
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 <          }
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 339 | 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      
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 –    
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