OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
openmdformat.cpp
1/**********************************************************************
2Copyright (C) 2000 by OpenEye Scientific Software, Inc.
3Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
4Some portions Copyright (C) 2004 by Chris Morley
5Some portions Copyright (C) 2004-present by J. Daniel Gezelter
6
7This program is free software; you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation version 2 of the License.
10
11This program is distributed in the hope that it will be useful,
12but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15***********************************************************************/
16
17#include <fstream>
18
19#include <openbabel/atom.h>
20#include <openbabel/babelconfig.h>
21#include <openbabel/bond.h>
22#include <openbabel/chains.h>
23#include <openbabel/data.h>
24#include <openbabel/mol.h>
25#include <openbabel/obiter.h>
26#include <openbabel/obmolecformat.h>
27#include <openbabel/obutil.h>
28#include <openbabel/stereo/tetrahedral.h>
29
30#include "utils/StringUtils.hpp"
31
32using namespace std;
33namespace OpenBabel {
34
35 class OpenMDFormat : public OBMoleculeFormat {
36 public:
37 // Register this format type ID
38 OpenMDFormat() { OBConversion::RegisterFormat("omd", this); }
39
40 virtual const char* Description() // required
41 {
42 return "OpenMD combined meta-data / cartesian coordinates format\n\
43 No comments yet\n";
44 };
45
46 virtual const char* SpecificationURL() {
47 return "http://openmd.org";
48 }; // optional
49
50 virtual const char* GetMIMEType() { return "chemical/x-omd"; };
51
52 virtual unsigned int Flags() { return NOTREADABLE | WRITEONEONLY; }
53
54 virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
55
56 private:
57 bool AreSameFragments(OBMol& mol, vector<int>& frag1, vector<int>& frag2);
58 OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
59 void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os,
60 OBMol& mol, vector<int>& indices);
61 void CalcBoundingBox(OBMol& mol, double& min_x, double& max_x,
62 double& min_y, double& max_y, double& min_z,
63 double& max_z);
64 };
65
66 // Make an instance of the format class
67 OpenMDFormat theOpenMDFormat;
68
69 bool OpenMDFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) {
70 OBMol* pmol = dynamic_cast<OBMol*>(pOb);
71 if (pmol == NULL) return false;
72
73 vector<vector<int>> fragmentLists;
74 pmol->ContigFragList(fragmentLists);
75 OBBitVec unused;
76 vector<bool> used(fragmentLists.size(), 0);
77 vector<vector<int>> molecules;
78 vector<int> indices;
79 for (std::size_t i = 0; i < used.size(); ++i) {
80 if (used[i]) continue;
81
82 used[i] = true;
83 vector<int> sameMolTypes;
84 sameMolTypes.push_back(i);
85 indices.insert(indices.end(), fragmentLists[i].begin(),
86 fragmentLists[i].end());
87 for (std::size_t j = i + 1; j < used.size(); ++j) {
88 if (used[j]) continue;
89
90 if (AreSameFragments(*pmol, fragmentLists[i], fragmentLists[j])) {
91 sameMolTypes.push_back(j);
92 indices.insert(indices.end(), fragmentLists[j].begin(),
93 fragmentLists[j].end());
94 used[j] = true;
95 }
96 }
97 molecules.push_back(sameMolTypes);
98 }
99
100 vector<OBMol*> mdMols;
101 vector<int> numMols;
102 for (vector<vector<int>>::iterator i = molecules.begin();
103 i != molecules.end(); ++i) {
104 mdMols.push_back(createMolFromFragment(*pmol, fragmentLists[i->front()]));
105 numMols.push_back((*i).size());
106 }
107
108 string OutputFileName = pConv->GetInFilename();
109 size_t pos = OutputFileName.rfind(".");
110 if (pos != string::npos)
111 OutputFileName = OutputFileName.substr(0, pos) + ".omd";
112 else
113 OutputFileName += ".omd";
114
115 ofstream ofs(OutputFileName.c_str());
116 if (!ofs) {
117 cerr << "Cannot write to " << OutputFileName << endl;
118 return false;
119 }
120
121 WriteMDFile(mdMols, numMols, ofs, *pmol, indices);
122
123 for (vector<OBMol*>::iterator i = mdMols.begin(); i != mdMols.end(); ++i) {
124 delete *i;
125 }
126
127 return (true);
128 }
129
130 bool OpenMDFormat::AreSameFragments(OBMol& mol, vector<int>& frag1,
131 vector<int>& frag2) {
132 if (frag1.size() != frag2.size()) return false;
133
134 // Exact graph matching is an NP complete problem.
135 // This just matches all of the atom atomic numbers and may falsely
136 // detect identical fragments which aren't really identical.
137
138 for (unsigned int i = 0; i < frag1.size(); ++i) {
139 OBAtom* atom1 = mol.GetAtom(frag1[i]);
140 OBAtom* atom2 = mol.GetAtom(frag2[i]);
141
142 if (atom1->GetAtomicNum() != atom2->GetAtomicNum()) return false;
143 }
144
145 OBMol* m1 = new OBMol();
146
147 m1->ReserveAtoms(frag1.size());
148 m1->BeginModify();
149 for (unsigned int i = 0; i < frag1.size(); ++i) {
150 OBAtom* atom1 = mol.GetAtom(frag1[i]);
151 OBAtom* newAtom = m1->NewAtom();
152 newAtom->Duplicate(atom1);
153 }
154
155 m1->EndModify();
156 m1->ConnectTheDots();
157 m1->PerceiveBondOrders();
158 OBStereoFacade facade1(m1, true);
159
160 OBMol* m2 = new OBMol();
161
162 m2->ReserveAtoms(frag2.size());
163 m2->BeginModify();
164 for (unsigned int i = 0; i < frag2.size(); ++i) {
165 OBAtom* atom2 = mol.GetAtom(frag2[i]);
166 OBAtom* newAtom = m2->NewAtom();
167 newAtom->Duplicate(atom2);
168 }
169
170 m2->EndModify();
171 m2->ConnectTheDots();
172 m2->PerceiveBondOrders();
173 OBStereoFacade facade2(m2, true);
174
175 OBTetrahedralStereo *stereo1;
176 OBTetrahedralStereo *stereo2;
177 OBTetrahedralStereo::Config config1;
178 OBTetrahedralStereo::Config config2;
179
180 OBAtomIterator a;
181 OBAtomIterator b;
182
183 for( a = m1->BeginAtoms(), b = m2->BeginAtoms(); a!= m1->EndAtoms(); ++a, ++b) {
184
185 std::cerr << (*a)->GetId() << " " << (*b)->GetId() << "\n";
186
187 if (facade1.HasTetrahedralStereo((*a)->GetId())) {
188 stereo1 = facade1.GetTetrahedralStereo((*a)->GetId());
189 config1 = stereo1->GetConfig();
190 std::cerr << *stereo1 << "\n";
191 }
192
193 if (facade2.HasTetrahedralStereo((*b)->GetId())) {
194 stereo2 = facade2.GetTetrahedralStereo((*b)->GetId());
195 config2 = stereo2->GetConfig();
196 std::cerr << *stereo2 << "\n";
197 }
198
199 if (config1 == config2) {
200 std::cerr << "configs are same";
201 } else {
202 std::cerr << "configs are different\n";
203 }
204 std::cerr << "\n";
205
206 }
207 return true;
208 }
209
210 struct SameAngle {
211 bool operator()(const triple<OBAtom*, OBAtom*, OBAtom*> t1,
212 const triple<OBAtom*, OBAtom*, OBAtom*> t2) const {
213 return (t1.second == t2.second) &&
214 ((t1.first == t2.first && t1.third == t2.third) ||
215 (t1.first == t2.third && t1.third == t2.first));
216 }
217 };
218
219 OBMol* OpenMDFormat::createMolFromFragment(OBMol& mol,
220 vector<int>& fragment) {
221 OBMol* newMol = new OBMol();
222
223 newMol->ReserveAtoms(fragment.size());
224 newMol->BeginModify();
225 for (vector<int>::iterator i = fragment.begin(); i != fragment.end(); ++i) {
226 OBAtom* newAtom = newMol->NewAtom();
227 *newAtom = *mol.GetAtom(*i);
228 }
229
230 newMol->EndModify();
231 newMol->ConnectTheDots();
232 newMol->PerceiveBondOrders();
233
234 return newMol;
235 }
236
237 void OpenMDFormat::WriteMDFile(vector<OBMol*> mols, vector<int> numMols,
238 ostream& os, OBMol& mol,
239 vector<int>& indices) {
240 std::string molPrefix("MolName");
241 std::string resName;
242 unsigned int i;
243 const int BUFFLEN = 1024;
244 char buffer[BUFFLEN];
245 string str, str1, str2, str3;
246 bool molIsWater = false;
247 OBResidue* r;
248 double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
249
250 os << "<OpenMD version=2>" << endl;
251 os << " <MetaData>" << endl << endl;
252
253 for (i = 0; i < mols.size(); ++i) {
254 OBMol* pmol = mols[i];
255 map<OBAtom*, int> atomMap;
256
257 molIsWater = false;
258 FOR_RESIDUES_OF_MOL(residue, *pmol) {
259 if (residue->GetName().compare("HOH") == 0) { molIsWater = true; }
260 }
261
262 if (molIsWater) {
263 // water include files define all of the known water types
264 os << "#include \"water.omd\";\n";
265 pmol->SetTitle("HOH");
266 } else {
267 os << "molecule {\n";
268 snprintf(buffer, BUFFLEN, "%u", i);
269 os << " name = \"" << molPrefix << buffer << "\";\n";
270
271 int ai = 0;
272 FOR_ATOMS_OF_MOL(atom, *pmol) {
273 str = atom->GetType();
274
275 r = atom->GetResidue();
276
277 if (r == NULL)
278 resName = "NULL";
279 else
280 resName = r->GetName();
281
282 if (resName.compare("NULL") == 0 || resName.compare("LIG") == 0 ||
283 resName.compare("UNL") == 0 || resName.compare("UNK") == 0) {
284 // Either couldn't find a residue at all or couldn't find a
285 // reasonable residue name to use. We'll punt and use
286 // OpenBabel's internal atom typing:
287 ttab.SetFromType("INT");
288 ttab.SetToType("INT");
289 ttab.Translate(str1, str);
290 } else {
291 // If we know what residue we've got, the specific atom name can
292 // be used to help specify partial charges.
293
294 // resdat.SetResName(resName);
295
296 // atom type from residue:
297 str = r->GetAtomID(&*atom);
298
299 // arginine has separate indices for chemically-identical
300 // nitrogen atoms:
301 if (resName.compare("ARG") == 0) {
302 if (str.compare("NH1") == 0 || str.compare("NH2") == 0) {
303 str = "NH";
304 }
305 }
306 if (resName.compare("VAL") == 0) {
307 if (str.compare("CG1") == 0 || str.compare("CG2") == 0) {
308 str = "CG";
309 }
310 }
311 if (resName.compare("LEU") == 0) {
312 if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
313 str = "CD";
314 }
315 }
316 if (resName.compare("ASP") == 0) {
317 if (str.compare("OD1") == 0 || str.compare("OD2") == 0) {
318 str = "OD";
319 }
320 }
321 if (resName.compare("GLU") == 0) {
322 if (str.compare("OE1") == 0 || str.compare("OE2") == 0) {
323 str = "OE";
324 }
325 }
326 if (resName.compare("TYR") == 0) {
327 if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
328 str = "CD";
329 }
330 if (str.compare("CE1") == 0 || str.compare("CE2") == 0) {
331 str = "CE";
332 }
333 }
334
335 if ((&*atom)->GetAtomicNum() == 1) {
336 FOR_NBORS_OF_ATOM(nbr, *atom) {
337 str2 = r->GetAtomID(&*nbr);
338 size_t startpos = str2.find_first_not_of(" ");
339 size_t endpos = str2.find_last_not_of(" ");
340 if ((endpos - startpos) < 1) {
341 // if the bonded atom type has only one character (i.e. N)
342 // then the hydrogen will be labeled "HN" to show what
343 // kind of proton it is:
344 str3 = str2;
345 } else {
346 if (str2.compare("OH") == 0) {
347 str3 = "O";
348 } else {
349 // When the bonded atom type is more specific, we drop
350 // the first character: i.e. H bonded to OG1 is HG1 type:
351 str3 = str2.substr(startpos + 1, endpos - startpos);
352 }
353 }
354 str = "H" + str3;
355 }
356 // same problem with arginine NH atoms, but now for connected
357 // hydrogens
358 if (resName.compare("ARG") == 0) {
359 if (str.compare("HH1") == 0 || str.compare("HH2") == 0) {
360 str = "HH";
361 }
362 }
363 if (resName.compare("VAL") == 0) {
364 if (str.compare("HG1") == 0 || str.compare("HG2") == 0) {
365 str = "HG";
366 }
367 }
368 if (resName.compare("LEU") == 0) {
369 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
370 str = "HD";
371 }
372 }
373 if (resName.compare("TYR") == 0) {
374 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
375 str = "HD";
376 }
377 if (str.compare("HE1") == 0 || str.compare("HE2") == 0) {
378 str = "HE";
379 }
380 }
381 }
382
383 // atom type from residue table:
384 // resdat.LookupType(str, str2, hyb);
385 size_t startpos = str.find_first_not_of(" ");
386 size_t endpos = str.find_last_not_of(" ");
387 str = str.substr(startpos, endpos - startpos + 1);
388 str1 = resName + "-" + str;
389 }
390 os << " atom[" << ai << "] { ";
391 os << "type = "
392 << "\"" << str1 << "\""
393 << "; ";
394 os << "position( " << (&*atom)->GetX() << ", " << (&*atom)->GetY()
395 << ", " << (&*atom)->GetZ() << ");";
396 os << "}\n";
397 atomMap[&(*atom)] = ai++;
398 }
399 os << "\n";
400
401 // bond
402
403 int bi;
404 double bo;
405 OBAtom *a, *b;
406
407 FOR_BONDS_OF_MOL(bond, *pmol) {
408 a = bond->GetBeginAtom();
409 b = bond->GetEndAtom();
410
411 ai = atomMap[a];
412 bi = atomMap[b];
413 bo = bond->GetBondOrder();
414 if (bond->IsAromatic()) bo = 1.5;
415 // e.g., in Cp rings, may not be "aromatic" by OB but check
416 // for explicit hydrogen counts (e.g., biphenyl inter-ring
417 // is not aromatic)
418 if ((a->GetType()[2] == 'R' && b->GetType()[2] == 'R') &&
419 (a->ExplicitHydrogenCount() == 1 &&
420 b->ExplicitHydrogenCount() == 1))
421 bo = 1.5;
422 if (bond->IsAmide()) bo = 1.41;
423
424 os << " bond { ";
425
426 if (ai < bi)
427 os << "members(" << ai << ", " << bi << "); ";
428 else
429 os << "members(" << bi << ", " << ai << "); ";
430
431 if (bo != 1) os << "bondOrder = " << bo << "; ";
432
433 os << "}" << endl;
434 }
435
436 os << endl;
437
438 os << "}" << endl;
439 os << endl;
440 }
441 }
442
443 os << endl;
444
445 for (i = 0; i < mols.size(); ++i) {
446 OBMol* pmol = mols[i];
447 os << "component{" << endl;
448 if (std::string(pmol->GetTitle()).compare("HOH") == 0) {
449 os << " type = "
450 << "\"HOH\""
451 << "; // change to appropriate water model" << endl;
452 } else {
453 snprintf(buffer, BUFFLEN, "%u", i);
454 os << " type = " << molPrefix << buffer << ";" << endl;
455 }
456 os << " nMol = " << numMols[i] << ";" << endl;
457 os << "}" << endl;
458 }
459
460 os << " </MetaData>" << endl;
461 os << " <Snapshot>" << endl;
462 os << " <FrameData>" << endl;
463
464 snprintf(buffer, BUFFLEN, " Time: %.10g", 0.0);
465
466 os << buffer << endl;
467
468 CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);
469
470 // still to do: should compute a bounding box here
471 snprintf(
472 buffer, BUFFLEN,
473 " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { "
474 "%.10g, "
475 "%.10g, %.10g }}",
476 max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0,
477 max_z - min_z);
478
479 os << buffer << endl;
480 os << " </FrameData>" << endl;
481 os << " <StuntDoubles>" << endl;
482
483 OBAtom* atom;
484
485 for (vector<int>::iterator i = indices.begin(); i != indices.end(); ++i) {
486 atom = mol.GetAtom(*i);
487 snprintf(buffer, BUFFLEN,
488 "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1, "pv",
489 atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
490 os << buffer << endl;
491 }
492 os << " </StuntDoubles>" << endl;
493 os << " </Snapshot>" << endl;
494 os << "</OpenMD>" << endl;
495 }
496
497 void OpenMDFormat::CalcBoundingBox(OBMol& mol, double& min_x, double& max_x,
498 double& min_y, double& max_y,
499 double& min_z, double& max_z) {
500 /* ---- Init bounding-box variables ---- */
501 min_x = (double)0.0;
502 max_x = (double)0.0;
503 min_y = (double)0.0;
504 max_y = (double)0.0;
505 min_z = (double)0.0;
506 max_z = (double)0.0;
507
508 /* ---- Check all atoms ---- */
509 for (unsigned int i = 1; i <= mol.NumAtoms(); ++i) {
510 /* ---- Get a pointer to ith atom ---- */
511 OBAtom* atom = mol.GetAtom(i);
512
513 /* ---- Check for minimal/maximal x-position ---- */
514 if (atom->GetX() < min_x) min_x = atom->GetX();
515 if (atom->GetX() > max_x) max_x = atom->GetX();
516
517 /* ---- Check for minimal/maximal y-position ---- */
518 if (atom->GetY() < min_y) min_y = atom->GetY();
519 if (atom->GetY() > max_y) max_y = atom->GetY();
520
521 /* ---- Check for minimal/maximal z-position ---- */
522 if (atom->GetZ() < min_z) min_z = atom->GetZ();
523 if (atom->GetZ() > max_z) max_z = atom->GetZ();
524 }
525 }
526} // namespace OpenBabel