1 /********************************************************************** 2 atom.cpp - Handle OBAtom class. 3 4 Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. 5 Some portions Copyright (C) 2001-2008 by Geoffrey R. Hutchison 6 Some portions Copyright (C) 2003 by Michael Banck 7 8 This file is part of the Open Babel project. 9 For more information, see <http://openbabel.org/> 10 11 This program is free software; you can redistribute it and/or modify 12 it under the terms of the GNU General Public License as published by 13 the Free Software Foundation version 2 of the License. 14 15 This program is distributed in the hope that it will be useful, 16 but WITHOUT ANY WARRANTY; without even the implied warranty of 17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 GNU General Public License for more details. 19 ***********************************************************************/ 20 #include <openbabel/babelconfig.h> 21 22 #include <openbabel/atom.h> 23 #include <openbabel/bond.h> 24 #include <openbabel/mol.h> 25 #include <openbabel/obiter.h> 26 #include <openbabel/generic.h> 27 #include <openbabel/molchrg.h> 28 #include <openbabel/ring.h> 29 #include <openbabel/phmodel.h> 30 #include <openbabel/builder.h> 31 #include <openbabel/elements.h> 32 #include <openbabel/chains.h> 33 #include <openbabel/obutil.h> 34 #include <openbabel/residue.h> 35 #include <openbabel/chains.h> 36 37 #include <openbabel/math/matrix3x3.h> 38 39 #if !HAVE_STRNCASECMP 40 extern "C" int strncasecmp(const char *s1, const char *s2, size_t n); 41 #endif 42 43 using namespace std; 44 45 46 namespace OpenBabel 47 { 48 OB_EXTERN OBChainsParser chainsparser; 49 /** \class OBAtom atom.h <openbabel/atom.h> 50 \brief Atom class 51 52 To understand the OBAtom class it is important to state a key 53 decision on which the design was based. The OBAtom class not only 54 holds data, but facilitates extraction of data perceived from both 55 the atom and the molecule. This prevents the OBMol class from 56 becoming overly large and complicated. 57 58 A number of data extraction methods perform what is called 59 <a href="http://en.wikipedia.org/wiki/Lazy_evaluation">`Lazy Evaluation,'</a> 60 which is essentially on-the-fly evaluation. 61 For example, when an atom is queried as to whether it is cyclic 62 or what it's hybridization state is the information is perceived 63 automatically. The perception of a particular trait is actually 64 performed on the entire molecule the first time it is requested of 65 an atom or bond, and stored for subsequent requests for the same 66 trait of additional atoms or bonds.The OBAtom class is similar to 67 OBMol and the whole of Open Babel in that data access and modification 68 is done through Get and Set methods. 69 70 The following code demonstrates how to print out the atom numbers, 71 element numbers, and coordinates of a molecule: 72 \code 73 OBMol mol; 74 75 FOR_ATOMS_OF_MOL(atom, mol) 76 { 77 cout << atom->GetIdx() << ` `; 78 cout << atom->GetAtomicNum() << ` `; 79 cout << atom->GetVector() << endl; 80 } 81 \endcode 82 A number of the property member functions indicate that atoms 83 have some knowledge of their covalently attached neighbor atoms. 84 Bonding information is partly redundant within a molecule in 85 that an OBMol has a complete list of bonds in a molecule, and 86 an OBAtom has a list bonds of which it is a member. The following 87 code demonstrates how an OBAtom uses its bond information to loop 88 over atoms attached to itself: 89 \code 90 OBMol mol; 91 OBAtom *atom; 92 93 atom = mol.GetAtom(1); 94 FOR_NBORS_OF_ATOM(nbr, atom) 95 { 96 cout << "atom #" << atom->GetIdx() << " is attached to atom #" 97 << nbr->GetIdx() << endl; 98 } 99 \endcode 100 101 should produce an output like 102 \code 103 atom #1 is attached to atom #2 104 \endcode 105 */ 106 107 extern THREAD_LOCAL OBAromaticTyper aromtyper; 108 extern THREAD_LOCAL OBAtomTyper atomtyper; 109 extern THREAD_LOCAL OBPhModel phmodel; 110 OB_EXTERN OBTypeTable ttab; 111 112 // 113 // OBAtom member functions 114 // 115 OBAtom()116 OBAtom::OBAtom() 117 { 118 _parent = nullptr; 119 Clear(); 120 } 121 ~OBAtom()122 OBAtom::~OBAtom() 123 { 124 if (_residue != nullptr) 125 { 126 _residue->RemoveAtom(this); 127 } 128 /* 129 if (!_vdata.empty()) 130 { 131 vector<OBGenericData*>::iterator m; 132 for (m = _vdata.begin();m != _vdata.end();++m) 133 delete *m; 134 _vdata.clear(); 135 } 136 */ 137 } 138 Clear()139 bool OBAtom::Clear() 140 { 141 _c = nullptr; 142 _cidx = 0; 143 _flags=0; 144 _idx = 0; 145 _hyb = 0; 146 _ele = (char)0; 147 _isotope = 0; 148 _spinmultiplicity=0; // CM 18 Sept 2003 149 _imph = 0; 150 _fcharge = 0; 151 _type[0] = '\0'; 152 _pcharge = 0.0; 153 _vbond.clear(); 154 _vbond.reserve(4); 155 _residue = nullptr; 156 _id = NoId; 157 158 return(OBBase::Clear()); 159 } 160 operator =(OBAtom & src)161 OBAtom &OBAtom::operator=(OBAtom &src) 162 //copy atom information 163 //bond info is not copied here as ptrs may be invalid 164 { 165 if (this != &src) { 166 _idx = src.GetIdx(); 167 Duplicate(&src); 168 } 169 return(*this); 170 } 171 Duplicate(OBAtom * src)172 void OBAtom::Duplicate(OBAtom *src) 173 { 174 if (!src) 175 return; 176 177 _hyb = src->_hyb; 178 _ele = src->_ele; 179 _imph = src->_imph; 180 _isotope = src->_isotope; 181 _fcharge = src->_fcharge; 182 _spinmultiplicity = src->_spinmultiplicity; 183 strncpy(_type,src->_type, sizeof(_type) - 1); 184 _type[sizeof(_type) - 1] = '\0'; 185 _pcharge = src->_pcharge; 186 _v = src->GetVector(); 187 _flags = src->_flags; 188 _residue = nullptr; 189 _id = src->_id; 190 191 _vdata.clear(); 192 //Copy all the OBGenericData, providing the new atom 193 vector<OBGenericData*>::iterator itr; 194 for(itr=src->BeginData();itr!=src->EndData();++itr) 195 { 196 OBGenericData* pCopiedData = (*itr)->Clone(this); 197 SetData(pCopiedData); 198 } 199 } 200 IsConnected(OBAtom * a1)201 bool OBAtom::IsConnected(OBAtom *a1) 202 { 203 OBBondIterator i; 204 OBBond *bond; 205 206 for (bond = BeginBond(i);bond;bond = NextBond(i)) 207 if (bond->GetBeginAtom() == a1 || bond->GetEndAtom() == a1) 208 return(true); 209 210 return(false); 211 } 212 IsOneThree(OBAtom * a1)213 bool OBAtom::IsOneThree(OBAtom *a1) 214 { 215 OBAtom *atom1,*atom2; 216 OBBond *bond1,*bond2; 217 OBBondIterator i,j; 218 atom1 = this; 219 atom2 = a1; 220 221 for (bond1 = atom1->BeginBond(i);bond1;bond1 = atom1->NextBond(i)) 222 for (bond2 = atom2->BeginBond(j);bond2;bond2 = atom2->NextBond(j)) 223 if (bond1->GetNbrAtom(atom1) == bond2->GetNbrAtom(atom2)) 224 return(true); 225 226 return(false); 227 } 228 IsOneFour(OBAtom * a1)229 bool OBAtom::IsOneFour(OBAtom *a1) 230 { 231 OBAtom *atom1,*atom2; 232 OBBond *bond1,*bond2; 233 OBBondIterator i,j; 234 atom1 = this; 235 atom2 = a1; 236 237 for (bond1 = atom1->BeginBond(i);bond1;bond1 = atom1->NextBond(i)) 238 for (bond2 = atom2->BeginBond(j);bond2;bond2 = atom2->NextBond(j)) 239 if ((bond1->GetNbrAtom(atom1))->IsConnected(bond2->GetNbrAtom(atom2))) 240 return(true); 241 242 return(false); 243 } 244 IsAxial()245 bool OBAtom::IsAxial() 246 { 247 double tor; 248 OBAtom *a,*b,*c; 249 OBBondIterator i,j,k; 250 251 for (a = BeginNbrAtom(i);a;a = NextNbrAtom(i)) 252 if (a->GetHyb() == 3 && a->IsInRing() && !(*i)->IsInRing()) 253 for (b = a->BeginNbrAtom(j);b;b = a->NextNbrAtom(j)) 254 if (b != this && b->IsInRing() && b->GetHyb() == 3) 255 for (c = b->BeginNbrAtom(k);c;c = b->NextNbrAtom(k)) 256 if (c != a && c->IsInRing()) 257 { 258 tor = fabs(((OBMol*)GetParent())->GetTorsion(this,a,b,c)); 259 return(tor > 55.0 && tor < 75.0); 260 } 261 262 return(false); 263 } 264 265 /** This can be sketched as follows 266 \code 267 '*' 268 \ 269 a=b 270 \endcode 271 where a and b are the 'apha' and 'beta' atoms, respectively and '*' 272 indicates the current atom. 273 **/ HasAlphaBetaUnsat(bool includePandS)274 bool OBAtom::HasAlphaBetaUnsat(bool includePandS) 275 { 276 OBAtom *a1,*a2; 277 OBBondIterator i,j; 278 279 for (a1 = BeginNbrAtom(i);a1;a1 = NextNbrAtom(i)) 280 if (includePandS || (a1->GetAtomicNum() != OBElements::Phosphorus && a1->GetAtomicNum() != OBElements::Sulfur)) 281 for (a2 = a1->BeginNbrAtom(j);a2;a2 = a1->NextNbrAtom(j)) 282 if (a2 != this && ((*j)->GetBondOrder() == 2 || (*j)->GetBondOrder() == 3 || (*j)->GetBondOrder() == 5)) 283 return(true); 284 285 return(false); 286 } 287 HasBondOfOrder(unsigned int order)288 bool OBAtom::HasBondOfOrder(unsigned int order) 289 { 290 OBBond *bond; 291 OBBondIterator i; 292 for (bond = BeginBond(i);bond;bond = NextBond(i)) 293 if (bond->GetBondOrder() == order) 294 return(true); 295 296 return(false); 297 } 298 CountBondsOfOrder(unsigned int order)299 int OBAtom::CountBondsOfOrder(unsigned int order) 300 { 301 int count = 0; 302 OBBond *bond; 303 OBBondIterator i; 304 for (bond = BeginBond(i);bond;bond = NextBond(i)) 305 if (bond->GetBondOrder() == order) 306 count++; 307 308 return(count); 309 } 310 HighestBondOrder()311 int OBAtom::HighestBondOrder() 312 { 313 unsigned int highest = 0; 314 OBBond *bond; 315 OBBondIterator i; 316 for(bond = BeginBond(i); bond; bond = NextBond(i)) 317 if(bond->GetBondOrder() > highest) 318 highest = bond->GetBondOrder(); 319 320 return(highest); 321 } 322 HasNonSingleBond()323 bool OBAtom::HasNonSingleBond() 324 { 325 OBBond *bond; 326 OBBondIterator i; 327 for (bond = BeginBond(i);bond;bond = NextBond(i)) 328 if (bond->GetBondOrder() != 1) 329 return(true); 330 331 return(false); 332 } 333 IsPolarHydrogen()334 bool OBAtom::IsPolarHydrogen() 335 { 336 if (GetAtomicNum() != OBElements::Hydrogen) 337 return(false); 338 339 OBAtom *atom; 340 OBBond *bond; 341 OBBondIterator i; 342 for (bond = BeginBond(i);bond;bond = NextBond(i)) 343 { 344 atom = bond->GetNbrAtom(this); 345 if (atom->GetAtomicNum() == 7) 346 return(true); 347 if (atom->GetAtomicNum() == 8) 348 return(true); 349 if (atom->GetAtomicNum() == 15) 350 return(true); 351 if (atom->GetAtomicNum() == 16) 352 return(true); 353 } 354 355 return(false); 356 } 357 IsNonPolarHydrogen()358 bool OBAtom::IsNonPolarHydrogen() 359 { 360 if (GetAtomicNum() != OBElements::Hydrogen) 361 return(false); 362 363 OBAtom *atom; 364 OBBond *bond; 365 OBBondIterator i; 366 for (bond = BeginBond(i);bond;bond = NextBond(i)) 367 { 368 atom = bond->GetNbrAtom(this); 369 if (atom->GetAtomicNum() == 6) 370 return(true); 371 } 372 373 return(false); 374 } 375 GetVector()376 vector3 &OBAtom::GetVector() 377 { 378 if (!_c) 379 return(_v); 380 381 _v.Set((*_c)[_cidx],(*_c)[_cidx+1],(*_c)[_cidx+2]); 382 return(_v); 383 } 384 GetVector() const385 const vector3 &OBAtom::GetVector() const 386 { 387 if (!_c) 388 return(_v); 389 390 _v.Set((*_c)[_cidx],(*_c)[_cidx+1],(*_c)[_cidx+2]); 391 return(_v); 392 } 393 SetVector()394 void OBAtom::SetVector() 395 { 396 // obAssert(_c); 397 if (_c) 398 _v.Set((*_c)[_cidx],(*_c)[_cidx+1],(*_c)[_cidx+2]); 399 } 400 SetVector(const vector3 & v)401 void OBAtom::SetVector(const vector3 &v) 402 { 403 if (!_c) 404 _v = v; 405 else 406 { 407 (*_c)[_cidx ] = v.x(); 408 (*_c)[_cidx+1] = v.y(); 409 (*_c)[_cidx+2] = v.z(); 410 } 411 } 412 SetVector(const double v_x,const double v_y,const double v_z)413 void OBAtom::SetVector(const double v_x,const double v_y,const double v_z) 414 { 415 if (!_c) 416 _v.Set(v_x,v_y,v_z); 417 else 418 { 419 (*_c)[_cidx ] = v_x; 420 (*_c)[_cidx+1] = v_y; 421 (*_c)[_cidx+2] = v_z; 422 } 423 } 424 SetType(const char * type)425 void OBAtom::SetType(const char *type) 426 { 427 strncpy(_type,type, sizeof(_type) - 1); 428 _type[sizeof(_type) - 1] = '\0'; 429 if (_ele == 1 && type[0] == 'D') 430 _isotope = 2; 431 } 432 SetType(const string & type)433 void OBAtom::SetType(const string &type) 434 { 435 strncpy(_type,type.c_str(), sizeof(_type) - 1); 436 _type[sizeof(_type) - 1] = '\0'; 437 if (_ele == 1 && type[0] == 'D') 438 _isotope = 2; 439 } 440 SetIsotope(unsigned int iso)441 void OBAtom::SetIsotope(unsigned int iso) 442 { 443 _isotope = iso; 444 } 445 GetResidue()446 OBResidue *OBAtom::GetResidue() 447 { 448 OBMol *mol = this->GetParent(); 449 if (!mol->HasChainsPerceived()) 450 chainsparser.PerceiveChains(*mol); 451 452 return _residue; 453 } 454 GetAtomicMass() const455 double OBAtom::GetAtomicMass() const 456 { 457 if (_isotope == 0) 458 return OBElements::GetMass(_ele); 459 else 460 return OBElements::GetExactMass(_ele, _isotope); 461 } 462 GetExactMass() const463 double OBAtom::GetExactMass() const 464 { 465 return OBElements::GetExactMass(_ele, _isotope); 466 } 467 GetType()468 char *OBAtom::GetType() 469 { 470 OBMol *mol = (OBMol*)GetParent(); 471 if (mol) 472 if (!mol->HasAtomTypesPerceived()) 473 atomtyper.AssignTypes(*((OBMol*)GetParent())); 474 475 if (strlen(_type) == 0) // Somehow we still don't have a type! 476 { 477 char num[6]; 478 string fromType = ttab.GetFromType(); // save previous types 479 string toType = ttab.GetToType(); 480 481 ttab.SetFromType("ATN"); 482 ttab.SetToType("INT"); 483 snprintf(num, 6, "%d", GetAtomicNum()); 484 ttab.Translate(_type, num); 485 486 ttab.SetFromType(fromType.c_str()); 487 ttab.SetToType(toType.c_str()); 488 } 489 if (_ele == 1 && _isotope == 2) 490 snprintf(_type, 6, "%s", "D"); 491 492 return(_type); 493 } 494 GetHyb() const495 unsigned int OBAtom::GetHyb() const 496 { 497 //hybridization is assigned when atoms are typed 498 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 499 if (mol && !mol->HasHybridizationPerceived()) 500 atomtyper.AssignHyb(*mol); 501 502 return(_hyb); 503 } 504 505 GetHvyDegree() const506 unsigned int OBAtom::GetHvyDegree() const 507 { 508 unsigned int count=0; 509 510 OBBond *bond; 511 OBBondIterator i; 512 for (bond = ((OBAtom*)this)->BeginBond(i); bond; bond = ((OBAtom*)this)->NextBond(i)) 513 if (bond->GetNbrAtom((OBAtom*)this)->GetAtomicNum() != OBElements::Hydrogen) 514 count++; 515 516 return(count); 517 } 518 GetHeteroDegree() const519 unsigned int OBAtom::GetHeteroDegree() const 520 { 521 unsigned int count=0; 522 OBBond *bond; 523 OBBondIterator i; 524 for (bond = ((OBAtom*)this)->BeginBond(i);bond;bond = ((OBAtom*)this)->NextBond(i)) 525 if (bond->GetNbrAtom((OBAtom*)this)->IsHeteroatom()) 526 count++; 527 528 return((unsigned int)count); 529 } 530 GetPartialCharge()531 double OBAtom::GetPartialCharge() 532 { 533 if (!GetParent()) 534 return(_pcharge); 535 if (!((OBMol*)GetParent())->AutomaticPartialCharge()) 536 return(_pcharge); 537 538 if (!((OBMol*)GetParent())->HasPartialChargesPerceived()) 539 { 540 //seed partial charges are set in the atom typing procedure 541 OBAtom *atom; 542 OBMol *mol = (OBMol*)GetParent(); 543 vector<OBAtom*>::iterator i; 544 for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i)) 545 atom->SetPartialCharge(0.0); 546 547 phmodel.AssignSeedPartialCharge(*((OBMol*)GetParent())); 548 OBGastChrg gc; 549 gc.AssignPartialCharges(*((OBMol*)GetParent())); 550 } 551 552 return(_pcharge); 553 } 554 555 //! Returns true if nitrogen is part of an amide IsAmideNitrogen()556 bool OBAtom::IsAmideNitrogen() 557 { 558 if (GetAtomicNum() != OBElements::Nitrogen) 559 return(false); 560 561 OBAtom *nbratom,*atom; 562 OBBond *abbond,*bond; 563 564 OBBondIterator i,j; 565 atom = this; 566 for (bond = BeginBond(i);bond;bond = NextBond(i)) 567 { 568 nbratom = bond->GetNbrAtom(atom); 569 for (abbond = nbratom->BeginBond(j);abbond;abbond = nbratom->NextBond(j)) 570 if (abbond->GetBondOrder() == 2 && 571 (((abbond->GetNbrAtom(nbratom))->GetAtomicNum() == 8) || 572 ((abbond->GetNbrAtom(nbratom))->GetAtomicNum() == 16))) 573 return(true); 574 } 575 576 return(false); 577 } 578 IsAromaticNOxide()579 bool OBAtom::IsAromaticNOxide() 580 { 581 if (GetAtomicNum() != OBElements::Nitrogen || !IsAromatic()) 582 return(false); 583 584 OBAtom *atom; 585 OBBondIterator i; 586 587 for (atom = BeginNbrAtom(i);atom;atom = NextNbrAtom(i)) 588 if (atom->GetAtomicNum() == OBElements::Oxygen && !(*i)->IsInRing() && (*i)->GetBondOrder() == 2) 589 return(true); 590 591 return(false); 592 } 593 IsCarboxylOxygen()594 bool OBAtom::IsCarboxylOxygen() 595 { 596 if (GetAtomicNum() != OBElements::Oxygen) 597 return(false); 598 if (GetHvyDegree() != 1) 599 return(false); 600 601 OBAtom *atom; 602 OBBond *bond; 603 OBBondIterator i; 604 605 atom = nullptr; 606 for (bond = BeginBond(i);bond;bond = NextBond(i)) 607 if ((bond->GetNbrAtom(this))->GetAtomicNum() == OBElements::Carbon) 608 { 609 atom = bond->GetNbrAtom(this); 610 break; 611 } 612 if (!atom) 613 return(false); 614 if (!(atom->CountFreeOxygens() == 2) 615 && !(atom->CountFreeOxygens() == 1 && atom->CountFreeSulfurs() == 1)) 616 return(false); 617 618 //atom is connected to a carbon that has a total 619 //of 2 attached free oxygens or 1 free oxygen and 1 free sulfur 620 return(true); 621 } 622 IsPhosphateOxygen()623 bool OBAtom::IsPhosphateOxygen() 624 { 625 if (GetAtomicNum() != OBElements::Oxygen) 626 return(false); 627 if (GetHvyDegree() != 1) 628 return(false); 629 OBAtom *atom; 630 OBBond *bond; 631 OBBondIterator i; 632 633 atom = nullptr; 634 for (bond = BeginBond(i);bond;bond = NextBond(i)) 635 if ((bond->GetNbrAtom(this))->GetAtomicNum() == OBElements::Phosphorus) 636 { 637 atom = bond->GetNbrAtom(this); 638 break; 639 } 640 if (!atom) 641 return(false); 642 if (atom->CountFreeOxygens() > 2) 643 return(true); 644 645 //atom is connected to a carbon that has a total 646 //of 2 attached free oxygens 647 return(false); 648 } 649 IsSulfateOxygen()650 bool OBAtom::IsSulfateOxygen() 651 { 652 if (GetAtomicNum() != OBElements::Oxygen) 653 return(false); 654 if (GetHvyDegree() != 1) 655 return(false); 656 657 OBAtom *atom; 658 OBBond *bond; 659 OBBondIterator i; 660 661 atom = nullptr; 662 for (bond = BeginBond(i);bond;bond = NextBond(i)) 663 if ((bond->GetNbrAtom(this))->GetAtomicNum() == OBElements::Sulfur) 664 { 665 atom = bond->GetNbrAtom(this); 666 break; 667 } 668 if (!atom) 669 return(false); 670 if (atom->CountFreeOxygens() < 3) 671 return(false); 672 673 //atom is connected to a carbon that has a total 674 //of 2 attached free oxygens 675 return(true); 676 } 677 678 // Helper function for IsHBondAcceptor IsSulfoneOxygen(OBAtom * atm)679 static bool IsSulfoneOxygen(OBAtom* atm) 680 // Stefano Forli 681 //atom is connected to a sulfur that has a total 682 //of 2 attached free oxygens, and it's not a sulfonamide 683 //e.g. C-SO2-C 684 // Is this atom an oxygen in a sulfone(R1 - SO2 - R2) group ? 685 { 686 if (atm->GetAtomicNum() != OBElements::Oxygen) 687 return(false); 688 if (atm->GetHvyDegree() != 1){ 689 //cerr << "sulfone> O valence is not 1\n"; 690 return(false); 691 } 692 693 OBAtom *nbr = nullptr; 694 OBBond *bond1,*bond2; 695 OBBondIterator i,j; 696 697 // searching for attached sulfur 698 for (bond1 = atm->BeginBond(i); bond1; bond1 = atm->NextBond(i)) 699 if ((bond1->GetNbrAtom(atm))->GetAtomicNum() == OBElements::Sulfur) 700 { nbr = bond1->GetNbrAtom(atm); 701 break; } 702 if (!nbr){ 703 //cerr << "sulfone> atom null\n" ; 704 return(false); } 705 706 // check for sulfate 707 //cerr << "sulfone> If we're here... " << atom->GetAtomicNum() <<"\n" << atom->GetAtomicNum() == OBElements::Sulfur << "\n"; 708 //cerr << "sulfone> number of free oxygens:" << atom->CountFreeOxygens() << "\n"; 709 if (nbr->CountFreeOxygens() != 2){ 710 //cerr << "sulfone> count of free oxygens not 2" << atom->CountFreeOxygens() << '\n' ; 711 return(false); } 712 713 // check for sulfonamide 714 for (bond2 = nbr->BeginBond(j);bond2;bond2 = nbr->NextBond(j)){ 715 //cerr<<"NEIGH: " << (bond2->GetNbrAtom(atom))->GetAtomicNum()<<"\n"; 716 if ((bond2->GetNbrAtom(nbr))->GetAtomicNum() == OBElements::Nitrogen){ 717 //cerr << "sulfone> sulfonamide null\n" ; 718 return(false);}} 719 //cerr << "sulfone> none of the above\n"; 720 return(true); // true sulfone 721 } 722 723 724 725 726 IsNitroOxygen()727 bool OBAtom::IsNitroOxygen() 728 { 729 if (GetAtomicNum() != OBElements::Oxygen) 730 return(false); 731 if (GetHvyDegree() != 1) 732 return(false); 733 734 OBAtom *atom; 735 OBBond *bond; 736 OBBondIterator i; 737 738 atom = nullptr; 739 for (bond = BeginBond(i);bond;bond = NextBond(i)) 740 if ((bond->GetNbrAtom(this))->GetAtomicNum() == OBElements::Nitrogen) 741 { 742 atom = bond->GetNbrAtom(this); 743 break; 744 } 745 if (!atom) 746 return(false); 747 if (atom->CountFreeOxygens() != 2) 748 return(false); 749 750 //atom is connected to a nitrogen that has a total 751 //of 2 attached free oxygens 752 return(true); 753 } 754 IsHeteroatom()755 bool OBAtom::IsHeteroatom() 756 { 757 switch(GetAtomicNum()) 758 { 759 case 7: 760 case 8: 761 case 15: 762 case 16: 763 case 33: 764 case 34: 765 case 51: 766 case 52: 767 case 83: 768 case 84: 769 return(true); 770 } 771 return(false); 772 } 773 IsAromatic() const774 bool OBAtom::IsAromatic() const 775 { 776 OBMol *mol = ((OBAtom*)this)->GetParent(); 777 if (!mol->HasAromaticPerceived()) 778 aromtyper.AssignAromaticFlags(*mol); 779 780 if (((OBAtom*)this)->HasFlag(OB_AROMATIC_ATOM)) 781 return true; 782 783 return false; 784 } 785 IsInRing() const786 bool OBAtom::IsInRing() const 787 { 788 OBMol *mol = ((OBAtom*)this)->GetParent(); 789 if (!mol->HasRingAtomsAndBondsPerceived()) 790 mol->FindRingAtomsAndBonds(); 791 792 if (((OBAtom*)this)->HasFlag(OB_RING_ATOM)) 793 return true; 794 795 return false; 796 } 797 798 //! @todo IsChiral()799 bool OBAtom::IsChiral() 800 { 801 OBMol *mol = (OBMol*) GetParent(); 802 OBStereoFacade stereoFacade(mol); 803 return stereoFacade.HasTetrahedralStereo(_id); 804 } 805 IsPeriodic() const806 bool OBAtom::IsPeriodic() const 807 { 808 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 809 return mol->IsPeriodic(); 810 } 811 IsInRingSize(int size) const812 bool OBAtom::IsInRingSize(int size) const 813 { 814 vector<OBRing*> rlist; 815 vector<OBRing*>::iterator i; 816 817 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 818 if (!mol->HasSSSRPerceived()) 819 mol->FindSSSR(); 820 821 if (!((OBAtom*)this)->HasFlag(OB_RING_ATOM)) 822 return(false); 823 824 rlist = mol->GetSSSR(); 825 for (i = rlist.begin();i != rlist.end();++i) 826 if ((*i)->IsInRing(GetIdx()) && static_cast<int>((*i)->PathSize()) == size) 827 return(true); 828 829 return(false); 830 } 831 MemberOfRingCount() const832 unsigned int OBAtom::MemberOfRingCount() const 833 { 834 vector<OBRing*> rlist; 835 vector<OBRing*>::iterator i; 836 unsigned int count=0; 837 838 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 839 840 if (!mol->HasSSSRPerceived()) 841 mol->FindSSSR(); 842 843 if (!((OBAtom*)this)->IsInRing()) 844 return(0); 845 846 rlist = mol->GetSSSR(); 847 848 for (i = rlist.begin();i != rlist.end();++i) 849 if ((*i)->IsInRing(GetIdx())) 850 count++; 851 852 return((unsigned int)count); 853 } 854 MemberOfRingSize() const855 unsigned int OBAtom::MemberOfRingSize() const 856 { 857 vector<OBRing*> rlist; 858 vector<OBRing*>::iterator i; 859 860 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 861 862 if (!mol->HasSSSRPerceived()) 863 mol->FindSSSR(); 864 865 if (!((OBAtom*)this)->IsInRing()) 866 return(0); 867 868 rlist = mol->GetSSSR(); 869 870 for (i = rlist.begin();i != rlist.end();++i) 871 if ((*i)->IsInRing(GetIdx())) 872 return((*i)->Size()); 873 874 return(0); 875 } 876 CountRingBonds() const877 unsigned int OBAtom::CountRingBonds() const 878 { 879 unsigned int count = 0; 880 OBBond *bond; 881 OBBondIterator i; 882 883 for (bond = ((OBAtom*)this)->BeginBond(i); 884 bond; 885 bond = ((OBAtom*)this)->NextBond(i)) 886 { 887 if (bond->IsInRing()) 888 count++; 889 } 890 return count; 891 } 892 SmallestBondAngle()893 double OBAtom::SmallestBondAngle() 894 { 895 OBAtom *b, *c; 896 vector3 v1, v2; 897 double degrees, minDegrees; 898 // vector<OBAtom*>::iterator i; 899 OBBondIterator j,k; 900 901 minDegrees = 360.0; 902 903 for (b = BeginNbrAtom(j); b; b = NextNbrAtom(j)) 904 { 905 k = j; 906 for (c = NextNbrAtom(k); c; c = NextNbrAtom(k)) 907 { 908 degrees = b->GetAngle((OBAtom*)this, c); 909 if (degrees < minDegrees) 910 minDegrees = degrees; 911 } 912 } 913 return minDegrees; 914 } 915 AverageBondAngle()916 double OBAtom::AverageBondAngle() 917 { 918 OBAtom *b, *c; 919 vector3 v1, v2; 920 double degrees, avgDegrees; 921 // vector<OBAtom*>::iterator i; 922 OBBondIterator j,k; 923 int n=0; 924 925 avgDegrees = 0.0; 926 927 for (b = BeginNbrAtom(j); b; b = NextNbrAtom(j)) 928 { 929 k = j; 930 for (c = NextNbrAtom(k); c; c = NextNbrAtom(k)) 931 { 932 degrees = b->GetAngle((OBAtom*)this, c); 933 avgDegrees += degrees; 934 n++; 935 } 936 } 937 938 if (n >=1) 939 avgDegrees /= n; 940 941 return avgDegrees; 942 } 943 CountFreeOxygens() const944 unsigned int OBAtom::CountFreeOxygens() const 945 { 946 unsigned int count = 0; 947 OBAtom *atom; 948 OBBond *bond; 949 OBBondIterator i; 950 951 for (bond = ((OBAtom*)this)->BeginBond(i);bond;bond = ((OBAtom*)this)->NextBond(i)) 952 { 953 atom = bond->GetNbrAtom((OBAtom*)this); 954 if (atom->GetAtomicNum() == OBElements::Oxygen && atom->GetHvyDegree() == 1) 955 count++; 956 } 957 958 return(count); 959 } 960 CountFreeSulfurs() const961 unsigned int OBAtom::CountFreeSulfurs() const 962 { 963 unsigned int count = 0; 964 OBAtom *atom; 965 OBBond *bond; 966 OBBondIterator i; 967 968 for (bond = ((OBAtom*)this)->BeginBond(i);bond;bond = ((OBAtom*)this)->NextBond(i)) 969 { 970 atom = bond->GetNbrAtom((OBAtom*)this); 971 if (atom->GetAtomicNum() == OBElements::Sulfur && atom->GetHvyDegree() == 1) 972 count++; 973 } 974 975 return(count); 976 } 977 GetExplicitValence() const978 unsigned int OBAtom::GetExplicitValence() const 979 { 980 unsigned int bosum = 0; 981 982 OBBondIterator i; 983 for (OBBond *bond = ((OBAtom*)this)->BeginBond(i); bond; bond = ((OBAtom*)this)->NextBond(i)) 984 bosum += bond->GetBondOrder(); 985 986 return bosum; 987 } 988 GetTotalValence() const989 unsigned int OBAtom::GetTotalValence() const 990 { 991 return GetExplicitValence() + _imph; 992 } 993 ExplicitHydrogenCount(bool ExcludeIsotopes) const994 unsigned int OBAtom::ExplicitHydrogenCount(bool ExcludeIsotopes) const 995 { 996 //If ExcludeIsotopes is true, H atoms with _isotope!=0 are not included. 997 //This excludes D, T and H when _isotope exlicitly set to 1 rather than the default 0. 998 int numH=0; 999 OBAtom *atom; 1000 OBBondIterator i; 1001 for (atom = ((OBAtom*)this)->BeginNbrAtom(i);atom;atom = ((OBAtom*)this)->NextNbrAtom(i)) 1002 if (atom->GetAtomicNum() == OBElements::Hydrogen && !(ExcludeIsotopes && atom->GetIsotope()!=0)) 1003 numH++; 1004 1005 return(numH); 1006 } 1007 1008 /** 1009 * The returned values count whole lone pairs, so the acid count is the number of 1010 * electron pairs desired and the base count is the number of electron pairs 1011 * available. 1012 * 1013 @verbatim 1014 Algorithm from: 1015 Clark, A. M. Accurate Specification of Molecular Structures: The Case for 1016 Zero-Order Bonds and Explicit Hydrogen Counting. Journal of Chemical Information 1017 and Modeling, 51, 3149-3157 (2011). http://pubs.acs.org/doi/abs/10.1021/ci200488k 1018 @endverbatim 1019 */ LewisAcidBaseCounts() const1020 pair<int, int> OBAtom::LewisAcidBaseCounts() const 1021 { 1022 // TODO: Is this data stored elsewhere? 1023 // The number of valence electrons in a free atom 1024 const int VALENCE[113] = {0,1,2,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,9,10, 1025 11,12,3,4,5,6,7,8,1,2,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,1,2, 1026 4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7, 1027 8,1,1,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,4,5,6,7,8,9,10,11,12}; 1028 // The number of electrons required to make up a full valence shell 1029 const int SHELL[113] = {0,2,2,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,18,18,18,18,18,18, 1030 18,18,18,18,8,8,8,8,8,8,8,8,18,18,18,18,18,18,18,18,18,18,8, 1031 8,8,8,8,8,8,8,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18, 1032 18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,8,8,18,18,18,18, 1033 18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18}; 1034 1035 pair<int, int> counts; 1036 int N = GetAtomicNum(); 1037 if (N == 0 || N > 112) { 1038 counts.first = 0; 1039 counts.second = 0; 1040 } else { 1041 int S = SHELL[N]; 1042 int V = VALENCE[N]; 1043 int C = GetFormalCharge(); 1044 int B = GetTotalValence(); 1045 // TODO: Do we actually want to divide by 2 here? (counting pairs instead of single) 1046 counts.first = (S - V - B + C) / 2; // Acid: Number of electrons pairs desired 1047 counts.second = (V - B - C) / 2; // Base: Number of electrons pairs available 1048 } 1049 return counts; 1050 } 1051 DeleteBond(OBBond * bond)1052 bool OBAtom::DeleteBond(OBBond *bond) 1053 { 1054 OBBondIterator i; 1055 for (i = _vbond.begin();i != _vbond.end();++i) 1056 if ((OBBond*)bond == *i) 1057 { 1058 _vbond.erase(i); 1059 return(true); 1060 } 1061 return(false); 1062 } 1063 MatchesSMARTS(const char * pattern)1064 bool OBAtom::MatchesSMARTS(const char *pattern) 1065 { 1066 OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent(); 1067 vector<vector<int> > mlist; 1068 vector<vector<int> >::iterator l; 1069 1070 OBSmartsPattern test; 1071 test.Init(pattern); 1072 if (test.Match(*mol)) 1073 { 1074 mlist = test.GetUMapList(); 1075 for (l = mlist.begin(); l != mlist.end(); ++l) 1076 if (GetIdx() == mol->GetAtom((*l)[0])->GetIdx()) 1077 return true; 1078 } 1079 return false; 1080 } 1081 BeginBond(OBBondIterator & i)1082 OBBond *OBAtom::BeginBond(OBBondIterator &i) 1083 { 1084 i = _vbond.begin(); 1085 return i == _vbond.end() ? nullptr : (OBBond*)*i; 1086 } 1087 NextBond(OBBondIterator & i)1088 OBBond *OBAtom::NextBond(OBBondIterator &i) 1089 { 1090 i++; 1091 return i == _vbond.end() ? nullptr : (OBBond*)*i; 1092 } 1093 BeginNbrAtom(OBBondIterator & i)1094 OBAtom *OBAtom::BeginNbrAtom(OBBondIterator &i) 1095 { 1096 i = _vbond.begin(); 1097 return i != _vbond.end() ? ((OBBond*) *i)->GetNbrAtom(this) : nullptr; 1098 } 1099 NextNbrAtom(OBBondIterator & i)1100 OBAtom *OBAtom::NextNbrAtom(OBBondIterator &i) 1101 { 1102 i++; 1103 return i != _vbond.end() ? ((OBBond*) *i)->GetNbrAtom(this) : nullptr; 1104 } 1105 GetDistance(OBAtom * b)1106 double OBAtom::GetDistance(OBAtom *b) 1107 { 1108 if (!IsPeriodic()) 1109 { 1110 return(( this->GetVector() - b->GetVector() ).length()); 1111 } 1112 else 1113 { 1114 OBUnitCell *box = (OBUnitCell*)GetParent()->GetData(OBGenericDataType::UnitCell); 1115 return (box->MinimumImageCartesian(this->GetVector() - b->GetVector())).length(); 1116 } 1117 } 1118 GetDistance(int b)1119 double OBAtom::GetDistance(int b) 1120 { 1121 OBMol *mol = (OBMol*)GetParent(); 1122 return( this->GetDistance(mol->GetAtom(b)) ); 1123 } 1124 GetDistance(vector3 * v)1125 double OBAtom::GetDistance(vector3 *v) 1126 { 1127 return(( this->GetVector() - *v ).length()); 1128 } 1129 GetAngle(OBAtom * b,OBAtom * c)1130 double OBAtom::GetAngle(OBAtom *b, OBAtom *c) 1131 { 1132 vector3 v1,v2; 1133 1134 v1 = this->GetVector() - b->GetVector(); 1135 v2 = c->GetVector() - b->GetVector(); 1136 if (IsPeriodic()) 1137 { 1138 OBUnitCell *box = (OBUnitCell*)GetParent()->GetData(OBGenericDataType::UnitCell); 1139 v1 = box->MinimumImageCartesian(v1); 1140 v2 = box->MinimumImageCartesian(v2); 1141 } 1142 1143 if (IsNearZero(v1.length(), 1.0e-3) 1144 || IsNearZero(v2.length(), 1.0e-3)) { 1145 return(0.0); 1146 } 1147 1148 return(vectorAngle(v1, v2)); 1149 } 1150 GetAngle(int b,int c)1151 double OBAtom::GetAngle(int b, int c) 1152 { 1153 OBMol *mol = (OBMol*)GetParent(); 1154 return(this->GetAngle(mol->GetAtom(b), mol->GetAtom(c))); 1155 } 1156 GetNewBondVector(vector3 & v,double length)1157 bool OBAtom::GetNewBondVector(vector3 &v,double length) 1158 { 1159 v = OBBuilder::GetNewBondVector(this, length); 1160 return(true); 1161 } 1162 1163 //! \return a "corrected" bonding radius based on the hybridization. 1164 //! Scales the covalent radius by 0.95 for sp2 and 0.90 for sp hybrids CorrectedBondRad(unsigned int elem,unsigned int hyb)1165 static double CorrectedBondRad(unsigned int elem, unsigned int hyb) 1166 { 1167 double rad = OBElements::GetCovalentRad(elem); 1168 switch (hyb) { 1169 case 2: 1170 return rad * 0.95; 1171 case 1: 1172 return rad * 0.90; 1173 default: 1174 return rad; 1175 } 1176 } 1177 HtoMethyl()1178 bool OBAtom::HtoMethyl() 1179 { 1180 if (GetAtomicNum() != OBElements::Hydrogen) 1181 return(false); 1182 1183 obErrorLog.ThrowError(__FUNCTION__, 1184 "Ran OpenBabel::HtoMethyl", obAuditMsg); 1185 1186 OBMol *mol = (OBMol*)GetParent(); 1187 1188 mol->BeginModify(); 1189 1190 SetAtomicNum(6); 1191 SetType("C3"); 1192 SetHyb(3); 1193 1194 OBAtom *atom; 1195 OBBond *bond; 1196 OBBondIterator i; 1197 atom = BeginNbrAtom(i); 1198 bond = (OBBond *)*i; 1199 if (!atom) 1200 { 1201 mol->EndModify(); 1202 return(false); 1203 } 1204 1205 double br1,br2; 1206 br1 = CorrectedBondRad(6,3); 1207 br2 = CorrectedBondRad(atom->GetAtomicNum(),atom->GetHyb()); 1208 bond->SetLength(atom,br1+br2); 1209 1210 OBAtom *hatom; 1211 br2 = CorrectedBondRad(1,0); 1212 vector3 v; 1213 1214 for (int j = 0;j < 3;++j) 1215 { 1216 hatom = mol->NewAtom(); 1217 hatom->SetAtomicNum(1); 1218 hatom->SetType("H"); 1219 1220 GetNewBondVector(v,br1+br2); 1221 hatom->SetVector(v); 1222 mol->AddBond(GetIdx(),mol->NumAtoms(),1); 1223 } 1224 1225 mol->EndModify(); 1226 return(true); 1227 } 1228 ApplyRotMatToBond(OBMol & mol,matrix3x3 & m,OBAtom * a1,OBAtom * a2)1229 static void ApplyRotMatToBond(OBMol &mol,matrix3x3 &m,OBAtom *a1,OBAtom *a2) 1230 { 1231 vector<int> children; 1232 mol.FindChildren(children,a1->GetIdx(),a2->GetIdx()); 1233 children.push_back(a2->GetIdx()); 1234 1235 vector3 v; 1236 vector<int>::iterator i; 1237 for (i = children.begin();i != children.end();++i) 1238 { 1239 v = mol.GetAtom(*i)->GetVector(); 1240 v -= a1->GetVector(); 1241 v *= m; 1242 v += a1->GetVector(); 1243 mol.GetAtom(*i)->SetVector(v); 1244 } 1245 1246 } 1247 1248 //! \deprecated This will be removed in future versions of Open Babel SetHybAndGeom(int hyb)1249 bool OBAtom::SetHybAndGeom(int hyb) 1250 { 1251 obErrorLog.ThrowError(__FUNCTION__, 1252 "Ran OpenBabel::SetHybridizationAndGeometry", 1253 obAuditMsg); 1254 1255 //if (hyb == GetHyb()) return(true); 1256 if (GetAtomicNum() == 1) 1257 return(false); 1258 if (hyb == 0 && GetHvyDegree() > 1) 1259 return(false); 1260 if (hyb == 1 && GetHvyDegree() > 2) 1261 return(false); 1262 if (hyb == 2 && GetHvyDegree() > 3) 1263 return(false); 1264 if (hyb == 3 && GetHvyDegree() > 4) 1265 return(false); 1266 1267 OBMol *mol = (OBMol*)GetParent(); 1268 1269 OBAtom *nbr; 1270 vector<OBAtom*> delatm; 1271 OBBondIterator i; 1272 1273 for (nbr = BeginNbrAtom(i);nbr;nbr = NextNbrAtom(i)) 1274 if (nbr->GetAtomicNum() == OBElements::Hydrogen) 1275 delatm.push_back(nbr); 1276 1277 //delete attached hydrogens 1278 mol->IncrementMod(); 1279 vector<OBAtom*>::iterator j; 1280 for (j = delatm.begin();j != delatm.end();++j) 1281 mol->DeleteAtom(*j); 1282 mol->DecrementMod(); 1283 1284 double targetAngle; 1285 if (hyb == 3) 1286 targetAngle = 109.5; 1287 else if (hyb == 2) 1288 targetAngle = 120.0; 1289 else if (hyb == 1) 1290 targetAngle = 180.0; 1291 else 1292 targetAngle = 0.0; 1293 1294 if (IsInRing()) 1295 targetAngle = 180.0 - (360.0 / MemberOfRingSize()); 1296 1297 //adjust attached acyclic bond lengths 1298 double br1,br2, length; 1299 br1 = CorrectedBondRad(GetAtomicNum(),hyb); 1300 for (nbr = BeginNbrAtom(i);nbr;nbr = NextNbrAtom(i)) 1301 if (!(*i)->IsInRing()) 1302 { 1303 br2 = CorrectedBondRad(nbr->GetAtomicNum(),nbr->GetHyb()); 1304 length = br1 + br2; 1305 if ((*i)->IsAromatic()) 1306 length *= 0.93; 1307 else if ((*i)->GetBondOrder() == 2) 1308 length *= 0.91; 1309 else if ((*i)->GetBondOrder() == 3) 1310 length *= 0.87; 1311 ((OBBond*) *i)->SetLength(this, length); 1312 } 1313 1314 if (GetExplicitDegree() > 1) 1315 { 1316 double angle; 1317 matrix3x3 m; 1318 vector3 v1,v2,v3,v4,n,s; 1319 OBAtom *r1,*r2,*r3,*a1,*a2,*a3,*a4; 1320 r1 = r2 = r3 = a1 = a2 = a3 = a4 = nullptr; 1321 1322 //find ring atoms first 1323 for (nbr = BeginNbrAtom(i);nbr;nbr = NextNbrAtom(i)) 1324 if ((*i)->IsInRing()) { 1325 if (!r1) { 1326 r1 = nbr; 1327 } 1328 else if (!r2) { 1329 r2 = nbr; 1330 } 1331 else if (!r3) { 1332 r3 = nbr; 1333 } 1334 } 1335 1336 //find non-ring atoms 1337 for (nbr = BeginNbrAtom(i);nbr;nbr = NextNbrAtom(i)) 1338 if (!(*i)->IsInRing()) { 1339 if (!a1) { 1340 a1 = nbr; 1341 } 1342 else if (!a2) { 1343 a2 = nbr; 1344 } 1345 else if (!a3) { 1346 a3 = nbr; 1347 } 1348 else if (!a4) { 1349 a4 = nbr; 1350 } 1351 } 1352 1353 //adjust geometries of heavy atoms according to hybridization 1354 if (hyb == 1) 1355 { 1356 if (a2) 1357 { 1358 v1 = a1->GetVector()-GetVector(); 1359 v1.normalize(); 1360 v2 = a2->GetVector()-GetVector(); 1361 v2.normalize(); 1362 n = cross(v1,v2); 1363 angle = vectorAngle(v1,v2)-targetAngle; 1364 m.RotAboutAxisByAngle(n,-angle); 1365 ApplyRotMatToBond(*mol,m,this,a1); 1366 } 1367 } 1368 else if (hyb == 2) 1369 { 1370 if (r1 && r2 && a1) 1371 { 1372 v1 = r1->GetVector()-GetVector(); 1373 v1.normalize(); 1374 v2 = r2->GetVector()-GetVector(); 1375 v2.normalize(); 1376 v3 = a1->GetVector()-GetVector(); 1377 s = v1+v2; 1378 s.normalize(); 1379 s *= -1.0; 1380 n = cross(s,v3); 1381 angle = vectorAngle(s,v3); 1382 m.RotAboutAxisByAngle(n,angle); 1383 ApplyRotMatToBond(*mol,m,this,a1); 1384 } 1385 else 1386 { 1387 if (a2) 1388 { 1389 v1 = a1->GetVector()-GetVector(); 1390 v1.normalize(); 1391 v2 = a2->GetVector()-GetVector(); 1392 v2.normalize(); 1393 n = cross(v1,v2); 1394 angle = vectorAngle(v1,v2)-targetAngle; 1395 m.RotAboutAxisByAngle(n,-angle); 1396 ApplyRotMatToBond(*mol,m,this,a1); 1397 } 1398 if (a3) 1399 { 1400 v1 = a1->GetVector()-GetVector(); 1401 v1.normalize(); 1402 v2 = a2->GetVector()-GetVector(); 1403 v2.normalize(); 1404 v3 = a3->GetVector()-GetVector(); 1405 s = v1+v2; 1406 s.normalize(); 1407 s *= -1.0; 1408 n = cross(s,v3); 1409 angle = vectorAngle(s,v3); 1410 m.RotAboutAxisByAngle(n,angle); 1411 ApplyRotMatToBond(*mol,m,this,a3); 1412 } 1413 } 1414 } 1415 else if (hyb == 3) 1416 { 1417 if (r1 && r2 && r3 && a1) 1418 { 1419 v1 = r1->GetVector()-GetVector(); 1420 v1.normalize(); 1421 v2 = r2->GetVector()-GetVector(); 1422 v2.normalize(); 1423 v3 = r3->GetVector()-GetVector(); 1424 v3.normalize(); 1425 v4 = a1->GetVector()-GetVector(); 1426 s = v1 + v2 + v3; 1427 s *= -1.0; 1428 s.normalize(); 1429 n = cross(s,v4); 1430 angle = vectorAngle(s,v4); 1431 m.RotAboutAxisByAngle(n,angle); 1432 ApplyRotMatToBond(*mol,m,this,a1); 1433 } 1434 else if (r1 && r2 && !r3 && a1) 1435 { 1436 v1 = r1->GetVector()-GetVector(); 1437 v1.normalize(); 1438 v2 = r2->GetVector()-GetVector(); 1439 v2.normalize(); 1440 v3 = a1->GetVector()-GetVector(); 1441 s = v1+v2; 1442 s *= -1.0; 1443 s.normalize(); 1444 n = cross(v1,v2); 1445 n.normalize(); 1446 #ifndef ONE_OVER_SQRT3 1447 #define ONE_OVER_SQRT3 0.57735026918962576451 1448 #endif //SQRT_TWO_THIRDS 1449 #ifndef SQRT_TWO_THIRDS 1450 #define SQRT_TWO_THIRDS 0.81649658092772603272 1451 #endif //ONE_OVER_SQRT3 1452 s *= ONE_OVER_SQRT3; //1/sqrt(3) 1453 n *= SQRT_TWO_THIRDS; //sqrt(2/3) 1454 s += n; 1455 s.normalize(); 1456 n = cross(s,v3); 1457 angle = vectorAngle(s,v3); 1458 m.RotAboutAxisByAngle(n,angle); 1459 ApplyRotMatToBond(*mol,m,this,a1); 1460 1461 if (a2) 1462 { 1463 v1 = r1->GetVector()-GetVector(); 1464 v1.normalize(); 1465 v2 = r2->GetVector()-GetVector(); 1466 v2.normalize(); 1467 v3 = a1->GetVector()-GetVector(); 1468 v3.normalize(); 1469 v4 = a2->GetVector()-GetVector(); 1470 s = v1 + v2 + v3; 1471 s *= -1.0; 1472 s.normalize(); 1473 n = cross(s,v4); 1474 angle = vectorAngle(s,v4); 1475 m.RotAboutAxisByAngle(n,angle); 1476 ApplyRotMatToBond(*mol,m,this,a2); 1477 } 1478 } 1479 else 1480 { 1481 if (a2) 1482 { 1483 v1 = a1->GetVector()-GetVector(); 1484 v1.normalize(); 1485 v2 = a2->GetVector()-GetVector(); 1486 v2.normalize(); 1487 n = cross(v1,v2); 1488 angle = vectorAngle(v1,v2)-targetAngle; 1489 m.RotAboutAxisByAngle(n,-angle); 1490 ApplyRotMatToBond(*mol,m,this,a1); 1491 } 1492 if (a3) 1493 { 1494 v1 = a1->GetVector()-GetVector(); 1495 v1.normalize(); 1496 v2 = a2->GetVector()-GetVector(); 1497 v2.normalize(); 1498 v3 = a3->GetVector()-GetVector(); 1499 s = v1+v2; 1500 s *= -1.0; 1501 s.normalize(); 1502 n = cross(v1,v2); 1503 n.normalize(); 1504 s *= ONE_OVER_SQRT3; //1/sqrt(3) 1505 n *= SQRT_TWO_THIRDS; //sqrt(2/3) 1506 s += n; 1507 s.normalize(); 1508 n = cross(s,v3); 1509 angle = vectorAngle(s,v3); 1510 m.RotAboutAxisByAngle(n,angle); 1511 ApplyRotMatToBond(*mol,m,this,a3); 1512 } 1513 } 1514 } 1515 } 1516 1517 //add hydrogens back to atom 1518 int impval=1; 1519 switch (GetAtomicNum()) 1520 { 1521 case 6: 1522 if (hyb == 3) 1523 impval = 4; 1524 if (hyb == 2) 1525 impval = 3; 1526 if (hyb == 1) 1527 impval = 2; 1528 break; 1529 case 7: 1530 if (hyb == 3) 1531 impval = 3; 1532 if (hyb == 2) 1533 impval = 2; 1534 if (hyb == 1) 1535 impval = 1; 1536 break; 1537 case 8: 1538 if (hyb == 3) 1539 impval = 2; 1540 if (hyb == 2) 1541 impval = 2; 1542 if (hyb == 1) 1543 impval = 2; 1544 break; 1545 case 16: 1546 if (hyb == 3) 1547 impval = 2; 1548 if (hyb == 2) 1549 impval = 2; 1550 if (hyb == 1) 1551 impval = 0; 1552 break; 1553 case 15: 1554 if (hyb == 3) 1555 impval = 4; 1556 if (hyb == 2) 1557 impval = 3; 1558 if (hyb == 1) 1559 impval = 2; 1560 break; 1561 default: 1562 impval = 1; 1563 } 1564 1565 int hcount = impval-GetHvyDegree(); 1566 if (hcount) 1567 { 1568 int k; 1569 vector3 v; 1570 OBAtom *atom; 1571 double brsum = CorrectedBondRad(1,0)+ 1572 CorrectedBondRad(GetAtomicNum(),GetHyb()); 1573 SetHyb(hyb); 1574 1575 mol->BeginModify(); 1576 for (k = 0;k < hcount;++k) 1577 { 1578 GetNewBondVector(v,brsum); 1579 atom = mol->NewAtom(); 1580 atom->SetAtomicNum(1); 1581 atom->SetType("H"); 1582 atom->SetVector(v); 1583 mol->AddBond(atom->GetIdx(),GetIdx(),1); 1584 } 1585 mol->EndModify(); 1586 } 1587 1588 1589 return(true); 1590 } 1591 GetBond(OBAtom * nbr)1592 OBBond *OBAtom::GetBond(OBAtom *nbr) 1593 { 1594 OBBond *bond; 1595 vector<OBBond *>::iterator i; 1596 for (bond = BeginBond(i) ; bond ; bond = NextBond(i)) 1597 if (bond->GetNbrAtom(this) == nbr) 1598 return bond; 1599 return nullptr; 1600 } 1601 1602 /*Now in OBBase 1603 // OBGenericData methods 1604 bool OBAtom::HasData(string &s) 1605 //returns true if the generic attribute/value pair exists 1606 { 1607 if (_vdata.empty()) 1608 return(false); 1609 1610 vector<OBGenericData*>::iterator i; 1611 1612 for (i = _vdata.begin();i != _vdata.end();++i) 1613 if ((*i)->GetAttribute() == s) 1614 return(true); 1615 1616 return(false); 1617 } 1618 1619 bool OBAtom::HasData(const char *s) 1620 //returns true if the generic attribute/value pair exists 1621 { 1622 if (_vdata.empty()) 1623 return(false); 1624 1625 vector<OBGenericData*>::iterator i; 1626 1627 for (i = _vdata.begin();i != _vdata.end();++i) 1628 if ((*i)->GetAttribute() == s) 1629 return(true); 1630 1631 return(false); 1632 } 1633 1634 bool OBAtom::HasData(unsigned int dt) 1635 //returns true if the generic attribute/value pair exists 1636 { 1637 if (_vdata.empty()) 1638 return(false); 1639 1640 vector<OBGenericData*>::iterator i; 1641 1642 for (i = _vdata.begin();i != _vdata.end();++i) 1643 if ((*i)->GetDataType() == dt) 1644 return(true); 1645 1646 return(false); 1647 } 1648 1649 OBGenericData *OBAtom::GetData(string &s) 1650 //returns the value given an attribute 1651 { 1652 vector<OBGenericData*>::iterator i; 1653 1654 for (i = _vdata.begin();i != _vdata.end();++i) 1655 if ((*i)->GetAttribute() == s) 1656 return(*i); 1657 1658 return(NULL); 1659 } 1660 1661 OBGenericData *OBAtom::GetData(const char *s) 1662 //returns the value given an attribute 1663 { 1664 vector<OBGenericData*>::iterator i; 1665 1666 for (i = _vdata.begin();i != _vdata.end();++i) 1667 if ((*i)->GetAttribute() == s) 1668 return(*i); 1669 1670 return(NULL); 1671 } 1672 1673 OBGenericData *OBAtom::GetData(unsigned int dt) 1674 { 1675 vector<OBGenericData*>::iterator i; 1676 for (i = _vdata.begin();i != _vdata.end();++i) 1677 if ((*i)->GetDataType() == dt) 1678 return(*i); 1679 return(NULL); 1680 } 1681 1682 void OBAtom::DeleteData(unsigned int dt) 1683 { 1684 vector<OBGenericData*> vdata; 1685 vector<OBGenericData*>::iterator i; 1686 for (i = _vdata.begin();i != _vdata.end();++i) 1687 if ((*i)->GetDataType() == dt) 1688 delete *i; 1689 else 1690 vdata.push_back(*i); 1691 _vdata = vdata; 1692 } 1693 1694 void OBAtom::DeleteData(vector<OBGenericData*> &vg) 1695 { 1696 vector<OBGenericData*> vdata; 1697 vector<OBGenericData*>::iterator i,j; 1698 1699 bool del; 1700 for (i = _vdata.begin();i != _vdata.end();++i) 1701 { 1702 del = false; 1703 for (j = vg.begin();j != vg.end();++j) 1704 if (*i == *j) 1705 { 1706 del = true; 1707 break; 1708 } 1709 if (del) 1710 delete *i; 1711 else 1712 vdata.push_back(*i); 1713 } 1714 _vdata = vdata; 1715 } 1716 1717 void OBAtom::DeleteData(OBGenericData *gd) 1718 { 1719 vector<OBGenericData*>::iterator i; 1720 for (i = _vdata.begin();i != _vdata.end();++i) 1721 if (*i == gd) 1722 { 1723 delete *i; 1724 _vdata.erase(i); 1725 } 1726 1727 } 1728 */ 1729 IsHbondAcceptorSimple()1730 bool OBAtom::IsHbondAcceptorSimple() 1731 { 1732 // Changes from Liu Zhiguo 1733 if (_ele == 8 || _ele == 9) 1734 return true; 1735 if (_ele == 7) { 1736 // N+ ions and sp2 hybrid N with 3 valences should not be Hbond acceptors 1737 if (!((GetExplicitDegree() == 4 && GetHyb() == 3) 1738 || (GetExplicitDegree() == 3 && GetHyb() == 2))) 1739 return true; 1740 } 1741 // Changes from Paolo Tosco 1742 if (_ele == 16 && GetFormalCharge() == -1) 1743 return true; 1744 return false; 1745 } 1746 1747 // new function, Stefano Forli 1748 // Incorporate ideas and data from Kubyni and others. 1749 // [1] Kubinyi, H. "Changing paradigms in drug discovery. 1750 // "SPECIAL PUBLICATION-ROYAL SOCIETY OF CHEMISTRY 304.1 (2006): 219-232. 1751 // 1752 // [2] Kingsbury, Charles A. "Why are the Nitro and Sulfone 1753 // Groups Poor Hydrogen Bonders?." (2015). 1754 // 1755 // [3] Per Restorp, Orion B. Berryman, Aaron C. Sather, Dariush Ajami 1756 // and Julius Rebek Jr., Chem. Commun., 2009, 5692 DOI: 10.1039/b914171e 1757 // 1758 // [4] Dunitz, Taylor. "Organic fluorine hardly ever accepts 1759 // hydrogen bonds." Chemistry-A European Journal 3.1 (1997): 83-92. 1760 // 1761 // This function has a finer grain than the original 1762 // implementation, checking also the neighbor atoms. 1763 // Accordingly to these rules, the function will return: 1764 // 1765 // aliph-O-aliph ether -> true [1] 1766 // hydroxy O-sp3 -> true [1] 1767 // aro-O-aliph ether -> true [1] 1768 // ester O-sp2 -> true [1] 1769 // sulfate O (R-SO3) -> true [2] 1770 // sulfoxyde O (R-SO-R) -> true [2] 1771 // organoboron-F (R-BF3) -> true [3] 1772 // ester O-sp3 -> false [1] 1773 // sulfone (R1-SO2-R2 ) -> false [2] 1774 // aro-O-aro -> false [1] 1775 // aromatic O -> false [1] 1776 // O-nitro -> false [2] 1777 // organic F (R-F) -> false [4] 1778 // IsHbondAcceptor()1779 bool OBAtom::IsHbondAcceptor() { 1780 if (_ele == 8) { 1781 // oxygen; this should likely be a separate function 1782 // something like IsHbondAcceptorOxygen() 1783 unsigned int aroCount = 0; 1784 1785 OBBond *bond; 1786 OBBondIterator i; 1787 if (IsNitroOxygen()){ // maybe could be a bool option in the function? 1788 return (false); 1789 } 1790 if (IsAromatic()){ // aromatic oxygen (furan) (NO) 1791 return(false); 1792 } 1793 if (IsSulfoneOxygen(this)){ // sulfone (NO) 1794 return(false); 1795 } 1796 FOR_NBORS_OF_ATOM(nbr, this){ 1797 if (nbr->IsAromatic()){ 1798 aroCount += 1; 1799 if (aroCount == 2){ // aromatic ether (aro-O-aro) (NO) 1800 return(false); 1801 } 1802 } 1803 else { 1804 if (nbr->GetAtomicNum() == OBElements::Hydrogen) { // hydroxyl (YES) 1805 return(true); 1806 } 1807 else { 1808 bond = nbr->GetBond(this); 1809 if ( (bond->IsEster()) && (!(IsCarboxylOxygen()))) 1810 return(false); 1811 } 1812 } 1813 } 1814 return(true); // any other oxygen 1815 } // oxygen END 1816 // fluorine 1817 if (_ele == 9 ) { 1818 OBBondIterator i; 1819 // organic fluorine (NO) 1820 for (OBAtom* nbr = BeginNbrAtom(i);nbr;nbr=NextNbrAtom(i)) 1821 if (nbr->GetAtomicNum() == 6) 1822 return (false); 1823 else 1824 return (true); 1825 }; 1826 if (_ele == 7) { 1827 // N+ ions and sp2 hybrid N with 3 valences should not be Hbond acceptors 1828 if (!((GetExplicitDegree() == 4 && GetHyb() == 3) 1829 || (GetExplicitDegree() == 3 && GetHyb() == 2))) 1830 return true; 1831 }; 1832 // Changes from Paolo Tosco 1833 if (_ele == 16 && GetFormalCharge() == -1){ 1834 return (true); } 1835 // everything else 1836 return (false); 1837 } 1838 IsHbondDonor()1839 bool OBAtom::IsHbondDonor() 1840 { 1841 // return MatchesSMARTS("[$([#8,#7H,#9;!H0])]"); 1842 if (!(_ele == 7 || _ele == 8 || _ele == 9)) 1843 return false; 1844 1845 OBAtom *nbr; 1846 OBBondIterator i; 1847 for (nbr = BeginNbrAtom(i);nbr;nbr = NextNbrAtom(i)) 1848 if (nbr->GetAtomicNum() == OBElements::Hydrogen) 1849 return true; 1850 1851 return false; 1852 } 1853 IsHbondDonorH()1854 bool OBAtom::IsHbondDonorH() 1855 { 1856 if (GetAtomicNum() != OBElements::Hydrogen) return(false); 1857 1858 // Recursive definition -- this atom is an H atom 1859 // are any neighbors HbondDonors? 1860 1861 OBAtom *atom; 1862 OBBond *bond; 1863 OBBondIterator i; 1864 for (bond = BeginBond(i);bond;bond = NextBond(i)) { 1865 atom = bond->GetNbrAtom(this); 1866 if (atom->IsHbondDonor()) return(true); 1867 } 1868 1869 return(false); 1870 } 1871 IsMetal()1872 bool OBAtom::IsMetal() 1873 { 1874 const unsigned NMETALS = 78; 1875 const int metals[NMETALS] = { 1876 3,4,11,12,13,19,20,21,22,23,24,25,26,27,28,29, 1877 30,31,37,38,39,40,41,42,43,44,45,46,47,48,49,50,55,56,57,58,59,60,61,62,63, 1878 64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,87,88,89,90,91, 1879 92,93,94,95,96,97,98,99,100,101,102,103}; 1880 return std::find(metals, metals+78, GetAtomicNum())!=metals+78; 1881 } 1882 1883 } // end namespace OpenBabel 1884 1885 //! \file atom.cpp 1886 //! \brief Handle OBAtom class 1887