1 |
/********************************************************************** |
2 |
ring.cpp - Deal with rings, find smallest set of smallest rings (SSSR). |
3 |
|
4 |
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. |
5 |
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison |
6 |
|
7 |
This file is part of the Open Babel project. |
8 |
For more information, see <http://openbabel.sourceforge.net/> |
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 "mol.hpp" |
21 |
#include <deque> |
22 |
|
23 |
#include "obutil.hpp" |
24 |
|
25 |
using namespace std; |
26 |
|
27 |
namespace OpenBabel |
28 |
{ |
29 |
|
30 |
/*! \class OBRing |
31 |
\brief Stores information on rings in a molecule from SSSR perception. |
32 |
|
33 |
Ring information beyond atom and bond membership is usually not |
34 |
necessary, but more information can be had about the rings in a |
35 |
molecule. The OBRing class is used to store the information from a |
36 |
Smallest Set of Smallest Rings (SSSR) perception of a molecule. The |
37 |
OBMol member function OBMol::GetSSSR() stores the information it perceives in |
38 |
a vector<OBRing*> inside the molecule. Perception is only done once |
39 |
for a molecule unless the connection table is modified. The following |
40 |
code demonstrates how to extract the SSSR information: |
41 |
\code |
42 |
OBMol mol; |
43 |
|
44 |
vector<OBRing*> vr; |
45 |
vr = mol.GetSSSR(); |
46 |
\endcode |
47 |
OBRings store the atom numbers of the atoms in each of the smallest |
48 |
set of smallest rings in both a vector<int> and an OBBitVec. |
49 |
An example of how to print out the atom numbers present in all SSSR |
50 |
rings is show below: |
51 |
\code |
52 |
vector<OBRing*>::iterator i; |
53 |
vector<int>::iterator j; |
54 |
vector<OBRing*> *rlist = (vector<OBRing*>*)mol.GetData("RingList"); |
55 |
for (i = rlist->begin();i != rlist->end();i++) |
56 |
{ |
57 |
for(j = (*i)->_path.begin();j != (*i)->_path.end();j++) |
58 |
cout << *j << ` `; |
59 |
cout << endl; |
60 |
} |
61 |
\endcode |
62 |
will produce something like the following output for benzene: |
63 |
\code |
64 |
1 2 3 4 5 6 |
65 |
\endcode |
66 |
Ring information is automatically deleted from an OBMol when it goes |
67 |
out of scope or the Clear() member function is called. |
68 |
|
69 |
Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#findSmallestSetOfSmallestRings">blue-obelisk:findSmallestSetOfSmallestRings</a>. |
70 |
*/ |
71 |
|
72 |
static int DetermineFRJ(OBMol &); |
73 |
static void BuildOBRTreeVector(OBAtom*,OBRTree*,vector<OBRTree*>&,OBBitVec&); |
74 |
|
75 |
|
76 |
void OBMol::FindSSSR() |
77 |
{ |
78 |
if (HasSSSRPerceived()) |
79 |
return; |
80 |
SetSSSRPerceived(); |
81 |
obErrorLog.ThrowError(__func__, |
82 |
"Ran OpenBabel::FindSSSR", obAuditMsg); |
83 |
|
84 |
OBRing *ring; |
85 |
vector<OBRing*>::iterator j; |
86 |
|
87 |
//get frerejaque taking int account multiple possible spanning graphs |
88 |
int frj = DetermineFRJ(*this); |
89 |
if (frj) |
90 |
{ |
91 |
vector<OBRing*> vr; |
92 |
FindRingAtomsAndBonds(); |
93 |
|
94 |
OBBond *bond; |
95 |
vector<OBEdgeBase*> cbonds; |
96 |
vector<OBEdgeBase*>::iterator k; |
97 |
|
98 |
//restrict search for rings around closure bonds |
99 |
for (bond = BeginBond(k);bond;bond = NextBond(k)) |
100 |
if (bond->IsClosure()) |
101 |
cbonds.push_back(bond); |
102 |
|
103 |
if (!cbonds.empty()) |
104 |
{ |
105 |
OBRingSearch rs; |
106 |
//search for all rings about closures |
107 |
vector<OBEdgeBase*>::iterator i; |
108 |
|
109 |
for (i = cbonds.begin();i != cbonds.end();i++) |
110 |
rs.AddRingFromClosure(*this,(OBBond*)*i); |
111 |
|
112 |
rs.SortRings(); |
113 |
rs.RemoveRedundant(frj); |
114 |
//store the SSSR set |
115 |
|
116 |
for (j = rs.BeginRings();j != rs.EndRings();j++) |
117 |
{ |
118 |
ring = new OBRing ((*j)->_path,NumAtoms()+1); |
119 |
ring->SetParent(this); |
120 |
vr.push_back(ring); |
121 |
} |
122 |
//rs.WriteRings(); |
123 |
} |
124 |
|
125 |
if (!HasData(OBGenericDataType::RingData)) |
126 |
SetData(new OBRingData); |
127 |
OBRingData *rd = (OBRingData*)GetData(OBGenericDataType::RingData); |
128 |
rd->SetData(vr); |
129 |
} |
130 |
} |
131 |
|
132 |
static int DetermineFRJ(OBMol &mol) |
133 |
{ |
134 |
vector<vector<int> >::iterator i; |
135 |
vector<vector<int> > cfl; |
136 |
//find all continuous graphs in the mol area |
137 |
mol.ContigFragList(cfl); |
138 |
|
139 |
if (cfl.empty()) |
140 |
return(0); |
141 |
if (cfl.size() == 1) |
142 |
return(mol.NumBonds() - mol.NumAtoms() + 1); |
143 |
|
144 |
//count up the atoms and bonds belonging to each graph |
145 |
OBBond *bond; |
146 |
vector<OBEdgeBase*>::iterator j; |
147 |
int numatoms,numbonds,frj=0; |
148 |
OBBitVec frag; |
149 |
for (i = cfl.begin();i != cfl.end();i++) |
150 |
{ |
151 |
frag.Clear(); |
152 |
frag.FromVecInt(*i); |
153 |
numatoms = (*i).size(); |
154 |
numbonds=0; |
155 |
for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j)) |
156 |
if (frag.BitIsOn(bond->GetBeginAtomIdx()) && |
157 |
frag.BitIsOn(bond->GetEndAtomIdx())) |
158 |
numbonds++; |
159 |
frj += numbonds - numatoms + 1; |
160 |
} |
161 |
|
162 |
return(frj); |
163 |
} |
164 |
|
165 |
void OBRingSearch::RemoveRedundant(int frj) |
166 |
{ |
167 |
OBBitVec tmp; |
168 |
int i,j; |
169 |
|
170 |
//remove identical rings |
171 |
for (i = _rlist.size()-1;i > 0;i--) |
172 |
for (j = i-1;j >= 0;j--) |
173 |
if ((_rlist[i])->_pathset == (_rlist[j])->_pathset) |
174 |
{ |
175 |
delete _rlist[i]; |
176 |
_rlist.erase(_rlist.begin()+i); |
177 |
break; |
178 |
} |
179 |
|
180 |
//make sure tmp is the same size as the rings |
181 |
for (j = 0;j < (signed)_rlist.size();j++) |
182 |
tmp = (_rlist[j])->_pathset; |
183 |
|
184 |
//remove larger rings that cover the same atoms as smaller rings |
185 |
for (i = _rlist.size()-1;i >= 0;i--) |
186 |
{ |
187 |
tmp.Clear(); |
188 |
for (j = 0;j < (signed)_rlist.size();j++) |
189 |
if ((_rlist[j])->_path.size() <= (_rlist[i])->_path.size() && i != j) |
190 |
tmp |= (_rlist[j])->_pathset; |
191 |
|
192 |
tmp = tmp & (_rlist[i])->_pathset; |
193 |
|
194 |
if (tmp == (_rlist[i])->_pathset) |
195 |
{ |
196 |
delete _rlist[i]; |
197 |
_rlist.erase(_rlist.begin()+i); |
198 |
} |
199 |
|
200 |
if (_rlist.size() == (unsigned)frj) |
201 |
break; |
202 |
} |
203 |
} |
204 |
|
205 |
|
206 |
void OBRingSearch::AddRingFromClosure(OBMol &mol,OBBond *cbond) |
207 |
{ |
208 |
vector<OBRTree*> t1(mol.NumAtoms()+1,(OBRTree*)NULL); |
209 |
vector<OBRTree*> t2(mol.NumAtoms()+1,(OBRTree*)NULL); |
210 |
OBBitVec bv1,bv2; |
211 |
|
212 |
bv1.SetBitOn(cbond->GetEndAtomIdx()); |
213 |
bv2.SetBitOn(cbond->GetBeginAtomIdx()); |
214 |
BuildOBRTreeVector(cbond->GetBeginAtom(),NULL,t1,bv1); |
215 |
BuildOBRTreeVector(cbond->GetEndAtom(),NULL,t2,bv2); |
216 |
|
217 |
bool pathok; |
218 |
deque<int> p1,p2; |
219 |
vector<OBNodeBase*> path1,path2; |
220 |
vector<OBNodeBase*>::iterator m,n; |
221 |
vector<OBRTree*>::iterator i; |
222 |
|
223 |
for (i = t1.begin();i != t1.end();i++) |
224 |
if (*i) |
225 |
{ |
226 |
path1.clear(); |
227 |
(*i)->PathToRoot(path1); |
228 |
|
229 |
if (t2[(*i)->GetAtomIdx()]) |
230 |
{ |
231 |
pathok = true; |
232 |
path2.clear(); |
233 |
t2[(*i)->GetAtomIdx()]->PathToRoot(path2); |
234 |
|
235 |
p1.clear(); |
236 |
m = path1.begin(); |
237 |
if (m != path1.end()) |
238 |
p1.push_back((*m)->GetIdx()); |
239 |
for (m = path1.begin(),m++;m != path1.end();m++) |
240 |
{ |
241 |
p1.push_back((*m)->GetIdx()); |
242 |
p2.clear(); |
243 |
for (n = path2.begin(),n++;n != path2.end();n++) |
244 |
{ |
245 |
p2.push_front((*n)->GetIdx()); |
246 |
if (*n == *m)//don't traverse across identical atoms |
247 |
{ |
248 |
p2.pop_front(); |
249 |
if (p1.size()+p2.size() > 2) |
250 |
SaveUniqueRing(p1,p2); |
251 |
pathok = false; |
252 |
break; |
253 |
} |
254 |
if ((*n)->IsConnected(*m) && p1.size()+p2.size() > 2) |
255 |
SaveUniqueRing(p1,p2); |
256 |
} |
257 |
if (!pathok) |
258 |
break; |
259 |
} |
260 |
} |
261 |
} |
262 |
|
263 |
//clean up OBRTree vectors |
264 |
for (i = t1.begin();i != t1.end();i++) |
265 |
if (*i) |
266 |
delete *i; |
267 |
|
268 |
for (i = t2.begin();i != t2.end();i++) |
269 |
if (*i) |
270 |
delete *i; |
271 |
} |
272 |
|
273 |
bool OBRingSearch::SaveUniqueRing(deque<int> &d1,deque<int> &d2) |
274 |
{ |
275 |
vector<int> path; |
276 |
OBBitVec bv; |
277 |
deque<int>::iterator i; |
278 |
|
279 |
for (i = d1.begin();i != d1.end();i++) |
280 |
{ |
281 |
bv.SetBitOn(*i); |
282 |
path.push_back(*i); |
283 |
} |
284 |
|
285 |
for (i = d2.begin();i != d2.end();i++) |
286 |
{ |
287 |
bv.SetBitOn(*i); |
288 |
path.push_back(*i); |
289 |
} |
290 |
|
291 |
vector<OBRing*>::iterator j; |
292 |
for (j = _rlist.begin();j != _rlist.end();j++) |
293 |
if (bv == (*j)->_pathset) |
294 |
return(false); |
295 |
|
296 |
OBRing *ring = new OBRing; |
297 |
ring->_path = path; |
298 |
ring->_pathset = bv; |
299 |
_rlist.push_back(ring); |
300 |
|
301 |
return(true); |
302 |
} |
303 |
|
304 |
OBRingSearch::~OBRingSearch() |
305 |
{ |
306 |
vector<OBRing*>::iterator i; |
307 |
for (i = _rlist.begin();i != _rlist.end();i++) |
308 |
delete *i; |
309 |
} |
310 |
|
311 |
bool CompareRingSize(const OBRing *a,const OBRing *b) |
312 |
{ |
313 |
return(a->PathSize() < b->PathSize()); |
314 |
} |
315 |
|
316 |
void OBRingSearch::WriteRings() |
317 |
{ |
318 |
vector<OBRing*>::iterator i; |
319 |
|
320 |
for (i = _rlist.begin();i != _rlist.end();i++) |
321 |
cout << (*i)->_pathset << endl; |
322 |
} |
323 |
|
324 |
static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit, |
325 |
OBBitVec &bvisit, int natom,int depth ); |
326 |
|
327 |
void OBMol::FindRingAtomsAndBonds() |
328 |
{ |
329 |
if (HasFlag(OB_RINGFLAGS_MOL)) |
330 |
return; |
331 |
SetFlag(OB_RINGFLAGS_MOL); |
332 |
|
333 |
obErrorLog.ThrowError(__func__, |
334 |
"Ran OpenBabel::FindRingAtomsAndBonds", obAuditMsg); |
335 |
|
336 |
OBBitVec avisit,bvisit; |
337 |
avisit.Resize(NumAtoms()+1); |
338 |
bvisit.Resize(NumAtoms()+1); |
339 |
vector<int> path; |
340 |
path.resize(NumAtoms()+1); |
341 |
|
342 |
for(unsigned int i=1; i<= NumAtoms(); i++ ) |
343 |
if(!avisit[i]) |
344 |
FindRings(*this,path,avisit,bvisit,i,0); |
345 |
} |
346 |
|
347 |
static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit, |
348 |
OBBitVec &bvisit, int natom,int depth ) |
349 |
{ |
350 |
OBAtom *atom; |
351 |
OBBond *bond; |
352 |
vector<OBEdgeBase*>::iterator k; |
353 |
|
354 |
if (avisit[natom]) |
355 |
{ |
356 |
int j = depth-1; |
357 |
bond=mol.GetBond(path[j--]); |
358 |
bond->SetInRing(); |
359 |
while( j >= 0 ) |
360 |
{ |
361 |
bond=mol.GetBond(path[j--]); |
362 |
bond->SetInRing(); |
363 |
(bond->GetBeginAtom())->SetInRing(); |
364 |
(bond->GetEndAtom())->SetInRing(); |
365 |
if(bond->GetBeginAtomIdx()==static_cast<unsigned int>(natom) || bond-> |
366 |
GetEndAtomIdx()==static_cast<unsigned int>(natom)) |
367 |
break; |
368 |
} |
369 |
} |
370 |
else |
371 |
{ |
372 |
avisit.SetBitOn(natom); |
373 |
atom = mol.GetAtom(natom); |
374 |
for(bond = atom->BeginBond(k);bond;bond=atom->NextBond(k)) |
375 |
if( !bvisit[bond->GetIdx()]) |
376 |
{ |
377 |
path[depth] = bond->GetIdx(); |
378 |
bvisit.SetBitOn(bond->GetIdx()); |
379 |
FindRings(mol,path,avisit,bvisit,bond->GetNbrAtomIdx(atom), |
380 |
depth+1); |
381 |
} |
382 |
} |
383 |
} |
384 |
|
385 |
bool OBRing::IsAromatic() |
386 |
{ |
387 |
OBMol *mol = _parent; |
388 |
vector<int>::iterator i; |
389 |
for (i = _path.begin();i != _path.end();i++) |
390 |
if (!(mol->GetAtom(*i))->IsAromatic()) |
391 |
return(false); |
392 |
|
393 |
return(true); |
394 |
} |
395 |
|
396 |
bool OBRing::IsMember(OBAtom *a) |
397 |
{ |
398 |
return(_pathset.BitIsOn(a->GetIdx())); |
399 |
} |
400 |
|
401 |
bool OBRing::IsMember(OBBond *b) |
402 |
{ |
403 |
return((_pathset.BitIsOn(b->GetBeginAtomIdx()))&&(_pathset.BitIsOn(b->GetEndAtomIdx()))); |
404 |
} |
405 |
|
406 |
OBRing::OBRing(vector<int> &path,int size) |
407 |
{ |
408 |
_path = path; |
409 |
_pathset.FromVecInt(_path); |
410 |
_pathset.Resize(size); |
411 |
} |
412 |
|
413 |
/*! |
414 |
**\brief OBRing copy constructor |
415 |
**\param src reference to original OBRing object (rhs) |
416 |
*/ |
417 |
OBRing::OBRing(const OBRing &src) |
418 |
//no base class |
419 |
: |
420 |
_path(src._path) |
421 |
, _pathset(src._pathset) //chain to member classes |
422 |
{ |
423 |
//member data |
424 |
_parent = src._parent; //this is messed up, but what can you do? |
425 |
} |
426 |
|
427 |
/*! |
428 |
**\brief OBRing assignment operator |
429 |
**\param src reference to original OBRing object (rhs) |
430 |
**\return reference to modified OBRing object (lhs) |
431 |
*/ |
432 |
OBRing& OBRing::operator =(const OBRing &src) |
433 |
{ |
434 |
//on identity, return |
435 |
if(this == &src) |
436 |
return(*this); |
437 |
|
438 |
//no base class |
439 |
|
440 |
//memeber classes & data |
441 |
_path = src._path; |
442 |
_pathset = src._pathset; |
443 |
_parent = src._parent; //note, this may not be what you want |
444 |
|
445 |
return(*this); |
446 |
} |
447 |
void BuildOBRTreeVector(OBAtom *atom,OBRTree *prv,vector<OBRTree*> &vt,OBBitVec &bv) |
448 |
{ |
449 |
vt[atom->GetIdx()] = new OBRTree (atom,prv); |
450 |
|
451 |
int i; |
452 |
OBAtom *nbr; |
453 |
OBMol *mol = (OBMol*)atom->GetParent(); |
454 |
OBBitVec curr,used,next; |
455 |
vector<OBEdgeBase*>::iterator j; |
456 |
curr |= atom->GetIdx(); |
457 |
used = bv|curr; |
458 |
|
459 |
#define OB_RTREE_CUTOFF 20 |
460 |
|
461 |
int level=0; |
462 |
for (;;) |
463 |
{ |
464 |
next.Clear(); |
465 |
for (i = curr.NextBit(0);i != bv.EndBit();i = curr.NextBit(i)) |
466 |
{ |
467 |
atom = mol->GetAtom(i); |
468 |
for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j)) |
469 |
if (!used[nbr->GetIdx()]) |
470 |
{ |
471 |
next |= nbr->GetIdx(); |
472 |
used |= nbr->GetIdx(); |
473 |
vt[nbr->GetIdx()] = new OBRTree (nbr,vt[atom->GetIdx()]); |
474 |
} |
475 |
} |
476 |
|
477 |
if (next.Empty()) |
478 |
break; |
479 |
curr = next; |
480 |
level++; |
481 |
if (level > OB_RTREE_CUTOFF) |
482 |
break; |
483 |
} |
484 |
#undef OB_RTREE_CUTOFF |
485 |
} |
486 |
|
487 |
OBRTree::OBRTree(OBAtom *atom,OBRTree *prv) |
488 |
{ |
489 |
_atom = atom; |
490 |
_prv = prv; |
491 |
} |
492 |
|
493 |
//! The supplied path is built up of OBAtom nodes, with the root atom |
494 |
//! the last item in the vector. |
495 |
void OBRTree::PathToRoot(vector<OBNodeBase*> &path) |
496 |
{ |
497 |
path.push_back(_atom); |
498 |
if (_prv) |
499 |
_prv->PathToRoot(path); |
500 |
} |
501 |
|
502 |
int OBRTree::GetAtomIdx() |
503 |
{ |
504 |
return(_atom->GetIdx()); |
505 |
} |
506 |
|
507 |
bool OBRing::findCenterAndNormal(vector3 & center, vector3 &norm1, vector3 &norm2) |
508 |
{ |
509 |
OBMol *mol= this->_parent; |
510 |
int j= 0; |
511 |
const int nA= this->_path.size(); |
512 |
vector3 tmp; |
513 |
|
514 |
center.Set(0.0,0.0,0.0); |
515 |
norm1.Set(0.0,0.0,0.0); |
516 |
norm2.Set(0.0,0.0,0.0); |
517 |
for (j = 0; j != nA; j++) |
518 |
{ |
519 |
center += (mol->GetAtom(_path[j]))->GetVector(); |
520 |
} |
521 |
center/= double(nA); |
522 |
|
523 |
for (j = 0; j != nA; j++) |
524 |
{ |
525 |
vector3 v1= (mol->GetAtom(_path[j]))->GetVector() - center; |
526 |
vector3 v2= (mol->GetAtom(_path[j+1==nA?0:j+1]))->GetVector() - center; |
527 |
tmp= cross(v1,v2); |
528 |
norm1+= tmp; |
529 |
} |
530 |
norm1/= double(nA); |
531 |
norm1.normalize(); |
532 |
norm2= norm1; |
533 |
norm2 *= -1.0; |
534 |
return(true); |
535 |
} |
536 |
|
537 |
} // end namespace OpenBabel |
538 |
|
539 |
//! \file ring.cpp |
540 |
//! \brief Deal with rings, find smallest set of smallest rings (SSSR). |