/src/openbabel/src/reactionfacade.cpp
Line  | Count  | Source  | 
1  |  | /**********************************************************************  | 
2  |  | reactionfacade.cpp - A facade class to help interrogate and manipulate  | 
3  |  |                      reactions  | 
4  |  |  | 
5  |  | Copyright (C) 2018 by Noel M. O'Boyle  | 
6  |  |  | 
7  |  | This file is part of the Open Babel project.  | 
8  |  | For more information, see <http://openbabel.org/>  | 
9  |  |  | 
10  |  | This program is free software; you can redistribute it and/or modify  | 
11  |  | it under the terms of the GNU General Public License as published by  | 
12  |  | the Free Software Foundation version 2 of the License.  | 
13  |  |  | 
14  |  | This program is distributed in the hope that it will be useful,  | 
15  |  | but WITHOUT ANY WARRANTY; without even the implied warranty of  | 
16  |  | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the  | 
17  |  | GNU General Public License for more details.  | 
18  |  | ***********************************************************************/  | 
19  |  |  | 
20  |  | #include <openbabel/mol.h>  | 
21  |  | #include <openbabel/atom.h>  | 
22  |  | #include <openbabel/obiter.h>  | 
23  |  | #include <openbabel/generic.h>  | 
24  |  | #include <openbabel/oberror.h>  | 
25  |  | #include <openbabel/reactionfacade.h>  | 
26  |  |  | 
27  |  | #include <set>  | 
28  |  |  | 
29  |  | namespace OpenBabel  | 
30  |  | { | 
31  |  |   extern OBMessageHandler obErrorLog;  | 
32  |  |  | 
33  |  |   class OBReactionFacadePrivate  | 
34  |  |   { | 
35  |  |   public:  | 
36  | 0  |     OBReactionFacadePrivate(OBMol *mol) : _mol(mol), _found_components(false)  | 
37  | 0  |     { | 
38  | 0  |     }  | 
39  |  |       | 
40  |  |     void AddComponent(OBMol* mol, OBReactionRole rxnrole);  | 
41  |  |     void AssignComponentIds(bool wipe);  | 
42  |  |     void ClearInternalState();  | 
43  |  |     bool GetComponent(OBMol* mol, OBReactionRole rxnrole, unsigned int num);  | 
44  |  |     unsigned int GetComponentId(OBAtom *atom);  | 
45  |  |     OBReactionRole GetRole(OBAtom *atom);  | 
46  |  |     bool IsValid();  | 
47  |  |     unsigned int NumComponents(OBReactionRole rxnrole);  | 
48  |  |     bool ReassignComponent(OBReactionRole oldrole, unsigned int num, OBReactionRole newrole);  | 
49  |  |     void SetComponentId(OBAtom* atom, unsigned int compid);  | 
50  |  |     void SetRole(OBAtom* atom, OBReactionRole rxnrole);  | 
51  |  |  | 
52  |  |   private:  | 
53  |  |     OBMol* _mol;  | 
54  |  |     bool _found_components;  | 
55  |  |     std::vector<unsigned int> _unassigned_components;  | 
56  |  |     std::vector<unsigned int> _reactant_components;  | 
57  |  |     std::vector<unsigned int> _product_components;  | 
58  |  |     std::vector<unsigned int> _agent_components;  | 
59  |  |  | 
60  |  |     int GetId(const char* idtype, OBAtom *atom);  | 
61  |  |     void SetId(const char* idtype, OBAtom* atom, int idval);  | 
62  |  |     void FindComponents();  | 
63  |  |     std::vector<unsigned int>* GetComponentIds(OBReactionRole rxnrole);  | 
64  |  |   };  | 
65  |  |  | 
66  |  |   void OBReactionFacadePrivate::AddComponent(OBMol* mol, OBReactionRole rxnrole)  | 
67  | 0  |   { | 
68  | 0  |     if (!_found_components)  | 
69  | 0  |       FindComponents();  | 
70  |  |  | 
71  |  |     // Get the current maximum component id  | 
72  | 0  |     unsigned int max_compid = 0;  | 
73  | 0  |     if (!_product_components.empty())  | 
74  | 0  |       max_compid = _product_components.back();  | 
75  | 0  |     if (!_agent_components.empty() && _agent_components.back()>max_compid)  | 
76  | 0  |       max_compid = _agent_components.back();  | 
77  | 0  |     if (!_reactant_components.empty() && _reactant_components.back()>max_compid)  | 
78  | 0  |       max_compid = _reactant_components.back();  | 
79  | 0  |     if (!_unassigned_components.empty() && _unassigned_components.back()>max_compid)  | 
80  | 0  |       max_compid = _unassigned_components.back();  | 
81  |  | 
  | 
82  | 0  |     int new_compid = max_compid + 1;  | 
83  | 0  |     if (new_compid == 0)  | 
84  | 0  |       new_compid = 1;  | 
85  |  | 
  | 
86  | 0  |     FOR_ATOMS_OF_MOL(atom, mol) { | 
87  | 0  |       SetRole(&*atom, rxnrole);  | 
88  | 0  |       SetComponentId(&*atom, new_compid);  | 
89  | 0  |     }  | 
90  | 0  |     *_mol += *mol;  | 
91  |  | 
  | 
92  | 0  |     GetComponentIds(rxnrole)->push_back(new_compid);  | 
93  | 0  |   }  | 
94  |  |  | 
95  |  |   void OBReactionFacadePrivate::AssignComponentIds(bool wipe)  | 
96  | 0  |   { | 
97  | 0  |     unsigned int compid = 1;  | 
98  |  | 
  | 
99  | 0  |     OBMolAtomDFSIter iter(_mol);  | 
100  | 0  |     while(iter) { // for each connected component | 
101  | 0  |       do { // for each atom in connected component | 
102  | 0  |         if (wipe || !(iter->HasData("rxncomp"))) | 
103  | 0  |           SetComponentId(&*iter, compid);  | 
104  | 0  |       } while ((iter++).next());  | 
105  | 0  |       compid++;  | 
106  | 0  |     }  | 
107  | 0  |   }  | 
108  |  |  | 
109  |  |   void OBReactionFacadePrivate::ClearInternalState()  | 
110  | 0  |   { | 
111  | 0  |     _found_components = false;  | 
112  | 0  |   }  | 
113  |  |  | 
114  |  |   void OBReactionFacadePrivate::FindComponents()  | 
115  | 0  |   { | 
116  | 0  |     std::set<unsigned int> reactant_components;  | 
117  | 0  |     std::set<unsigned int> product_components;  | 
118  | 0  |     std::set<unsigned int> agent_components;  | 
119  | 0  |     std::set<unsigned int> unassigned_components;  | 
120  |  | 
  | 
121  | 0  |     FOR_ATOMS_OF_MOL(atom, _mol) { | 
122  | 0  |       OBAtom* atm = &*atom;  | 
123  | 0  |       unsigned int component = GetComponentId(&*atom);  | 
124  | 0  |       switch (GetRole(&*atom)) { | 
125  | 0  |       case REACTANT:  | 
126  | 0  |         reactant_components.insert(component);  | 
127  | 0  |         break;  | 
128  | 0  |       case PRODUCT:  | 
129  | 0  |         product_components.insert(component);  | 
130  | 0  |         break;  | 
131  | 0  |       case AGENT:  | 
132  | 0  |         agent_components.insert(component);  | 
133  | 0  |         break;  | 
134  | 0  |       default:  | 
135  | 0  |         unassigned_components.insert(component);  | 
136  | 0  |       }  | 
137  | 0  |     }  | 
138  |  |     // Convert to vector so we have random access - note: these will be in sorted order  | 
139  | 0  |     for (std::set<unsigned int>::iterator sit = reactant_components.begin(); sit != reactant_components.end(); ++sit)  | 
140  | 0  |       _reactant_components.push_back(*sit);  | 
141  | 0  |     for (std::set<unsigned int>::iterator sit = product_components.begin(); sit != product_components.end(); ++sit)  | 
142  | 0  |       _product_components.push_back(*sit);  | 
143  | 0  |     for (std::set<unsigned int>::iterator sit = agent_components.begin(); sit != agent_components.end(); ++sit)  | 
144  | 0  |       _agent_components.push_back(*sit);  | 
145  | 0  |     for (std::set<unsigned int>::iterator sit = unassigned_components.begin(); sit != unassigned_components.end(); ++sit)  | 
146  | 0  |       _unassigned_components.push_back(*sit);  | 
147  |  | 
  | 
148  | 0  |     _found_components = true;  | 
149  | 0  |   }  | 
150  |  |  | 
151  |  |   bool OBReactionFacadePrivate::GetComponent(OBMol* mol, OBReactionRole rxnrole, unsigned int num)  | 
152  | 0  |   { | 
153  | 0  |     std::vector<unsigned int> *component_ids = GetComponentIds(rxnrole);  | 
154  | 0  |     if (num >= component_ids->size())  | 
155  | 0  |       return false;  | 
156  | 0  |     OBBitVec atoms;  | 
157  | 0  |     unsigned int componentId = (*component_ids)[num];  | 
158  | 0  |     FOR_ATOMS_OF_MOL(atom, _mol) { | 
159  | 0  |       if (GetRole(&*atom) == rxnrole && GetComponentId(&*atom) == componentId)  | 
160  | 0  |         atoms.SetBitOn(atom->GetIdx());  | 
161  | 0  |     }  | 
162  | 0  |     bool ok = _mol->CopySubstructure(*mol, &atoms);  | 
163  | 0  |     return ok;  | 
164  | 0  |   }  | 
165  |  |  | 
166  |  |   unsigned int OBReactionFacadePrivate::GetComponentId(OBAtom *atom)  | 
167  | 0  |   { | 
168  | 0  |     return GetId("rxncomp", atom); | 
169  | 0  |   }  | 
170  |  |  | 
171  |  |   std::vector<unsigned int>* OBReactionFacadePrivate::GetComponentIds(OBReactionRole rxnrole)  | 
172  | 0  |   { | 
173  | 0  |     if (!_found_components)  | 
174  | 0  |       FindComponents();  | 
175  |  | 
  | 
176  | 0  |     switch (rxnrole) { | 
177  | 0  |     case REACTANT:  | 
178  | 0  |       return &_reactant_components;  | 
179  | 0  |     case PRODUCT:  | 
180  | 0  |       return &_product_components;  | 
181  | 0  |     case AGENT:  | 
182  | 0  |       return &_agent_components;  | 
183  | 0  |     case NO_REACTIONROLE:  | 
184  | 0  |       return &_unassigned_components;  | 
185  | 0  |     }  | 
186  | 0  |     return nullptr;  | 
187  | 0  |   }  | 
188  |  |  | 
189  |  |   int OBReactionFacadePrivate::GetId(const char* idtype, OBAtom *atom)  | 
190  | 0  |   { | 
191  | 0  |     int idval = 0;  | 
192  | 0  |     OBGenericData *data = atom->GetData(idtype);  | 
193  | 0  |     if (data) { | 
194  | 0  |       OBPairInteger *pi = (OBPairInteger*)data;  | 
195  | 0  |       idval = pi->GetGenericValue();  | 
196  | 0  |     }  | 
197  | 0  |     return idval;  | 
198  | 0  |   }  | 
199  |  |  | 
200  |  |   OBReactionRole OBReactionFacadePrivate::GetRole(OBAtom *atom)  | 
201  | 0  |   { | 
202  | 0  |     int rxnrole = GetId("rxnrole", atom); | 
203  | 0  |     switch (rxnrole) { | 
204  | 0  |     default: case 0: return NO_REACTIONROLE;  | 
205  | 0  |     case 1: return REACTANT;  | 
206  | 0  |     case 2: return AGENT;  | 
207  | 0  |     case 3: return PRODUCT;  | 
208  | 0  |     }  | 
209  | 0  |   }  | 
210  |  |  | 
211  |  |   bool OBReactionFacadePrivate::IsValid()  | 
212  | 0  |   { | 
213  | 0  |     if (!_mol->IsReaction()) { | 
214  | 0  |       obErrorLog.ThrowError(__FUNCTION__, "The molecule is not marked as a reaction. Use SetIsReaction().", obWarning);  | 
215  | 0  |       return false;  | 
216  | 0  |     }  | 
217  | 0  |     FOR_ATOMS_OF_MOL(atom, _mol) { | 
218  | 0  |       OBGenericData *data;  | 
219  | 0  |       OBPairInteger *pi;  | 
220  | 0  |       int val;  | 
221  |  | 
  | 
222  | 0  |       data = atom->GetData("rxncomp"); | 
223  | 0  |       if (!data) { | 
224  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "The molecule contains an atom that is missing a reaction component Id. Use SetComponentId().", obWarning);  | 
225  | 0  |         return false;  | 
226  | 0  |       }  | 
227  | 0  |       pi = dynamic_cast<OBPairInteger*>(data);  | 
228  | 0  |       if (!pi) { | 
229  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "A reaction component Id has been stored using a data type that is not an OBPairInteger.", obWarning);  | 
230  | 0  |         return false;  | 
231  | 0  |       }  | 
232  | 0  |       val = pi->GetGenericValue();  | 
233  | 0  |       if (val <= 0) { | 
234  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "Reaction component Ids should all be non-zero positive integers.", obWarning);  | 
235  | 0  |         return false;  | 
236  | 0  |       }  | 
237  |  |  | 
238  | 0  |       data = atom->GetData("rxnrole"); | 
239  | 0  |       if (!data) { | 
240  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "The molecule contains an atom that is missing reaction role information. Use SetRole().", obWarning);  | 
241  | 0  |         return false;  | 
242  | 0  |       }  | 
243  | 0  |       pi = dynamic_cast<OBPairInteger*>(data);  | 
244  | 0  |       if (!pi) { | 
245  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "Reaction role information has been stored using a data type that is not an OBPairInteger.", obWarning);  | 
246  | 0  |         return false;  | 
247  | 0  |       }  | 
248  | 0  |       val = pi->GetGenericValue();  | 
249  | 0  |       if (val < 0 || val > 3) { | 
250  | 0  |         obErrorLog.ThrowError(__FUNCTION__, "Reaction roles should be in the range 0 to 3 inclusive.", obWarning);  | 
251  | 0  |         return false;  | 
252  | 0  |       }  | 
253  | 0  |     }  | 
254  |  |  | 
255  |  |     // Ensure that every atom in a particular connected component has the same component Id and rxn role  | 
256  | 0  |     OBMolAtomDFSIter iter(_mol);  | 
257  | 0  |     while (iter) { // for each connected component | 
258  | 0  |       unsigned int rxncomp = GetComponentId(&*iter);  | 
259  | 0  |       OBReactionRole rxnrole = GetRole(&*iter);  | 
260  | 0  |       do { // for each atom in connected component | 
261  | 0  |         unsigned my_rxncomp = GetComponentId(&*iter);  | 
262  | 0  |         if (my_rxncomp != rxncomp) { | 
263  | 0  |           obErrorLog.ThrowError(__FUNCTION__, "The molecule contains a connected component that contains atoms with different reaction component Ids. All atoms in a particular connected component should have the same value.", obWarning);  | 
264  | 0  |           return false;  | 
265  | 0  |         }  | 
266  | 0  |         unsigned int my_rxnrole = GetRole(&*iter);  | 
267  | 0  |         if (my_rxnrole != rxnrole) { | 
268  | 0  |           obErrorLog.ThrowError(__FUNCTION__, "The molecule contains a connected component that contains atoms with different reaction roles. All atoms in a particular connected component should have the same role.", obWarning);  | 
269  | 0  |           return false;  | 
270  | 0  |         }  | 
271  | 0  |       } while ((iter++).next());  | 
272  | 0  |     }  | 
273  |  |  | 
274  | 0  |     return true;  | 
275  | 0  |   }  | 
276  |  |  | 
277  |  |   unsigned int OBReactionFacadePrivate::NumComponents(OBReactionRole rxnrole)  | 
278  | 0  |   { | 
279  | 0  |     std::vector<unsigned int> *component_ids = GetComponentIds(rxnrole);  | 
280  | 0  |     return (unsigned int)component_ids->size();  | 
281  | 0  |   }  | 
282  |  |  | 
283  |  |   bool OBReactionFacadePrivate::ReassignComponent(OBReactionRole oldrole, unsigned int num, OBReactionRole newrole)  | 
284  | 0  |   { | 
285  | 0  |     std::vector<unsigned int> *component_ids = GetComponentIds(oldrole);  | 
286  | 0  |     if (num >= component_ids->size())  | 
287  | 0  |       return false;  | 
288  | 0  |     unsigned int componentId = (*component_ids)[num];  | 
289  | 0  |     FOR_ATOMS_OF_MOL(atom, _mol) { | 
290  | 0  |       if (GetRole(&*atom) == oldrole && GetComponentId(&*atom) == componentId) { | 
291  | 0  |         SetRole(&*atom, newrole);  | 
292  | 0  |       }  | 
293  | 0  |     }  | 
294  |  |     // remove the entry from the original component_ids and update the new one  | 
295  | 0  |     component_ids->erase(component_ids->begin() + num);  | 
296  | 0  |     GetComponentIds(newrole)->push_back(componentId);  | 
297  | 0  |     return true;  | 
298  | 0  |   }  | 
299  |  |  | 
300  |  |   void OBReactionFacadePrivate::SetId(const char* idtype, OBAtom* atom, int idval)  | 
301  | 0  |   { | 
302  |  |     // Note: this replaces any existing data  | 
303  | 0  |     OBGenericData *data = atom->GetData(idtype);  | 
304  | 0  |     OBPairInteger *pi;  | 
305  | 0  |     if (data) { | 
306  | 0  |       pi = (OBPairInteger*)data;  | 
307  | 0  |       pi->SetValue(idval);  | 
308  | 0  |     }  | 
309  | 0  |     else { | 
310  | 0  |       pi = new OBPairInteger();  | 
311  | 0  |       pi->SetAttribute(idtype);  | 
312  | 0  |       pi->SetValue(idval);  | 
313  | 0  |       atom->SetData(pi);  | 
314  | 0  |     }  | 
315  | 0  |   }  | 
316  |  |  | 
317  |  |   void OBReactionFacadePrivate::SetComponentId(OBAtom* atom, unsigned int compid)  | 
318  | 0  |   { | 
319  | 0  |     SetId("rxncomp", atom, compid); | 
320  | 0  |   }  | 
321  |  |  | 
322  |  |   void OBReactionFacadePrivate::SetRole(OBAtom* atom, OBReactionRole rxnrole)  | 
323  | 0  |   { | 
324  | 0  |     SetId("rxnrole", atom, rxnrole); | 
325  | 0  |   }  | 
326  |  |  | 
327  |  |     | 
328  |  |     | 
329  |  |   // OBReactionFacade methods  | 
330  |  |   // These are all forwarded to the private implementation (pimpl idiom)  | 
331  | 0  |   OBReactionFacade::OBReactionFacade(OBMol *mol) : d(new OBReactionFacadePrivate(mol))  | 
332  | 0  |   { | 
333  | 0  |   };  | 
334  |  |   OBReactionFacade::~OBReactionFacade()  | 
335  | 0  |   { | 
336  | 0  |     delete d;  | 
337  | 0  |   }  | 
338  |  |   void OBReactionFacade::AddComponent(OBMol* mol, OBReactionRole rxnrole)  | 
339  | 0  |   { | 
340  | 0  |     d->AddComponent(mol, rxnrole);  | 
341  | 0  |   }  | 
342  |  |   void OBReactionFacade::AssignComponentIds(bool wipe)  | 
343  | 0  |   { | 
344  | 0  |     d->AssignComponentIds(wipe);  | 
345  | 0  |   }  | 
346  |  |   void OBReactionFacade::ClearInternalState()  | 
347  | 0  |   { | 
348  | 0  |     d->ClearInternalState();  | 
349  | 0  |   }  | 
350  |  |   bool OBReactionFacade::GetComponent(OBMol* mol, OBReactionRole rxnrole, unsigned int num)  | 
351  | 0  |   { | 
352  | 0  |     return d->GetComponent(mol, rxnrole, num);  | 
353  | 0  |   }  | 
354  |  |   unsigned int OBReactionFacade::GetComponentId(OBAtom *atom)  | 
355  | 0  |   { | 
356  | 0  |     return d->GetComponentId(atom);  | 
357  | 0  |   }  | 
358  |  |   OBReactionRole OBReactionFacade::GetRole(OBAtom *atom)  | 
359  | 0  |   { | 
360  | 0  |     return d->GetRole(atom);  | 
361  | 0  |   }  | 
362  |  |   bool OBReactionFacade::IsValid()  | 
363  | 0  |   { | 
364  | 0  |     return d->IsValid();  | 
365  | 0  |   }  | 
366  |  |   unsigned int OBReactionFacade::NumComponents(OBReactionRole rxnrole)  | 
367  | 0  |   { | 
368  | 0  |     return d->NumComponents(rxnrole);  | 
369  | 0  |   }  | 
370  |  |   bool OBReactionFacade::ReassignComponent(OBReactionRole oldrole, unsigned int num, OBReactionRole newrole)  | 
371  | 0  |   { | 
372  | 0  |     return d->ReassignComponent(oldrole, num, newrole);  | 
373  | 0  |   }  | 
374  |  |   void OBReactionFacade::SetComponentId(OBAtom* atom, unsigned int compid)  | 
375  | 0  |   { | 
376  | 0  |     d->SetComponentId(atom, compid);  | 
377  | 0  |   }  | 
378  |  |   void OBReactionFacade::SetRole(OBAtom* atom, OBReactionRole rxnrole)  | 
379  | 0  |   { | 
380  | 0  |     d->SetRole(atom, rxnrole);  | 
381  | 0  |   }  | 
382  |  |  | 
383  |  | } // namespace OpenBabel  | 
384  |  |  | 
385  |  |  | 
386  |  | //! \file reactionfacade.cpp  | 
387  |  | //! \brief A facade class to help interrogate and manipulate reactions  |