ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/openbabel/bondtyper.cpp
Revision: 2450
Committed: Thu Nov 17 20:38:45 2005 UTC (18 years, 8 months ago) by gezelter
File size: 6729 byte(s)
Log Message:
Unifying config.h stuff and making sure the OpenBabel codes can find
our default (and environment variable) Force Field directories.

File Contents

# Content
1 /**********************************************************************
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 STR_DEFINE(_dir, FRC_PATH );
45 _envvar = "FORCE_PARAM_PATH";
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