1 /********************************************************************** 2 typer.cpp - Open Babel atom typer. 3 4 Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. 5 Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison 6 7 This file is part of the Open Babel project. 8 For more information, see <http://openbabel.org/> 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 #include <openbabel/babelconfig.h> 20 21 #include <openbabel/mol.h> 22 #include <openbabel/atom.h> 23 #include <openbabel/bond.h> 24 #include <openbabel/ring.h> 25 #include <openbabel/obiter.h> 26 #include <openbabel/oberror.h> 27 #include <openbabel/typer.h> 28 #include <openbabel/elements.h> 29 30 // private data headers with default parameters 31 #include "atomtyp.h" 32 33 #ifdef WIN32 34 #pragma warning (disable : 4786) 35 #endif 36 37 using namespace std; 38 39 namespace OpenBabel 40 { 41 // Initialize globals declared in typer.h 42 THREAD_LOCAL OBAromaticTyper aromtyper; 43 THREAD_LOCAL OBAtomTyper atomtyper; 44 45 /*! \class OBAtomTyper typer.h <openbabel/typer.h> 46 \brief Assigns atom types, hybridization, and formal charges 47 48 The OBAtomTyper class is designed to read in a list of atom typing 49 rules and apply them to molecules. The code that performs atom 50 typing is not usually used directly as atom typing, hybridization 51 assignment, and charge are all done 52 automatically when their corresponding values are requested of 53 atoms: 54 \code 55 atom->GetType(); 56 atom->GetFormalCharge(); 57 atom->GetHyb(); 58 \endcode 59 */ OBAtomTyper()60 OBAtomTyper::OBAtomTyper() 61 { 62 _init = false; 63 _dir = BABEL_DATADIR; 64 _envvar = "BABEL_DATADIR"; 65 _filename = "atomtyp.txt"; 66 _subdir = "data"; 67 _dataptr = AtomTypeData; 68 } 69 ParseLine(const char * buffer)70 void OBAtomTyper::ParseLine(const char *buffer) 71 { 72 vector<string> vs; 73 OBSmartsPattern *sp; 74 75 if (EQn(buffer,"INTHYB",6)) 76 { 77 tokenize(vs,buffer); 78 if (vs.size() < 3) 79 { 80 obErrorLog.ThrowError(__FUNCTION__, " Could not parse INTHYB line in atom type table from atomtyp.txt", obInfo); 81 return; 82 } 83 84 sp = new OBSmartsPattern; 85 if (sp->Init(vs[1])) 86 _vinthyb.push_back(pair<OBSmartsPattern*,int> (sp,atoi((char*)vs[2].c_str()))); 87 else 88 { 89 delete sp; 90 sp = nullptr; 91 obErrorLog.ThrowError(__FUNCTION__, " Could not parse INTHYB line in atom type table from atomtyp.txt", obInfo); 92 return; 93 } 94 } 95 else if (EQn(buffer,"EXTTYP",6)) 96 { 97 tokenize(vs,buffer); 98 if (vs.size() < 3) 99 { 100 obErrorLog.ThrowError(__FUNCTION__, " Could not parse EXTTYP line in atom type table from atomtyp.txt", obInfo); 101 return; 102 } 103 sp = new OBSmartsPattern; 104 if (sp->Init(vs[1])) 105 _vexttyp.push_back(pair<OBSmartsPattern*,string> (sp,vs[2])); 106 else 107 { 108 delete sp; 109 sp = nullptr; 110 obErrorLog.ThrowError(__FUNCTION__, " Could not parse EXTTYP line in atom type table from atomtyp.txt", obInfo); 111 return; 112 } 113 } 114 } 115 ~OBAtomTyper()116 OBAtomTyper::~OBAtomTyper() 117 { 118 vector<pair<OBSmartsPattern*,int> >::iterator i; 119 for (i = _vinthyb.begin();i != _vinthyb.end();++i) 120 { 121 delete i->first; 122 i->first = nullptr; 123 } 124 125 vector<pair<OBSmartsPattern*,string> >::iterator j; 126 for (j = _vexttyp.begin();j != _vexttyp.end();++j) 127 { 128 delete j->first; 129 j->first = nullptr; 130 } 131 132 } 133 AssignTypes(OBMol & mol)134 void OBAtomTyper::AssignTypes(OBMol &mol) 135 { 136 if (!_init) 137 Init(); 138 139 obErrorLog.ThrowError(__FUNCTION__, 140 "Ran OpenBabel::AssignTypes", obAuditMsg); 141 142 mol.SetAtomTypesPerceived(); 143 144 vector<vector<int> >::iterator j; 145 vector<pair<OBSmartsPattern*,string> >::iterator i; 146 147 for (i = _vexttyp.begin(); i != _vexttyp.end(); ++i) { 148 std::vector<std::vector<int> > mlist; 149 if (i->first->Match(mol, mlist)) 150 { 151 for (j = mlist.begin(); j != mlist.end(); ++j) 152 mol.GetAtom((*j)[0])->SetType(i->second); 153 } 154 } 155 156 // Special cases 157 vector<OBAtom*>::iterator a; 158 OBAtom* atom; 159 for (atom = mol.BeginAtom(a); atom; atom = mol.NextAtom(a)) { 160 // guanidinium. Fixes PR#1800964 161 if (strncasecmp(atom->GetType(),"C2", 2) == 0) { 162 int guanidineN = 0; 163 OBAtom *nbr; 164 vector<OBBond*>::iterator k; 165 for (nbr = atom->BeginNbrAtom(k);nbr;nbr = atom->NextNbrAtom(k)) { 166 if (strncasecmp(nbr->GetType(),"Npl", 3) == 0 || 167 strncasecmp(nbr->GetType(),"N2", 2) == 0 || 168 strncasecmp(nbr->GetType(),"Ng+", 3) == 0) 169 ++guanidineN; 170 } 171 if (guanidineN == 3) 172 atom->SetType("C+"); 173 174 } // end C2 carbon for guanidinium 175 176 } // end special cases 177 } 178 AssignHyb(OBMol & mol)179 void OBAtomTyper::AssignHyb(OBMol &mol) 180 { 181 if (!_init) 182 Init(); 183 184 aromtyper.AssignAromaticFlags(mol); 185 186 mol.SetHybridizationPerceived(); 187 obErrorLog.ThrowError(__FUNCTION__, 188 "Ran OpenBabel::AssignHybridization", obAuditMsg); 189 190 OBAtom *atom; 191 vector<OBAtom*>::iterator k; 192 for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k)) 193 atom->SetHyb(0); 194 195 vector<vector<int> >::iterator j; 196 vector<pair<OBSmartsPattern*,int> >::iterator i; 197 198 for (i = _vinthyb.begin(); i != _vinthyb.end(); ++i) { 199 std::vector<std::vector<int> > mlist; 200 if (i->first->Match(mol, mlist)) 201 { 202 for (j = mlist.begin(); j != mlist.end(); ++j) 203 mol.GetAtom((*j)[0])->SetHyb(i->second); 204 } 205 } 206 207 // check all atoms to make sure *some* hybridization is assigned 208 for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k)) 209 if (atom->GetHyb() == 0) { 210 switch (atom->GetExplicitDegree()) { 211 case 0: 212 case 1: 213 case 2: 214 atom->SetHyb(1); 215 break; 216 case 3: 217 atom->SetHyb(2); 218 break; 219 case 4: 220 atom->SetHyb(3); 221 break; 222 default: 223 atom->SetHyb(atom->GetExplicitDegree()); 224 } 225 } 226 } 227 228 /*! \class OBRingTyper typer.h <openbabel/typer.h> 229 \brief Assigns ring types 230 231 The OBRingTyper class is designed to read in a list of ring typing 232 rules and apply them to molecules. The code that performs ring 233 typing is not usually used directly as ring typing is done 234 automatically when the ring type is requested of rings: 235 \code 236 vector<OBRing*>::iterator i; 237 vector<OBRing*> rlist = mol.GetSSSR(); 238 239 for (i = rlist.begin();i != rlist.end();++i) 240 cout << "ring type = " << (*i)->GetType() << endl; 241 \endcode 242 */ OBRingTyper()243 OBRingTyper::OBRingTyper() 244 { 245 _init = false; 246 _dir = BABEL_DATADIR; 247 _envvar = "BABEL_DATADIR"; 248 _filename = "ringtyp.txt"; 249 _subdir = "data"; 250 //_dataptr = RingTypeData; 251 } 252 ParseLine(const char * buffer)253 void OBRingTyper::ParseLine(const char *buffer) 254 { 255 vector<string> vs; 256 OBSmartsPattern *sp; 257 258 if (EQn(buffer,"RINGTYP",7)) { 259 tokenize(vs,buffer); 260 if (vs.size() < 3) { 261 obErrorLog.ThrowError(__FUNCTION__, " Could not parse RING line in ring type table from ringtyp.txt", obInfo); 262 return; 263 } 264 sp = new OBSmartsPattern; 265 if (sp->Init(vs[2])) 266 _ringtyp.push_back(pair<OBSmartsPattern*,string> (sp,vs[1])); 267 else { 268 delete sp; 269 sp = nullptr; 270 obErrorLog.ThrowError(__FUNCTION__, " Could not parse RING line in ring type table from ringtyp.txt", obInfo); 271 return; 272 } 273 } 274 } 275 ~OBRingTyper()276 OBRingTyper::~OBRingTyper() 277 { 278 vector<pair<OBSmartsPattern*,string> >::iterator i; 279 for (i = _ringtyp.begin();i != _ringtyp.end();++i) { 280 delete i->first; 281 i->first = nullptr; 282 } 283 } 284 AssignTypes(OBMol & mol)285 void OBRingTyper::AssignTypes(OBMol &mol) 286 { 287 if (!_init) 288 Init(); 289 290 obErrorLog.ThrowError(__FUNCTION__, 291 "Ran OBRing::AssignTypes", obAuditMsg); 292 293 mol.SetRingTypesPerceived(); 294 295 vector<vector<int> >::iterator j2; 296 vector<pair<OBSmartsPattern*,string> >::iterator i2; 297 298 vector<OBRing*>::iterator i; 299 vector<int>::iterator j; 300 vector<OBRing*> rlist = mol.GetSSSR(); 301 302 unsigned int member_count; 303 for (i2 = _ringtyp.begin();i2 != _ringtyp.end();++i2) { // for each ring type 304 std::vector<std::vector<int> > mlist; 305 if (i2->first->Match(mol, mlist)) { 306 for (j2 = mlist.begin();j2 != mlist.end();++j2) { // for each found match 307 308 for (i = rlist.begin();i != rlist.end();++i) { // for each ring 309 member_count = 0; 310 311 for(j = j2->begin(); j != j2->end(); ++j) { // for each atom in the match 312 if ((*i)->IsMember(mol.GetAtom(*j))) 313 member_count++; 314 } 315 316 if ((*i)->Size() == member_count) 317 (*i)->SetType(i2->second); 318 } 319 } 320 } 321 } 322 } 323 324 // Start of helper functions for AssignOBAromaticityModel 325 enum ExocyclicAtom { NO_EXOCYCLIC_ATOM, EXO_OXYGEN, EXO_NONOXYGEN }; 326 FindExocyclicAtom(OBAtom * atm)327 static ExocyclicAtom FindExocyclicAtom(OBAtom *atm) 328 { 329 FOR_BONDS_OF_ATOM(bond, atm) { 330 if (bond->GetBondOrder() == 2 && !bond->IsInRing()) { 331 unsigned int atomicnum = bond->GetNbrAtom(atm)->GetAtomicNum(); 332 switch (atomicnum) { 333 case OBElements::Oxygen: 334 return EXO_OXYGEN; 335 default: 336 return EXO_NONOXYGEN; 337 } 338 } 339 } 340 return NO_EXOCYCLIC_ATOM; 341 } 342 HasExocyclicBondToOxygenMinus(OBAtom * atm)343 static bool HasExocyclicBondToOxygenMinus(OBAtom *atm) 344 { 345 FOR_BONDS_OF_ATOM(bond, atm) { 346 if (bond->GetBondOrder() == 1 && !bond->IsInRing()) { 347 OBAtom* nbr = bond->GetNbrAtom(atm); 348 if (nbr->GetAtomicNum() == OBElements::Oxygen && nbr->GetFormalCharge() == -1) 349 return true; 350 } 351 } 352 return false; 353 } 354 HasExocyclicDblBondToOxygen(OBAtom * atm)355 static bool HasExocyclicDblBondToOxygen(OBAtom *atm) 356 { 357 FOR_BONDS_OF_ATOM(bond, atm) { 358 if (bond->GetBondOrder() == 2 && !bond->IsInRing() && 359 bond->GetNbrAtom(atm)->GetAtomicNum() == OBElements::Oxygen) 360 return true; 361 } 362 return false; 363 } 364 HasExocyclicDblBondToHet(OBAtom * atm)365 static bool HasExocyclicDblBondToHet(OBAtom *atm) 366 { 367 FOR_BONDS_OF_ATOM(bond, atm) { 368 if (bond->GetBondOrder() == 2 && !bond->IsInRing()) { 369 unsigned int atomicnum = bond->GetNbrAtom(atm)->GetAtomicNum(); 370 switch (atomicnum) { 371 case OBElements::Carbon: 372 case OBElements::Hydrogen: 373 break; 374 default: 375 return true; 376 } 377 } 378 } 379 return false; 380 } 381 // End of predicates for AssignOBAromaticityModel 382 AssignOBAromaticityModel(OBAtom * atm,int & min,int & max)383 static bool AssignOBAromaticityModel(OBAtom *atm, int &min, int &max) 384 { 385 // The Open Babel aromaticity model 386 // 387 // Return minimum and maximum pi-electrons contributed to an aromatic system 388 // The return value indicates a potentially aromatic atom (i.e. any patterns matched) 389 // 390 // The original code used SMARTS patterns organised in increasing order of 391 // prioirity (i.e. later matches overrode earlier ones). These SMARTS patterns 392 // are now implemented in the code below, but are included in the comments 393 // for reference (Case 1->22). 394 395 if (!atm->IsInRing()) { 396 min = 0; max = 0; return false; 397 } 398 399 unsigned int elem = atm->GetAtomicNum(); 400 int chg = atm->GetFormalCharge(); 401 unsigned int deg = atm->GetExplicitDegree() + atm->GetImplicitHCount(); 402 unsigned int val = atm->GetExplicitValence() + atm->GetImplicitHCount(); 403 404 switch (elem) { 405 case OBElements::Carbon: 406 switch (chg) { 407 case 0: 408 if (val == 4 && deg == 3) { 409 if (HasExocyclicDblBondToHet(atm)) { 410 min = 0; max = 0; return true; 411 } 412 else { 413 min = 1; max = 1; return true; 414 } 415 } 416 break; 417 case 1: 418 if (val == 3) { 419 switch (deg) { 420 case 3: 421 min = 0; max = 0; return true; 422 case 2: 423 min = 1; max = 1; return true; 424 } 425 } 426 break; 427 case -1: 428 if (val == 3) { 429 switch (deg) { 430 case 3: 431 min = 2; max = 2; return true; 432 case 2: 433 min = 1; max = 1; return true; 434 } 435 } 436 break; 437 } 438 break; 439 440 case OBElements::Nitrogen: // nitrogen 441 case OBElements::Phosphorus: // phosphorus 442 switch (chg) { 443 case 0: 444 switch (val) { 445 case 3: 446 switch (deg) { 447 case 3: 448 min = 2; max = 2; return true; 449 case 2: 450 min = 1; max = 1; return true; 451 } 452 break; 453 case 5: 454 if (deg == 3) { 455 ExocyclicAtom exoatom = FindExocyclicAtom(atm); 456 switch (exoatom) { 457 case EXO_OXYGEN: 458 min = 1; max = 1; return true; 459 case EXO_NONOXYGEN: 460 min = 2; max = 2; return true; 461 default: 462 ; 463 } 464 } 465 } 466 break; 467 case 1: 468 if (val == 4 && deg == 3) { 469 min = 1; max = 1; return true; 470 } 471 break; 472 case -1: 473 if (val == 2 && deg == 2) { 474 min = 2; max = 2; return true; 475 } 476 } 477 break; 478 479 case OBElements::Oxygen: 480 case OBElements::Selenium: 481 switch (chg) { 482 case 0: 483 if (val == 2 && deg == 2) { 484 min = 2; max = 2; return true; 485 } 486 case 1: 487 if (val == 3 && deg == 2) { 488 min = 1; max = 1; return true; 489 } 490 } 491 break; 492 493 case OBElements::Sulfur: 494 switch (chg) { 495 case 0: 496 switch (val) { 497 case 2: 498 if (deg == 2) { 499 min = 2; max = 2; return true; 500 } 501 break; 502 case 4: 503 if (deg == 3 && HasExocyclicDblBondToOxygen(atm)) { 504 min = 2; max = 2; return true; 505 } 506 } 507 break; 508 case 1: 509 if (val == 3) { 510 switch (deg) { 511 case 2: 512 min = 1; max = 1; return true; 513 case 3: 514 if (HasExocyclicBondToOxygenMinus(atm)) { 515 min = 2; max = 2; return true; 516 } 517 } 518 } 519 } 520 break; 521 522 case OBElements::Boron: 523 if (chg == 0 && val == 3) { 524 switch (deg) { 525 case 2: 526 min = 1; max = 1; return true; 527 case 3: 528 min = 0; max = 0; return true; 529 } 530 } 531 break; 532 533 case OBElements::Arsenic: 534 switch (chg) { 535 case 0: 536 if (val == 3) { 537 switch (deg) { 538 case 2: 539 min = 1; max = 1; return true; 540 case 3: 541 min = 2; max = 2; return true; 542 } 543 } 544 break; 545 case 1: 546 if (val == 4 && deg == 3) { 547 min = 1; max = 1; return true; 548 } 549 } 550 break; 551 552 553 case 0: // Asterisk 554 if (chg == 0) { 555 switch (val) { 556 case 2: 557 switch (deg) { 558 case 2: 559 case 3: 560 min = 0; max = 2; return true; 561 } 562 break; 563 case 3: 564 switch (deg) { 565 case 2: 566 case 3: 567 min = 0; max = 1; return true; 568 } 569 break; 570 } 571 } 572 break; 573 } 574 575 min = 0; max = 0; // nothing matched 576 return false; 577 } 578 579 class OBAromaticTyperMolState 580 { 581 public: OBAromaticTyperMolState(OBMol & mol)582 OBAromaticTyperMolState(OBMol &mol) : mol(mol) 583 { 584 _vpa.resize(mol.NumAtoms() + 1); 585 _velec.resize(mol.NumAtoms() + 1); 586 _root.resize(mol.NumAtoms() + 1); 587 _visit.resize(mol.NumAtoms() + 1); 588 } 589 void AssignAromaticFlags(); 590 private: 591 OBMol &mol; 592 std::vector<bool> _vpa; //!< potentially aromatic atoms 593 std::vector<bool> _visit; 594 std::vector<bool> _root; 595 std::vector<std::pair<int, int> > _velec; //!< # electrons an atom contributes 596 597 //! "Anti-alias" potentially aromatic flags around a molecule 598 //! (aromatic atoms need to have >= 2 neighboring ring atoms) 599 void PropagatePotentialAromatic(OBAtom*); 600 // Documentation in typer.cpp 601 void SelectRootAtoms(bool avoidInnerRingAtoms = true); 602 //! Remove 3-member rings from consideration 603 void ExcludeSmallRing(); 604 //! Check aromaticity starting from the root atom, up to a specified depth 605 void CheckAromaticity(OBAtom *root, int searchDepth); 606 // Documentation in typer.cpp 607 bool TraverseCycle(OBAtom *root, OBAtom *atom, OBBond *prev, 608 std::pair<int, int> &er, int depth); 609 }; 610 AssignAromaticFlags()611 void OBAromaticTyperMolState::AssignAromaticFlags() 612 { 613 OBBond *bond; 614 OBAtom *atom; 615 vector<OBAtom*>::iterator i; 616 vector<OBBond*>::iterator j; 617 618 //unset all aromatic flags 619 for (atom = mol.BeginAtom(i); atom; atom = mol.NextAtom(i)) 620 atom->SetAromatic(false); 621 for (bond = mol.BeginBond(j); bond; bond = mol.NextBond(j)) 622 bond->SetAromatic(false); 623 624 // New code using lookups instead of SMARTS patterns 625 FOR_ATOMS_OF_MOL(atom, mol) { 626 unsigned int idx = atom->GetIdx(); 627 _vpa[idx] = AssignOBAromaticityModel(&(*atom), _velec[idx].first, _velec[idx].second); 628 } 629 630 //propagate potentially aromatic atoms 631 for (atom = mol.BeginAtom(i); atom; atom = mol.NextAtom(i)) 632 if (_vpa[atom->GetIdx()]) 633 PropagatePotentialAromatic(atom); 634 635 //select root atoms 636 SelectRootAtoms(); 637 638 ExcludeSmallRing(); //remove 3 membered rings from consideration 639 640 //loop over root atoms and look for aromatic rings 641 642 for (atom = mol.BeginAtom(i); atom; atom = mol.NextAtom(i)) 643 if (_root[atom->GetIdx()]) 644 CheckAromaticity(atom, 14); 645 646 //for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) 647 // if (atom->IsAromatic()) 648 // cerr << "aro = " <<atom->GetIdx() << endl; 649 650 //for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j)) 651 //if (bond->IsAromatic()) 652 //cerr << bond->GetIdx() << ' ' << bond->IsAromatic() << endl; 653 } 654 655 /*! \class OBAromaticTyper typer.h <openbabel/typer.h> 656 \brief Assigns aromatic typing to atoms and bonds 657 658 The OBAromaticTyper class applies a set of 659 aromatic perception rules to molecules. The code 660 that performs typing is not usually used directly -- it is usually 661 done automatically when their corresponding values are requested of atoms 662 or bonds. 663 \code 664 atom->IsAromatic(); 665 bond->IsAromatic(); 666 \endcode 667 */ 668 AssignAromaticFlags(OBMol & mol)669 void OBAromaticTyper::AssignAromaticFlags(OBMol &mol) 670 { 671 if (mol.HasAromaticPerceived()) 672 return; 673 mol.SetAromaticPerceived(); 674 obErrorLog.ThrowError(__FUNCTION__, 675 "Ran OpenBabel::AssignAromaticFlags", obAuditMsg); 676 677 OBAromaticTyperMolState molstate(mol); 678 molstate.AssignAromaticFlags(); 679 } 680 681 /** \brief Traverse a potentially aromatic cycle starting at @p root. 682 \return True if the cycle is likely aromatic 683 \param root The initial, "root" atom in traversing this ring 684 \param atom The current atom to visit and check 685 \param prev The bond traversed in moving to this @p atom 686 \param er The min and max number of pi electrons for this ring 687 \param depth The maximum number of atoms to visit in a ring (e.g., 6) 688 689 This method traverses a potentially aromatic ring, adding up the possible 690 pi electrons for each atom. At the end (e.g., when @p atom == @p root) 691 the Huekel 4n+2 rule is checked to see if there is a possible electronic 692 configuration which corresponds to aromaticity. 693 **/ TraverseCycle(OBAtom * root,OBAtom * atom,OBBond * prev,std::pair<int,int> & er,int depth)694 bool OBAromaticTyperMolState::TraverseCycle(OBAtom *root, OBAtom *atom, OBBond *prev, 695 std::pair<int,int> &er,int depth) 696 { 697 if (atom == root) 698 { 699 int i; 700 for (i = er.first;i <= er.second;++i) 701 if (i%4 == 2 && i > 2) 702 return(true); 703 704 return(false); 705 } 706 707 if (!depth || !_vpa[atom->GetIdx()] || _visit[atom->GetIdx()]) 708 return(false); 709 710 bool result = false; 711 712 depth--; 713 er.first += _velec[atom->GetIdx()].first; 714 er.second += _velec[atom->GetIdx()].second; 715 716 _visit[atom->GetIdx()] = true; 717 OBAtom *nbr; 718 vector<OBBond*>::iterator i; 719 for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i)) 720 if (*i != prev && (*i)->IsInRing() && _vpa[nbr->GetIdx()]) 721 { 722 if (TraverseCycle(root,nbr,(OBBond*)(*i),er,depth)) 723 { 724 result = true; 725 ((OBBond*) *i)->SetAromatic(); 726 } 727 } 728 729 _visit[atom->GetIdx()] = false; 730 if (result) 731 atom->SetAromatic(); 732 733 er.first -= _velec[atom->GetIdx()].first; 734 er.second -= _velec[atom->GetIdx()].second; 735 736 return(result); 737 } 738 CheckAromaticity(OBAtom * atom,int depth)739 void OBAromaticTyperMolState::CheckAromaticity(OBAtom *atom,int depth) 740 { 741 OBAtom *nbr; 742 vector<OBBond*>::iterator i; 743 744 pair<int,int> erange; 745 for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i)) 746 if ((*i)->IsInRing()) // check all rings, regardless of assumed aromaticity 747 { 748 erange = _velec[atom->GetIdx()]; 749 750 if (TraverseCycle(atom,nbr,(OBBond *)*i,erange,depth-1)) 751 { 752 atom->SetAromatic(); 753 ((OBBond*) *i)->SetAromatic(); 754 } 755 } 756 } 757 PropagatePotentialAromatic(OBAtom * atom)758 void OBAromaticTyperMolState::PropagatePotentialAromatic(OBAtom *atom) 759 { 760 int count = 0; 761 OBAtom *nbr; 762 vector<OBBond*>::iterator i; 763 764 for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i)) 765 if ((*i)->IsInRing() && _vpa[nbr->GetIdx()]) 766 count++; 767 768 if (count < 2) 769 { 770 _vpa[atom->GetIdx()] = false; 771 if (count == 1) 772 for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i)) 773 if ((*i)->IsInRing() && _vpa[nbr->GetIdx()]) 774 PropagatePotentialAromatic(nbr); 775 } 776 } 777 778 /** 779 * \brief Select the root atoms for traversing atoms in rings. 780 * 781 * Picking only the begin atom of a closure bond can cause 782 * difficulties when the selected atom is an inner atom 783 * with three neighbour ring atoms. Why ? Because this atom 784 * can get trapped by the other atoms when determining aromaticity, 785 * because a simple visited flag is used in the 786 * OBAromaticTyper::TraverseCycle() method. 787 * 788 * Ported from JOELib, copyright Joerg Wegner, 2003 under the GPL version 2 789 * Improved by Fabian (fab5) in 2009 -- PR#2889708 790 * 791 * @param mol the molecule 792 * @param avoidInnerRingAtoms inner closure ring atoms with more than 2 neighbours will be avoided 793 * 794 */ SelectRootAtoms(bool avoidInnerRingAtoms)795 void OBAromaticTyperMolState::SelectRootAtoms(bool avoidInnerRingAtoms) 796 { 797 OBBond *bond; 798 OBAtom *atom, *nbr, *nbr2; 799 OBRing *ring; 800 // vector<OBAtom*>::iterator i; 801 vector<OBBond*>::iterator j, l, nbr2Iter; 802 vector<OBRing*> sssRings = mol.GetSSSR(); 803 vector<OBRing*>::iterator k; 804 805 int rootAtom; 806 int ringNbrs; 807 int heavyNbrs; 808 int newRoot = -1; 809 vector<int> tmpRootAtoms; 810 vector<int> tmp; 811 812 vector<OBBond*> cbonds; 813 vector< vector<OBRing*> > ringAtoms; // store ring pointers on an atom basis 814 815 //generate list of closure bonds 816 for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j)) 817 { 818 if( bond->IsClosure() ) 819 { 820 cbonds.push_back(bond); 821 if(avoidInnerRingAtoms) 822 { 823 tmpRootAtoms.push_back(bond->GetBeginAtomIdx()); 824 } 825 } 826 } 827 828 if(avoidInnerRingAtoms) 829 { 830 //for every atom fill vector with ring pointer it's associated with 831 ringAtoms.resize(mol.NumAtoms()+1); 832 for (k = sssRings.begin();k != sssRings.end();++k) 833 { 834 tmp = (*k)->_path; 835 for (unsigned int j (0),j_end(tmp.size()); j < j_end; ++j) 836 { 837 ringAtoms[tmp[j]].push_back(*k); 838 } 839 } 840 } 841 842 843 //loop over closure bonds 844 for(OBBondIterator bd(cbonds.begin()),bd_end(cbonds.end());bd!=bd_end;++bd) 845 { 846 bond = *bd; 847 848 // BASIC APPROACH 849 // pick beginning atom at closure bond 850 // this is really ready, isn't it ! ;-) 851 rootAtom = bond->GetBeginAtomIdx(); 852 _root[rootAtom] = true; 853 854 // EXTENDED APPROACH 855 if (avoidInnerRingAtoms) 856 { 857 // count the number of neighbor ring atoms 858 atom = mol.GetAtom(rootAtom); 859 ringNbrs = heavyNbrs = 0; 860 861 for (nbr = atom->BeginNbrAtom(l);nbr;nbr = atom->NextNbrAtom(l)) 862 { 863 // we can get this from atom->GetHvyDegree() 864 // but we need to find neighbors in rings too 865 // so let's save some time 866 if (nbr->GetAtomicNum() != OBElements::Hydrogen) 867 { 868 heavyNbrs++; 869 if (nbr->IsInRing()) 870 ringNbrs++; 871 } 872 873 // if this atom has more than 2 neighbor ring atoms 874 // we could get trapped later when traversing cycles 875 // which can cause aromaticity false detection 876 newRoot = -1; 877 878 if (ringNbrs > 2) 879 { 880 // try to find another root atom 881 // only loop over rings which contain rootAtom 882 for(k = ringAtoms[rootAtom].begin() ; k != ringAtoms[rootAtom].end(); ++k) 883 { 884 ring = (*k); 885 tmp = ring->_path; 886 887 bool checkThisRing = false; 888 int rootAtomNumber=0; 889 int idx=0; 890 // avoiding two root atoms in one ring ! 891 for (unsigned int j = 0; j < tmpRootAtoms.size(); ++j) 892 { 893 idx= tmpRootAtoms[j]; 894 if(ring->IsInRing(idx)) 895 { 896 rootAtomNumber++; 897 if(rootAtomNumber>=2) 898 break; 899 } 900 } 901 if(rootAtomNumber<2) 902 { 903 for (unsigned int j = 0; j < tmp.size(); ++j) 904 { 905 // find critical ring 906 if (tmp[j] == rootAtom) 907 { 908 checkThisRing = true; 909 } 910 else 911 { 912 // second root atom in this ring ? 913 if (_root[tmp[j]] == true) 914 { 915 // when there is a second root 916 // atom this ring can not be 917 // used for getting an other 918 // root atom 919 checkThisRing = false; 920 921 break; 922 } 923 } 924 } 925 } 926 927 // check ring for getting another 928 // root atom to avoid aromaticity typer problems 929 if (checkThisRing) 930 { 931 // check if we can find another root atom 932 for (unsigned int m = 0; m < tmp.size(); ++m) 933 { 934 ringNbrs = heavyNbrs = 0; 935 for (nbr2 = (mol.GetAtom(tmp[m]))->BeginNbrAtom(nbr2Iter); 936 nbr2;nbr2 = (mol.GetAtom(tmp[m]))->NextNbrAtom(nbr2Iter)) 937 { 938 if (nbr2->GetAtomicNum() != OBElements::Hydrogen) 939 { 940 heavyNbrs++; 941 942 if (nbr2->IsInRing()) 943 ringNbrs++; 944 } 945 } 946 947 // if the number of neighboured heavy atoms is also 948 // the number of neighboured ring atoms, the aromaticity 949 // typer could be stuck in a local traversing trap 950 if (ringNbrs <= 2 && ring->IsInRing((mol.GetAtom(tmp[m])->GetIdx()))) 951 { 952 newRoot = tmp[m]; 953 } 954 } 955 } 956 } 957 958 if ((newRoot != -1) && (rootAtom != newRoot)) 959 { 960 // unset root atom 961 _root[rootAtom] = false; 962 963 // pick new root atom 964 _root[newRoot] = true; 965 } 966 } // if (ringNbrs > 2) 967 968 } // end for 969 } // if (avoid) 970 } // end for(closure bonds) 971 } 972 ExcludeSmallRing()973 void OBAromaticTyperMolState::ExcludeSmallRing() 974 { 975 OBAtom *atom,*nbr1,*nbr2; 976 vector<OBAtom*>::iterator i; 977 vector<OBBond*>::iterator j,k; 978 979 for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) 980 if (_root[atom->GetIdx()]) 981 for (nbr1 = atom->BeginNbrAtom(j);nbr1;nbr1 = atom->NextNbrAtom(j)) 982 if ((*j)->IsInRing() && _vpa[nbr1->GetIdx()]) 983 for (nbr2 = nbr1->BeginNbrAtom(k);nbr2;nbr2 = nbr1->NextNbrAtom(k)) 984 if (nbr2 != atom && (*k)->IsInRing() && _vpa[nbr2->GetIdx()]) 985 if (atom->IsConnected(nbr2)) 986 _root[atom->GetIdx()] = false; 987 } 988 989 } //namespace OpenBabel; 990 991 //! \file typer.cpp 992 //! \brief Open Babel atom and aromaticity typer. 993