src/openbabel/bondtyper.cpp

/**********************************************************************
bondtyper.cpp - Bond typer to perceive connectivity and bond orders/types.

Copyright (C) 2003-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 "bondtyper.hpp"
#include "bondtyp.hpp"

#ifdef WIN32
#pragma warning (disable : 4786)
#endif

using namespace std;

namespace OpenBabel
{

  //! Global OBBondTyper for perception of bond order assignment.
OBBondTyper  bondtyper;

/*! \class OBBondTyper
    \brief Assigns bond types for file formats without bond information

The OBBondTyper class is designed to read in a list of bond typing
rules and apply them to molecules.
*/
OBBondTyper::OBBondTyper()
{
    _init = false;
    _dir = BABEL_DATADIR;
    _envvar = "BABEL_DATADIR";
    _filename = "bondtyp.txt";
    _subdir = "data";
    _dataptr = BondTypeData;
}

void OBBondTyper::ParseLine(const char *buffer)
{
    vector<string> vs;
    vector<int>    bovector;
    OBSmartsPattern *sp;

    if (buffer[0] != '#')
    {
        tokenize(vs,buffer);
        // Make sure we actually have a SMARTS pattern plus at least one triple
        // and make sure we have the correct number of integers
        if (vs.empty() || vs.size() < 4)
          return; // just ignore empty (or short lines)
        else if (!vs.empty() && vs.size() >= 4 && (vs.size() % 3 != 1))
        {
#ifdef HAVE_SSTREAM
          stringstream errorMsg;
#else
          strstream errorMsg;
#endif
          errorMsg << " Error in OBBondTyper. Pattern is incorrect, found "
                   << vs.size() << " tokens." << endl;
          errorMsg << " Buffer is: " << buffer << endl;
          obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obInfo);
          return;
        }

        sp = new OBSmartsPattern;
        if (sp->Init(vs[0]))
        {
            for (unsigned int i = 1; i < vs.size() ; i++)
            {
                bovector.push_back( atoi((char *)vs[i].c_str()) );
            }
            _fgbonds.push_back(pair<OBSmartsPattern*,vector<int> >
                               (sp, bovector));
        }
        else
        {
            delete sp;
            sp = NULL;
        }
    }
}

OBBondTyper::~OBBondTyper()
{
    vector<pair<OBSmartsPattern*, vector<int> > >::iterator i;
    for (i = _fgbonds.begin();i != _fgbonds.end();i++)
    {
        delete i->first;
        i->first = NULL;
    }
}

