1 |
tim |
741 |
/********************************************************************* |
2 |
|
|
chiral.cpp - Detect chiral atoms and molecules. |
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 <list> |
21 |
|
|
|
22 |
|
|
#include "mol.hpp" |
23 |
|
|
#include "obutil.hpp" |
24 |
|
|
#include "matrix.hpp" |
25 |
|
|
#include "chiral.hpp" |
26 |
|
|
#include "matrix3x3.hpp" |
27 |
|
|
|
28 |
|
|
using namespace std; |
29 |
|
|
namespace OpenBabel |
30 |
|
|
{ |
31 |
|
|
|
32 |
|
|
//! Sets atom->IsChiral() to true for chiral atoms |
33 |
|
|
void OBMol::FindChiralCenters() |
34 |
|
|
{ |
35 |
|
|
if (HasChiralityPerceived()) |
36 |
|
|
return; |
37 |
|
|
SetChiralityPerceived(); |
38 |
|
|
|
39 |
tim |
819 |
obErrorLog.ThrowError(__func__, |
40 |
tim |
741 |
"Ran OpenBabel::FindChiralCenters", obAuditMsg); |
41 |
|
|
|
42 |
|
|
//do quick test to see if there are any possible chiral centers |
43 |
|
|
bool mayHaveChiralCenter=false; |
44 |
|
|
OBAtom *atom,*nbr; |
45 |
|
|
vector<OBNodeBase*>::iterator i; |
46 |
|
|
for (atom = BeginAtom(i);atom;atom = NextAtom(i)) |
47 |
|
|
if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3) |
48 |
|
|
{ |
49 |
|
|
mayHaveChiralCenter=true; |
50 |
|
|
break; |
51 |
|
|
} |
52 |
|
|
|
53 |
|
|
if (!mayHaveChiralCenter) |
54 |
|
|
return; |
55 |
|
|
|
56 |
|
|
OBBond *bond; |
57 |
|
|
vector<OBEdgeBase*>::iterator j; |
58 |
|
|
for (bond = BeginBond(j);bond;bond = NextBond(j)) |
59 |
|
|
if (bond->IsWedge() || bond->IsHash()) |
60 |
|
|
(bond->GetBeginAtom())->SetChiral(); |
61 |
|
|
|
62 |
|
|
#define INTMETHOD |
63 |
|
|
|
64 |
|
|
#ifdef INTMETHOD |
65 |
|
|
|
66 |
|
|
vector<unsigned int> vgid; |
67 |
|
|
GetGIDVector(vgid); |
68 |
|
|
vector<unsigned int> tlist; |
69 |
|
|
vector<unsigned int>::iterator k; |
70 |
|
|
#else //use Golender floating point method |
71 |
|
|
|
72 |
|
|
vector<double> gp; |
73 |
|
|
GraphPotentials(*this,gp); |
74 |
|
|
vector<double> tlist; |
75 |
|
|
vector<double>::iterator k; |
76 |
|
|
#endif |
77 |
|
|
|
78 |
|
|
bool ischiral; |
79 |
|
|
for (atom = BeginAtom(i);atom;atom = NextAtom(i)) |
80 |
|
|
if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral()) |
81 |
|
|
{ |
82 |
|
|
tlist.clear(); |
83 |
|
|
ischiral = true; |
84 |
|
|
|
85 |
|
|
for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j)) |
86 |
|
|
{ |
87 |
|
|
for (k = tlist.begin();k != tlist.end();k++) |
88 |
|
|
#ifdef INTMETHOD |
89 |
|
|
|
90 |
|
|
if (vgid[nbr->GetIdx()-1] == *k) |
91 |
|
|
#else |
92 |
|
|
|
93 |
|
|
if (fabs(gp[nbr->GetIdx()]-*k) < 0.001) |
94 |
|
|
#endif |
95 |
|
|
|
96 |
|
|
ischiral = false; |
97 |
|
|
|
98 |
|
|
#ifdef INTMETHOD |
99 |
|
|
|
100 |
|
|
if (ischiral) |
101 |
|
|
tlist.push_back(vgid[nbr->GetIdx()-1]); |
102 |
|
|
#else |
103 |
|
|
|
104 |
|
|
if (ischiral) |
105 |
|
|
tlist.push_back(gp[nbr->GetIdx()]); |
106 |
|
|
#endif |
107 |
|
|
|
108 |
|
|
else |
109 |
|
|
break; |
110 |
|
|
} |
111 |
|
|
|
112 |
|
|
if (ischiral) |
113 |
|
|
atom->SetChiral(); |
114 |
|
|
} |
115 |
|
|
} |
116 |
|
|
|
117 |
|
|
// Seems to make a vector chirality become filled with array of +/- 1 for chiral atoms. |
118 |
|
|
void GetChirality(OBMol &mol, std::vector<int> &chirality) |
119 |
|
|
{ |
120 |
|
|
chirality.resize(mol.NumAtoms()+1); |
121 |
|
|
fill(chirality.begin(),chirality.end(),0); |
122 |
|
|
|
123 |
|
|
OBAtom *atom; |
124 |
|
|
vector<OBNodeBase*>::iterator i; |
125 |
|
|
for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) |
126 |
|
|
if (atom->IsChiral()) |
127 |
|
|
{ |
128 |
|
|
if (!atom->HasChiralVolume()) |
129 |
|
|
{ |
130 |
|
|
double sv = CalcSignedVolume(mol,atom); |
131 |
|
|
if (sv < 0.0) |
132 |
|
|
{ |
133 |
|
|
chirality[atom->GetIdx()-1] = -1; |
134 |
|
|
atom->SetNegativeStereo(); |
135 |
|
|
} |
136 |
|
|
else if (sv > 0.0) |
137 |
|
|
{ |
138 |
|
|
chirality[atom->GetIdx()-1] = 1; |
139 |
|
|
atom->SetPositiveStereo(); |
140 |
|
|
} |
141 |
|
|
} |
142 |
|
|
else // already calculated signed volume (e.g., imported from somewhere) |
143 |
|
|
{ |
144 |
|
|
if (atom ->IsPositiveStereo()) |
145 |
|
|
chirality[atom->GetIdx()-1] = 1; |
146 |
|
|
else |
147 |
|
|
chirality[atom->GetIdx()-1] = -1; |
148 |
|
|
} |
149 |
|
|
} |
150 |
|
|
} |
151 |
|
|
|
152 |
|
|
int GetParity4Ref(vector<unsigned int> pref) |
153 |
|
|
{ |
154 |
|
|
if(pref.size()!=4)return(-1); // should be given a vector of size 4. |
155 |
|
|
int parity=0; |
156 |
|
|
for (int i=0;i<3;i++) // do the bubble sort this many times |
157 |
|
|
{ |
158 |
|
|
for(int j=0;j<3;j++) // iterate across the array 4th element has no |
159 |
|
|
{ // right hand neighbour so no need to sort |
160 |
|
|
if (pref[j+1] < pref[j]) // compare the two neighbors |
161 |
|
|
{ |
162 |
|
|
unsigned int tmp = pref[j]; // swap a[j] and a[j+1] |
163 |
|
|
pref[j] = pref[j+1]; |
164 |
|
|
pref[j+1] = tmp; |
165 |
|
|
parity++; // parity odd will invert stereochem |
166 |
|
|
} |
167 |
|
|
} |
168 |
|
|
} // End Bubble Sort |
169 |
|
|
return(parity%2); |
170 |
|
|
} |
171 |
|
|
|
172 |
|
|
bool CorrectChirality(OBMol &mol, OBAtom *atm, atomreftype i, atomreftype o) |
173 |
|
|
{ |
174 |
|
|
if (!atm->HasChiralitySpecified()) // if no chirality defined can't do any more for 0D |
175 |
|
|
return(false); |
176 |
|
|
|
177 |
|
|
int parityI=0,parityO=0; |
178 |
|
|
OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); |
179 |
|
|
if ((cd->GetAtom4Refs(input)).size()!=4)return(false); // must have 4 refs |
180 |
|
|
parityI=GetParity4Ref(cd->GetAtom4Refs(i)); // Gets Atom4Refs used to define the chirality |
181 |
|
|
parityO=GetParity4Ref(cd->GetAtom4Refs(o));//GetsOutput parity. |
182 |
|
|
/* switch (CHTYPE) |
183 |
|
|
{ |
184 |
|
|
case SMILES: // SMILES always uses 1234 atom refs |
185 |
|
|
parityO=0; // if parityO==parityI then clockwise in input = clockwise in output |
186 |
|
|
break; |
187 |
|
|
case MOLV3000: // MOLV3000 uses 1234 unless an H then 123H |
188 |
|
|
if (atm->GetHvyValence()==3) |
189 |
|
|
{ |
190 |
|
|
OBAtom *nbr; |
191 |
|
|
int Hid=1000;// max Atom ID +1 should be used here |
192 |
|
|
vector<int> nbr_atms; |
193 |
|
|
vector<OBEdgeBase*>::iterator i; |
194 |
|
|
for (nbr = atm->BeginNbrAtom(i);nbr;nbr = atm->NextNbrAtom(i)) |
195 |
|
|
{ |
196 |
|
|
if (nbr->IsHydrogen()){Hid=nbr->GetIdx();continue;} |
197 |
|
|
nbr_atms.push_back(nbr->GetIdx()); |
198 |
|
|
} |
199 |
|
|
sort(nbr_atms.begin(),nbr_atms.end()); |
200 |
|
|
nbr_atms.push_back(Hid); |
201 |
|
|
int tmp[4]; |
202 |
|
|
for(int i=0;i<4;i++){tmp[i]=nbr_atms[i];} |
203 |
|
|
parityO=GetParity4Ref(tmp); |
204 |
|
|
} |
205 |
|
|
else if (atm->GetHvyValence()==4) |
206 |
|
|
parityO=0; |
207 |
|
|
break; |
208 |
|
|
default: |
209 |
|
|
parityO=0; |
210 |
|
|
}*/ |
211 |
|
|
if (parityO==parityI) |
212 |
|
|
{//cout << "Parity is the same"<<endl; |
213 |
|
|
return(true); |
214 |
|
|
} |
215 |
|
|
else if(parityO!=parityI) // Need to invert the Chirality which has been set |
216 |
|
|
{ //cout << "Parity is Opposite"<<endl; |
217 |
|
|
if (atm->IsClockwise()) |
218 |
|
|
{atm->UnsetStereo();atm->SetAntiClockwiseStereo();} |
219 |
|
|
else if (atm->IsAntiClockwise()) |
220 |
|
|
{atm->UnsetStereo();atm->SetClockwiseStereo();} |
221 |
|
|
else |
222 |
|
|
return(false); |
223 |
|
|
return(true); |
224 |
|
|
} |
225 |
|
|
return false; |
226 |
|
|
} |
227 |
|
|
|
228 |
|
|
//! Calculate the signed volume for an atom. If the atom has a valence of 3 |
229 |
|
|
//! the coordinates of an attached hydrogen are calculated |
230 |
|
|
//! Puts attached Hydrogen last at the moment, like mol V3000 format. |
231 |
|
|
double CalcSignedVolume(OBMol &mol,OBAtom *atm) |
232 |
|
|
{ |
233 |
|
|
vector3 tmp_crd; |
234 |
|
|
vector<unsigned int> nbr_atms; |
235 |
|
|
vector<vector3> nbr_crds; |
236 |
|
|
bool use_central_atom = false,is2D=false; |
237 |
|
|
double hbrad = etab.CorrectedBondRad(1,0); |
238 |
|
|
|
239 |
|
|
if (!mol.Has3D()) //give peudo Z coords if mol is 2D |
240 |
|
|
{ |
241 |
|
|
vector3 v,vz(0.0,0.0,1.0); |
242 |
|
|
is2D = true; |
243 |
|
|
OBAtom *nbr; |
244 |
|
|
OBBond *bond; |
245 |
|
|
vector<OBEdgeBase*>::iterator i; |
246 |
|
|
for (bond = atm->BeginBond(i);bond;bond = atm->NextBond(i)) |
247 |
|
|
{ |
248 |
|
|
nbr = bond->GetEndAtom(); |
249 |
|
|
if (nbr != atm) |
250 |
|
|
{ |
251 |
|
|
v = nbr->GetVector(); |
252 |
|
|
if (bond->IsWedge()) |
253 |
|
|
v += vz; |
254 |
|
|
else |
255 |
|
|
if (bond->IsHash()) |
256 |
|
|
v -= vz; |
257 |
|
|
|
258 |
|
|
nbr->SetVector(v); |
259 |
|
|
} |
260 |
|
|
else |
261 |
|
|
{ |
262 |
|
|
nbr = bond->GetBeginAtom(); |
263 |
|
|
v = nbr->GetVector(); |
264 |
|
|
if (bond->IsWedge()) |
265 |
|
|
v -= vz; |
266 |
|
|
else |
267 |
|
|
if (bond->IsHash()) |
268 |
|
|
v += vz; |
269 |
|
|
|
270 |
|
|
nbr->SetVector(v); |
271 |
|
|
} |
272 |
|
|
} |
273 |
|
|
} |
274 |
|
|
|
275 |
|
|
if (atm->GetHvyValence() < 3) |
276 |
|
|
{ |
277 |
|
|
#ifdef HAVE_SSTREAM |
278 |
|
|
stringstream errorMsg; |
279 |
|
|
#else |
280 |
|
|
strstream errorMsg; |
281 |
|
|
#endif |
282 |
|
|
errorMsg << "Cannot calculate a signed volume for an atom with a heavy atom valence of " << atm->GetHvyValence() << endl; |
283 |
tim |
819 |
obErrorLog.ThrowError(__func__, errorMsg.str(), obInfo); |
284 |
tim |
741 |
return(0.0); |
285 |
|
|
} |
286 |
|
|
|
287 |
|
|
// Create a vector with the coordinates of the neighbor atoms |
288 |
|
|
// Also make a vector with Atom IDs |
289 |
|
|
OBAtom *nbr; |
290 |
|
|
vector<OBEdgeBase*>::iterator bint; |
291 |
|
|
for (nbr = atm->BeginNbrAtom(bint);nbr;nbr = atm->NextNbrAtom(bint)) |
292 |
|
|
{ |
293 |
|
|
nbr_atms.push_back(nbr->GetIdx()); |
294 |
|
|
} |
295 |
|
|
// sort the neighbor atoms to insure a consistent ordering |
296 |
|
|
sort(nbr_atms.begin(),nbr_atms.end()); |
297 |
|
|
for (unsigned int i = 0; i < nbr_atms.size(); i++) |
298 |
|
|
{ |
299 |
|
|
OBAtom *tmp_atm = mol.GetAtom(nbr_atms[i]); |
300 |
|
|
nbr_crds.push_back(tmp_atm->GetVector()); |
301 |
|
|
} |
302 |
|
|
/* |
303 |
|
|
// If we have three heavy atoms we need to calculate the position of the fourth |
304 |
|
|
if (atm->GetHvyValence() == 3) |
305 |
|
|
{ |
306 |
|
|
double bondlen = hbrad+etab.CorrectedBondRad(atm->GetAtomicNum(),atm->GetHyb()); |
307 |
|
|
atm->GetNewBondVector(tmp_crd,bondlen); |
308 |
|
|
nbr_crds.push_back(tmp_crd); |
309 |
|
|
} |
310 |
|
|
*/ |
311 |
|
|
for(int j=0;j < nbr_crds.size();j++) // Checks for a neighbour having 0 co-ords (added hydrogen etc) |
312 |
|
|
{ |
313 |
|
|
if (nbr_crds[j]==0 && use_central_atom==false)use_central_atom=true; |
314 |
|
|
else if (nbr_crds[j]==0) |
315 |
|
|
{ |
316 |
tim |
819 |
obErrorLog.ThrowError(__func__, "More than 2 neighbours have 0 co-ords when attempting 3D chiral calculation", obInfo); |
317 |
tim |
741 |
} |
318 |
|
|
} |
319 |
|
|
|
320 |
|
|
// If we have three heavy atoms we can use the chiral center atom itself for the fourth |
321 |
|
|
// will always give same sign (for tetrahedron), magnitude will be smaller. |
322 |
|
|
if(nbr_atms.size()==3 || use_central_atom==true) |
323 |
|
|
{ |
324 |
|
|
nbr_crds.push_back(atm->GetVector()); |
325 |
|
|
nbr_atms.push_back(mol.NumAtoms()+1); // meed to add largest number on end to work |
326 |
|
|
} |
327 |
|
|
OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); //Set the output atom4refs to the ones used |
328 |
|
|
if(cd==NULL) |
329 |
|
|
{ |
330 |
|
|
cd = new OBChiralData; |
331 |
|
|
atm->SetData(cd); |
332 |
|
|
} |
333 |
|
|
cd->SetAtom4Refs(nbr_atms,calcvolume); |
334 |
|
|
|
335 |
|
|
//re-zero psuedo-coords |
336 |
|
|
if (is2D) |
337 |
|
|
{ |
338 |
|
|
vector3 v; |
339 |
|
|
OBAtom *atom; |
340 |
|
|
vector<OBNodeBase*>::iterator k; |
341 |
|
|
for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k)) |
342 |
|
|
{ |
343 |
|
|
v = atom->GetVector(); |
344 |
|
|
v.SetZ(0.0); |
345 |
|
|
atom->SetVector(v); |
346 |
|
|
} |
347 |
|
|
} |
348 |
|
|
|
349 |
|
|
return(signed_volume(nbr_crds[0],nbr_crds[1],nbr_crds[2],nbr_crds[3])); |
350 |
|
|
} |
351 |
|
|
|
352 |
|
|
//! Calculate a signed volume given a set of 4 coordinates |
353 |
|
|
double signed_volume(const vector3 &a, const vector3 &b, const vector3 &c, const vector3 &d) |
354 |
|
|
{ |
355 |
|
|
vector3 A,B,C; |
356 |
|
|
A = b-a; |
357 |
|
|
B = c-a; |
358 |
|
|
C = d-a; |
359 |
|
|
matrix3x3 m(A,B,C); |
360 |
|
|
return(m.determinant()); |
361 |
|
|
} |
362 |
|
|
|
363 |
|
|
//! \brief Calculate the Graph Potentials of a molecule |
364 |
|
|
//! |
365 |
|
|
//! based on |
366 |
|
|
//! V.E. and Rozenblit, A.B. Golender |
367 |
|
|
//! <em>Logical and Combinatorial Algorithms for Drug Design</em>. \n |
368 |
|
|
//! For an example see: |
369 |
|
|
//! Walters, W. P., Yalkowsky, S. H., \em JCICS, 1996, 36(5), 1015-1017. |
370 |
|
|
//! <a href="http://dx.doi.org/10.1021/ci950278o">DOI: 10.1021/ci950278o</a> |
371 |
|
|
void GraphPotentials(OBMol &mol, std::vector<double> &pot) |
372 |
|
|
{ |
373 |
|
|
double det; |
374 |
|
|
|
375 |
|
|
vector<vector<double> > g,c,h; |
376 |
|
|
construct_g_matrix(mol,g); |
377 |
|
|
invert_matrix(g,det); |
378 |
|
|
construct_c_matrix(mol,c); |
379 |
|
|
mult_matrix(h,g,c); |
380 |
|
|
pot.resize(mol.NumAtoms()+1); |
381 |
|
|
|
382 |
|
|
for (unsigned int i = 0; i < mol.NumAtoms();i++) |
383 |
|
|
pot[i+1] = h[i][0]; |
384 |
|
|
} |
385 |
|
|
|
386 |
|
|
|
387 |
|
|
//! Construct the matrix G, which puts each atoms valence+1 |
388 |
|
|
//! on the diagonal and and -1 on the off diagonal if two |
389 |
|
|
//! atoms are connected. |
390 |
|
|
void construct_g_matrix(OBMol &mol, std::vector<std::vector<double> > &m) |
391 |
|
|
{ |
392 |
|
|
unsigned int i,j; |
393 |
|
|
|
394 |
|
|
OBAtom *atm1,*atm2; |
395 |
|
|
vector<OBNodeBase*>::iterator aint,bint; |
396 |
|
|
|
397 |
|
|
m.resize(mol.NumAtoms()); |
398 |
|
|
for (i = 0; i < m.size(); i++) |
399 |
|
|
m[i].resize(mol.NumAtoms()); |
400 |
|
|
|
401 |
|
|
for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++) |
402 |
|
|
for (atm2 = mol.BeginAtom(bint),j=0;atm2;atm2 = mol.NextAtom(bint),j++) |
403 |
|
|
{ |
404 |
|
|
if (i == j) |
405 |
|
|
{ |
406 |
|
|
m[i][j] = atm1->GetValence() + 1; |
407 |
|
|
m[i][j] += (double)atm1->GetAtomicNum()/10.0; |
408 |
|
|
m[i][j] += (double)atm1->GetHyb()/100.0; |
409 |
|
|
} |
410 |
|
|
else |
411 |
|
|
{ |
412 |
|
|
if (atm1->IsConnected(atm2)) |
413 |
|
|
m[i][j] = -1; |
414 |
|
|
else |
415 |
|
|
m[i][j] = 0; |
416 |
|
|
} |
417 |
|
|
} |
418 |
|
|
} |
419 |
|
|
|
420 |
|
|
//! Construct the matrix C, which is simply a column vector |
421 |
|
|
//! consisting of the valence for each atom |
422 |
|
|
void construct_c_matrix(OBMol &mol,std::vector<std::vector<double > > &m) |
423 |
|
|
{ |
424 |
|
|
unsigned int i; |
425 |
|
|
OBAtom *atm1; |
426 |
|
|
vector<OBNodeBase*>::iterator aint; |
427 |
|
|
|
428 |
|
|
m.resize(mol.NumAtoms()); |
429 |
|
|
for (i = 0; i < m.size(); i++) |
430 |
|
|
m[i].resize(1); |
431 |
|
|
for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++) |
432 |
|
|
{ |
433 |
|
|
m[i][0] = atm1->GetValence(); |
434 |
|
|
} |
435 |
|
|
} |
436 |
|
|
|
437 |
|
|
} // namespace OpenBabel |
438 |
|
|
|
439 |
|
|
//! \file chiral.cpp |
440 |
|
|
//! \brief Detect chiral atoms and molecules. |