src/openbabel/ring.cpp

/**********************************************************************
ring.cpp - Deal with rings, find smallest set of smallest rings (SSSR).
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#include "mol.hpp"
#include <deque>

#include "obutil.hpp"

using namespace std;

namespace OpenBabel
{

/*! \class OBRing
    \brief Stores information on rings in a molecule from SSSR perception.

Ring information beyond atom and bond membership is usually not
necessary, but more information can be had about the rings in a
molecule. The OBRing class is used to store the information from a
Smallest Set of Smallest Rings (SSSR) perception of a molecule. The
OBMol member function OBMol::GetSSSR() stores the information it perceives in
a vector<OBRing*> inside the molecule. Perception is only done once
for a molecule unless the connection table is modified. The following
code demonstrates how to extract the SSSR information:
\code
OBMol mol;

vector<OBRing*> vr;
vr = mol.GetSSSR();
\endcode
OBRings store the atom numbers of the atoms in each of the smallest
set of smallest rings in both a vector<int> and an OBBitVec.
An example of how to print out the atom numbers present in all SSSR
rings is show below:
\code
vector<OBRing*>::iterator i;
vector<int>::iterator j;
vector<OBRing*> *rlist = (vector<OBRing*>*)mol.GetData("RingList");
for (i = rlist->begin();i != rlist->end();i++)
{
for(j = (*i)->_path.begin();j != (*i)->_path.end();j++)
    cout << *j << ` `;
cout << endl;
}
\endcode
will produce something like the following output for benzene:
\code
1 2 3 4 5 6
\endcode
Ring information is automatically deleted from an OBMol when it goes
out of scope or the Clear() member function is called.

Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#findSmallestSetOfSmallestRings">blue-obelisk:findSmallestSetOfSmallestRings</a>.
 */

static int DetermineFRJ(OBMol &);
static void BuildOBRTreeVector(OBAtom*,OBRTree*,vector<OBRTree*>&,OBBitVec&);


void OBMol::FindSSSR()
{
    if (HasSSSRPerceived())
        return;
    SetSSSRPerceived();
    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::FindSSSR", obAuditMsg);

    OBRing *ring;
    vector<OBRing*>::iterator j;

    //get frerejaque taking int account multiple possible spanning graphs
    int frj = DetermineFRJ(*this);
    if (frj)
    {
        vector<OBRing*> vr;
        FindRingAtomsAndBonds();

        OBBond *bond;
        vector<OBEdgeBase*> cbonds;
        vector<OBEdgeBase*>::iterator k;

        //restrict search for rings around closure bonds
        for (bond = BeginBond(k);bond;bond = NextBond(k))
            if (bond->IsClosure())
                cbonds.push_back(bond);

        if (!cbonds.empty())
        {
            OBRingSearch rs;
            //search for all rings about closures
            vector<OBEdgeBase*>::iterator i;

            for (i = cbonds.begin();i != cbonds.end();i++)
                rs.AddRingFromClosure(*this,(OBBond*)*i);

            rs.SortRings();
            rs.RemoveRedundant(frj);
            //store the SSSR set

            for (j = rs.BeginRings();j != rs.EndRings();j++)
            {
                ring = new OBRing ((*j)->_path,NumAtoms()+1);
                ring->SetParent(this);
                vr.push_back(ring);
            }
            //rs.WriteRings();
        }

        if (!HasData(OBGenericDataType::RingData))
            SetData(new OBRingData);
        OBRingData *rd = (OBRingData*)GetData(OBGenericDataType::RingData);
        rd->SetData(vr);
    }
}

