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 3429 by gezelter, Thu Jul 3 12:45:30 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    };
64 <
64 >  
65    //Make an instance of the format class
66    OOPSEFormat theOOPSEFormat;
67 <  
68 <  /////////////////////////////////////////////////////////////////
67 <
68 <
69 <  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
70 <  {
67 >
68 >  bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
69      OBMol* pmol = dynamic_cast<OBMol*>(pOb);
70      if(pmol==NULL)
71        return false;
# Line 79 | Line 77 | namespace OpenBabel
77      vector<vector<int> > molecules;
78      vector<int> indices;
79      for(int i =0; i < used.size(); ++i) {
80 <      if (used[i])
81 <        {
82 <          continue;
85 <        }
80 >      if (used[i])
81 >        continue;
82 >
83        used[i] = true;
84        vector<int> sameMolTypes;
85        sameMolTypes.push_back(i);
86        indices.insert(indices.end(), fragmentLists[i].begin(),
87                       fragmentLists[i].end());
88 <      for (int j = i + 1;j < used.size(); ++j)
89 <        {
90 <          if (used[j])
91 <            {
92 <              continue;
93 <            }
94 <          
95 <          if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
96 <            {
100 <              sameMolTypes.push_back(j);
101 <              indices.insert(indices.end(), fragmentLists[j].begin(),
102 <                             fragmentLists[j].end());
103 <              used[j]=true;
104 <            }
88 >      for (int j = i + 1;j < used.size(); ++j) {
89 >        if (used[j])
90 >          continue;
91 >                
92 >        if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j])) {
93 >          sameMolTypes.push_back(j);
94 >          indices.insert(indices.end(), fragmentLists[j].begin(),
95 >                         fragmentLists[j].end());
96 >          used[j]=true;
97          }
98 <      molecules.push_back(sameMolTypes);
99 <      
98 >      }
99 >      molecules.push_back(sameMolTypes);      
100      }
101      
110 –    //
102      vector<OBMol*> mdMols;    
103      vector<int> numMols;
104      for(vector<vector<int> >::iterator  i = molecules.begin();
105 <        i != molecules.end(); ++i)
106 <      {
107 <        mdMols.push_back(createMolFromFragment(*pmol,
108 <                                               fragmentLists[i->front()]));
109 <        numMols.push_back((*i).size());
110 <      }
105 >        i != molecules.end(); ++i) {
106 >      
107 >      mdMols.push_back(createMolFromFragment(*pmol,
108 >                                             fragmentLists[i->front()]));
109 >      numMols.push_back((*i).size());
110 >    }
111      
112      string OutputFileName = pConv->GetInFilename();
113 <    unsigned int pos = OutputFileName.rfind(".");
114 <    if(pos==string::npos)
124 <      OutputFileName += ".md";
125 <    else
113 >    size_t pos = OutputFileName.rfind(".");
114 >    if(pos!=string::npos)
115        OutputFileName = OutputFileName.substr(0, pos) + ".md";      
116 +    else
117 +      OutputFileName += ".md";
118 +    
119      ofstream ofs(OutputFileName.c_str());
120 <    if(!ofs)
129 <      {
120 >    if(!ofs) {
121          cerr << "Cannot write to " << OutputFileName <<endl;
122          return false;
123 <      }
123 >    }
124      
125      WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
126      
127 <    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i)
128 <      {
129 <        delete *i;
139 <      }
127 >    for(vector<OBMol*>::iterator  i = mdMols.begin(); i != mdMols.end(); ++i) {
128 >      delete *i;
129 >    }
130      
141 –    //    
142 –    
131      return(true);
132    }
133    
134    bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
135 <                                     vector<int>& frag2)
148 <  {
135 >                                     vector<int>& frag2) {
136      if (frag1.size() != frag2.size())
137 <      {
151 <        return false;
152 <      }
137 >      return false;
138      
139 <    //exact graph matching is a NP complete problem
140 <    /** @todo using sparse matrix to store the connectivities*/
141 <    for (unsigned int i =0 ; i < frag1.size(); ++i)
142 <      {
139 >    // Exact graph matching is an NP complete problem.
140 >    // This just matches all of the atom atomic numbers and may falsely
141 >    // detect identical fragments which aren't really identical.
142 >    // @todo using sparse matrix to store the connectivities
143 >
144 >    for (unsigned int i =0 ; i < frag1.size(); ++i) {
145          OBAtom* atom1 = mol.GetAtom(frag1[i]);
146          OBAtom* atom2 = mol.GetAtom(frag2[i]);
147          
148 <        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
149 <          {
150 <            return false;
151 <          }
165 <      }
148 >        if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
149 >          return false;
150 >        
151 >    }
152      return true;
153    }
154    
155 <  struct SameAngle
170 <  {
155 >  struct SameAngle {
156      bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
157 <                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
173 <    {
157 >                    const triple<OBAtom*,OBAtom*,OBAtom*> t2) const {
158        return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) || (t1.first == t2.third && t1.third == t2.first));
159      }
160    };
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;
161  
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 –    }
162  
163 <    //get new data and attach it to molecule
164 <    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) {
163 >  OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol,
164 >                                            vector<int>& fragment) {
165      
166      OBMol* newMol = new OBMol();
167      newMol->ReserveAtoms(fragment.size());
# Line 225 | Line 170 | namespace OpenBabel
170        OBAtom* newAtom = newMol->NewAtom();
171        *newAtom = *mol.GetAtom(*i);
172      }
173 +
174      newMol->EndModify();
175      newMol->ConnectTheDots();
176 <    findAngles(*newMol);
177 <    newMol->FindTorsions();
176 >    newMol->PerceiveBondOrders();
177 >
178      return newMol;
179    }
180    
181 <  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices) {
182 <    std::string indentLevel1("  ");
183 <    std::string indentLevel2("    ");
181 >  void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols,
182 >                                ostream& os, OBMol& mol,
183 >                                vector<int>& indices) {
184 >    
185      std::string molPrefix("MolName");
186 +    std::string resName;
187      unsigned int i;
188      const int BUFFLEN = 1024;
189      char buffer[BUFFLEN];
190      string str, str1;
191 <    
191 >    OBAtom *a, *b, *c, *d;    
192 >    bool molIsWater = false;
193 >    OBResidue *r;
194 >    int resKey;
195 >    char type_name[10];
196 >    char *element_name;
197 >    int res_num;
198 >    OBChainsParser* chainParser = new OBChainsParser();
199 >
200      
201      os << "<OOPSE version=4>" << endl;
202      os << "  <MetaData>" << endl << endl;
203      
204      for(i = 0; i < mols.size(); ++i) {
205        OBMol* pmol = mols[i];
250 –      
251 –      pmol->ConnectTheDots();
252 –      pmol->PerceiveBondOrders();
253 –      //pmol->FindSSSR();
254 –      //pmol->SetAromaticPerceived();
255 –      //pmol->Kekulize();
256 –      
206        map<OBAtom*, int> atomMap;
207 <      os << "molecule {\n";
208 <      sprintf(buffer, "%d", i);
209 <      os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
207 >
208 >      chainParser->PerceiveChains(*pmol, false);
209 >      molIsWater = false;
210 >      FOR_RESIDUES_OF_MOL(residue, *pmol) {
211 >        std::cerr << "residue = " << residue->GetName() << "\n";
212 >        if (residue->GetName().compare("HOH") == 0) {
213 >          molIsWater = true;
214 >        }
215 >      }
216        
217 <      //atom
218 <      int ai = 0;
219 <      FOR_ATOMS_OF_MOL(atom, *pmol ) {
220 <        str = atom->GetType();
221 <        ttab.SetFromType("INT");
222 <        ttab.SetToType("INT");
223 <        ttab.Translate(str1,str);
224 <        os << indentLevel1 << "atom[" << ai << "] {\n";
225 <        // os << indentLevel2 << "type = " << "\"" << etab.GetSymbol(atom->GetAtomicNum()) << "\"" << ";\n";
226 <        os << indentLevel2 << "type = " << "\"" << str1 << "\"" << ";\n";
227 <        os << indentLevel1 << "}\n";
228 <        atomMap[&(*atom)] = ai++;
229 <      }        
230 <      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 <          }
217 >      if (molIsWater) {
218 >        // water include files define all of the known water types
219 >        os << "#include \"water.md\";\n";
220 >        pmol->SetTitle("HOH");
221 >      } else {
222 >
223 >        os << "molecule {\n";
224 >        sprintf(buffer, "%d", i);
225 >        os << "  name = \"" << molPrefix << buffer << "\";\n";
226 >        
227 >        int ai = 0;
228 >        FOR_ATOMS_OF_MOL(atom, *pmol ) {
229 >          str = atom->GetType();
230 >          r = atom->GetResidue();
231            
232 <          //torsion
233 <          pGenericData = pmol->GetData(OBGenericDataType::TorsionData);
234 <          OBTorsionData* pTorsionData = dynamic_cast<OBTorsionData*>(pGenericData);
235 <          vector<OBTorsion> torsions = pTorsionData->GetData();
236 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > torsionArray;
237 <          for (vector<OBTorsion>::iterator ti = torsions.begin(); ti != torsions.end(); ++ti)
238 <          {
239 <          vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> > tmpTorsions = ti->getTorsions();
240 <          torsionArray.insert(torsionArray.end(), tmpTorsions.begin(), tmpTorsions.end());            
232 >          if (r == NULL)
233 >            resName = "NULL";
234 >          else
235 >            resName = r->GetName();
236 >          
237 >          if (resName.compare("NULL") ==0 ||
238 >              resName.compare("LIG") == 0 ||
239 >              resName.compare("UNK") == 0) {
240 >            // Either couldn't find a residue at all or couldn't find a
241 >            // reasonable residue name to use.  We'll punt and use
242 >            // OpenBabel's internal atom typing:
243 >            ttab.SetFromType("INT");
244 >            ttab.SetToType("INT");
245 >            ttab.Translate(str1, str);
246 >          } else {                        
247 >            
248 >            // If we know what residue we've got, the specific atom name can
249 >            // be used to help specify partial charges.
250 >            
251 >            str = r->GetAtomID(&*atom);
252 >            size_t startpos = str.find_first_not_of(" ");
253 >            size_t endpos = str.find_last_not_of(" ");
254 >            str = str.substr( startpos, endpos-startpos+1 );
255 >            str1 = resName + "-" + str;
256 >          }      
257 >          os << "  atom[" << ai << "] { ";
258 >          os << "type = " << "\"" << str1 << "\"" << "; ";
259 >          os << "}\n";
260 >          atomMap[&(*atom)] = ai++;
261 >        }        
262 >        os << "\n";
263 >        
264 >        //bond
265 >        
266 >        int b1, b2;
267 >        FOR_BONDS_OF_MOL(bond, *pmol ) {
268 >          b1 = atomMap[bond->GetBeginAtom()];
269 >          b2 = atomMap[bond->GetEndAtom()];
270 >          
271 >          os << "  bond { ";
272 >          
273 >          if (b1 < b2)
274 >            os << "members(" << b1 <<  ", " << b2 << "); ";
275 >          else
276 >            os << "members(" << b2 <<  ", " << b1 << "); ";
277 >          
278 >          os << "}\n";
279 >        }
280 >        
281 >        os << "\n";
282 >
283 >        std::vector<int> possibleInversion;
284 >        FOR_ATOMS_OF_MOL(atom, *pmol) {
285 >          possibleInversion.clear();
286 >          
287 >          FOR_NBORS_OF_ATOM(nbor, &*atom) {
288 >            possibleInversion.push_back(atomMap[&(*nbor)]);          
289            }
290            
291 <          os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
292 <          int torsionIndex = 0;
293 <          for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
294 <          {
295 <          os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
315 <          os << indentLevel2 << "member(" << atomMap[ti->first] <<  ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
316 <          os << indentLevel1 << "}\n";          
291 >          if (possibleInversion.size() == 3) {
292 >
293 >            os << "  inversion { ";
294 >            os << "center(" << atomMap[&(*atom)] << "); ";
295 >            os << "}\n";
296            }
297 <      */
298 <      os << "}" << endl;
299 <      os << endl;
300 <      
297 >          
298 >        }
299 >        os << "}" << endl;
300 >        os << endl;
301 >      }
302      }
303      
304      os << endl;
305      
306      
307 <    for(i=0; i < mols.size(); ++i) {      
307 >    for(i=0; i < mols.size(); ++i) {
308 >      OBMol* pmol = mols[i];      
309        os << "component{" << endl;
310 <      sprintf(buffer, "%d", i);
311 <      os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
312 <      os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
310 >      if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
311 >        os << "  type = " << "HOH" << ";" << endl;
312 >      } else {
313 >        sprintf(buffer, "%d", i);
314 >        os << "  type = " << molPrefix << buffer << ";" << endl;
315 >      }
316 >      os << "  nMol = " << numMols[i]<< ";" << endl;
317        os << "}" << endl;
318      }
319      
# Line 340 | Line 325 | namespace OpenBabel
325      
326      os << buffer << endl;
327      
328 +
329 +    // should really compute a bounding box here:
330      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);
331      
332      os << buffer << endl;
# Line 348 | Line 335 | namespace OpenBabel
335      
336      OBAtom *atom;
337      
338 +    // still to do: intercept waters and recompute pvqj lines
339 +
340      for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {    
341        atom = mol.GetAtom(*i);
342        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);

Diff Legend

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