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

# Content
1 /**********************************************************************
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 #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 {
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 //void findAngles(OBMol& mol);
64 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 size_t pos = OutputFileName.rfind(".");
128 if(pos!=string::npos)
129 OutputFileName = OutputFileName.substr(0, pos) + ".md";
130 else
131 OutputFileName += ".md";
132
133 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
184 /*
185 void OOPSEFormat::findAngles(OBMol& mol) {
186
187 //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
224 }
225 */
226
227 OBMol* OOPSEFormat::createMolFromFragment(OBMol& mol, vector<int>& fragment) {
228
229 OBMol* newMol = new OBMol();
230 OBChainsParser* chainParser = new OBChainsParser();
231 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 newMol->PerceiveBondOrders();
240 newMol->FindAngles();
241 newMol->FindTorsions();
242 //chainParser->PerceiveChains(*newMol, false);
243
244 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 std::string resName;
252 unsigned int i;
253 const int BUFFLEN = 1024;
254 char buffer[BUFFLEN];
255 string str, str1;
256 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
264 std::cerr << "yo\n";
265 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 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
278 map<OBAtom*, int> atomMap;
279 os << "molecule {\n";
280 sprintf(buffer, "%d", i);
281 os << indentLevel1 << "name = " << "\"" << molPrefix << buffer << "\"" << ";\n";
282
283
284 //atom
285
286 int ai = 0;
287 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 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 }
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