void OBBondTyper::AssignFunctionalGroupBonds(OBMol &mol)
{
  if (!_init)
    Init();

  OBSmartsPattern *currentPattern;
  OBBond *b1, *b2;
  OBAtom *a1,*a2, *a3;
  double angle, dist1, dist2;
  vector<int> assignments;
  vector<vector<int> > mlist;
  vector<vector<int> >::iterator matches, l;
  vector<pair<OBSmartsPattern*, vector<int> > >::iterator i;  
  unsigned int j;

  // Loop through for all the functional groups and assign bond orders
  for (i = _fgbonds.begin();i != _fgbonds.end();i++)
    {
      currentPattern = i->first;
      assignments = i->second;

      if (currentPattern && currentPattern->Match(mol))
        {
          mlist = currentPattern->GetUMapList();
          for (matches = mlist.begin(); matches != mlist.end(); matches++)
            {
              // Now loop through the bonds to assign from _fgbonds
              for (j = 0; j < assignments.size(); j += 3)
                {
                  // along the assignments vector: atomID1 atomID2 bondOrder
                  a1 = mol.GetAtom((*matches)[ assignments[j] ]);
                  a2 = mol.GetAtom((*matches)[ assignments[j+1 ] ]);
                  if (!a1 || !a2) continue;

                  b1 = a1->GetBond(a2);

                  if (!b1) continue;
                  b1->SetBO(assignments[j+2]);
                } // bond order assignments
            } // each match
        } // current pattern matches

    } // for(functional groups)

  // FG with distance and/or bond criteria
  // Carbonyl oxygen C=O
  OBSmartsPattern carbo; carbo.Init("[#8D1][#6](*)(*)");
  
  if (carbo.Match(mol))
    {
      mlist = carbo.GetUMapList();
      for (l = mlist.begin(); l != mlist.end(); l++)
        {
          a1 = mol.GetAtom((*l)[0]);
          a2 = mol.GetAtom((*l)[1]);
           
          angle = a2->AverageBondAngle();
          dist1 = a1->GetDistance(a2);
            
          // carbonyl geometries ?
          if (angle > 115 && angle < 150 && dist1 < 1.28) { 
                
            if ( !a1->HasDoubleBond() ) {// no double bond already assigned
              b1 = a1->GetBond(a2); 

              if (!b1 ) continue;
              b1->SetBO(2);
            }
          }
        }
    } // Carbonyl oxygen

  // thione C=S
  OBSmartsPattern thione; thione.Init("[#16D1][#6](*)(*)");
  
  if (thione.Match(mol))
    {
      mlist = thione.GetUMapList();
      for (l = mlist.begin(); l != mlist.end(); l++)
        {
          a1 = mol.GetAtom((*l)[0]);
          a2 = mol.GetAtom((*l)[1]);
           
          angle = a2->AverageBondAngle();
          dist1 = a1->GetDistance(a2);
            
          // thione geometries ?
          if (angle > 115 && angle < 150 && dist1 < 1.72) { 
                
            if ( !a1->HasDoubleBond() ) {// no double bond already assigned
              b1 = a1->GetBond(a2); 

              if (!b1 ) continue;
              b1->SetBO(2);
            }
          }
        }
    } // thione

  // Isocyanate N=C=O or Isothiocyanate
  bool dist1OK;
  OBSmartsPattern isocyanate; isocyanate.Init("[#8,#16;D1][#6D2][#7D2]");
  if (isocyanate.Match(mol))
    {
      mlist = isocyanate.GetUMapList();
      for (l = mlist.begin(); l != mlist.end(); l++)
        {
          a1 = mol.GetAtom((*l)[0]);
          a2 = mol.GetAtom((*l)[1]);
          a3 = mol.GetAtom((*l)[2]);

          angle = a2->AverageBondAngle();
          dist1 = a1->GetDistance(a2);
          dist2 = a2->GetDistance(a3);
            
          // isocyanate geometry or Isotiocyanate geometry ?
          if (a1->IsOxygen()) 
            dist1OK =  dist1 < 1.28;
          else
            dist1OK =  dist1 < 1.72;
            
          if (angle > 150 && dist1OK && dist2 < 1.34) { 
                  
            b1 = a1->GetBond(a2); 
            b2 = a2->GetBond(a3);
            if (!b1 || !b2) continue;
            b1->SetBO(2);
            b2->SetBO(2);
                
          }
            
        }      
    } // Isocyanate

}

} //namespace OpenBabel;