static int DetermineFRJ(OBMol &mol)
{
    vector<vector<int> >::iterator i;
    vector<vector<int> > cfl;
    //find all continuous graphs in the mol area
    mol.ContigFragList(cfl);

    if (cfl.empty())
        return(0);
    if (cfl.size() == 1)
        return(mol.NumBonds() - mol.NumAtoms() + 1);

    //count up the atoms and bonds belonging to each graph
    OBBond *bond;
    vector<OBEdgeBase*>::iterator j;
    int numatoms,numbonds,frj=0;
    OBBitVec frag;
    for (i = cfl.begin();i != cfl.end();i++)
    {
        frag.Clear();
        frag.FromVecInt(*i);
        numatoms = (*i).size();
        numbonds=0;
        for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
            if (frag.BitIsOn(bond->GetBeginAtomIdx()) &&
                    frag.BitIsOn(bond->GetEndAtomIdx()))
                numbonds++;
        frj += numbonds - numatoms + 1;
    }

    return(frj);
}

void OBRingSearch::RemoveRedundant(int frj)
{
    OBBitVec tmp;
    int i,j;

    //remove identical rings
    for (i = _rlist.size()-1;i > 0;i--)
        for (j = i-1;j >= 0;j--)
            if ((_rlist[i])->_pathset == (_rlist[j])->_pathset)
            {
                delete _rlist[i];
                _rlist.erase(_rlist.begin()+i);
                break;
            }

    //make sure tmp is the same size as the rings
    for (j = 0;j < (signed)_rlist.size();j++)
        tmp = (_rlist[j])->_pathset;

    //remove larger rings that cover the same atoms as smaller rings
    for (i = _rlist.size()-1;i >= 0;i--)
    {
        tmp.Clear();
        for (j = 0;j < (signed)_rlist.size();j++)
            if ((_rlist[j])->_path.size() <= (_rlist[i])->_path.size() && i != j)
                tmp |= (_rlist[j])->_pathset;

        tmp = tmp & (_rlist[i])->_pathset;

        if (tmp == (_rlist[i])->_pathset)
        {
            delete _rlist[i];
            _rlist.erase(_rlist.begin()+i);
        }

        if (_rlist.size() == (unsigned)frj)
            break;
    }
}


void OBRingSearch::AddRingFromClosure(OBMol &mol,OBBond *cbond)
{
    vector<OBRTree*> t1(mol.NumAtoms()+1,(OBRTree*)NULL);
    vector<OBRTree*> t2(mol.NumAtoms()+1,(OBRTree*)NULL);
    OBBitVec bv1,bv2;

    bv1.SetBitOn(cbond->GetEndAtomIdx());
    bv2.SetBitOn(cbond->GetBeginAtomIdx());
    BuildOBRTreeVector(cbond->GetBeginAtom(),NULL,t1,bv1);
    BuildOBRTreeVector(cbond->GetEndAtom(),NULL,t2,bv2);

    bool pathok;
    deque<int> p1,p2;
    vector<OBNodeBase*> path1,path2;
    vector<OBNodeBase*>::iterator m,n;
    vector<OBRTree*>::iterator i;

    for (i = t1.begin();i != t1.end();i++)
        if (*i)
        {
            path1.clear();
            (*i)->PathToRoot(path1);

            if (t2[(*i)->GetAtomIdx()])
            {
                pathok = true;
                path2.clear();
                t2[(*i)->GetAtomIdx()]->PathToRoot(path2);

                p1.clear();
                m = path1.begin();
                if (m != path1.end())
                    p1.push_back((*m)->GetIdx());
                for (m = path1.begin(),m++;m != path1.end();m++)
                {
                    p1.push_back((*m)->GetIdx());
                    p2.clear();
                    for (n = path2.begin(),n++;n != path2.end();n++)
                    {
                        p2.push_front((*n)->GetIdx());
                        if (*n == *m)//don't traverse across identical atoms
                        {
                            p2.pop_front();
                            if (p1.size()+p2.size() > 2)
                                SaveUniqueRing(p1,p2);
                            pathok = false;
                            break;
                        }
                        if ((*n)->IsConnected(*m) && p1.size()+p2.size() > 2)
                            SaveUniqueRing(p1,p2);
                    }
                    if (!pathok)
                        break;
                }
            }
        }

    //clean up OBRTree vectors
    for (i = t1.begin();i != t1.end();i++)
        if (*i)
            delete *i;

    for (i = t2.begin();i != t2.end();i++)
        if (*i)
            delete *i;
}

