ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/atom2md/oopseformat.cpp
Revision: 1269
Committed: Tue Jul 1 13:28:23 2008 UTC (16 years, 10 months ago) by gezelter
File size: 16323 byte(s)
Log Message:
Adding infrastructure for Amber force field.

File Contents

# User Rev Content
1 gezelter 1210 /**********************************************************************
2     Copyright (C) 2000 by OpenEye Scientific Software, Inc.
3     Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
4     Some portions Copyright (C) 2004 by Chris Morley
5     Some portions Copyright (C) 2008 by J. Daniel Gezelter
6    
7     This program is free software; you can redistribute it and/or modify
8     it under the terms of the GNU General Public License as published by
9     the Free Software Foundation version 2 of the License.
10    
11     This program is distributed in the hope that it will be useful,
12     but WITHOUT ANY WARRANTY; without even the implied warranty of
13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14     GNU General Public License for more details.
15     ***********************************************************************/
16    
17     #include <openbabel/babelconfig.h>
18     #include <openbabel/obmolecformat.h>
19 gezelter 1269 #include <openbabel/obiter.h>
20     #include <openbabel/mol.h>
21     #include <openbabel/chains.h>
22 gezelter 1210 #include <fstream>
23    
24 gezelter 1269 #include "utils/StringUtils.hpp"
25    
26 gezelter 1210 using namespace std;
27     namespace OpenBabel
28     {
29    
30     class OOPSEFormat : public OBMoleculeFormat
31     {
32     public:
33     //Register this format type ID
34     OOPSEFormat()
35     {
36     OBConversion::RegisterFormat("md",this);
37     }
38    
39     virtual const char* Description() //required
40     {
41     return
42     "OOPSE combined meta-data / cartesian coordinates format\nNo comments yet\n";
43     };
44    
45     virtual const char* SpecificationURL()
46     {return "http://www.oopse.org";}; //optional
47    
48     virtual const char* GetMIMEType()
49     {return "chemical/x-md"; };
50    
51     virtual unsigned int Flags()
52     {
53     return NOTREADABLE | WRITEONEONLY;
54     }
55    
56     //*** This section identical for most OBMol conversions ***
57     ////////////////////////////////////////////////////
58     /// The "API" interface functions
59     virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
60    
61     private:
62     bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
63 gezelter 1269 //void findAngles(OBMol& mol);
64 gezelter 1210 OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
65     void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices);
66     };
67    
68     //Make an instance of the format class
69     OOPSEFormat theOOPSEFormat;
70    
71     /////////////////////////////////////////////////////////////////
72    
73    
74     bool OOPSEFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv)
75     {
76     OBMol* pmol = dynamic_cast<OBMol*>(pOb);
77     if(pmol==NULL)
78     return false;
79    
80     vector<vector<int> > fragmentLists;
81     pmol->ContigFragList(fragmentLists);
82     OBBitVec unused;
83     vector<bool> used(fragmentLists.size(), 0);
84     vector<vector<int> > molecules;
85     vector<int> indices;
86     for(int i =0; i < used.size(); ++i) {
87     if (used[i])
88     {
89     continue;
90     }
91     used[i] = true;
92     vector<int> sameMolTypes;
93     sameMolTypes.push_back(i);
94     indices.insert(indices.end(), fragmentLists[i].begin(),
95     fragmentLists[i].end());
96     for (int j = i + 1;j < used.size(); ++j)
97     {
98     if (used[j])
99     {
100     continue;
101     }
102    
103     if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j]))
104     {
105     sameMolTypes.push_back(j);
106     indices.insert(indices.end(), fragmentLists[j].begin(),
107     fragmentLists[j].end());
108     used[j]=true;
109     }
110     }
111     molecules.push_back(sameMolTypes);
112    
113     }
114    
115     //
116     vector<OBMol*> mdMols;
117     vector<int> numMols;
118     for(vector<vector<int> >::iterator i = molecules.begin();
119     i != molecules.end(); ++i)
120     {
121     mdMols.push_back(createMolFromFragment(*pmol,
122     fragmentLists[i->front()]));
123     numMols.push_back((*i).size());
124     }
125    
126     string OutputFileName = pConv->GetInFilename();
127 gezelter 1269 size_t pos = OutputFileName.rfind(".");
128     if(pos!=string::npos)
129     OutputFileName = OutputFileName.substr(0, pos) + ".md";
130     else
131 gezelter 1210 OutputFileName += ".md";
132 gezelter 1269
133 gezelter 1210 ofstream ofs(OutputFileName.c_str());
134     if(!ofs)
135     {
136     cerr << "Cannot write to " << OutputFileName <<endl;
137     return false;
138     }
139    
140     WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
141    
142     for(vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i)
143     {
144     delete *i;
145     }
146    
147     //
148    
149     return(true);
150     }
151    
152     bool OOPSEFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
153     vector<int>& frag2)
154     {
155     if (frag1.size() != frag2.size())
156     {
157     return false;
158     }
159    
160     //exact graph matching is a NP complete problem
161     /** @todo using sparse matrix to store the connectivities*/
162     for (unsigned int i =0 ; i < frag1.size(); ++i)
163     {
164     OBAtom* atom1 = mol.GetAtom(frag1[i]);
165     OBAtom* atom2 = mol.GetAtom(frag2[i]);
166    
167     if (atom1->GetAtomicNum() != atom2->GetAtomicNum())
168     {
169     return false;
170     }
171     }
172     return true;
173     }
174    
175     struct SameAngle
176     {
177     bool operator()(const triple<OBAtom*,OBAtom*,OBAtom*> t1,
178     const triple<OBAtom*,OBAtom*,OBAtom*> t2) const
179     {
180     return (t1.second == t2.second) && ( (t1.first == t2.first && t1.third == t2.third) || (t1.first == t2.third && t1.third == t2.first));
181     }
182     };
183 gezelter 1269
184     /*
185     void OOPSEFormat::findAngles(OBMol& mol) {
186    
187 gezelter 1210 //if already has data return
188     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 gezelter 1269
224 gezelter 1210 }
225 gezelter 1269 */
226 gezelter 1210
227     OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment) {
228    
229     OBMol* newMol = new OBMol();
230 gezelter 1269 OBChainsParser* chainParser = new OBChainsParser();
231 gezelter 1210 newMol->ReserveAtoms(fragment.size());
232     newMol->BeginModify();
233     for(vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
234     OBAtom* newAtom = newMol->NewAtom();
235     *newAtom = *mol.GetAtom(*i);
236     }
237     newMol->EndModify();
238     newMol->ConnectTheDots();
239 gezelter 1269 newMol->PerceiveBondOrders();
240     newMol->FindAngles();
241 gezelter 1210 newMol->FindTorsions();
242 gezelter 1269 //chainParser->PerceiveChains(*newMol, false);
243    
244 gezelter 1210 return newMol;
245     }
246    
247     void OOPSEFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os, OBMol& mol, vector<int>& indices) {
248     std::string indentLevel1(" ");
249     std::string indentLevel2(" ");
250     std::string molPrefix("MolName");
251 gezelter 1269 std::string resName;
252 gezelter 1210 unsigned int i;
253     const int BUFFLEN = 1024;
254     char buffer[BUFFLEN];
255     string str, str1;
256 gezelter 1269 OBAtom *a, *b, *c, *d;
257     OBChainsParser* chainParser = new OBChainsParser();
258     OBResidue *r;
259     int resKey;
260     char type_name[10];
261     char *element_name;
262     int res_num;
263 gezelter 1210
264 gezelter 1269 std::cerr << "yo\n";
265 gezelter 1210 os << "<OOPSE version=4>" << endl;
266     os << " <MetaData>" << endl << endl;
267    
268     for(i = 0; i < mols.size(); ++i) {
269     OBMol* pmol = mols[i];
270 gezelter 1269 std::cerr << "yo1\n";
271 gezelter 1210 pmol->ConnectTheDots();
272     pmol->PerceiveBondOrders();
273 gezelter 1269 chainParser->PerceiveChains(*pmol, false);
274 gezelter 1210 //pmol->FindSSSR();
275     //pmol->SetAromaticPerceived();
276     //pmol->Kekulize();
277    
278     map<OBAtom*, int> atomMap;
279     os << "molecule {\n";
280     sprintf(buffer, "%d", i);
281     os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
282    
283 gezelter 1269
284 gezelter 1210 //atom
285 gezelter 1269
286 gezelter 1210 int ai = 0;
287 gezelter 1269 FOR_RESIDUES_OF_MOL(res, *pmol) {
288    
289     std::cerr << "found residue" << res->GetName() << "\n";
290     }
291    
292     // FOR_RESIDUES_OF_MOL(res, *pmol) {
293    
294     // resName = res->GetName();
295     // resKey = res->GetResKey();
296    
297    
298    
299     // std::cerr << "found residue " << resName << "\n";
300     // std::cerr << " num = " << res->GetNum() << "\n";
301     // std::cerr << " numAtoms = " << res->GetNumAtoms() << "\n";
302     // std::cerr << " num = " << res->GetNum() << "\n";
303     // std::cerr << " chain = " << res->GetChain() << "\n";
304     // 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 gezelter 1210 FOR_ATOMS_OF_MOL(atom, *pmol ) {
325     str = atom->GetType();
326 gezelter 1269 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 gezelter 1210 atomMap[&(*atom)] = ai++;
372     }
373     os << "\n";
374    
375     //bond
376     FOR_BONDS_OF_MOL(bond, *pmol ) {
377 gezelter 1269 // 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 gezelter 1210 /*
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     }
438    
439     //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     }
449    
450     os << indentLevel1 << "nTorsions = " << torsionArray.size() << ";\n";
451     int torsionIndex = 0;
452     for (vector<quad<OBAtom*,OBAtom*,OBAtom*,OBAtom*> >::iterator ti = torsionArray.begin(); ti != torsionArray.end(); ++ti)
453     {
454     os << indentLevel1 << "torsion[" << torsionIndex++ << "] {\n";
455     os << indentLevel2 << "member(" << atomMap[ti->first] << ", " << atomMap[ti->second] <<", " << atomMap[ti->third] <<", " << atomMap[ti->forth] << ");\n";
456     os << indentLevel1 << "}\n";
457     }
458     */
459     os << "}" << endl;
460     os << endl;
461    
462     }
463    
464     os << endl;
465    
466    
467     for(i=0; i < mols.size(); ++i) {
468     os << "component{" << endl;
469     sprintf(buffer, "%d", i);
470     os << indentLevel1 << "type = " << molPrefix << buffer << ";" << endl;
471     os << indentLevel1 << "nMol = " << numMols[i]<< ";" << endl;
472     os << "}" << endl;
473     }
474    
475     os << " </MetaData>" << endl;
476     os << " <Snapshot>" << endl;
477     os << " <FrameData>" << endl;
478    
479     sprintf(buffer, " Time: %.10g", 0.0);
480    
481     os << buffer << endl;
482    
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    
485     os << buffer << endl;
486     os << " </FrameData>" << endl;
487     os << " <StuntDoubles>" << endl;
488    
489     OBAtom *atom;
490    
491     for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {
492     atom = mol.GetAtom(*i);
493     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);
494     os << buffer << endl;
495     }
496     os << " </StuntDoubles>" << endl;
497     os << " </Snapshot>" << endl;
498     os << "</OOPSE>" << endl;
499     }
500    
501     } //namespace OpenBabel
502