//! \file bondtyper.cpp
//! \brief Bond typer to perceive connectivity and bond orders/types.
//! \todo Needs to add aromatic ring bond order assignment.
//!   Eventually need to migrate OBMol::PerceiveBondOrders(),
//!   OBMol::ConnectTheDots(), and possibly some of the Kekulize routines
Revision:	2440
Committed:	Wed Nov 16 19:42:11 2005 UTC (18 years, 8 months ago) by tim
File size:	6719 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	2440	/**********************************************************************
2			bondtyper.cpp - Bond typer to perceive connectivity and bond orders/types.
3
4			Copyright (C) 2003-2005 by Geoffrey R. Hutchison
5
6			This file is part of the Open Babel project.
7			For more information, see <http://openbabel.sourceforge.net/>
8
9			This program is free software; you can redistribute it and/or modify
10			it under the terms of the GNU General Public License as published by
11			the Free Software Foundation version 2 of the License.
12
13			This program is distributed in the hope that it will be useful,
14			but WITHOUT ANY WARRANTY; without even the implied warranty of
15			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16			GNU General Public License for more details.
17			***********************************************************************/
18
19			#include "mol.hpp"
20			#include "bondtyper.hpp"
21			#include "bondtyp.hpp"
22
23			#ifdef WIN32
24			#pragma warning (disable : 4786)
25			#endif
26
27			using namespace std;
28
29			namespace OpenBabel
30			{
31
32			//! Global OBBondTyper for perception of bond order assignment.
33			OBBondTyper bondtyper;
34
35			/*! \class OBBondTyper
36			\brief Assigns bond types for file formats without bond information
37
38			The OBBondTyper class is designed to read in a list of bond typing
39			rules and apply them to molecules.
40			*/
41			OBBondTyper::OBBondTyper()
42			{
43			_init = false;
44			_dir = BABEL_DATADIR;
45			_envvar = "BABEL_DATADIR";
46			_filename = "bondtyp.txt";
47			_subdir = "data";
48			_dataptr = BondTypeData;
49			}
50
51			void OBBondTyper::ParseLine(const char *buffer)
52			{
53			vector<string> vs;
54			vector<int> bovector;
55			OBSmartsPattern *sp;
56
57			if (buffer[0] != '#')
58			{
59			tokenize(vs,buffer);
60			// Make sure we actually have a SMARTS pattern plus at least one triple
61			// and make sure we have the correct number of integers
62			if (vs.empty() \|\| vs.size() < 4)
63			return; // just ignore empty (or short lines)
64			else if (!vs.empty() && vs.size() >= 4 && (vs.size() % 3 != 1))
65			{
66			#ifdef HAVE_SSTREAM
67			stringstream errorMsg;
68			#else
69			strstream errorMsg;
70			#endif
71			errorMsg << " Error in OBBondTyper. Pattern is incorrect, found "
72			<< vs.size() << " tokens." << endl;
73			errorMsg << " Buffer is: " << buffer << endl;
74			obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obInfo);
75			return;
76			}
77
78			sp = new OBSmartsPattern;
79			if (sp->Init(vs[0]))
80			{
81			for (unsigned int i = 1; i < vs.size() ; i++)
82			{
83			bovector.push_back( atoi((char *)vs[i].c_str()) );
84			}
85			_fgbonds.push_back(pair<OBSmartsPattern*,vector<int> >
86			(sp, bovector));
87			}
88			else
89			{
90			delete sp;
91			sp = NULL;
92			}
93			}
94			}
95
96			OBBondTyper::~OBBondTyper()
97			{
98			vector<pair<OBSmartsPattern*, vector<int> > >::iterator i;
99			for (i = _fgbonds.begin();i != _fgbonds.end();i++)
100			{
101			delete i->first;
102			i->first = NULL;
103			}
104			}
105
106			void OBBondTyper::AssignFunctionalGroupBonds(OBMol &mol)
107			{
108			if (!_init)
109			Init();
110
111			OBSmartsPattern *currentPattern;
112			OBBond b1, b2;
113			OBAtom a1,a2, *a3;
114			double angle, dist1, dist2;
115			vector<int> assignments;
116			vector<vector<int> > mlist;
117			vector<vector<int> >::iterator matches, l;
118			vector<pair<OBSmartsPattern*, vector<int> > >::iterator i;
119			unsigned int j;
120
121			// Loop through for all the functional groups and assign bond orders
122			for (i = _fgbonds.begin();i != _fgbonds.end();i++)
123			{
124			currentPattern = i->first;
125			assignments = i->second;
126
127			if (currentPattern && currentPattern->Match(mol))
128			{
129			mlist = currentPattern->GetUMapList();
130			for (matches = mlist.begin(); matches != mlist.end(); matches++)
131			{
132			// Now loop through the bonds to assign from _fgbonds
133			for (j = 0; j < assignments.size(); j += 3)
134			{
135			// along the assignments vector: atomID1 atomID2 bondOrder
136			a1 = mol.GetAtom((*matches)[ assignments[j] ]);
137			a2 = mol.GetAtom((*matches)[ assignments[j+1 ] ]);
138			if (!a1 \|\| !a2) continue;
139
140			b1 = a1->GetBond(a2);
141
142			if (!b1) continue;
143			b1->SetBO(assignments[j+2]);
144			} // bond order assignments
145			} // each match
146			} // current pattern matches
147
148			} // for(functional groups)
149
150			// FG with distance and/or bond criteria
151			// Carbonyl oxygen C=O
152			OBSmartsPattern carbo; carbo.Init("[#8D1][#6]()()");
153
154			if (carbo.Match(mol))
155			{
156			mlist = carbo.GetUMapList();
157			for (l = mlist.begin(); l != mlist.end(); l++)
158			{
159			a1 = mol.GetAtom((*l)[0]);
160			a2 = mol.GetAtom((*l)[1]);
161
162			angle = a2->AverageBondAngle();
163			dist1 = a1->GetDistance(a2);
164
165			// carbonyl geometries ?
166			if (angle > 115 && angle < 150 && dist1 < 1.28) {
167
168			if ( !a1->HasDoubleBond() ) {// no double bond already assigned
169			b1 = a1->GetBond(a2);
170
171			if (!b1 ) continue;
172			b1->SetBO(2);
173			}
174			}
175			}
176			} // Carbonyl oxygen
177
178			// thione C=S
179			OBSmartsPattern thione; thione.Init("[#16D1][#6]()()");
180
181			if (thione.Match(mol))
182			{
183			mlist = thione.GetUMapList();
184			for (l = mlist.begin(); l != mlist.end(); l++)
185			{
186			a1 = mol.GetAtom((*l)[0]);
187			a2 = mol.GetAtom((*l)[1]);
188
189			angle = a2->AverageBondAngle();
190			dist1 = a1->GetDistance(a2);
191
192			// thione geometries ?
193			if (angle > 115 && angle < 150 && dist1 < 1.72) {
194
195			if ( !a1->HasDoubleBond() ) {// no double bond already assigned
196			b1 = a1->GetBond(a2);
197
198			if (!b1 ) continue;
199			b1->SetBO(2);
200			}
201			}
202			}
203			} // thione
204
205			// Isocyanate N=C=O or Isothiocyanate
206			bool dist1OK;
207			OBSmartsPattern isocyanate; isocyanate.Init("[#8,#16;D1][#6D2][#7D2]");
208			if (isocyanate.Match(mol))
209			{
210			mlist = isocyanate.GetUMapList();
211			for (l = mlist.begin(); l != mlist.end(); l++)
212			{
213			a1 = mol.GetAtom((*l)[0]);
214			a2 = mol.GetAtom((*l)[1]);
215			a3 = mol.GetAtom((*l)[2]);
216
217			angle = a2->AverageBondAngle();
218			dist1 = a1->GetDistance(a2);
219			dist2 = a2->GetDistance(a3);
220
221			// isocyanate geometry or Isotiocyanate geometry ?
222			if (a1->IsOxygen())
223			dist1OK = dist1 < 1.28;
224			else
225			dist1OK = dist1 < 1.72;
226
227			if (angle > 150 && dist1OK && dist2 < 1.34) {
228
229			b1 = a1->GetBond(a2);
230			b2 = a2->GetBond(a3);
231			if (!b1 \|\| !b2) continue;
232			b1->SetBO(2);
233			b2->SetBO(2);
234
235			}
236
237			}
238			} // Isocyanate
239
240			}
241
242			} //namespace OpenBabel;
243
244			//! \file bondtyper.cpp
245			//! \brief Bond typer to perceive connectivity and bond orders/types.
246			//! \todo Needs to add aromatic ring bond order assignment.
247			//! Eventually need to migrate OBMol::PerceiveBondOrders(),
248			//! OBMol::ConnectTheDots(), and possibly some of the Kekulize routines