bool OBRingSearch::SaveUniqueRing(deque<int> &d1,deque<int> &d2)
{
    vector<int> path;
    OBBitVec bv;
    deque<int>::iterator i;

    for (i = d1.begin();i != d1.end();i++)
    {
        bv.SetBitOn(*i);
        path.push_back(*i);
    }

    for (i = d2.begin();i != d2.end();i++)
    {
        bv.SetBitOn(*i);
        path.push_back(*i);
    }

    vector<OBRing*>::iterator j;
    for (j = _rlist.begin();j != _rlist.end();j++)
        if (bv == (*j)->_pathset)
            return(false);

    OBRing *ring = new OBRing;
    ring->_path = path;
    ring->_pathset = bv;
    _rlist.push_back(ring);

    return(true);
}

OBRingSearch::~OBRingSearch()
{
    vector<OBRing*>::iterator i;
    for (i = _rlist.begin();i != _rlist.end();i++)
        delete *i;
}

bool CompareRingSize(const OBRing *a,const OBRing *b)
{
    return(a->PathSize() < b->PathSize());
}

void OBRingSearch::WriteRings()
{
    vector<OBRing*>::iterator i;

    for (i = _rlist.begin();i != _rlist.end();i++)
        cout << (*i)->_pathset << endl;
}

static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit,
                      OBBitVec &bvisit, int natom,int depth );

void OBMol::FindRingAtomsAndBonds()
{
    if (HasFlag(OB_RINGFLAGS_MOL))
        return;
    SetFlag(OB_RINGFLAGS_MOL);

    obErrorLog.ThrowError(__FUNCTION__,
                          "Ran OpenBabel::FindRingAtomsAndBonds", obAuditMsg);

    OBBitVec avisit,bvisit;
    avisit.Resize(NumAtoms()+1);
    bvisit.Resize(NumAtoms()+1);
    vector<int> path;
    path.resize(NumAtoms()+1);

    for(unsigned int i=1; i<= NumAtoms(); i++ )
        if(!avisit[i])
            FindRings(*this,path,avisit,bvisit,i,0);
}

static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit,
                      OBBitVec &bvisit, int natom,int depth )
{
    OBAtom *atom;
    OBBond *bond;
    vector<OBEdgeBase*>::iterator k;

    if (avisit[natom])
    {
        int j = depth-1;
        bond=mol.GetBond(path[j--]);
        bond->SetInRing();
        while( j >= 0 )
        {
            bond=mol.GetBond(path[j--]);
            bond->SetInRing();
            (bond->GetBeginAtom())->SetInRing();
            (bond->GetEndAtom())->SetInRing();
            if(bond->GetBeginAtomIdx()==static_cast<unsigned int>(natom) || bond->
                    GetEndAtomIdx()==static_cast<unsigned int>(natom))
                break;
        }
    }
    else
    {
        avisit.SetBitOn(natom);
        atom = mol.GetAtom(natom);
        for(bond = atom->BeginBond(k);bond;bond=atom->NextBond(k))
            if( !bvisit[bond->GetIdx()])
            {
                path[depth] = bond->GetIdx();
                bvisit.SetBitOn(bond->GetIdx());
                FindRings(mol,path,avisit,bvisit,bond->GetNbrAtomIdx(atom),
                          depth+1);
            }
    }
}

bool OBRing::IsAromatic()
{
    OBMol *mol = _parent;
    vector<int>::iterator i;
    for (i = _path.begin();i != _path.end();i++)
        if (!(mol->GetAtom(*i))->IsAromatic())
            return(false);

    return(true);
}

bool OBRing::IsMember(OBAtom *a)
{
    return(_pathset.BitIsOn(a->GetIdx()));
}

bool OBRing::IsMember(OBBond *b)
{
    return((_pathset.BitIsOn(b->GetBeginAtomIdx()))&&(_pathset.BitIsOn(b->GetEndAtomIdx())));
}

