| 1 |
/********************************************************************* |
| 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 |
obErrorLog.ThrowError(__func__, |
| 40 |
"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 |
obErrorLog.ThrowError(__func__, errorMsg.str(), obInfo); |
| 284 |
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 |
obErrorLog.ThrowError(__func__, "More than 2 neighbours have 0 co-ords when attempting 3D chiral calculation", obInfo); |
| 317 |
} |
| 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. |
