ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/openbabel/ring.cpp
Revision: 2440
Committed: Wed Nov 16 19:42:11 2005 UTC (18 years, 8 months ago) by tim
File size: 14843 byte(s)
Log Message:
adding openbabel

File Contents

# Content
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(__FUNCTION__,
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(__FUNCTION__,
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).