OBRing::OBRing(vector<int> &path,int size)
{
    _path = path;
    _pathset.FromVecInt(_path);
    _pathset.Resize(size);
}

/*!
**\brief OBRing copy constructor
**\param src reference to original OBRing object (rhs)
*/
OBRing::OBRing(const OBRing &src)
//no base class
        :
        _path(src._path)
        ,       _pathset(src._pathset)  //chain to member classes
{
    //member data
    _parent = src._parent;  //this is messed up, but what can you do?
}

/*!
**\brief OBRing assignment operator
**\param src reference to original OBRing object (rhs)
**\return reference to modified OBRing object (lhs)
*/
OBRing& OBRing::operator =(const OBRing &src)
{
    //on identity, return
    if(this == &src)
        return(*this);

    //no base class

    //memeber classes & data
    _path = src._path;
    _pathset = src._pathset;
    _parent = src._parent; //note, this may not be what you want

    return(*this);
}
void BuildOBRTreeVector(OBAtom *atom,OBRTree *prv,vector<OBRTree*> &vt,OBBitVec &bv)
{
    vt[atom->GetIdx()] = new OBRTree (atom,prv);

    int i;
    OBAtom *nbr;
    OBMol *mol = (OBMol*)atom->GetParent();
    OBBitVec curr,used,next;
    vector<OBEdgeBase*>::iterator j;
    curr |= atom->GetIdx();
    used = bv|curr;

#define OB_RTREE_CUTOFF 20

    int level=0;
    for (;;)
    {
        next.Clear();
        for (i = curr.NextBit(0);i != bv.EndBit();i = curr.NextBit(i))
        {
            atom = mol->GetAtom(i);
            for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
                if (!used[nbr->GetIdx()])
                {
                    next |= nbr->GetIdx();
                    used |= nbr->GetIdx();
                    vt[nbr->GetIdx()] = new OBRTree (nbr,vt[atom->GetIdx()]);
                }
        }

        if (next.Empty())
            break;
        curr = next;
        level++;
        if (level > OB_RTREE_CUTOFF)
            break;
    }
#undef OB_RTREE_CUTOFF
}

OBRTree::OBRTree(OBAtom *atom,OBRTree *prv)
{
    _atom = atom;
    _prv = prv;
}

//! The supplied path is built up of OBAtom nodes, with the root atom 
//! the last item in the vector.
void OBRTree::PathToRoot(vector<OBNodeBase*> &path)
{
    path.push_back(_atom);
    if (_prv)
        _prv->PathToRoot(path);
}

int OBRTree::GetAtomIdx()
{
    return(_atom->GetIdx());
}

bool OBRing::findCenterAndNormal(vector3 & center, vector3 &norm1, vector3 &norm2)
{
    OBMol *mol= this->_parent;
    int j= 0;
    const int nA= this->_path.size();
    vector3 tmp;

    center.Set(0.0,0.0,0.0);
    norm1.Set(0.0,0.0,0.0);
    norm2.Set(0.0,0.0,0.0);
    for (j = 0; j != nA; j++)
    {
        center += (mol->GetAtom(_path[j]))->GetVector();
    }
    center/= double(nA);

    for (j = 0; j != nA; j++)
    {
        vector3 v1= (mol->GetAtom(_path[j]))->GetVector() - center;
        vector3 v2= (mol->GetAtom(_path[j+1==nA?0:j+1]))->GetVector() - center;
        tmp= cross(v1,v2);
        norm1+= tmp;
    }
    norm1/= double(nA);
    norm1.normalize();
    norm2= norm1;
    norm2 *= -1.0;
    return(true);
}

} // end namespace OpenBabel

//! \file ring.cpp
//! \brief Deal with rings, find smallest set of smallest rings (SSSR).
Revision:	741
Committed:	Wed Nov 16 19:42:11 2005 UTC (19 years, 11 months ago) by tim
File size:	14843 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	741	/**********************************************************************
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).