1 package Molecule;
   2 #
   3 # File: Molecule.pm
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2024 Manish Sud. All rights reserved.
   7 #
   8 # This file is part of MayaChemTools.
   9 #
  10 # MayaChemTools is free software; you can redistribute it and/or modify it under
  11 # the terms of the GNU Lesser General Public License as published by the Free
  12 # Software Foundation; either version 3 of the License, or (at your option) any
  13 # later version.
  14 #
  15 # MayaChemTools is distributed in the hope that it will be useful, but without
  16 # any warranty; without even the implied warranty of merchantability of fitness
  17 # for a particular purpose.  See the GNU Lesser General Public License for more
  18 # details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public License
  21 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  22 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  23 # Boston, MA, 02111-1307, USA.
  24 #
  25 
  26 use strict;
  27 use Carp;
  28 use Exporter;
  29 use Storable ();
  30 use Scalar::Util ();
  31 use ObjectProperty;
  32 use MathUtil;
  33 use PeriodicTable;
  34 use Text::ParseWords;
  35 use TextUtil;
  36 use FileUtil;
  37 use Graph;
  38 use Atom;
  39 use Bond;
  40 use MolecularFormula;
  41 
  42 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  43 
  44 @ISA = qw(Graph ObjectProperty Exporter);
  45 @EXPORT = qw(IsMolecule);
  46 @EXPORT_OK = qw(FormatElementalCompositionInformation GetSupportedAromaticityModels IsSupportedAromaticityModel);
  47 
  48 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
  49 
  50 # Setup class variables...
  51 my($ClassName, $ObjectID, %AromaticityModelsDataMap, %CanonicalAromaticityModelNamesMap);
  52 _InitializeClass();
  53 
  54 # Overload Perl functions...
  55 use overload '""' => 'StringifyMolecule';
  56 
  57 # Class constructor...
  58 sub new {
  59   my($Class, %NamesAndValues) = @_;
  60 
  61   # Initialize object...
  62   my $This = $Class->SUPER::new();
  63   bless $This, ref($Class) || $Class;
  64   $This->_InitializeMolecule();
  65 
  66   if (keys %NamesAndValues) { $This->_InitializeMoleculeProperties(%NamesAndValues); }
  67 
  68   return $This;
  69 }
  70 
  71 # Initialize object data...
  72 sub _InitializeMolecule {
  73   my($This) = @_;
  74   my($ObjectID) = _GetNewObjectID();
  75 
  76   # All other property names and values along with all Set/Get<PropertyName> methods
  77   # are implemented on-demand using ObjectProperty class.
  78   $This->{ID} = $ObjectID;
  79   $This->{Name} = "Molecule ${ObjectID}";
  80 }
  81 
  82 # Initialize molecule properties...
  83 sub _InitializeMoleculeProperties {
  84   my($This, %NamesAndValues) = @_;
  85 
  86   my($Name, $Value, $MethodName);
  87   while (($Name, $Value) = each  %NamesAndValues) {
  88     $MethodName = "Set${Name}";
  89     $This->$MethodName($Value);
  90   }
  91 }
  92 
  93 # Initialize class ...
  94 sub _InitializeClass {
  95   #Class name...
  96   $ClassName = __PACKAGE__;
  97 
  98   # ID to keep track of objects...
  99   $ObjectID = 0;
 100 
 101   # Load molecule class data...
 102   _LoadMoleculeClassData();
 103 }
 104 
 105 # Setup an explicit SetID method to block setting of ID by AUTOLOAD function...
 106 sub SetID {
 107   my($This, $Value) = @_;
 108 
 109   carp "Warning: ${ClassName}->SetID: Object ID can't be changed: it's used for internal tracking...";
 110 
 111   return $This;
 112 }
 113 
 114 # Add an atom...
 115 sub AddAtom {
 116   my($This, $Atom) = @_;
 117 
 118   if (!defined $Atom) {
 119     carp "Warning: ${ClassName}->AddAtom: No atom added: Atom must be specified...";
 120     return undef;
 121   }
 122   if ($This->HasAtom($Atom)) {
 123     carp "Warning: ${ClassName}->AddAtom: No atom added: Atom already exists...";
 124     return undef;
 125   }
 126   return $This->_AddAtom($Atom);
 127 }
 128 
 129 # Add an atom...
 130 sub _AddAtom {
 131   my($This, $Atom) = @_;
 132 
 133   # Assign atom to this molecule...
 134   $Atom->_SetMolecule($This);
 135 
 136   # Add it to the graph as a vertex...
 137   my($AtomID);
 138   $AtomID = $Atom->GetID();
 139   $This->AddVertex($AtomID);
 140   $This->SetVertexProperty('Atom', $Atom, $AtomID);
 141 
 142   return $This;
 143 }
 144 
 145 # Add atoms...
 146 sub AddAtoms {
 147   my($This, @Atoms) = @_;
 148 
 149   if (!@Atoms) {
 150     carp "Warning: ${ClassName}->AddAtoms: No atoms added: Atoms list must be specified...";
 151     return undef;
 152   }
 153   my($Atom);
 154   for $Atom (@Atoms) {
 155     $This->AddAtom($Atom);
 156   }
 157   return $This;
 158 }
 159 
 160 # Create an atom and add it to molecule...
 161 sub NewAtom {
 162   my($This, %NamesAndValues) = @_;
 163   my($Atom);
 164 
 165   $Atom = new Atom(%NamesAndValues);
 166   $This->AddAtom($Atom);
 167 
 168   return $Atom;
 169 }
 170 
 171 # Delete an atom...
 172 sub DeleteAtom {
 173   my($This, $Atom) = @_;
 174 
 175   if (!defined $Atom) {
 176     carp "Warning: ${ClassName}->DeleteAtom: No atom deleted: Atom must be specified...";
 177     return undef;
 178   }
 179   # Does the atom exist in  molecule?
 180   if (!$This->HasAtom($Atom)) {
 181     carp "Warning: ${ClassName}->DeleteAtom: No atom deleted: Atom doesn't exist...";
 182     return undef;
 183   }
 184   return $This->_DeleteAtom($Atom);
 185 }
 186 
 187 # Delete atom...
 188 sub _DeleteAtom {
 189   my($This, $Atom) = @_;
 190 
 191   my($AtomID);
 192   $AtomID = $Atom->GetID();
 193   $This->DeleteVertex($AtomID);
 194 
 195   return $This;
 196 }
 197 
 198 # Delete atoms...
 199 sub DeleteAtoms {
 200   my($This, @Atoms) = @_;
 201 
 202   if (!@Atoms) {
 203     carp "Warning: ${ClassName}->DeleteAtoms: No atoms added: Atoms list must be specified...";
 204     return undef;
 205   }
 206   my($Atom);
 207   for $Atom (@Atoms) {
 208     $This->DeleteAtom($Atom);
 209   }
 210 
 211   return $This;
 212 }
 213 
 214 # Is this atom present?
 215 sub HasAtom {
 216   my($This, $Atom) = @_;
 217 
 218   if (!defined $Atom) {
 219     return 0;
 220   }
 221   if (!$Atom->HasProperty('Molecule')) {
 222     # It's not in molecule...
 223     return 0;
 224   }
 225   my($AtomID);
 226   $AtomID = $Atom->GetID();
 227   if (!$This->HasVertex($AtomID)) {
 228     # It's not in molecule...
 229     return 0;
 230   }
 231   my($Molecule);
 232   $Molecule = $Atom->GetProperty('Molecule');
 233 
 234   return ($This->HasVertex($AtomID) && $This->GetID() == $Molecule->GetID()) ? 1 : 0;
 235 }
 236 
 237 # Return an array of atoms. In scalar context,  it returns number of atoms. Additionally,
 238 # atoms array can be filtered by any user specifiable Atom class method...
 239 #
 240 sub GetAtoms {
 241   my($This, $AtomCheckMethodName, $NegateMethodResult) = @_;
 242   my(@Atoms, @AtomIDs);
 243 
 244   @Atoms = (); @AtomIDs = ();
 245 
 246   @AtomIDs = $This->GetVertices();
 247   if (!@AtomIDs) {
 248     return wantarray ? @Atoms : scalar @Atoms;
 249   }
 250 
 251   @Atoms = $This->GetVerticesProperty('Atom', @AtomIDs);
 252 
 253   if (!defined $AtomCheckMethodName) {
 254     return wantarray ? @Atoms : scalar @Atoms;
 255   }
 256   $NegateMethodResult = (defined($NegateMethodResult) &&  $NegateMethodResult) ? 1 : 0;
 257 
 258   # Filter out atoms...
 259   my($Atom, $KeepAtom, @FilteredAtoms);
 260   @FilteredAtoms = ();
 261   for $Atom (@Atoms) {
 262     $KeepAtom = $NegateMethodResult ? (!$Atom->$AtomCheckMethodName()) : $Atom->$AtomCheckMethodName();
 263     if ($KeepAtom) {
 264       push @FilteredAtoms, $Atom;
 265     }
 266   }
 267   return wantarray ? @FilteredAtoms : scalar @FilteredAtoms;
 268 }
 269 
 270 # Return an array of bonds. In scalar context, it returns number of bonds...
 271 sub GetBonds {
 272   my($This) = @_;
 273   my(@Bonds, @EdgesAtomsIDs);
 274 
 275   @Bonds = (); @EdgesAtomsIDs = ();
 276 
 277   @EdgesAtomsIDs = $This->GetEdges();
 278   if (@EdgesAtomsIDs) {
 279     @Bonds = $This->GetEdgesProperty('Bond', @EdgesAtomsIDs);
 280   }
 281   return wantarray ? @Bonds : scalar @Bonds;
 282 }
 283 
 284 # Get number of atoms in molecule...
 285 sub GetNumOfAtoms {
 286   my($This) = @_;
 287   my($NumOfAtoms);
 288 
 289   $NumOfAtoms = $This->GetAtoms();
 290 
 291   return $NumOfAtoms;
 292 }
 293 
 294 # Get number of bonds in molecule...
 295 sub GetNumOfBonds {
 296   my($This) = @_;
 297   my($NumOfBonds);
 298 
 299   $NumOfBonds = $This->GetBonds();
 300 
 301   return $NumOfBonds;
 302 }
 303 
 304 # Get number of heavy atoms in molecule...
 305 sub GetNumOfHeavyAtoms {
 306   my($This) = @_;
 307 
 308   return $This->GetNumOfNonHydrogenAtoms();
 309 }
 310 
 311 # Get number of non-hydrogen atoms in molecule...
 312 sub GetNumOfNonHydrogenAtoms {
 313   my($This) = @_;
 314   my($NumOfNonHydrogenAtoms, $Atom, @Atoms);
 315 
 316   @Atoms = $This->GetAtoms();
 317   $NumOfNonHydrogenAtoms = 0;
 318   for $Atom (@Atoms) {
 319     if (!$Atom->IsHydrogen()) {
 320       $NumOfNonHydrogenAtoms++;
 321     }
 322   }
 323   return $NumOfNonHydrogenAtoms;
 324 }
 325 
 326 # Get number of hydrogen atoms in molecule...
 327 sub GetNumOfHydrogenAtoms {
 328   my($This) = @_;
 329   my($NumOfHydrogenAtoms, $Atom, @Atoms);
 330 
 331   @Atoms = $This->GetAtoms();
 332   $NumOfHydrogenAtoms = 0;
 333   for $Atom (@Atoms) {
 334     if ($Atom->IsHydrogen()) {
 335       $NumOfHydrogenAtoms++;
 336     }
 337   }
 338   return $NumOfHydrogenAtoms;
 339 }
 340 
 341 # Get number of missing hydrogen atoms in molecule...
 342 sub GetNumOfMissingHydrogenAtoms {
 343   my($This) = @_;
 344   my($NumOfMissingHydrogenAtoms, $Atom, @Atoms);
 345 
 346   @Atoms = $This->GetAtoms();
 347   $NumOfMissingHydrogenAtoms = 0;
 348   for $Atom (@Atoms) {
 349     if (!$Atom->IsHydrogen()) {
 350       $NumOfMissingHydrogenAtoms += $Atom->GetNumOfMissingHydrogens();
 351     }
 352   }
 353   return $NumOfMissingHydrogenAtoms;
 354 }
 355 
 356 # Add bond...
 357 sub AddBond {
 358   my($This, $Bond) = @_;
 359   my($Atom1, $Atom2);
 360 
 361   ($Atom1, $Atom2) = $Bond->GetAtoms();
 362   if (!(defined($Atom1) && defined($Atom2))) {
 363     carp "Warning: ${ClassName}->AddBond: No bond added: Both atoms must be specified...";
 364     return undef;
 365   }
 366   if (!($This->HasAtom($Atom1) && $This->HasAtom($Atom2))) {
 367     carp "Warning: ${ClassName}->AddBond: No bond added: Both atoms must be present...";
 368     return undef;
 369   }
 370   if ($This->HasBond($Bond)) {
 371     carp "Warning: ${ClassName}->AddBond: No bond added: Bond already exists...";
 372     return undef;
 373   }
 374   return $This->_AddBond($Bond);
 375 }
 376 
 377 # Add bond...
 378 sub _AddBond {
 379   my($This, $Bond) = @_;
 380 
 381   # Assign bond to this molecule...
 382   $Bond->_SetMolecule($This);
 383 
 384   # Add it to the graph as an edge...
 385   my($Atom1, $Atom2, $AtomID1, $AtomID2);
 386   ($Atom1, $Atom2) = $Bond->GetAtoms();
 387   $AtomID1 = $Atom1->GetID(); $AtomID2 = $Atom2->GetID();
 388   $This->AddEdge($AtomID1, $AtomID2);
 389   $This->SetEdgeProperty('Bond', $Bond, $AtomID1, $AtomID2);
 390 
 391   return $This;
 392 }
 393 
 394 # Add Bonds...
 395 sub AddBonds {
 396   my($This, @Bonds) = @_;
 397 
 398   if (!@Bonds) {
 399     carp "Warning: ${ClassName}->AddBonds: No bonds added: Bonds list must be specified...";
 400     return undef;
 401   }
 402   my($Bond);
 403   for $Bond (@Bonds) {
 404     $This->AddBond($Bond);
 405   }
 406   return $This;
 407 }
 408 
 409 # Create a bond and add it to molecule...
 410 sub NewBond {
 411   my($This, %NamesAndValues) = @_;
 412   my($Bond);
 413 
 414   $Bond = new Bond(%NamesAndValues);
 415   $This->AddBond($Bond);
 416 
 417   return $Bond;
 418 }
 419 
 420 # Delete a bond...
 421 sub DeleteBond {
 422   my($This, $Bond) = @_;
 423 
 424   if (!defined $Bond) {
 425     carp "Warning: ${ClassName}->DeleteBond: No bond deleted: Bond must be specified...";
 426     return undef;
 427   }
 428   # Does the bond exist in molecule?
 429   if (!$This->HasBond($Bond)) {
 430     carp "Warning: ${ClassName}->DeleteBond: No bond deleted: Bond doesn't exist...";
 431     return undef;
 432   }
 433   return $This->_DeleteBond($Bond);
 434 }
 435 
 436 # Delete bond...
 437 sub _DeleteBond {
 438   my($This, $Bond) = @_;
 439 
 440   my($Atom1, $Atom2, $AtomID1, $AtomID2);
 441   ($Atom1, $Atom2) = $Bond->GetAtoms();
 442   $AtomID1 = $Atom1->GetID(); $AtomID2 = $Atom2->GetID();
 443   $This->DeleteEdge($AtomID1, $AtomID2);
 444 
 445   return $This;
 446 }
 447 
 448 # Delete bonds...
 449 sub DeleteBonds {
 450   my($This, @Bonds) = @_;
 451 
 452   if (!@Bonds) {
 453     carp "Warning: ${ClassName}->DeleteBonds: No bonds added: Bonds list must be specified...";
 454     return undef;
 455   }
 456   my($Bond);
 457   for $Bond (@Bonds) {
 458     $This->DeleteBond($Bond);
 459   }
 460 
 461   return $This;
 462 }
 463 
 464 # Has bond...
 465 sub HasBond {
 466   my($This, $Bond) = @_;
 467   my($Atom1, $Atom2);
 468 
 469   ($Atom1, $Atom2) = $Bond->GetAtoms();
 470   if (!($This->HasAtom($Atom1) && $This->HasAtom($Atom2))) {
 471     return 0;
 472   }
 473   if (!$Bond->HasProperty('Molecule')) {
 474     # It's not in molecule...
 475     return 0;
 476   }
 477   my($AtomID1, $AtomID2, $Molecule);
 478   $AtomID1 = $Atom1->GetID(); $AtomID2 = $Atom2->GetID();
 479   $Molecule = $Bond->GetMolecule();
 480 
 481   return ($This->HasEdge($AtomID1, $AtomID2) && $This->GetID() == $Molecule->GetID()) ? 1 : 0;
 482 }
 483 
 484 # Get atom neighbors...
 485 sub _GetAtomNeighbors {
 486   my($This, $Atom) = @_;
 487 
 488   my($AtomID, @Atoms, @AtomIDs);
 489 
 490   @Atoms = (); @AtomIDs = ();
 491   $AtomID = $Atom->GetID();
 492   @AtomIDs = $This->GetNeighbors($AtomID);
 493   if (@AtomIDs) {
 494     @Atoms = $This->GetVerticesProperty('Atom', @AtomIDs);
 495   }
 496   return wantarray ? @Atoms : scalar @Atoms;
 497 }
 498 
 499 # Get atom bonds...
 500 sub _GetAtomBonds {
 501   my($This, $Atom) = @_;
 502   my($AtomID, @AtomIDs, @Bonds);
 503 
 504   @Bonds = (); @AtomIDs = ();
 505   $AtomID = $Atom->GetID();
 506   @AtomIDs = $This->GetEdges($AtomID);
 507   if (@AtomIDs) {
 508     @Bonds = $This->GetEdgesProperty('Bond', @AtomIDs);
 509   }
 510   return wantarray ? @Bonds : scalar @Bonds;
 511 }
 512 
 513 # Get bond to specified atom...
 514 sub _GetBondToAtom {
 515   my($This, $Atom1, $Atom2) = @_;
 516   my($AtomID1, $AtomID2);
 517 
 518   $AtomID1 = $Atom1->GetID();
 519   $AtomID2 = $Atom2->GetID();
 520 
 521   return $This->GetEdgeProperty('Bond', $AtomID1, $AtomID2);
 522 }
 523 
 524 # Are two specified atoms bonded?
 525 sub _IsBondedToAtom {
 526   my($This, $Atom1, $Atom2) = @_;
 527   my($AtomID1, $AtomID2);
 528 
 529   $AtomID1 = $Atom1->GetID();
 530   $AtomID2 = $Atom2->GetID();
 531 
 532   return $This->HasEdgeProperty('Bond', $AtomID1, $AtomID2);
 533 }
 534 
 535 # Add hydrogens to each atoms in molecule and return total number of hydrogens added...
 536 sub AddHydrogens {
 537   my($This) = @_;
 538 
 539   return $This->_AddHydrogens();
 540 }
 541 
 542 # Add hydrogens to polar atoms (N, O, P, S) in molecule and return total number of hydrogens added...
 543 sub AddPolarHydrogens {
 544   my($This) = @_;
 545   my($PolarHydrogensOnly) = 1;
 546 
 547   return $This->_AddHydrogens($PolarHydrogensOnly);
 548 }
 549 
 550 # Add all the hydrogens or hydrogens for polar atoms only...
 551 #
 552 # Note:
 553 #   . The current release of MayaChemTools doesn't assign any hydrogen positions.
 554 #
 555 sub _AddHydrogens {
 556   my($This, $PolarHydrogensOnly) = @_;
 557   my($Atom, $NumOfHydrogensAdded, $HydrogenPositionsWarning, @Atoms);
 558 
 559   if (! defined $PolarHydrogensOnly) {
 560     $PolarHydrogensOnly = 0;
 561   }
 562 
 563   $NumOfHydrogensAdded = 0;
 564   @Atoms = $This->GetAtoms();
 565   $HydrogenPositionsWarning = 0;
 566 
 567   ATOM: for $Atom (@Atoms) {
 568     if ($PolarHydrogensOnly) {
 569       if (!$Atom->IsPolarAtom()) {
 570         next ATOM;
 571       }
 572     }
 573     $NumOfHydrogensAdded += $Atom->AddHydrogens($HydrogenPositionsWarning);
 574   }
 575   return $NumOfHydrogensAdded;
 576 }
 577 
 578 # Delete all hydrogens atoms in molecule and return total number of hydrogens removed...
 579 sub DeleteHydrogens {
 580   my($This) = @_;
 581 
 582   return $This->_DeleteHydrogens();
 583 }
 584 
 585 # Delete hydrogens to polar atoms (N, O, P, S) in molecule and return total number of hydrogens removed...
 586 sub DeletePolarHydrogens {
 587   my($This) = @_;
 588   my($PolarHydrogensOnly) = 1;
 589 
 590   return $This->_DeleteHydrogens($PolarHydrogensOnly);
 591 }
 592 
 593 # Delete all hydrogens atoms in molecule and return total number of hydrogens removed...
 594 sub _DeleteHydrogens {
 595   my($This, $PolarHydrogensOnly) = @_;
 596   my($Atom, $NumOfHydrogensRemoved, @Atoms);
 597 
 598   if (! defined $PolarHydrogensOnly) {
 599     $PolarHydrogensOnly = 0;
 600   }
 601 
 602   $NumOfHydrogensRemoved = 0;
 603   @Atoms = $This->GetAtoms();
 604 
 605   ATOM: for $Atom (@Atoms) {
 606     if ($PolarHydrogensOnly) {
 607       if (!$Atom->IsPolarHydrogen()) {
 608         next ATOM;
 609       }
 610     }
 611     elsif (!$Atom->IsHydrogen()) {
 612       next ATOM;
 613     }
 614     $This->_DeleteAtom($Atom);
 615     $NumOfHydrogensRemoved++;
 616   }
 617   return $NumOfHydrogensRemoved;
 618 }
 619 
 620 # Get molecular weight by summing up atomic weights of all the atoms...
 621 sub GetMolecularWeight {
 622   my($This, $IncludeMissingHydrogens) = @_;
 623   my($MolecularWeight, $AtomicWeight, @Atoms, $Atom);
 624 
 625   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 1;
 626 
 627   $MolecularWeight = 0;
 628   @Atoms = $This->GetAtoms();
 629   for $Atom (@Atoms) {
 630     $AtomicWeight = $Atom->GetAtomicWeight();
 631     if (defined $AtomicWeight) {
 632       $MolecularWeight += $AtomicWeight;
 633     }
 634   }
 635 
 636   if (!$IncludeMissingHydrogens) {
 637     return $MolecularWeight;
 638   }
 639 
 640   # Account for missing hydrogen atoms...
 641   my($NumOfMissingHydrogenAtoms);
 642 
 643   $NumOfMissingHydrogenAtoms = $This->GetNumOfMissingHydrogenAtoms();
 644   if ($NumOfMissingHydrogenAtoms) {
 645     $MolecularWeight += $NumOfMissingHydrogenAtoms * PeriodicTable::GetElementAtomicWeight('H');
 646   }
 647 
 648   return $MolecularWeight;
 649 }
 650 
 651 # Get exact mass by summing up exact masses of all the atoms...
 652 sub GetExactMass {
 653   my($This, $IncludeMissingHydrogens) = @_;
 654   my($ExactMass, $AtomicMass, @Atoms, $Atom);
 655 
 656   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 1;
 657 
 658   $ExactMass = 0;
 659   @Atoms = $This->GetAtoms();
 660   for $Atom (@Atoms) {
 661     $AtomicMass = $Atom->GetExactMass();
 662     if (defined $AtomicMass) {
 663       $ExactMass += $AtomicMass;
 664     }
 665   }
 666 
 667   if (!$IncludeMissingHydrogens) {
 668     return $ExactMass;
 669   }
 670 
 671   # Account for missing hydrogen atoms...
 672   my($NumOfMissingHydrogenAtoms);
 673 
 674   $NumOfMissingHydrogenAtoms = $This->GetNumOfMissingHydrogenAtoms();
 675   if ($NumOfMissingHydrogenAtoms) {
 676     $ExactMass += $NumOfMissingHydrogenAtoms * PeriodicTable::GetElementMostAbundantNaturalIsotopeMass('H');
 677   }
 678 
 679   return $ExactMass;
 680 }
 681 
 682 # Get net formal charge on the molecule using one of the following two methods:
 683 #   . Using explicitly set FormalCharge property
 684 #   . Adding up formal charge on each atom in the molecule
 685 #
 686 # Caveats:
 687 #   . FormalCharge is different from Charge property of the molecule which corresponds to
 688 #     sum of partial atomic charges explicitly set for each atom using a specific methodology.
 689 #
 690 sub GetFormalCharge {
 691   my($This) = @_;
 692 
 693   # Is FormalCharge property explicitly set?
 694   if ($This->HasProperty('FormalCharge')) {
 695     return $This->GetProperty('FormalCharge');
 696   }
 697   my($FormalCharge, $AtomicFormalCharge, @Atoms, $Atom);
 698 
 699   $FormalCharge = 0;
 700   @Atoms = $This->GetAtoms();
 701   for $Atom (@Atoms) {
 702     $AtomicFormalCharge = $Atom->GetFormalCharge();
 703     if (defined $AtomicFormalCharge) {
 704       $FormalCharge += $AtomicFormalCharge;
 705     }
 706   }
 707   return $FormalCharge;
 708 }
 709 
 710 # Get net charge on the molecule using one of the following two methods:
 711 #   . Using explicitly set Charge property
 712 #   . Adding up charge on each atom in the molecule
 713 #
 714 # Caveats:
 715 #   . FormalCharge is different from Charge property of the molecule which corresponds to
 716 #     sum of partial atomic charges explicitly set for each atom using a specific methodology.
 717 #
 718 sub GetCharge {
 719   my($This) = @_;
 720 
 721   # Is Charge property explicitly set?
 722   if ($This->HasProperty('Charge')) {
 723     return $This->GetProperty('Charge');
 724   }
 725   my($Charge, $AtomicCharge, @Atoms, $Atom);
 726 
 727   $Charge = 0;
 728   @Atoms = $This->GetAtoms();
 729   for $Atom (@Atoms) {
 730     $AtomicCharge = $Atom->GetCharge();
 731     if (defined $AtomicCharge) {
 732       $Charge += $AtomicCharge;
 733     }
 734   }
 735   return $Charge;
 736 }
 737 
 738 # Get total SpinMultiplicity for the molecule using one of the following two methods:
 739 #   . Using explicitly set SpinMultiplicity property
 740 #   . Adding up SpinMultiplicity on each atom in the molecule
 741 #
 742 #
 743 sub GetSpinMultiplicity {
 744   my($This) = @_;
 745 
 746   # Is SpinMultiplicity property explicitly set?
 747   if ($This->HasProperty('SpinMultiplicity')) {
 748     return $This->GetProperty('SpinMultiplicity');
 749   }
 750   my($AtomicSpinMultiplicity, $SpinMultiplicity, @Atoms, $Atom);
 751 
 752   $SpinMultiplicity = 0;
 753   @Atoms = $This->GetAtoms();
 754   for $Atom (@Atoms) {
 755     $AtomicSpinMultiplicity = $Atom->GetSpinMultiplicity();
 756     if (defined $AtomicSpinMultiplicity) {
 757       $SpinMultiplicity += $AtomicSpinMultiplicity;
 758     }
 759   }
 760   return $SpinMultiplicity;
 761 }
 762 
 763 # Get total FreeRadicalElectrons for the molecule using one of the following two methods:
 764 #   . Using explicitly set FreeRadicalElectrons property
 765 #   . Adding up FreeRadicalElectrons on each atom in the molecule
 766 #
 767 #
 768 sub GetFreeRadicalElectrons {
 769   my($This) = @_;
 770 
 771   # Is FreeRadicalElectrons property explicitly set?
 772   if ($This->HasProperty('FreeRadicalElectrons')) {
 773     return $This->GetProperty('FreeRadicalElectrons');
 774   }
 775   my($AtomicFreeRadicalElectrons, $FreeRadicalElectrons, @Atoms, $Atom);
 776 
 777   $FreeRadicalElectrons = 0;
 778   @Atoms = $This->GetAtoms();
 779   for $Atom (@Atoms) {
 780     $AtomicFreeRadicalElectrons = $Atom->GetFreeRadicalElectrons();
 781     if (defined $AtomicFreeRadicalElectrons) {
 782       $FreeRadicalElectrons += $AtomicFreeRadicalElectrons;
 783     }
 784   }
 785   return $FreeRadicalElectrons;
 786 }
 787 
 788 # Set valence model...
 789 #
 790 sub SetValenceModel {
 791   my($This, $ValenceModel) = @_;
 792 
 793   if ($ValenceModel !~ /^(MDLValenceModel|DaylightValenceModel|InternalValenceModel|MayaChemToolsValenceModel)$/i) {
 794     carp "Warning: ${ClassName}->SetValenceModel: The current release of MayaChemTools doesn't support the specified valence model $ValenceModel. Supported valence models: MDLValenceModel, DaylightValenceModel, InternalValenceModel or MayaChemToolsValenceModel. Using internal valence model...";
 795     $ValenceModel = 'InternalValenceModel';
 796   }
 797 
 798   $This->SetProperty('ValenceModel', $ValenceModel);
 799 
 800   return $This;
 801 }
 802 
 803 # Get valence model...
 804 #
 805 sub GetValenceModel {
 806   my($This) = @_;
 807 
 808   # Is ValenceModel property explicitly set?
 809   if ($This->HasProperty('ValenceModel')) {
 810     return $This->GetProperty('ValenceModel');
 811   }
 812 
 813   # Used internal valence model as default model...
 814   return 'InternalValenceModel';
 815 }
 816 
 817 # Get molecular formula by collecting information about all atoms in the molecule and
 818 # composing the formula using Hills ordering system:
 819 #   . C shows up first and H follows assuming C is present
 820 #   . All other standard elements are sorted alphanumerically
 821 #   . All other non-stanard atom symbols are also sorted alphanumerically and follow standard elements
 822 #
 823 # Caveats:
 824 #   . By default, missing hydrogens and nonelements are also included
 825 #   . Elements for disconnected fragments are combined into the same formula
 826 #
 827 # Handle formula generation for disconnected structures. e.g: molecule generated by
 828 # [Na+].[O-]c1ccccc1
 829 #
 830 sub GetMolecularFormula {
 831   my($This, $IncludeMissingHydrogens, $IncludeNonElements) = @_;
 832   my($MolecularFormula, $AtomSymbol, $ElementsCountRef, $NonElementsCountRef);
 833 
 834   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 1;
 835   $IncludeNonElements = defined($IncludeNonElements) ? $IncludeNonElements : 1;
 836 
 837   # Get elements count and setup molecular formula...
 838   ($ElementsCountRef, $NonElementsCountRef) = $This->GetElementsAndNonElements($IncludeMissingHydrogens);
 839   $MolecularFormula = '';
 840 
 841   # Count C and H first...
 842   if (exists $ElementsCountRef->{C} ) {
 843     $MolecularFormula .= 'C' . ($ElementsCountRef->{C} > 1 ? $ElementsCountRef->{C} : '');
 844     delete $ElementsCountRef->{C};
 845 
 846     if (exists $ElementsCountRef->{H} ) {
 847       $MolecularFormula .= 'H' . ($ElementsCountRef->{H} > 1 ? $ElementsCountRef->{H} : '');
 848       delete $ElementsCountRef->{H};
 849     }
 850   }
 851 
 852   # Sort elements...
 853   for $AtomSymbol (sort {$a cmp $b} keys %{$ElementsCountRef}) {
 854     $MolecularFormula .= $AtomSymbol . ($ElementsCountRef->{$AtomSymbol} > 1 ? $ElementsCountRef->{$AtomSymbol} : '');
 855   }
 856 
 857   # Sort non-elements...
 858   if ($IncludeNonElements) {
 859     for $AtomSymbol (sort {$a cmp $b} keys %{$NonElementsCountRef}) {
 860       $MolecularFormula .= $AtomSymbol . ($NonElementsCountRef->{$AtomSymbol} > 1 ? $NonElementsCountRef->{$AtomSymbol} : '');
 861     }
 862   }
 863 
 864   # Check formal charge...
 865   my($FormalCharge);
 866   $FormalCharge = $This->GetFormalCharge();
 867   if ($FormalCharge) {
 868     # Setup formal charge string...
 869     my($FormalChargeString);
 870     if ($FormalCharge == 1 ) {
 871       $FormalChargeString =  "+";
 872     }
 873     elsif ($FormalCharge == -1 ) {
 874       $FormalChargeString =  "-";
 875     }
 876     else {
 877       $FormalChargeString = ($FormalCharge > 0) ? ("+" . abs($FormalCharge)) : ("-" . abs($FormalCharge));
 878     }
 879     $MolecularFormula = "${MolecularFormula}${FormalChargeString}";
 880   }
 881 
 882   return $MolecularFormula;
 883 }
 884 
 885 # Count elements and non-elements in molecule and return references to hashes
 886 # containing count of elements and non-elements respectively. By default, missing
 887 # hydrogens are not added to the element hash.
 888 #
 889 #
 890 sub GetElementsAndNonElements {
 891   my($This, $IncludeMissingHydrogens) = @_;
 892   my($Atom, $AtomicNumber, $AtomSymbol, $NumOfMissingHydrogens, @Atoms, %ElementsCount, %NonElementsCount);
 893 
 894   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 0;
 895 
 896   %ElementsCount = (); %NonElementsCount = ();
 897   $NumOfMissingHydrogens = 0;
 898 
 899   # Count elements and non elements...
 900   @Atoms = $This->GetAtoms();
 901   for $Atom (@Atoms) {
 902     $AtomicNumber = $Atom->GetAtomicNumber();
 903     $AtomSymbol = $Atom->GetAtomSymbol();
 904     if ($AtomicNumber) {
 905       if (exists $ElementsCount{$AtomSymbol}) {
 906         $ElementsCount{$AtomSymbol} += 1;
 907       }
 908       else {
 909         $ElementsCount{$AtomSymbol} = 1;
 910       }
 911       if ($IncludeMissingHydrogens) {
 912         $NumOfMissingHydrogens += $Atom->GetNumOfMissingHydrogens();
 913       }
 914     }
 915     else {
 916       if (exists $NonElementsCount{$AtomSymbol}) {
 917         $NonElementsCount{$AtomSymbol} += 1;
 918       }
 919       else {
 920         $NonElementsCount{$AtomSymbol} = 1;
 921       }
 922     }
 923   }
 924   if ($IncludeMissingHydrogens && $NumOfMissingHydrogens) {
 925     $AtomSymbol = 'H';
 926     if (exists $ElementsCount{$AtomSymbol}) {
 927       $ElementsCount{$AtomSymbol} += $NumOfMissingHydrogens;
 928     }
 929     else {
 930       $ElementsCount{$AtomSymbol} = $NumOfMissingHydrogens;
 931     }
 932   }
 933 
 934   return (\%ElementsCount, \%NonElementsCount);
 935 }
 936 
 937 # Get number of element and non-elements in molecule. By default, missing
 938 # hydrogens are not added to element count.
 939 #
 940 sub GetNumOfElementsAndNonElements {
 941   my($This, $IncludeMissingHydrogens) = @_;
 942   my($ElementCount, $NonElementCount, $Atom);
 943 
 944   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 0;
 945 
 946   ($ElementCount, $NonElementCount) = (0, 0);
 947 
 948   ATOM: for $Atom ($This->GetAtoms()) {
 949     if (!$Atom->GetAtomicNumber()) {
 950       $NonElementCount++;
 951       next ATOM;
 952     }
 953     # Process element...
 954     $ElementCount++;
 955     if ($IncludeMissingHydrogens) {
 956       if (!$Atom->IsHydrogen()) {
 957         $ElementCount += $Atom->GetNumOfMissingHydrogens();
 958       }
 959     }
 960   }
 961 
 962   return ($ElementCount, $NonElementCount);
 963 }
 964 
 965 # Calculate elemental composition and return reference to arrays
 966 # containing elements and their percent composition.
 967 #
 968 # Caveats:
 969 #   . By default, missing hydrogens are included
 970 #   . Non elemnents are ignored
 971 #   . Mass number are ignored
 972 #
 973 sub GetElementalComposition {
 974   my($This, $IncludeMissingHydrogens) = @_;
 975   my($MolecularFormula, $IncludeNonElements, $ElementsCountRef, $NonElementsCountRef, $ElementsRef, $ElementsCompositionRef);
 976 
 977   $IncludeMissingHydrogens = defined($IncludeMissingHydrogens) ? $IncludeMissingHydrogens : 1;
 978 
 979   $IncludeNonElements = 0;
 980   ($ElementsCountRef, $NonElementsCountRef) = $This->GetElementsAndNonElements($IncludeMissingHydrogens);
 981 
 982   $MolecularFormula = $This->GetMolecularFormula($IncludeMissingHydrogens, $IncludeNonElements);
 983 
 984   ($ElementsRef, $ElementsCompositionRef) = MolecularFormula::CalculateElementalComposition($MolecularFormula);
 985 
 986   return ($ElementsRef, $ElementsCompositionRef);
 987 }
 988 
 989 # Using refernece to element and its composition arrays, format composition information
 990 # as: Element: Composition;...
 991 #
 992 sub FormatElementalCompositionInformation {
 993   my($FirstParameter, $SecondParameter, $ThirdParameter, $FourthParameter) = @_;
 994   my($This, $ElementsRef, $ElementCompositionRef, $Precision);
 995 
 996   if (_IsMolecule($FirstParameter)) {
 997     ($This, $ElementsRef, $ElementCompositionRef, $Precision) = ($FirstParameter, $SecondParameter, $ThirdParameter, $FourthParameter);
 998   }
 999   else {
1000     ($ElementsRef, $ElementCompositionRef, $Precision) = ($FirstParameter, $SecondParameter, $ThirdParameter);
1001   }
1002   my($FormattedInfo) = '';
1003 
1004   if (!(defined($ElementsRef) && @{$ElementsRef})) {
1005     carp "Warning: ${ClassName}->FormatElementalCompositionInformation: Elements list is not defined or empty...";
1006     return undef;
1007   }
1008   if (!(defined($ElementCompositionRef) && @{$ElementCompositionRef})) {
1009     carp "Warning: ${ClassName}->FormatElementalCompositionInformation: Elements composition list is not defined or empty...";
1010     return undef;
1011   }
1012 
1013   if (!defined $Precision) {
1014     $Precision = 2;
1015   }
1016 
1017   $FormattedInfo = MolecularFormula::FormatCompositionInfomation($ElementsRef, $ElementCompositionRef, $Precision);
1018 
1019   return $FormattedInfo;
1020 }
1021 
1022 # Copy molecule and its associated data using Storable::dclone and update:
1023 #
1024 #   o Atom references corresponding atoms and bonds objects
1025 #   o Bond object references
1026 #
1027 # Object IDs for Atoms and Bonds don't get changed. So there is no need to clear
1028 # up any exisiting ring data attached to molecule via graph as vertex IDs.
1029 #
1030 sub Copy {
1031   my($This) = @_;
1032   my($NewMolecule, $Atom, $NewAtom, $AtomID, @Atoms, @AtomIDs, %AtomsIDsToNewAtoms);
1033 
1034   $NewMolecule = Storable::dclone($This);
1035 
1036   # Update atom references stored as vertex property...
1037 
1038   @Atoms = (); @AtomIDs = ();
1039   %AtomsIDsToNewAtoms = ();
1040 
1041   @AtomIDs = $This->GetVertices();
1042   if (@AtomIDs) {
1043     @Atoms = $This->GetVerticesProperty('Atom', @AtomIDs);
1044   }
1045 
1046   for $Atom (@Atoms) {
1047     $AtomID = $Atom->GetID();
1048 
1049     # Setup a reference to copied atom object...
1050     $NewAtom = $Atom->Copy();
1051     $AtomsIDsToNewAtoms{$AtomID} = $NewAtom;
1052 
1053     # Update atom reference to new atom object...
1054     $NewMolecule->UpdateVertexProperty('Atom', $NewAtom, $AtomID);
1055   }
1056 
1057   # Update bond object and bond atom references stored as edge property...
1058   my($Index, $AtomID1, $AtomID2, $Bond, $NewBond, $NewAtom1, $NewAtom2, @EdgesAtomsIDs);
1059   @EdgesAtomsIDs = ();
1060   @EdgesAtomsIDs = $This->GetEdges();
1061   for ($Index = 0; $Index < $#EdgesAtomsIDs; $Index += 2) {
1062     $AtomID1 = $EdgesAtomsIDs[$Index]; $AtomID2 = $EdgesAtomsIDs[$Index + 1];
1063 
1064     # Get reference to bond object...
1065     $Bond = $This->GetEdgeProperty('Bond', $AtomID1, $AtomID2);
1066 
1067     # Make a new bond object and update its atom references...
1068     $NewBond = $Bond->Copy();
1069     $NewAtom1 = $AtomsIDsToNewAtoms{$AtomID1};
1070     $NewAtom2 = $AtomsIDsToNewAtoms{$AtomID2};
1071     $NewBond->SetAtoms($NewAtom1, $NewAtom2);
1072 
1073     # Update bond object reference in the new molecule...
1074     $NewMolecule->UpdateEdgeProperty('Bond', $NewBond, $AtomID1, $AtomID2);
1075   }
1076 
1077   return $NewMolecule;
1078 }
1079 
1080 # Get number of connected components...
1081 #
1082 sub GetNumOfConnectedComponents {
1083   my($This) = @_;
1084   my($NumOfComponents);
1085 
1086   $NumOfComponents = $This->GetConnectedComponentsVertices();
1087 
1088   return $NumOfComponents;
1089 }
1090 
1091 # Return a reference to an array containing molecules corresponding
1092 # to connected components sorted in decreasing order of component size...
1093 #
1094 sub GetConnectedComponents {
1095   my($This) = @_;
1096   my($Index, @ComponentMolecules, @ConnectedComponents);
1097 
1098   @ConnectedComponents = ();
1099   @ConnectedComponents = $This->GetConnectedComponentsVertices();
1100   @ComponentMolecules = ();
1101 
1102   for $Index (0 .. $#ConnectedComponents) {
1103     push @ComponentMolecules, $This->_GetConnectedComponent(\@ConnectedComponents, $Index);
1104   }
1105   return @ComponentMolecules;
1106 }
1107 
1108 # Return a reference to largest connected component as a molecule object...
1109 #
1110 sub GetLargestConnectedComponent {
1111   my($This) = @_;
1112   my($LargestComponentIndex, @ConnectedComponents);
1113 
1114   $LargestComponentIndex = 0;
1115   @ConnectedComponents = ();
1116   @ConnectedComponents = $This->GetConnectedComponentsVertices();
1117 
1118   return $This->_GetConnectedComponent(\@ConnectedComponents, $LargestComponentIndex);
1119 }
1120 
1121 # Return connected component as a molecule...
1122 #
1123 sub _GetConnectedComponent {
1124   my($This, $ConnectedComponentsRef, $ComponentIndex) = @_;
1125   my($ComponentMolecule);
1126 
1127   # Copy existing molecule...
1128   $ComponentMolecule = $This->Copy();
1129 
1130   # Delete all atoms besides the ones in specified component...
1131   $ComponentMolecule->_DeleteConnectedComponents($ConnectedComponentsRef, $ComponentIndex);
1132 
1133   # Clear any deteced rings...
1134   if ($ComponentMolecule->HasRings()) {
1135     $ComponentMolecule->ClearRings();
1136   }
1137   return $ComponentMolecule;
1138 }
1139 
1140 # Delete atoms corresponding to all connected components except the one specified...
1141 #
1142 sub _DeleteConnectedComponents {
1143   my($This, $ConnectedComponentsRef, $KeepComponentIndex) = @_;
1144   my($Index, $AtomID);
1145 
1146   INDEX: for $Index (0 .. $#{$ConnectedComponentsRef}) {
1147     if ($Index == $KeepComponentIndex) {
1148       next INDEX;
1149     }
1150     for $AtomID (@{$ConnectedComponentsRef->[$Index]}) {
1151       $This->DeleteVertex($AtomID);
1152     }
1153   }
1154   return $This;
1155 }
1156 
1157 # Return an array containing references to atom arrays corresponding to atoms of
1158 # connected components sorted in order of their decreasing size...
1159 #
1160 sub GetConnectedComponentsAtoms {
1161   my($This) = @_;
1162   my($Index, @ComponentsAtoms, @ConnectedComponents);
1163 
1164   @ConnectedComponents = ();
1165   @ConnectedComponents = $This->GetConnectedComponentsVertices();
1166 
1167   @ComponentsAtoms = ();
1168   for $Index (0 .. $#ConnectedComponents) {
1169     my(@ComponentAtoms);
1170 
1171     @ComponentAtoms = ();
1172     @ComponentAtoms = $This->_GetConnectedComponentAtoms(\@ConnectedComponents, $Index);
1173     push @ComponentsAtoms, \@ComponentAtoms;
1174   }
1175   return @ComponentsAtoms;
1176 }
1177 
1178 # Return an array containing atoms correspondig to largest connected component...
1179 #
1180 sub GetLargestConnectedComponentAtoms {
1181   my($This) = @_;
1182   my($LargestComponentIndex, @ConnectedComponents);
1183 
1184   $LargestComponentIndex = 0;
1185   @ConnectedComponents = ();
1186   @ConnectedComponents = $This->GetConnectedComponentsVertices();
1187 
1188   return $This->_GetConnectedComponentAtoms(\@ConnectedComponents, $LargestComponentIndex);
1189 }
1190 
1191 # Return an array containing atoms corresponding to specified connected component...
1192 #
1193 sub _GetConnectedComponentAtoms {
1194   my($This, $ConnectedComponentsRef, $ComponentIndex) = @_;
1195   my($AtomID, @AtomIDs, @ComponentAtoms);
1196 
1197   @ComponentAtoms = ();
1198   @AtomIDs = ();
1199 
1200   for $AtomID (@{$ConnectedComponentsRef->[$ComponentIndex]}) {
1201     push @AtomIDs, $AtomID;
1202   }
1203   @ComponentAtoms = $This->_GetAtomsFromAtomIDs(@AtomIDs);
1204 
1205   return @ComponentAtoms;
1206 }
1207 
1208 # Except for the largest connected component, delete atoms corresponding to all other
1209 # connected components...
1210 #
1211 sub KeepLargestComponent {
1212   my($This) = @_;
1213   my($LargestComponentIndex, @ConnectedComponents);
1214 
1215   @ConnectedComponents = ();
1216   @ConnectedComponents = $This->GetConnectedComponentsVertices();
1217   if (@ConnectedComponents == 1) {
1218     return $This;
1219   }
1220   $LargestComponentIndex = 0;
1221   $This->_DeleteConnectedComponents(\@ConnectedComponents, $LargestComponentIndex);
1222 
1223   # Clear any deteced rings...
1224   if ($This->HasRings()) {
1225     $This->ClearRings();
1226   }
1227 
1228   return $This;
1229 }
1230 
1231 # Get an array of topologically sorted atoms starting from a specified atom or
1232 # an arbitrary atom in the molecule...
1233 #
1234 sub GetTopologicallySortedAtoms {
1235   my($This, $StartAtom) = @_;
1236   my(@SortedAtoms);
1237 
1238   @SortedAtoms = ();
1239   if (defined($StartAtom) && !$This->HasAtom($StartAtom)) {
1240     carp "Warning: ${ClassName}->_GetTopologicallySortedAtoms: No atoms retrieved: Start atom doesn't exist...";
1241     return @SortedAtoms;
1242   }
1243   my($StartAtomID, @AtomIDs);
1244 
1245   @AtomIDs = ();
1246   $StartAtomID = defined($StartAtom) ? $StartAtom->GetID() : undef;
1247 
1248   @AtomIDs = $This->GetTopologicallySortedVertices($StartAtomID);
1249   @SortedAtoms = $This->_GetAtomsFromAtomIDs(@AtomIDs);
1250 
1251   return @SortedAtoms;
1252 }
1253 
1254 # Detect rings in molecule...
1255 #
1256 sub DetectRings {
1257   my($This) = @_;
1258 
1259   # Use graph method to detect all cycles and associate 'em to graph as graph
1260   # and vertex properties...
1261   return $This->DetectCycles();
1262 }
1263 
1264 # Clear rings in molecule...
1265 #
1266 sub ClearRings {
1267   my($This) = @_;
1268 
1269   # Use graph method to clear all cycles...
1270   $This->ClearCycles();
1271 
1272   return $This;
1273 }
1274 
1275 # Setup rings type paths to use during all ring related methods. Possible values:
1276 # Independent or All. Default is to use Independent rings.
1277 #
1278 sub SetActiveRings {
1279   my($This, $RingsType) = @_;
1280 
1281   if (!defined $This->SetActiveCyclicPaths($RingsType)) {
1282     return undef;
1283   }
1284   return $This;
1285 }
1286 
1287 # Is it a supported aromaticity model?
1288 #
1289 sub IsSupportedAromaticityModel {
1290   my($FirstParameter, $SecondParameter) = @_;
1291   my($This, $AromaticityModel);
1292 
1293   if (_IsMolecule($FirstParameter)) {
1294     ($This, $AromaticityModel) = ($FirstParameter, $SecondParameter);
1295   }
1296   else {
1297     ($This, $AromaticityModel) = (undef, $FirstParameter);
1298   }
1299 
1300   return exists $CanonicalAromaticityModelNamesMap{lc($AromaticityModel)} ? 1 : 0;
1301 }
1302 
1303 # Get a list of supported aromaticity model names...
1304 #
1305 sub GetSupportedAromaticityModels {
1306   return (sort values %CanonicalAromaticityModelNamesMap);
1307 }
1308 
1309 # Set aromaticity model...
1310 #
1311 sub SetAromaticityModel {
1312   my($This, $AromaticityModel) = @_;
1313 
1314   if (!$This->IsSupportedAromaticityModel($AromaticityModel)) {
1315     my(@SupportedModels) = $This->GetSupportedAromaticityModels();
1316 
1317     carp "Warning: ${ClassName}->SetAromaticityModel: The current release of MayaChemTools doesn't support the specified aromaticity model $AromaticityModel. Supported aromaticity models defined in AromaticityModelsData.csv file are: @SupportedModels . Using MayaChemToolsAromaticityModel...";
1318     $AromaticityModel = 'MayaChemToolsAromaticityModel';
1319   }
1320 
1321   $This->SetProperty('AromaticityModel', $AromaticityModel);
1322 
1323   return $This;
1324 }
1325 
1326 # Get aromaticity model...
1327 #
1328 sub GetAromaticityModel {
1329   my($This) = @_;
1330 
1331   # Is ValenceModel property explicitly set?
1332   if ($This->HasProperty('AromaticityModel')) {
1333     return $This->GetProperty('AromaticityModel');
1334   }
1335 
1336   # Used internal aromaticity model as default model...
1337   return 'MayaChemToolsAromaticityModel';
1338 }
1339 
1340 # Identify aromatic rings and ring systems in a molecule and set aromaticity for
1341 # corresponding atoms and bonds.
1342 #
1343 # What is aromaticity? [ Ref 124 ] It's in the code of the implementer, did you
1344 # say? Agree. The implementation of aromaticity varies widely across different
1345 # packages [ Ref 125 ]; additionally, the implementation details are not always
1346 # completely available, and it's not possible to figure out the exact implementation
1347 # of aromaticity across various packages. Using the publicly available information,
1348 # however, one can try to reproduce the available results to the extent possible,
1349 # along with parameterizing all the control parameters used to implement different
1350 # aromaticity models, and that's exactly what the current release of MayaChemTools
1351 # does.
1352 #
1353 # The implementation of aromaticity corresponding to various aromaticity models in
1354 # MayaChemTools package is driven by an external CSV file AromaticityModelsData.csv,
1355 # which is distributed with the package and is available in lib/data directory. The CSV
1356 # files contains names of supported aromaticity models, along with various control
1357 # parameters and their values. This file is loaded and processed during instantiation
1358 # of Molecule class and data corresponding to specific aromaticity model are used
1359 # to detect aromaticity for that model. Any new aromaticity model added to the
1360 # aromaticity data file, using different combinations of values for existing control
1361 # parameters would work without any changes to the code; the addition of any new
1362 # control parameters, however, requires its implementation in the code used to
1363 # calculate number of pi electrons available towards delocalization in a ring or ring
1364 # systems.
1365 #
1366 # The current release of MayaChemTools package supports these aromaticity
1367 # models: MDLAromaticityModel, TriposAromaticityModel, MMFFAromaticityModel,
1368 # ChemAxonBasicAromaticityModel, ChemAxonGeneralAromaticityModel,
1369 # DaylightAromaticityModel, MayaChemToolsAromaticityModel.
1370 #
1371 # The current list of control parameters available to detect aromaticity corresponding
1372 # to different aromaticity models are: AllowHeteroRingAtoms, HeteroRingAtomsList,
1373 # AllowExocyclicDoubleBonds, AllowHomoNuclearExocyclicDoubleBonds,
1374 # AllowElectronegativeRingAtomExocyclicDoubleBonds, AllowRingAtomFormalCharge,
1375 # AllowHeteroRingAtomFormalCharge, MinimumRingSize. The values for these control
1376 # parameters are specified in AromaticityModelsData.csv file.
1377 #
1378 # Although definition of aromaticity differs across various aromaticity models, a ring
1379 # or a ring system containing 4n + 2 pi electrons (Huckel's rule) corresponding to
1380 # alternate single and double bonds, in general, is considered aromatic.
1381 #
1382 # The available valence free electrons on heterocyclic ring atoms, involved in two single
1383 # ring bonds, are also allowed to participate in pi electron delocalizaiton for most of
1384 # the supported aromaticity models.
1385 #
1386 # The presence of exocyclic terminal double bond on ring atoms involved in pi electron
1387 # delocalization is only allowed for some of the aromaticity models. Additionally, the type
1388 # atoms involved in exocyclic terminal double bonds may result in making a ring or ring
1389 # system non-aromatic.
1390 #
1391 sub DetectAromaticity {
1392   my($This) = @_;
1393 
1394   # Delete aromaticity property for atoms and bonds...
1395   $This->DeleteAromaticity();
1396 
1397   # Any ring out there...
1398   if (!$This->HasRings()) {
1399     return $This;
1400   }
1401 
1402   if ($This->HasFusedRings()) {
1403     $This->_DetectAromaticityUsingFusedRingSets();
1404   }
1405   else {
1406     $This->_DetectAromaticityUsingIndividualRings();
1407   }
1408   return $This;
1409 }
1410 
1411 # Go over all rings and set aromaticity property for corresponding ring atoms
1412 # and bonds involved in aromatic rings...
1413 #
1414 sub _DetectAromaticityUsingIndividualRings {
1415   my($This) = @_;
1416 
1417   return $This->_DetectRingsAromaticity($This->GetRings());
1418 }
1419 
1420 # For each fused ring set, detect aromaticity by treating all of its ring as one aromatic
1421 # system for counting pi electrons to satisfy Huckel's rule; In case of a failure, rings in
1422 # fused set are treated individually for aromaticity detection. Additionally, non-fused
1423 # rings are handled on their own during aromaticity detection.
1424 #
1425 # Note:
1426 #   . pi electrons in common bonds involved in fused ring sets are only counted once.
1427 #
1428 #
1429 sub _DetectAromaticityUsingFusedRingSets {
1430   my($This) = @_;
1431   my($Index, $RingAtomsRef, $FusedRingSetRef, $FusedRingSetsRef, $NonFusedRingsRef, @FusedRingSetIsAromatic);
1432 
1433   ($FusedRingSetsRef, $NonFusedRingsRef) = $This->GetFusedAndNonFusedRings();
1434 
1435   @FusedRingSetIsAromatic = ();
1436   RINGSET: for $Index (0 .. $#{$FusedRingSetsRef}) {
1437     $FusedRingSetRef = $FusedRingSetsRef->[$Index];
1438     $FusedRingSetIsAromatic[$Index] = 0;
1439 
1440     my($NumOfPiElectronsInRingSet, $NumOfPiElectronsInRing, $Bond, $BondID, %FusedRingSetsBondsMap, %FusedRingSetsBondsVisitedMap, %FusedRingBondsMap);
1441 
1442     $NumOfPiElectronsInRingSet = 0;
1443 
1444     %FusedRingSetsBondsMap = ();
1445     %FusedRingSetsBondsVisitedMap = ();
1446     %FusedRingBondsMap = ();
1447 
1448     # Setup a bond ID map for all bonds in fused ring set and another one
1449     # for bonds involved in more than one ring...
1450     #
1451     for $RingAtomsRef (@{$FusedRingSetRef}) {
1452       for $Bond ($This->GetRingBonds(@{$RingAtomsRef})) {
1453         $BondID = $Bond->GetID();
1454         $FusedRingSetsBondsMap{$BondID} = $BondID;
1455 
1456         if ($Bond->GetNumOfRings() == 2) {
1457           $FusedRingBondsMap{$BondID} = $BondID;
1458         }
1459       }
1460     }
1461 
1462     for $RingAtomsRef (@{$FusedRingSetRef}) {
1463       my(@RingBonds);
1464 
1465       @RingBonds = ();
1466       @RingBonds = $This->GetRingBonds(@{$RingAtomsRef});
1467       $NumOfPiElectronsInRing = $This->_GetNumOfPiElectronsAvailableForDelocalization($RingAtomsRef, \@RingBonds, \%FusedRingSetsBondsMap, \%FusedRingSetsBondsVisitedMap, \%FusedRingBondsMap);
1468 
1469       if (!$NumOfPiElectronsInRing) {
1470         next RINGSET;
1471       }
1472       $NumOfPiElectronsInRingSet += $NumOfPiElectronsInRing;
1473     }
1474     if ($This->_DoPiElectronSatifyHuckelsRule($NumOfPiElectronsInRingSet)) {
1475       $FusedRingSetIsAromatic[$Index] = 1;
1476     }
1477   }
1478 
1479   # Set atom and bond aromatic flags for ring sets whose pi electrons satisfy Huckel's rule; otherwise,
1480   # treat rings in a ring set as individual rings for detecting aromaticity...
1481   for $Index (0 .. $#{$FusedRingSetsRef}) {
1482     $FusedRingSetRef = $FusedRingSetsRef->[$Index];
1483     if ($FusedRingSetIsAromatic[$Index]) {
1484       $This->_SetRingsAromaticity(@{$FusedRingSetRef});
1485     }
1486     else {
1487       $This->_DetectRingsAromaticity(@{$FusedRingSetRef});
1488     }
1489   }
1490 
1491   $This->_DetectRingsAromaticity(@{$NonFusedRingsRef});
1492 
1493   return $This;
1494 }
1495 
1496 # Detect and set aromaticity for rings...
1497 #
1498 sub _DetectRingsAromaticity {
1499   my($This, @Rings) = @_;
1500   my($RingAtom, $RingBond, $RingAtomsRef);
1501 
1502   RING: for $RingAtomsRef (@Rings) {
1503     if (!$This->_CheckRingAromaticity(@{$RingAtomsRef})) {
1504       next RING;
1505     }
1506     $This->_SetRingAromaticity(@{$RingAtomsRef});
1507   }
1508   return $This;
1509 }
1510 
1511 # Set aromatic property for all all atoms and bonds involved in all specified rings..
1512 #
1513 sub _SetRingsAromaticity {
1514   my($This, @Rings) = @_;
1515   my($RingAtomsRef );
1516 
1517   for $RingAtomsRef (@Rings) {
1518     $This->_SetRingAromaticity(@{$RingAtomsRef});
1519   }
1520   return $This;
1521 }
1522 
1523 # Set aromatic property for all all atoms and bonds involved in ring..
1524 #
1525 sub _SetRingAromaticity {
1526   my($This, @RingAtoms) = @_;
1527   my($RingAtom, $RingBond);
1528 
1529   for $RingAtom (@RingAtoms) {
1530     $RingAtom->SetAromatic(1);
1531   }
1532   for $RingBond ($This->GetRingBonds(@RingAtoms)) {
1533     $RingBond->SetAromatic(1);
1534   }
1535   return $This;
1536 }
1537 
1538 
1539 # For a ring to be an aromatic ring, all of its atoms must have aromatic property
1540 # set.
1541 #
1542 sub IsRingAromatic {
1543   my($This, @RingAtoms) = @_;
1544   my($RingAtom);
1545 
1546   for $RingAtom (@RingAtoms) {
1547     if (!$RingAtom->IsAromatic()) {
1548       return 0;
1549     }
1550   }
1551   return 1;
1552 }
1553 
1554 # Delete aromatic property for all atoms and bonds...
1555 #
1556 sub DeleteAromaticity {
1557   my($This) = @_;
1558 
1559   return $This->_DeleteAtomsAndBondsAromaticity();
1560 }
1561 
1562 # Check ring aromaticity...
1563 #
1564 sub _CheckRingAromaticity {
1565   my($This, @RingAtoms) = @_;
1566   my($NumOfPiElectrons, $BondID, @RingBonds);
1567 
1568   @RingBonds = ();
1569   @RingBonds = $This->GetRingBonds(@RingAtoms);
1570 
1571   $NumOfPiElectrons = $This->_GetNumOfPiElectronsAvailableForDelocalization(\@RingAtoms, \@RingBonds);
1572 
1573   return $This->_DoPiElectronSatifyHuckelsRule($NumOfPiElectrons);
1574 }
1575 
1576 # Get number of pi electrons available for delocalizaiton in a ring or ring system...
1577 #
1578 sub _GetNumOfPiElectronsAvailableForDelocalization {
1579   my($This, $RingAtomsRef, $RingBondsRef, $FusedRingSetsBondsMapRef, $FusedRingSetsBondsVisitedMapRef, $FusedRingBondsMapRef) = @_;
1580   my($AromaticityModelName, $AromaticityModelDataRef, $ExocyclicDoubleBondsDataMapRef, $NumOfConjugatedDoubleBonds, $NumOfExplicitAromaticBonds, $NumOfRingAtomElectronPairs, $NumOfRingBondsProcessed, $NumOfPiElectrons, $Index, $RingBond, $RingAtom, $RingAtomSymbol, $BondOrder, $RingAtomID, $RingBondID, $PreviousIndex, $PreviousRingBond, $ExcludeFreeRadicalElectrons, %ElectronPairContributionProcessedMap);
1581 
1582   $AromaticityModelName = $CanonicalAromaticityModelNamesMap{lc($This->GetAromaticityModel())};
1583   $AromaticityModelDataRef = \%{$AromaticityModelsDataMap{$AromaticityModelName}};
1584 
1585   # Perform an intial check for potential atomatic ring atoms..
1586   if (!$This->_CheckRingAtomsForPotentialAromaticity($RingAtomsRef, $RingBondsRef, $AromaticityModelDataRef)) {
1587     return 0;
1588   }
1589 
1590   $ExocyclicDoubleBondsDataMapRef = undef;
1591   $ExcludeFreeRadicalElectrons = 1;
1592 
1593   %ElectronPairContributionProcessedMap = ();
1594 
1595   ($NumOfPiElectrons, $NumOfRingBondsProcessed, $NumOfConjugatedDoubleBonds, $NumOfExplicitAromaticBonds, $NumOfRingAtomElectronPairs) = ('0') x 5;
1596 
1597   # Go over ring atoms and bonds to check their participation in aromaticity and count
1598   # pi electrons available for delocalization corresponding to various aromaticity models...
1599   #
1600   RINGBOND: for $Index (0 .. $#{$RingBondsRef}) {
1601     $RingBond = $RingBondsRef->[$Index];
1602     $RingAtom = $RingAtomsRef->[$Index];
1603     $BondOrder = $RingBond->GetBondOrder();
1604 
1605     # Is this ring bond part of a fused ring system which has been already processed?
1606     if (defined($FusedRingSetsBondsVisitedMapRef) && $RingBond->GetNumOfRings() == 2) {
1607       $RingBondID = $RingBond->GetID();
1608       if (exists $FusedRingSetsBondsVisitedMapRef->{$RingBondID}) {
1609         next RINGBOND;
1610       }
1611       $FusedRingSetsBondsVisitedMapRef->{$RingBondID} = $RingBondID;
1612     }
1613     $NumOfRingBondsProcessed++;
1614 
1615     # For first ring, previous ring bond corrresponds to last ring bond...
1616     $PreviousIndex = $Index ? ($Index -1) : $#{$RingBondsRef};
1617     $PreviousRingBond = $RingBondsRef->[$PreviousIndex];
1618 
1619     # Check for presence of alternate single/double bond configuration, and pesence of
1620     # hetero atoms with two single ring bonds along with any exocyclic double bonds...
1621     #
1622     BONDORDER: {
1623       # Is current ring double bond in an alternate single/double bond configuration?
1624       if ($BondOrder == 2) {
1625         if ($PreviousRingBond->GetBondOrder() != 1) {
1626           return 0;
1627         }
1628         $NumOfConjugatedDoubleBonds += 1;
1629         last BONDORDER;
1630       }
1631 
1632       # Is current ring bond order correspond to an explicit aromatic bond?
1633       if ($BondOrder == 1.5) {
1634         if ($PreviousRingBond->GetBondOrder() != 1.5) {
1635           return 0;
1636         }
1637         $NumOfExplicitAromaticBonds += 1;
1638         last BONDORDER;
1639       }
1640 
1641       # Check for potential hetero atoms involved in two single ring bonds along
1642       # with any terminal exocyclic bonds...
1643       if ($BondOrder == 1) {
1644         if ($PreviousRingBond->GetBondOrder() != 1) {
1645           # Part of a conjugated system...
1646           last BONDORDER;
1647         }
1648 
1649         # Identify any exocylic bonds on rings atoms...
1650         if (!defined $ExocyclicDoubleBondsDataMapRef) {
1651           $ExocyclicDoubleBondsDataMapRef = $This->_IdentifyRingAtomsInvolvedInExocyclicDoubleBonds($RingAtomsRef, $RingBondsRef, $FusedRingSetsBondsMapRef);
1652         }
1653 
1654         # Is current ring atom part of an allowed exocyclic terminal bond?
1655         if (!$This->_CheckPotentialAromaticRingAtomForExocylicDoubleBonds($RingAtom, $AromaticityModelDataRef, $ExocyclicDoubleBondsDataMapRef)) {
1656           return 0;
1657         }
1658 
1659         # Is it allowed to have any formal charge?
1660         if (!$This->_CheckPotentialAromaticRingAtomForFormalCharge($RingAtom, $AromaticityModelDataRef)) {
1661           return 0;
1662         }
1663 
1664         # It it an allowed hetero ring atom or a carbon atom?
1665         if (!$This->_CheckPotentialAromaticRingAtomForAllowedHeteroAtoms($RingAtom, $AromaticityModelDataRef)) {
1666           return 0;
1667         }
1668 
1669         $RingAtomID = $RingAtom->GetID();
1670         $ElectronPairContributionProcessedMap{$RingAtomID} = $RingAtomID;
1671 
1672         # Is it able to donate a pair for electrons towards pi electron delocalization?
1673         if ($RingAtom->GetValenceFreeElectrons($ExcludeFreeRadicalElectrons) >= 2) {
1674           # Possibilites:
1675           #   . Hetero atom with or without formal charge and an available electron pair
1676           #   . Carbon atom with -ve formal charge and with an available electron pair
1677           #
1678           $NumOfRingAtomElectronPairs += 1;
1679         }
1680         else {
1681           # Is ring atom involved in two single bonds without any electron pair allowed?
1682           if (!$This->_AllowRingAtomInTwoSingleBondsWithoutElectronPair($RingAtom, $RingBond, $PreviousRingBond, $ExocyclicDoubleBondsDataMapRef, $FusedRingBondsMapRef)) {
1683             return 0;
1684           }
1685         }
1686         last BONDORDER;
1687       }
1688 
1689       # Any other type of ring atom/bond is not allowed to contribute towards pi electron count
1690       # and caused loss of aromaticity...
1691       return 0;
1692     }
1693   }
1694 
1695   # Check for any electron pair contributions towards pi electron delocalization due to
1696   # -ve formal charge on ring atoms which haven't been already processed and part of
1697   # conjugated single/double bond system...
1698   #
1699   $NumOfRingAtomElectronPairs += $This->_GetElectronPairsContributionFromConjugatedRingAtoms($RingAtomsRef, $RingBondsRef, $ExcludeFreeRadicalElectrons, $AromaticityModelDataRef, \%ElectronPairContributionProcessedMap);
1700 
1701   # Setup pi electron count available for delocalization...
1702   COUNT: {
1703     if ($NumOfExplicitAromaticBonds == $NumOfRingBondsProcessed) {
1704       # Each aromatic bond contribute one electron towards pi electron delocalization...
1705       $NumOfPiElectrons = $NumOfExplicitAromaticBonds;
1706       last COUNT;
1707     }
1708 
1709     # Each conjugated double bond contribute two electrons towards pi electron delocalization...
1710     $NumOfPiElectrons = 2*$NumOfConjugatedDoubleBonds + 2*$NumOfRingAtomElectronPairs;
1711   }
1712 
1713   return $NumOfPiElectrons;
1714 }
1715 
1716 # Check ring atoms for their potential participation in aromatic systems..
1717 #
1718 sub _CheckRingAtomsForPotentialAromaticity {
1719   my($This, $RingAtomsRef, $RingBondsRef, $AromaticityModelDataRef) = @_;
1720   my($Index, $RingBond, $RingAtom);
1721 
1722   # Check availability of ring atoms and bonds...
1723   if (!(defined($RingAtomsRef) && @{$RingBondsRef})) {
1724     return 0;
1725   }
1726 
1727   # Is there any minimum ring size limit?
1728   if ($AromaticityModelDataRef->{MinimumRingSize}) {
1729     if (@{$RingAtomsRef} < $AromaticityModelDataRef->{MinimumRingSize}) {
1730       return 0;
1731     }
1732   }
1733 
1734   # Make sure ring bond order is not greater than 2 and ring atom is not connected to more
1735   # than 3 other atoms to eliminate any non sp2 carbon atoms and still allow for hetero atoms
1736   # to contrbute towards electron delocalization...
1737   #
1738   for $Index (0 .. $#{$RingBondsRef}) {
1739     $RingBond = $RingBondsRef->[$Index];
1740     $RingAtom = $RingAtomsRef->[$Index];
1741 
1742     if (($RingBond->GetBondOrder() > 2) || ($RingAtom->GetNumOfBonds() + $RingAtom->GetNumOfMissingHydrogens()) > 3) {
1743       return 0;
1744     }
1745   }
1746 
1747   return 1;
1748 }
1749 
1750 # Identify any exocylic double bonds on ring atoms...
1751 #
1752 sub _IdentifyRingAtomsInvolvedInExocyclicDoubleBonds {
1753   my($This, $RingAtomsRef, $RingBondsRef, $FusedRingSetsBondsMapRef) = @_;
1754   my($Index, $RingAtom, $RingBond, $RingAtomID, $Bond, $BondID, $BondedAtom, $RingBondsMapRef, %RingBondsMap, %ExocyclicDoubleBondsDataMap);
1755 
1756   # Setup a ring bond map to process exocyclic bonds...
1757   $RingBondsMapRef = undef;
1758   %RingBondsMap = ();
1759 
1760   if (defined $FusedRingSetsBondsMapRef) {
1761     $RingBondsMapRef = $FusedRingSetsBondsMapRef;
1762   }
1763   else {
1764     for $BondID (map { $_->GetID() } @{$RingBondsRef}) {
1765       $RingBondsMap{$BondID} = $BondID;
1766     }
1767     $RingBondsMapRef = \%RingBondsMap;
1768   }
1769 
1770   # Intialize exocyclic terminal double bond data...
1771   %ExocyclicDoubleBondsDataMap = ();
1772   %{$ExocyclicDoubleBondsDataMap{RingAtomID}} = ();
1773 
1774   for $Index (0 .. $#{$RingBondsRef}) {
1775     $RingBond = $RingBondsRef->[$Index];
1776     $RingAtom = $RingAtomsRef->[$Index];
1777 
1778     $RingAtomID = $RingAtom->GetID();
1779 
1780     BOND: for $Bond ($RingAtom->GetBonds()) {
1781       if ($Bond->GetBondOrder != 2) {
1782         next BOND;
1783       }
1784 
1785       # Is it part of ring or ring system under consideration?
1786       if (exists $RingBondsMapRef->{$Bond->GetID()}) {
1787         next BOND;
1788       }
1789 
1790       # Is bonded atom in a ring or a non-terminal atom?
1791       $BondedAtom = $Bond->GetBondedAtom($RingAtom);
1792       if ($BondedAtom->IsInRing() || !$BondedAtom->IsTerminal() ) {
1793         next BOND;
1794       }
1795 
1796       # Track exocyclic terminal double bond information...
1797       if (!exists $ExocyclicDoubleBondsDataMap{RingAtomID}{$RingAtomID}) {
1798         @{$ExocyclicDoubleBondsDataMap{RingAtomID}{$RingAtomID}} = ();
1799       }
1800       push @{$ExocyclicDoubleBondsDataMap{RingAtomID}{$RingAtomID}}, $BondedAtom;
1801     }
1802   }
1803 
1804   return \%ExocyclicDoubleBondsDataMap;
1805 }
1806 
1807 # Check to see whether ring atoms are allowed to participate in exocyclic terminal double
1808 # bonds...
1809 #
1810 sub _CheckPotentialAromaticRingAtomForExocylicDoubleBonds {
1811   my($This, $RingAtom, $AromaticityModelDataRef, $ExocyclicDoubleBondsDataMapRef) = @_;
1812   my($RingAtomID, $ExocyclicTerminalAtom, $RingAtomElectronegativity, $TerminalAtomElectronagativity);
1813 
1814   $RingAtomID = $RingAtom->GetID();
1815 
1816   # Is it part of an exocyclic terminal double bond?
1817   if (!exists $ExocyclicDoubleBondsDataMapRef->{RingAtomID}{$RingAtomID}) {
1818     return 1;
1819   }
1820 
1821   # Are exocyclic terminal double bonds allowed?
1822   if (!$AromaticityModelDataRef->{AllowExocyclicDoubleBonds}) {
1823     return 0;
1824   }
1825 
1826   # Are there multiple exocyclic double bonds?
1827   if (@{$ExocyclicDoubleBondsDataMapRef->{RingAtomID}{$RingAtomID}} > 1) {
1828     return 0;
1829   }
1830   ($ExocyclicTerminalAtom) = @{$ExocyclicDoubleBondsDataMapRef->{RingAtomID}{$RingAtomID}};
1831 
1832   # Are homo nuclear exocyclic terminal double bonds allowed?
1833   if (!$AromaticityModelDataRef->{AllowHomoNuclearExocyclicDoubleBonds}) {
1834     if ($RingAtom->GetAtomicNumber() == $ExocyclicTerminalAtom->GetAtomicNumber()) {
1835       return 0;
1836     }
1837   }
1838 
1839   # Are ring atoms with higher electronegativity allowed in exocyclic double bonds?
1840   if (!$AromaticityModelDataRef->{AllowElectronegativeRingAtomExocyclicDoubleBonds}) {
1841     $RingAtomElectronegativity = PeriodicTable::GetElementPaulingElectronegativity($RingAtom->GetAtomicNumber());
1842     $TerminalAtomElectronagativity = PeriodicTable::GetElementPaulingElectronegativity($ExocyclicTerminalAtom->GetAtomicNumber());
1843 
1844     if ($RingAtomElectronegativity && $TerminalAtomElectronagativity) {
1845       if ($RingAtomElectronegativity > $TerminalAtomElectronagativity) {
1846         return 0;
1847       }
1848     }
1849   }
1850 
1851   return 1;
1852 }
1853 
1854 #
1855 # Check for any formal charge participation into electron delocalization...
1856 #
1857 sub _CheckPotentialAromaticRingAtomForFormalCharge {
1858   my($This, $RingAtom, $AromaticityModelDataRef) = @_;
1859   my($FormalCharge);
1860 
1861   # Does atom has any formal charge?
1862   $FormalCharge = $RingAtom->GetFormalCharge();
1863   if (!$FormalCharge) {
1864     return 1;
1865   }
1866 
1867   # Are ring atoms with formal charge allowed to participate in electron delocalization?
1868   if (!$AromaticityModelDataRef->{AllowRingAtomFormalCharge}) {
1869     return 0;
1870   }
1871 
1872   # Are hetero ring atoms with formal charge allowed to participate in electron delocalization?
1873   if (!$RingAtom->IsCarbon()) {
1874     if (!$AromaticityModelDataRef->{AllowHeteroRingAtomFormalCharge}) {
1875       return 0;
1876     }
1877   }
1878 
1879   return 1;
1880 }
1881 
1882 #
1883 # Check ring atoms for allowed hetero atoms...
1884 #
1885 sub _CheckPotentialAromaticRingAtomForAllowedHeteroAtoms {
1886   my($This, $RingAtom, $AromaticityModelDataRef) = @_;
1887   my($RingAtomSymbol);
1888 
1889   # Is it a Carbon atom?
1890   if ($RingAtom->IsCarbon()) {
1891     return 1;
1892   }
1893 
1894   # Are heteroatoms allowed?
1895   if (!$AromaticityModelDataRef->{AllowHeteroRingAtoms}) {
1896     return 0;
1897   }
1898 
1899   # Is it an allowed hetero atom?
1900   $RingAtomSymbol = $RingAtom->GetAtomSymbol();
1901   if (!exists $AromaticityModelDataRef->{HeteroRingAtomsListMapRef}->{$RingAtomSymbol}) {
1902     return 0;
1903   }
1904 
1905   return 1;
1906 }
1907 
1908 # Check for any electron pair contributions toward pi electron delocalization due to
1909 # -ve formal charge on ring atoms which haven't been already processed and part of
1910 # conjugated single/double bond system...
1911 #
1912 sub _GetElectronPairsContributionFromConjugatedRingAtoms {
1913   my($This, $RingAtomsRef, $RingBondsRef, $ExcludeFreeRadicalElectrons, $AromaticityModelDataRef, $ElectronPairContributionProcessedMapRef) = @_;
1914   my($Index, $RingBond, $RingAtom, $NumOfRingAtomElectronPairs, $RingAtomID);
1915 
1916   # Is formal charge allowed on ring atoms?
1917   if (!$AromaticityModelDataRef->{AllowRingAtomFormalCharge}) {
1918     return 0;
1919   }
1920 
1921   $NumOfRingAtomElectronPairs = 0;
1922 
1923   # Process ring atoms...
1924   RINGBOND: for $Index (0 .. $#{$RingBondsRef}) {
1925     $RingBond = $RingBondsRef->[$Index];
1926     $RingAtom = $RingAtomsRef->[$Index];
1927     $RingAtomID = $RingAtom->GetID();
1928 
1929     # Is is already processed?
1930     if (exists $ElectronPairContributionProcessedMapRef->{$RingAtomID}) {
1931       next RINGBOND;
1932     }
1933     $ElectronPairContributionProcessedMapRef->{$RingAtomID} = $RingAtomID;
1934 
1935     # Is it allowed to have any formal charge?
1936     if (!$This->_CheckPotentialAromaticRingAtomForFormalCharge($RingAtom, $AromaticityModelDataRef)) {
1937       next RINGBOND;
1938     }
1939 
1940     # It it an allowed hetero ring atom or a carbon atom?
1941     if (!$This->_CheckPotentialAromaticRingAtomForAllowedHeteroAtoms($RingAtom, $AromaticityModelDataRef)) {
1942       next RINGBOND;
1943     }
1944 
1945     # It is an atom with -ve formal charge?
1946     if ($RingAtom->GetFormalCharge() >= 0) {
1947       next RINGBOND;
1948     }
1949 
1950     # Is it able to donate a pair for electrons towards pi electron delocalization?
1951     if ($RingAtom->GetValenceFreeElectrons($ExcludeFreeRadicalElectrons) < 2) {
1952       next RINGBOND;
1953     }
1954     $NumOfRingAtomElectronPairs += 1;
1955   }
1956 
1957   return $NumOfRingAtomElectronPairs;
1958 }
1959 
1960 # Check for ring atoms involved in two single ring bonds without any available electron
1961 # pair which are allowed to participate in aromatic system, after all other checks
1962 # corresponding to specified aromaticity models have already been performed...
1963 #
1964 sub _AllowRingAtomInTwoSingleBondsWithoutElectronPair {
1965   my($This, $RingAtom, $RingBond, $PreviousRingBond, $ExocyclicDoubleBondsDataMapRef, $FusedRingBondsMapRef) = @_;
1966 
1967   ALLOWRINGATOM: {
1968     if (exists $ExocyclicDoubleBondsDataMapRef->{RingAtomID}{$RingAtom->GetID()}) {
1969       # Ring atom in an exocylic terminal double bond without any available electron pair...
1970       last ALLOWRINGATOM;
1971     }
1972 
1973     if ($RingAtom->GetFormalCharge() > 0) {
1974       # Ring atom with positive formal charge without any available electron pair...
1975       last ALLOWRINGATOM;
1976     }
1977 
1978     if (defined $FusedRingBondsMapRef && (exists $FusedRingBondsMapRef->{$RingBond->GetID()} || exists $FusedRingBondsMapRef->{$PreviousRingBond->GetID()})) {
1979       # Ring atom involved in fused ring bond, which might end up being part of a conjugated
1980       # system in another fused ring...
1981       last ALLOWRINGATOM;
1982     }
1983 
1984     # Ring atom in any other environment is not allowed...
1985     return 0;
1986   }
1987 
1988   return 1;
1989 }
1990 
1991 # Do pi electrons satify huckel's rule: Number of pi electrons correspond to 4n + 2 where
1992 # n is a positive integer...
1993 #
1994 sub _DoPiElectronSatifyHuckelsRule {
1995   my($This, $NumOfPiElectrons) = @_;
1996 
1997   $NumOfPiElectrons = $NumOfPiElectrons - 2;
1998 
1999   return ($NumOfPiElectrons > 0) ? (($NumOfPiElectrons % 4) ? 0 : 1) : 0;
2000 }
2001 
2002 # Delete aromatic property for all atoms and bonds...
2003 #
2004 sub _DeleteAtomsAndBondsAromaticity {
2005   my($This) = @_;
2006   my($Atom, $Bond);
2007 
2008   for $Atom ($This->GetAtoms()) {
2009     $Atom->DeleteAromatic();
2010   }
2011   for $Bond ($This->GetBonds()) {
2012     $Bond->DeleteAromatic();
2013   }
2014   return $This;
2015 }
2016 
2017 # Kekulize marked ring and non-ring aromatic atoms in a molecule...
2018 #
2019 sub KekulizeAromaticAtoms {
2020   my($This) = @_;
2021 
2022   if (!$This->_KekulizeAromaticAtomsInRings()) {
2023     return 0;
2024   }
2025 
2026   if (!$This->_KekulizeAromaticAtomsNotInRings()) {
2027     return 0;
2028   }
2029 
2030   return 1;
2031 }
2032 
2033 # Kekulize marked aromatic atoms in rings and fused ring sets...
2034 #
2035 sub _KekulizeAromaticAtomsInRings {
2036   my($This) = @_;
2037 
2038   if (!$This->HasRings()) {
2039     # Nothing to do...
2040     return 1;
2041   }
2042 
2043   if (!$This->HasAromaticAtomsInRings()) {
2044     # Nothing to do...
2045     return 1;
2046   }
2047 
2048   # Identify fully aromatic fused and individual rings along with any partially aromatic ring components
2049   # using marked aromatic atoms in a molecule and kekulize them as individual stes...
2050   #
2051   my($AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef) = (undef) x 3;
2052   if ($This->HasFusedRings()) {
2053     ($AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef) = $This->_GetFusedAndNonFusedRingsContainingAromaticAtoms();
2054   }
2055   else {
2056     ($AromaticRingsRef, $PartiallyAromaticRingComponentsRef) = $This->_GetIndividualRingsContainingAromaticAtoms();
2057   }
2058 
2059   return $This->_KekulizeCompleteAndPartialAromaticRings($AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef);
2060 }
2061 
2062 # Identify fully aromatic fused and individual rings along with any partially aromatic ring components
2063 # using marked aromatic atoms in a molecule...
2064 #
2065 sub _GetFusedAndNonFusedRingsContainingAromaticAtoms {
2066   my($This)  = @_;
2067   my($Index, $SetAtomsCount, $SetAromaticAtomsCount, $FusedRingSetRef, $FusedRingSetsRef, $NonFusedRingsRef, $IndividualRingsRef, $RingAtomsRef, $RingAtomsCount, $AromaticAtomsCount, $RingAtom, $NonFusedFullyAromaticRingsRef, $NonFusedPartiallyAromaticRingComponentsRef, $PartiallyAromaticRingComponentsRef, @FullyAromaticFusedRingSets, @PotentialFullyAromaticRings, @FullyAromaticRings, @PotentialPartiallyAromaticRings, @PartiallyAromaticRingComponents);
2068 
2069   @FullyAromaticFusedRingSets = ();
2070 
2071   @PotentialFullyAromaticRings = ();
2072   @FullyAromaticRings = ();
2073 
2074   @PotentialPartiallyAromaticRings = ();
2075   @PartiallyAromaticRingComponents = ();
2076 
2077   ($FusedRingSetsRef, $NonFusedRingsRef) = $This->GetFusedAndNonFusedRings();
2078 
2079   # Go over fused ring sets...
2080   RINGSET: for $Index (0 .. $#{$FusedRingSetsRef}) {
2081     $FusedRingSetRef = $FusedRingSetsRef->[$Index];
2082 
2083     $SetAtomsCount = 0;
2084     $SetAromaticAtomsCount = 0;
2085 
2086     for $RingAtomsRef (@{$FusedRingSetRef}) {
2087       $SetAtomsCount += scalar @{$RingAtomsRef};
2088 
2089       for $RingAtom (@{$RingAtomsRef}) {
2090         if ($RingAtom->IsAromatic()) {
2091           $SetAromaticAtomsCount += 1;
2092         }
2093       }
2094     }
2095 
2096     if (!($SetAtomsCount && $SetAromaticAtomsCount)) {
2097       next RINGSET;
2098     }
2099 
2100     if ($SetAromaticAtomsCount == $SetAtomsCount) {
2101       push @FullyAromaticFusedRingSets, $FusedRingSetRef;
2102     }
2103     else {
2104       # Identify any individual rings in partial aromatic fused ring sets which might be
2105       # fully or partially aromatic...
2106       #
2107       RING: for $RingAtomsRef (@{$FusedRingSetRef}) {
2108         $RingAtomsCount = scalar @{$RingAtomsRef};
2109         $AromaticAtomsCount = 0;
2110 
2111         RINGATOM: for $RingAtom (@{$RingAtomsRef}) {
2112           if (!$RingAtom->IsAromatic()) {
2113             next RINGATOM;
2114           }
2115           $AromaticAtomsCount += 1;
2116         }
2117 
2118         if (!($RingAtomsCount && $AromaticAtomsCount)) {
2119           next RING;
2120         }
2121 
2122         if ($RingAtomsCount == $AromaticAtomsCount) {
2123           push @PotentialFullyAromaticRings, $RingAtomsRef;
2124         }
2125         else {
2126           #  Track partially aromatic rings in an different list before removing them for
2127           #  any overlap with other rings and then add to fully aromatic rings...
2128           push @PotentialPartiallyAromaticRings, $RingAtomsRef;
2129         }
2130       }
2131     }
2132   }
2133 
2134   if (@PotentialFullyAromaticRings > 1) {
2135     # Get any fully aromatic fused ring subsets from potentially fully aromatic rings...
2136     my($FullyAromaticFusedRingSetsRefs, $FullyAromaticNonFusedRingsRef);
2137     ($FullyAromaticFusedRingSetsRefs, $FullyAromaticNonFusedRingsRef) = $This->_GetFullyAromaticFusedAndNonFusedRingsInFusedSubset(\@PotentialFullyAromaticRings);
2138 
2139     if (@{$FullyAromaticFusedRingSetsRefs}) {
2140       push @FullyAromaticFusedRingSets, @{$FullyAromaticFusedRingSetsRefs};
2141     }
2142     if (@{$FullyAromaticNonFusedRingsRef}) {
2143       push @FullyAromaticRings, @{$FullyAromaticNonFusedRingsRef};
2144     }
2145   }
2146   else {
2147     push @FullyAromaticRings, @PotentialFullyAromaticRings;
2148   }
2149 
2150   # Go over partial aromatic ring components...
2151   if (@PotentialPartiallyAromaticRings) {
2152     $PartiallyAromaticRingComponentsRef = $This->_GetPartiallyAromaticRingComponents(\@PotentialPartiallyAromaticRings, \@PotentialFullyAromaticRings);
2153     if (@{$PartiallyAromaticRingComponentsRef}) {
2154       push @PartiallyAromaticRingComponents, @{$PartiallyAromaticRingComponentsRef};
2155     }
2156   }
2157 
2158   # Go over non-fused rings...
2159   if (@{$NonFusedRingsRef}) {
2160     ($NonFusedFullyAromaticRingsRef, $NonFusedPartiallyAromaticRingComponentsRef) = $This->_GetRingsContainingAromaticAtoms(@{$NonFusedRingsRef});
2161 
2162     if (@{$NonFusedFullyAromaticRingsRef}) {
2163       push @FullyAromaticRings, @{$NonFusedFullyAromaticRingsRef};
2164     }
2165     if (@{$NonFusedPartiallyAromaticRingComponentsRef}) {
2166       push @PartiallyAromaticRingComponents, @{$NonFusedPartiallyAromaticRingComponentsRef};
2167     }
2168   }
2169 
2170   return (\@FullyAromaticFusedRingSets, \@FullyAromaticRings, \@PartiallyAromaticRingComponents);
2171 }
2172 
2173 # Identify fully aromatic fused sets and non-fused rings in potentially fully aromatic
2174 # rings in fused ring sets...
2175 #
2176 # Fully aromatic rings in fused ring sets might contain fully aromatic fused subsets. These
2177 # fused subets need to be tracked and treated as fused sets.
2178 #
2179 # Note:
2180 #   . Fused ring sets share at least one common bond, which could be used to identify
2181 #     any multiple fully aromatic fused rings sets; absence of a shared ring bond implies
2182 #     there are no fused ring sets.
2183 #
2184 #
2185 sub _GetFullyAromaticFusedAndNonFusedRingsInFusedSubset {
2186   my($This, $PotentialFullyAromaticFusedRingsRef) = @_;
2187   my($RingIndex, $RingIndex1, $RingIndex2, $RingAtom, $RingAtomID, $RingIsFuesd, $RingIndicesGraph, $FusedRingSetIndicesRef, @RingIndices, @FusedRingPairIndices, @FusedRingSetIndicesRefs, @FullyAromaticFusedRingSets, @FullyAromaticRings, %RingIndexToAtomIDMap, %FullyAromaticFusedRingIndexMap);
2188 
2189   @FullyAromaticFusedRingSets = ();
2190   @FullyAromaticRings = ();
2191 
2192   # Setup a ring index map for ring atoms...
2193   #
2194   %RingIndexToAtomIDMap = ();
2195   for $RingIndex (0 .. $#{$PotentialFullyAromaticFusedRingsRef}) {
2196     %{$RingIndexToAtomIDMap{$RingIndex}} = ();
2197     for $RingAtom (@{$PotentialFullyAromaticFusedRingsRef->[$RingIndex]}) {
2198       $RingAtomID = $RingAtom->GetID();
2199       $RingIndexToAtomIDMap{$RingIndex}{$RingAtomID} = $RingAtomID;
2200     }
2201   }
2202 
2203   # Identify fused ring pairs...
2204   #
2205   @RingIndices = ();
2206   @FusedRingPairIndices = ();
2207 
2208   for $RingIndex1 (0 .. $#{$PotentialFullyAromaticFusedRingsRef}) {
2209     push @RingIndices, $RingIndex1;
2210     for $RingIndex2 (($RingIndex1 + 1) .. $#{$PotentialFullyAromaticFusedRingsRef}) {
2211       $RingIsFuesd = 0;
2212       RINGATOM: for $RingAtom (@{$PotentialFullyAromaticFusedRingsRef->[$RingIndex2]}) {
2213         $RingAtomID = $RingAtom->GetID();
2214         if (exists $RingIndexToAtomIDMap{$RingIndex1}{$RingAtomID}) {
2215           $RingIsFuesd = 1;
2216           last RINGATOM;
2217         }
2218       }
2219       if ($RingIsFuesd) {
2220         push @FusedRingPairIndices, ($RingIndex1, $RingIndex2);
2221       }
2222     }
2223   }
2224 
2225   if (!@FusedRingPairIndices) {
2226     # No fused ring subset out there...
2227     push @FullyAromaticRings, @{$PotentialFullyAromaticFusedRingsRef};
2228 
2229     return (\@FullyAromaticFusedRingSets, \@FullyAromaticRings);
2230   }
2231 
2232   # Identify fused ring sets...
2233   #
2234   $RingIndicesGraph = new Graph(@RingIndices);
2235   $RingIndicesGraph->AddEdges(@FusedRingPairIndices);
2236   @FusedRingSetIndicesRefs = $RingIndicesGraph->GetConnectedComponentsVertices();
2237 
2238   # Collect fully aromatic fused ring sets...
2239   #
2240   %FullyAromaticFusedRingIndexMap = ();
2241   for $FusedRingSetIndicesRef (@FusedRingSetIndicesRefs) {
2242     my(@FullyAromaticFusedRingSet) = ();
2243     for $RingIndex (@{$FusedRingSetIndicesRef}) {
2244       $FullyAromaticFusedRingIndexMap{$RingIndex} = $RingIndex;
2245       push @FullyAromaticFusedRingSet, $PotentialFullyAromaticFusedRingsRef->[$RingIndex];
2246     }
2247     if (@FullyAromaticFusedRingSet) {
2248       # Sort rings by size with in the fused ring set...
2249       @FullyAromaticFusedRingSet = sort { scalar @$a <=> scalar @$b } @FullyAromaticFusedRingSet;
2250       push @FullyAromaticFusedRingSets, \@FullyAromaticFusedRingSet;
2251     }
2252   }
2253 
2254   # Collect fully aromatic non-fused rings...
2255   #
2256   RINGINDEX: for $RingIndex (0 .. $#{$PotentialFullyAromaticFusedRingsRef}) {
2257     if (exists $FullyAromaticFusedRingIndexMap{$RingIndex}) {
2258       next RINGINDEX;
2259     }
2260     push @FullyAromaticRings, $PotentialFullyAromaticFusedRingsRef->[$RingIndex];
2261   }
2262 
2263   return (\@FullyAromaticFusedRingSets, \@FullyAromaticRings);
2264 }
2265 
2266 # Identify individual non-fused rings containing aromatic atoms...
2267 #
2268 sub _GetIndividualRingsContainingAromaticAtoms {
2269   my($This)  = @_;
2270 
2271   return $This->_GetRingsContainingAromaticAtoms($This->GetRings());
2272 }
2273 
2274 # Identify individual non-fused rings containing aromatic atoms...
2275 #
2276 sub _GetRingsContainingAromaticAtoms {
2277   my($This, @Rings)  = @_;
2278   my($RingAtom, $RingAtomsRef, $RingAtomsCount, $AromaticAtomsCount, $PartiallyAromaticRingComponentsRef, @FullyAromaticRings, @PartiallyAromaticRings);
2279 
2280   @FullyAromaticRings = ();
2281   @PartiallyAromaticRings = ();
2282 
2283   RING: for $RingAtomsRef (@Rings) {
2284     $RingAtomsCount = scalar @{$RingAtomsRef};
2285     $AromaticAtomsCount = 0;
2286 
2287     for $RingAtom (@{$RingAtomsRef}) {
2288       if ($RingAtom->IsAromatic()) {
2289         $AromaticAtomsCount += 1;
2290       }
2291     }
2292 
2293     if (!($AromaticAtomsCount && $RingAtomsCount)) {
2294       next RING;
2295     }
2296 
2297     if ($AromaticAtomsCount == $RingAtomsCount) {
2298       push @FullyAromaticRings, $RingAtomsRef;
2299     }
2300     else {
2301       push @PartiallyAromaticRings, $RingAtomsRef;
2302     }
2303   }
2304 
2305   $PartiallyAromaticRingComponentsRef = $This->_GetPartiallyAromaticRingComponents(\@PartiallyAromaticRings);
2306 
2307   return (\@FullyAromaticRings, $PartiallyAromaticRingComponentsRef);
2308 }
2309 
2310 # Get connected aromatic components with in partially aromatic rings...
2311 #
2312 sub _GetPartiallyAromaticRingComponents {
2313   my($This, $PotentialPartiallyAromaticRingsRef, $FullyAromaticRingsRef) = @_;
2314   my($RingAtomsRef, $RingAtom, $RingAtomID, $Index, @PartiallyAromaticRingComponents, %FullyAromaticRingAtomsMap);
2315 
2316   @PartiallyAromaticRingComponents = ();
2317 
2318   # Setup a map for atoms involve in fully aromatic rings to remove remove partial rings
2319   # containing only those atoms which are already part of some other fully aromatic ring
2320   # in fused ring scenarios or some other partially aromatic ring...
2321   #
2322   %FullyAromaticRingAtomsMap = ();
2323   if (defined $FullyAromaticRingsRef) {
2324     for $RingAtomsRef (@{$FullyAromaticRingsRef}) {
2325       for $RingAtom (@{$RingAtomsRef}) {
2326         $RingAtomID = $RingAtom->GetID();
2327         $FullyAromaticRingAtomsMap{$RingAtomID} = $RingAtomID;
2328       }
2329     }
2330   }
2331 
2332   # . Identify any connected components with in each partially aromatic ring.
2333   # . Use ring atom indices to figure out connnected components in rings: All ring atoms
2334   #   in a connected component have sequential indices and a difference by more than
2335   #   1 indicates a new component in the list.
2336   #
2337   RING: for $RingAtomsRef (@{$PotentialPartiallyAromaticRingsRef}) {
2338     my(@AromaticRingAtoms, @AromaticRingAtomsIndices);
2339 
2340     @AromaticRingAtoms = ();
2341     @AromaticRingAtomsIndices = ();
2342 
2343     RINGATOM: for $Index (0 .. $#{$RingAtomsRef}) {
2344       $RingAtom = $RingAtomsRef->[$Index];
2345       $RingAtomID = $RingAtom->GetID();
2346 
2347       if (defined $FullyAromaticRingsRef && exists $FullyAromaticRingAtomsMap{$RingAtomID}) {
2348         next RINGATOM;
2349       }
2350       if (!$RingAtom->IsAromatic()) {
2351         next RINGATOM;
2352       }
2353       push @AromaticRingAtoms, $RingAtom;
2354       push @AromaticRingAtomsIndices, $Index;
2355 
2356     }
2357     if (!@AromaticRingAtoms) {
2358       next RING;
2359     }
2360 
2361     # Start off with a new connected component...
2362     #
2363     my($ComponentNum);
2364     $ComponentNum = scalar @PartiallyAromaticRingComponents;
2365     @{$PartiallyAromaticRingComponents[$ComponentNum]} = ();
2366 
2367     $Index = 0;
2368     push @{$PartiallyAromaticRingComponents[$ComponentNum]}, $AromaticRingAtoms[$Index];
2369 
2370     for $Index (1 .. $#AromaticRingAtoms) {
2371       if (($AromaticRingAtomsIndices[$Index] - $AromaticRingAtomsIndices[$Index -1]) > 1) {
2372         # New connected component...
2373         $ComponentNum += 1;
2374         @{$PartiallyAromaticRingComponents[$ComponentNum]} = ();
2375       }
2376       push @{$PartiallyAromaticRingComponents[$ComponentNum]}, $AromaticRingAtoms[$Index];
2377     }
2378   }
2379 
2380   return (\@PartiallyAromaticRingComponents);
2381 }
2382 
2383 # Kekulize fully aromatic fused and individual rings along with any partially aromatic ring
2384 # components...
2385 #
2386 sub _KekulizeCompleteAndPartialAromaticRings {
2387   my($This, $AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef) = @_;
2388   my($Status, $ConnectedPathsAtomsSetsRef, $ConnectedPathsBondsSetsRef, $ConnectdPathsSetsTypesRef, $PathSetIndex, $PathAtom, $AtomID, $BondID, $PathBondsRef, $DeleteAtomsAromaticity, $DeleteBondsAromaticity, %PathAtomsProcessingStatusMap, %PathBondsProcessingStatusMap);
2389 
2390   ($ConnectedPathsAtomsSetsRef, $ConnectedPathsBondsSetsRef, $ConnectdPathsSetsTypesRef) = $This->_SetupCompleteAndPartialAromaticRingsForKekulizaiton($AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef);
2391 
2392   if (!@{$ConnectedPathsAtomsSetsRef}) {
2393     # Nothing to do...
2394     return 1;
2395   }
2396 
2397   # Delete any aromaticity property set for non-ring bonds connected any two ring
2398   # aromatic atoms...
2399   #
2400   $This->_ProcessNonRingAromaticBondsBetweenAromaticRingAtoms();
2401 
2402   %PathAtomsProcessingStatusMap = ();
2403   %PathBondsProcessingStatusMap = ();
2404 
2405   $Status = 1;
2406 
2407   PATHSET: for $PathSetIndex (0 .. $#{$ConnectedPathsAtomsSetsRef}) {
2408     my($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathSetProcessingStatusRef);
2409 
2410     ($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathSetProcessingStatusRef) = (undef) x 3;
2411 
2412     if ($ConnectdPathsSetsTypesRef->[$PathSetIndex] =~ /^FusedAromatic$/i) {
2413       # Fused set of connected paths...
2414       #
2415       my($FusedConnectedPathAtomsSetRef, $FusedConnectedPathBondsSetRef);
2416 
2417       $FusedConnectedPathAtomsSetRef = $ConnectedPathsAtomsSetsRef->[$PathSetIndex];
2418       $FusedConnectedPathBondsSetRef = $ConnectedPathsBondsSetsRef->[$PathSetIndex];
2419 
2420       # Prepare for kekulization...
2421       ($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathSetProcessingStatusRef) = $This->_SetupConnectedPathSetsForKekulization($FusedConnectedPathAtomsSetRef, $FusedConnectedPathBondsSetRef);
2422 
2423       # Perform kekulization starting with the first path set...
2424       $PathSetProcessingStatusRef->[0] = 'Processed';
2425       if (!$This->_KekulizeConnectedPathSets($FusedConnectedPathAtomsSetRef->[0],  $FusedConnectedPathBondsSetRef->[0], $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $FusedConnectedPathAtomsSetRef, $FusedConnectedPathBondsSetRef, $PathSetProcessingStatusRef)) {
2426         # Kekulization failed for the current fused paths set...
2427         $Status = 0;
2428       }
2429     }
2430     else {
2431       # An individual connected path...
2432       #
2433       my(@ConnectedPathAtomsSet, @ConnectedPathBondsSet);
2434 
2435       @ConnectedPathAtomsSet = ($ConnectedPathsAtomsSetsRef->[$PathSetIndex]);
2436       @ConnectedPathBondsSet = ($ConnectedPathsBondsSetsRef->[$PathSetIndex]);
2437 
2438       # Prepare for kekulization...
2439       ($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef) = $This->_SetupConnectedPathSetsForKekulization(\@ConnectedPathAtomsSet, \@ConnectedPathBondsSet);
2440 
2441       # Perform kekulization...
2442       if (!$This->_KekulizeConnectedPathSets($ConnectedPathsAtomsSetsRef->[$PathSetIndex], $ConnectedPathsBondsSetsRef->[$PathSetIndex], $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef)) {
2443         # Kekulization failed for the current path...
2444         $Status = 0;
2445       }
2446     }
2447 
2448     # Did kekulization succeed for the current path or path set?
2449     if (!$Status) {
2450       last PATHSET;
2451     }
2452 
2453     # Track atom and bond processing state for final assignment after kekulization
2454     # is successfully completed for all the paths and fused path sets...
2455     #
2456     for $AtomID (keys %{$AtomProcessingStatusMapRef}) {
2457       $PathAtomsProcessingStatusMap{$AtomID} = $AtomProcessingStatusMapRef->{$AtomID};
2458     }
2459 
2460     for $BondID (keys %{$BondProcessingStatusMapRef}) {
2461       $PathBondsProcessingStatusMap{$BondID} = $BondProcessingStatusMapRef->{$BondID};
2462     }
2463   }
2464 
2465   if (!$Status) {
2466     carp "Warning: ${ClassName}->_KekulizeCompleteAndPartialAromaticRings: Couldn't perform kekulization for marked ring aromatic atoms...";
2467     return 0;
2468   }
2469 
2470   # Use PathAtomsProcessingStatusMap and PathBondsProcessingStatusMap to set
2471   # single/double bonds in the molecule after successful kekulization along with modification of
2472   # any aromatic flags...
2473 
2474   for $PathSetIndex (0 .. $#{$ConnectedPathsAtomsSetsRef}) {
2475     $DeleteAtomsAromaticity = 0; $DeleteBondsAromaticity = 0;
2476 
2477     if ($ConnectdPathsSetsTypesRef->[$PathSetIndex] =~ /^FusedAromatic$/i) {
2478       for $PathBondsRef (@{$ConnectedPathsBondsSetsRef->[$PathSetIndex]}) {
2479         $This->_ProcessBondOrdersAssignedDuringSuccessfulKekulization($PathBondsRef, \%PathBondsProcessingStatusMap, $DeleteBondsAromaticity);
2480       }
2481     }
2482     else {
2483       if ($ConnectdPathsSetsTypesRef->[$PathSetIndex] =~ /^PartiallyAromatic$/i ) {
2484         $DeleteBondsAromaticity = 1; $DeleteAtomsAromaticity = 1;
2485       }
2486 
2487       if ($DeleteAtomsAromaticity) {
2488         for $PathAtom (@{$ConnectedPathsAtomsSetsRef->[$PathSetIndex]}) {
2489           $PathAtom->DeleteAromatic();
2490         }
2491       }
2492 
2493       $This->_ProcessBondOrdersAssignedDuringSuccessfulKekulization($ConnectedPathsBondsSetsRef->[$PathSetIndex], \%PathBondsProcessingStatusMap, $DeleteBondsAromaticity);
2494     }
2495   }
2496 
2497   return 1;
2498 }
2499 
2500 # Look for any aromatic bonds outside the rings between two ring aromatic atoms
2501 # and turn them into single non-aromatic bonds before kekulization; otherwise, kekulization
2502 # fails.
2503 #
2504 # Note:
2505 #   . Two atoms marked as aromatic atoms in two different rings, such as two rings
2506 #     connected through a single bond, are still aromatic, but the bond is outside
2507 #     the ring and shouldn't be marked as aromatic. It should be set to single bond without
2508 #     any aromatic property for kekulization to succeed.
2509 #
2510 #     For example, the molecule  generated by SMILES parser for biphenyl SMILES string
2511 #     "c1ccccc1c2ccccc2" sets up an aromatic bond between the two phenyl rings, as
2512 #     it's connected to two aromatic atoms.
2513 #
2514 sub _ProcessNonRingAromaticBondsBetweenAromaticRingAtoms {
2515   my($This) = @_;
2516   my($Bond, $Atom1, $Atom2);
2517 
2518   BOND: for $Bond ($This->GetBonds()) {
2519     if (!($Bond->IsAromatic() && $Bond->IsNotInRing())) {
2520       next BOND;
2521     }
2522 
2523     ($Atom1, $Atom2) = $Bond->GetAtoms();
2524     if (!($Atom1->IsAromatic() && $Atom2->IsAromatic() && $Atom1->IsInRing() && $Atom2->IsInRing())) {
2525       next BOND;
2526     }
2527 
2528     $Bond->SetBondOrder(1);
2529     $Bond->DeleteAromatic();
2530   }
2531 
2532   return $This;
2533 }
2534 
2535 # Setup completelty aromatic fused and individual rings along with partially aromatic ring
2536 # components as sets of connected paths...
2537 #
2538 sub _SetupCompleteAndPartialAromaticRingsForKekulizaiton {
2539   my($This, $AromaticFusedRingSetsRef, $AromaticRingsRef, $PartiallyAromaticRingComponentsRef) = @_;
2540   my(@ConnectedPathsSets, @ConnectedPathsBondsSets, @ConnectdPathsSetsTypes);
2541 
2542   @ConnectedPathsSets = ();
2543   @ConnectedPathsBondsSets = ();
2544   @ConnectdPathsSetsTypes = ();
2545 
2546   # Setup atoms and bonds for connected paths in fused aromatic ring sets...
2547   #
2548   if (defined $AromaticFusedRingSetsRef && @{$AromaticFusedRingSetsRef}) {
2549     my($RingSetIndex);
2550 
2551     push @ConnectdPathsSetsTypes, ('FusedAromatic') x scalar @{$AromaticFusedRingSetsRef};
2552     push @ConnectedPathsSets, @{$AromaticFusedRingSetsRef};
2553 
2554     for $RingSetIndex (0 .. $#{$AromaticFusedRingSetsRef}) {
2555       my(@AromaticFusedRingBondsSet);
2556 
2557       # Get ring bonds for each ring set...
2558       #
2559       @AromaticFusedRingBondsSet = $This->GetRingBondsFromRings(@{$AromaticFusedRingSetsRef->[$RingSetIndex]});
2560       push @ConnectedPathsBondsSets, \@AromaticFusedRingBondsSet;
2561     }
2562   }
2563 
2564   # Set up atoms and bonds for connected paths in aromatic rings...
2565   #
2566   if (defined $AromaticRingsRef && @{$AromaticRingsRef}) {
2567     my(@AromaticRingBondsSets);
2568 
2569     push @ConnectdPathsSetsTypes, ('Aromatic') x scalar @{$AromaticRingsRef};
2570     push @ConnectedPathsSets, @{$AromaticRingsRef};
2571 
2572     # Get ring bonds for each ring...
2573     @AromaticRingBondsSets = $This->GetRingBondsFromRings(@{$AromaticRingsRef});
2574     push @ConnectedPathsBondsSets, @AromaticRingBondsSets;
2575   }
2576 
2577   # Set up atoms and bonds for connected paths in partially aromatic rings...
2578   #
2579   if (defined $PartiallyAromaticRingComponentsRef && @{$PartiallyAromaticRingComponentsRef}) {
2580     my($ComponentIndex);
2581 
2582     push @ConnectedPathsSets, @{$PartiallyAromaticRingComponentsRef};
2583     push @ConnectdPathsSetsTypes, ('PartiallyAromatic') x scalar @{$PartiallyAromaticRingComponentsRef};
2584 
2585     for $ComponentIndex (0 .. $#{$PartiallyAromaticRingComponentsRef}) {
2586       my(@ComponentBonds);
2587       @ComponentBonds = $This->_GetPathBonds($This->_GetAtomsIDsFromAtoms(@{$PartiallyAromaticRingComponentsRef->[$ComponentIndex]}));
2588       push @ConnectedPathsBondsSets, \@ComponentBonds;
2589     }
2590   }
2591 
2592   return (\@ConnectedPathsSets, \@ConnectedPathsBondsSets, \@ConnectdPathsSetsTypes);
2593 }
2594 
2595 # Process non-ring connected atoms which are marked aromatic and set connected
2596 # bonds as alternate single/double bonds...
2597 #
2598 # Notes:
2599 #   . Atom and bond aromaticity is deleted during kekulization of non-ring atoms.
2600 #
2601 sub _KekulizeAromaticAtomsNotInRings {
2602   my($This) = @_;
2603   my($Status, $PathIndex, $PathAtom, $PathAtomID, $PathBondID, $ConnectedPathsAtomsRef, $ConnectedPathsBondsRef, $DeleteAtomsAromaticity, $DeleteBondsAromaticity, %PathAtomsProcessingStatusMap, %PathBondsProcessingStatusMap);
2604 
2605   if (!$This->HasAromaticAtomsNotInRings()) {
2606     # Nothing to do...
2607     return 1;
2608   }
2609 
2610   # Identify paths for connected components containing non-ring aromatic atoms...
2611   ($ConnectedPathsAtomsRef, $ConnectedPathsBondsRef) = $This->_GetConnectedComponentsPathsForNonRingAromaticAtoms();
2612 
2613   if (!@{$ConnectedPathsAtomsRef}) {
2614     carp "Warning: ${ClassName}->_KekulizeAromaticAtomsNotInRings: Couldn't perform kekulization for marked non-ring aromatic atoms...";
2615     return 0;
2616   }
2617 
2618   %PathAtomsProcessingStatusMap = ();
2619   %PathBondsProcessingStatusMap = ();
2620 
2621   $Status = 1;
2622 
2623   PATH: for $PathIndex (0 .. $#{$ConnectedPathsAtomsRef}) {
2624     my($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, @ConnectedPathAtomsSet, @ConnectedPathBondsSet);
2625 
2626     ($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef) = (undef) x 2;
2627 
2628     @ConnectedPathAtomsSet = ($ConnectedPathsAtomsRef->[$PathIndex]);
2629     @ConnectedPathBondsSet = ($ConnectedPathsBondsRef->[$PathIndex]);
2630 
2631     # Prepare for kekulization...
2632     ($AtomProcessingStatusMapRef, $BondProcessingStatusMapRef) = $This->_SetupConnectedPathSetsForKekulization(\@ConnectedPathAtomsSet, \@ConnectedPathBondsSet);
2633 
2634     # Perform kekulization...
2635     if (!$This->_KekulizeConnectedPathSets($ConnectedPathsAtomsRef->[$PathIndex], $ConnectedPathsBondsRef->[$PathIndex], $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef)) {
2636       # Kekulization failed for the current path...
2637       $Status = 0;
2638       last PATH;
2639     }
2640 
2641     # Track atom and bond processing state for final assignment after kekulization
2642     # is successfully completed for all the paths and fused path sets...
2643     #
2644     for $PathAtomID (keys %{$AtomProcessingStatusMapRef}) {
2645       $PathAtomsProcessingStatusMap{$PathAtomID} = $AtomProcessingStatusMapRef->{$PathAtomID};
2646     }
2647 
2648     for $PathBondID (keys %{$BondProcessingStatusMapRef}) {
2649       $PathBondsProcessingStatusMap{$PathBondID} = $BondProcessingStatusMapRef->{$PathBondID};
2650     }
2651   }
2652 
2653   if (!$Status) {
2654     carp "Warning: ${ClassName}->_KekulizeAromaticAtomsNotInRings: Couldn't perform kekulization for marked non-ring aromatic atoms...";
2655     return 0;
2656   }
2657 
2658   $DeleteAtomsAromaticity = 1; $DeleteBondsAromaticity = 1;
2659   for $PathIndex (0 .. $#{$ConnectedPathsAtomsRef}) {
2660     if ($DeleteAtomsAromaticity) {
2661       for $PathAtom (@{$ConnectedPathsAtomsRef->[$PathIndex]}) {
2662         $PathAtom->DeleteAromatic();
2663       }
2664     }
2665     $This->_ProcessBondOrdersAssignedDuringSuccessfulKekulization($ConnectedPathsBondsRef->[$PathIndex], \%PathBondsProcessingStatusMap, $DeleteBondsAromaticity);
2666   }
2667 
2668   return 1;
2669 }
2670 
2671 # Collect path atoms for connected components paths containing non-ring aromatic atoms...
2672 #
2673 sub _GetConnectedComponentsPathsForNonRingAromaticAtoms {
2674   my($This) = @_;
2675   my($ComponentRef, $AtomIDsRef, $AtomIDsMapRef, $ConnectedComponentsAtomIDsRef, $ConnectedComponentsAtomIDsMapRef, $ConnectedComponentsPathsAtomIDsRef, $ConnectedComponentsPathsAtomsRef, $ConnectedComponentsPathsBondsRef);
2676 
2677   # Retrieve information for marked aromatic atoms not in the rings...
2678   ($AtomIDsRef, $AtomIDsMapRef) = $This->_GetNonRingAromaticAtomIDs();
2679 
2680   # Identify connected components containing marked aromatic atoms not in the rings...
2681   ($ConnectedComponentsAtomIDsRef, $ConnectedComponentsAtomIDsMapRef) = $This->_GetConnectedComponentsForNonRingAromaticAtoms($AtomIDsRef);
2682 
2683   # Identify paths for connected components containing non-ring aromatic atoms...
2684   ($ConnectedComponentsPathsAtomsRef, $ConnectedComponentsPathsBondsRef) = $This->_GetConnectedComponentsPathsAtomsAndBondsForNonRingAromaticAtoms($AtomIDsMapRef, $ConnectedComponentsAtomIDsRef, $ConnectedComponentsAtomIDsMapRef);
2685 
2686   return ($ConnectedComponentsPathsAtomsRef, $ConnectedComponentsPathsBondsRef);
2687 }
2688 
2689 # Collect information for marked aromatic atoms not in the rings...
2690 #
2691 sub _GetNonRingAromaticAtomIDs {
2692   my($This) = @_;
2693   my($Atom, $AtomID, @AtomIDs, %AtomIDsMap);
2694 
2695   @AtomIDs = ();
2696   %AtomIDsMap = ();
2697 
2698   ATOM: for $Atom ($This->GetAtoms()) {
2699     if (!$Atom->IsAromatic()) {
2700       next ATOM;
2701     }
2702     if ($Atom->IsInRing()) {
2703       next ATOM;
2704     }
2705     $AtomID = $Atom->GetID();
2706 
2707     push @AtomIDs, $AtomID;
2708     $AtomIDsMap{$AtomID} = $Atom;
2709   }
2710 
2711   return (\@AtomIDs, \%AtomIDsMap);
2712 }
2713 
2714 # Retrieve connected non-ring atom components as a reference to an array of references
2715 # containing atom IDs of connecnted components...
2716 #
2717 sub _GetConnectedComponentsForNonRingAromaticAtoms {
2718   my($This, $AtomIDsRef) = @_;
2719   my($Index, $AtomID, $AtomIDsGraph, @BondedAtomPairIDs, @ComponentsAtomIDsRefs, @ComponentsAtomIDsMapRefs);
2720 
2721   @ComponentsAtomIDsRefs = ();
2722   @ComponentsAtomIDsMapRefs = ();
2723 
2724   # Get bonded atom pair IDs...
2725   @BondedAtomPairIDs = $This->_GetBondedAtomPairAtomIDsFromAtomIDs(@{$AtomIDsRef});
2726 
2727   if (!@BondedAtomPairIDs) {
2728     return (\@ComponentsAtomIDsRefs, \@ComponentsAtomIDsMapRefs);
2729   }
2730 
2731   $AtomIDsGraph = new Graph(@{$AtomIDsRef});
2732   $AtomIDsGraph->AddEdges(@BondedAtomPairIDs);
2733 
2734   @ComponentsAtomIDsRefs = $AtomIDsGraph->GetConnectedComponentsVertices();
2735 
2736   # Setup atom IDs map for each component...
2737   for $Index (0 .. $#ComponentsAtomIDsRefs) {
2738     %{$ComponentsAtomIDsMapRefs[$Index]} = ();
2739 
2740     for $AtomID (@{$ComponentsAtomIDsRefs[$Index]}) {
2741       $ComponentsAtomIDsMapRefs[$Index]{$AtomID} = $AtomID;
2742     }
2743   }
2744 
2745   return (\@ComponentsAtomIDsRefs, \@ComponentsAtomIDsMapRefs);
2746 }
2747 
2748 # Get linear paths for connected components starting and ending at terminal aromatic atoms,
2749 # which are connected to only one other aromatic atom in the connected component..
2750 #
2751 sub _GetConnectedComponentsPathsAtomsAndBondsForNonRingAromaticAtoms {
2752   my($This, $AtomIDsMapRef, $ComponentsAtomIDsRef, $ComponentsAtomIDsMapRef) = @_;
2753   my($Index, $AtomID, $Atom, $AtomNbr, $AtomNbrID, $NumOfNonRingAromaticNbrs, $AtomIndex1, $AtomIndex2, $AtomID1, $AtomID2, $Atom1, $Atom2, $AtomIDsGraph, $StartTerminalAtomID, $EndTerminalAtomID, @Paths, @PathAtomIDs, @PathsAtoms, @PathsBonds, @TerminalAtomIDs, @AtomIDs, @BondedAtomPairIDs);
2754 
2755   @PathsAtoms = ();
2756   @PathsBonds = ();
2757 
2758   @TerminalAtomIDs = ();
2759 
2760   $Index = 0;
2761   COMPONENT: for $Index (0 .. $#{$ComponentsAtomIDsRef}) {
2762     @{$TerminalAtomIDs[$Index]} = ();
2763 
2764     # Identify terminal atoms for connected components...
2765     #
2766     # Notes:
2767     #   . Terminal atoms are defined as atoms connected to only one marked
2768     #     aromatic atom.
2769     #   . Linear connected compoents contain only two terminal atoms.
2770     #
2771     ATOM: for $AtomID (@{$ComponentsAtomIDsRef->[$Index]}) {
2772       $Atom = $AtomIDsMapRef->{$AtomID};
2773       $NumOfNonRingAromaticNbrs = 0;
2774 
2775       ATOMNBRID: for $AtomNbr ($Atom->GetNeighbors()) {
2776         $AtomNbrID = $AtomNbr->GetID();
2777 
2778         # Is neighbor in the same connected components containing aromatic atoms?
2779         if (!exists $ComponentsAtomIDsMapRef->[$Index]{$AtomNbrID}) {
2780           next ATOMNBRID;
2781         }
2782         $NumOfNonRingAromaticNbrs++;
2783       }
2784 
2785       # Is it a terminal atom?
2786       if ($NumOfNonRingAromaticNbrs != 1) {
2787         next ATOM;
2788       }
2789       push @{$TerminalAtomIDs[$Index]}, $AtomID;
2790     }
2791 
2792     if (@{$TerminalAtomIDs[$Index]} != 2) {
2793       next COMPONENT;
2794     }
2795 
2796     # Setup bonded atom pair IDs for connected component...
2797     #
2798     @AtomIDs = @{$ComponentsAtomIDsRef->[$Index]};
2799     @BondedAtomPairIDs = ();
2800 
2801     for $AtomIndex1 ( 0 .. $#AtomIDs) {
2802       $AtomID1 = $AtomIDs[$AtomIndex1];
2803       $Atom1 = $AtomIDsMapRef->{$AtomID1};
2804 
2805       for $AtomIndex2 ( ($AtomIndex1 + 1) .. $#AtomIDs) {
2806         $AtomID2 = $AtomIDs[$AtomIndex2];
2807         $Atom2 = $AtomIDsMapRef->{$AtomID2};
2808 
2809         if ($Atom1->IsBondedToAtom($Atom2)) {
2810           push @BondedAtomPairIDs, ($AtomID1, $AtomID2);
2811         }
2812       }
2813     }
2814 
2815     if (!@BondedAtomPairIDs) {
2816       next COMPONENT;
2817     }
2818 
2819     # Get path for connected component...
2820     $AtomIDsGraph = new Graph(@AtomIDs);
2821     $AtomIDsGraph->AddEdges(@BondedAtomPairIDs);
2822 
2823     ($StartTerminalAtomID, $EndTerminalAtomID) = sort { $a <=> $b }  @{$TerminalAtomIDs[$Index]};
2824     @Paths = $AtomIDsGraph->GetPathsBetween($StartTerminalAtomID, $EndTerminalAtomID);
2825 
2826     if (@Paths != 1) {
2827       next COMPONENT;
2828     }
2829 
2830     @PathAtomIDs = $Paths[0]->GetVertices();
2831 
2832     my(@PathAtoms);
2833     @PathAtoms = $This->_GetAtomsFromAtomIDs(@PathAtomIDs);
2834     push @PathsAtoms, \@PathAtoms;
2835 
2836     my(@PathBonds);
2837     @PathBonds = $This->_GetPathBonds(@PathAtomIDs);
2838     push @PathsBonds, \@PathBonds;
2839 
2840   }
2841 
2842   return (\@PathsAtoms, \@PathsBonds);
2843 }
2844 
2845 # Setup initial processing status of atoms and bonds involved in connected paths
2846 # before starting kekulization...
2847 #
2848 # Possible atom processing status: DoubleBondPossible, DoubleBondAssigned, DoubleBondNotPossible
2849 # Initial status: DoubleBondPossible or DoubleBondNotPossible
2850 #
2851 # Possible bond processing status: DoubleBondAssigned, SingleBondAssigned, NotProcessed
2852 #
2853 # Possible paths processing status: Processed, NotProcessed
2854 # Initial status: NotProcessed
2855 #
2856 sub _SetupConnectedPathSetsForKekulization {
2857   my($This, $PathAtomsSetsRef, $PathBondsSetsRef) = @_;
2858   my($PathIndex, $PathAtomsRef, $PathBondsRef, $Atom, $AtomID, $Bond, $BondID, %AtomProcessingStatusMap, %BondProcessingStatusMap, @PathsProcessingStatus, %InitialPathBondOrderMap);
2859 
2860   # Possible path set status values: Processed, NotProcessed
2861   # Initial value: NotProcessed
2862   #
2863   @PathsProcessingStatus = ('NotProcessed') x scalar @{$PathAtomsSetsRef};
2864 
2865   # Collect initial bond order of path bonds before setting bond orders to 1
2866   # and use it to set the bond order back to intial value after it has been processed for
2867   # availability of double bonds...
2868   #
2869   %InitialPathBondOrderMap = ();
2870   for $PathBondsRef (@{$PathBondsSetsRef}) {
2871     BOND: for $Bond (@{$PathBondsRef}) {
2872       $BondID = $Bond->GetID();
2873       if (exists $InitialPathBondOrderMap{$BondID}) {
2874         next BOND;
2875       }
2876       $InitialPathBondOrderMap{$BondID} = $Bond->GetBondOrder();
2877       $Bond->SetBondOrder(1);
2878     }
2879   }
2880 
2881   %AtomProcessingStatusMap = ();
2882   %BondProcessingStatusMap = ();
2883 
2884   for $PathIndex (0 .. $#{$PathAtomsSetsRef}) {
2885 
2886     $PathAtomsRef = $PathAtomsSetsRef->[$PathIndex];
2887     ATOM: for $Atom (@{$PathAtomsRef}) {
2888       $AtomID = $Atom->GetID();
2889       if (exists $AtomProcessingStatusMap{$AtomID}) {
2890         next ATOM;
2891       }
2892       $AtomProcessingStatusMap{$AtomID} = ($Atom->GetNumOfBondsAvailableForNonHydrogenAtoms() >= 1) ? 'DoubleBondPossible' : 'DoubleBondNotPossible';
2893     }
2894 
2895     $PathBondsRef = $PathBondsSetsRef->[$PathIndex];
2896     BOND: for $Bond (@{$PathBondsRef}) {
2897       $BondID = $Bond->GetID();
2898       if (exists $BondProcessingStatusMap{$BondID}) {
2899         next BOND;
2900       }
2901       $BondProcessingStatusMap{$BondID} = 'NotProcessed';
2902     }
2903   }
2904 
2905   # Set bond orders back to initial bond orders...
2906   for $PathIndex (0 .. $#{$PathAtomsSetsRef}) {
2907     $PathBondsRef = $PathBondsSetsRef->[$PathIndex];
2908 
2909     for $Bond (@{$PathBondsRef}) {
2910       $BondID = $Bond->GetID();
2911       if (exists $InitialPathBondOrderMap{$BondID}) {
2912         $Bond->SetBondOrder($InitialPathBondOrderMap{$BondID});
2913       }
2914     }
2915   }
2916 
2917   return (\%AtomProcessingStatusMap, \%BondProcessingStatusMap, \@PathsProcessingStatus);
2918 }
2919 
2920 # Kekulize connected path sets corresponding to fused rings, individual rings, or any other
2921 # connected path...
2922 #
2923 # Note:
2924 #   . PathAtomsRef and PathBondsRef contain paths and bonds corresponding to path
2925 #     under consideration for kekulization
2926 #   . PathAtomsSetsRef and PathBondsSetsRef contain any other available paths fused
2927 #     to the path being kekulized
2928 #   . _KekulizeConnectedPathSets is invoked recursively to kekulize all available paths
2929 #
2930 sub _KekulizeConnectedPathSets {
2931   my($This, $PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef) = @_;
2932   my($PathBond);
2933 
2934   # Get next available path bond...
2935   $PathBond = $This->_GetNextAvailablePathBondForKekulization($PathBondsRef, $BondProcessingStatusMapRef);
2936 
2937   if ($PathBond) {
2938     return $This->_ProcessNextAvailablePathBondForKekulization($PathBond, $PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef);
2939   }
2940 
2941   # Did kekulization succeed for the current path bonds?
2942   if (!$This->_DidKekulizationSucceedForPathBonds($PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef)) {
2943     return 0;
2944   }
2945 
2946   # Is there any other path available for kekulization?
2947   ($PathAtomsRef, $PathBondsRef) = $This->_GetNextAvailablePathForKekulization($PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef);
2948 
2949   if ($PathAtomsRef && $PathBondsRef) {
2950     # Recursively call itself to kekulize next path, which could either be a new path or part
2951     # of a fused paths corresponding to fused ring sets...
2952     #
2953     return $This->_KekulizeConnectedPathSets($PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef);
2954   }
2955 
2956   return 1;
2957 }
2958 
2959 # Get next available path bond in a list of path bonds...
2960 #
2961 sub _GetNextAvailablePathBondForKekulization {
2962   my($This, $PathBondsRef, $BondProcessingStatusMapRef) = @_;
2963   my($AvailablePathBond, $PathBond, $PathBondID);
2964 
2965   $AvailablePathBond = undef;
2966 
2967   BOND: for $PathBond (@{$PathBondsRef}) {
2968     $PathBondID = $PathBond->GetID();
2969     if (!exists $BondProcessingStatusMapRef->{$PathBondID}) {
2970       next BOND;
2971     }
2972     if ($BondProcessingStatusMapRef->{$PathBondID} =~ /^NotProcessed$/i) {
2973       $AvailablePathBond = $PathBond;
2974       last BOND;
2975     }
2976   }
2977 
2978   return ($AvailablePathBond);
2979 }
2980 
2981 # Process next available path bond for kekulizaiton...
2982 #
2983 sub _ProcessNextAvailablePathBondForKekulization {
2984   my($This, $PathBond, $PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef) = @_;
2985   my($PathBondID, $PathAtom1, $PathAtom2, $PathAtomID1, $PathAtomID2, %CurrentAtomProcessingStatusMap, %CurrentBondProcessingStatusMap);
2986 
2987   $PathBondID = $PathBond->GetID();
2988 
2989   ($PathAtom1, $PathAtom2) = $PathBond->GetAtoms();
2990   ($PathAtomID1, $PathAtomID2) = ($PathAtom1->GetID(), $PathAtom2->GetID());
2991 
2992   %CurrentAtomProcessingStatusMap = %{$AtomProcessingStatusMapRef};
2993   %CurrentBondProcessingStatusMap = %{$BondProcessingStatusMapRef};
2994 
2995   # Is it possible to assign a double bond to the current path bond?
2996   if ($AtomProcessingStatusMapRef->{$PathAtomID1} =~ /^DoubleBondPossible$/i && $AtomProcessingStatusMapRef->{$PathAtomID2} =~ /^DoubleBondPossible$/i ) {
2997     # Set current bond to double bond by appropriately marking atom and bond process status...
2998     $AtomProcessingStatusMapRef->{$PathAtomID1} = 'DoubleBondAssigned';
2999     $AtomProcessingStatusMapRef->{$PathAtomID2} = 'DoubleBondAssigned';
3000 
3001     $BondProcessingStatusMapRef->{$PathBondID} = 'DoubleBondAssigned';
3002 
3003     # Recursively call  _KekulizeConnectedPathSets to kekulize next available bond...
3004     if ($This->_KekulizeConnectedPathSets($PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef)) {
3005       return 1;
3006     }
3007 
3008     # Double bond at the current ring bond position didn't lead to successful kekulization...
3009     %{$AtomProcessingStatusMapRef} = %CurrentAtomProcessingStatusMap;
3010     %{$BondProcessingStatusMapRef} = %CurrentBondProcessingStatusMap;
3011   }
3012 
3013   # Try single bond at the current ring bond position and recursively call _KekulizeConnectedPathSets to kekulize
3014   # rest of the ring bonds...
3015   #
3016   $BondProcessingStatusMapRef->{$PathBondID} = 'SingleBondAssigned';
3017 
3018   if ($This->_KekulizeConnectedPathSets($PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef)) {
3019     return 1;
3020   }
3021 
3022   %{$AtomProcessingStatusMapRef} = %CurrentAtomProcessingStatusMap;
3023   %{$BondProcessingStatusMapRef} = %CurrentBondProcessingStatusMap;
3024 
3025   # Kekulization didn't work out for path bonds...
3026 
3027   return 0;
3028 
3029 }
3030 
3031 # Get next available path for kekulization from a set of fused ring paths...
3032 #
3033 sub _GetNextAvailablePathForKekulization {
3034   my($This, $PathAtomsSetsRef, $PathBondsSetsRef, $PathsProcessingStatusRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef) = @_;
3035   my($PathIndex, $AvailablePathIndex, $PathAtomsRef, $PathBondsRef, $PathBond, $PathBondID, $MaxNumOfPathBondsProcessed, $NumOfPathBondsProcessed);
3036 
3037   ($PathAtomsRef, $PathBondsRef, $AvailablePathIndex) = (undef) x 3;
3038 
3039   if (!(defined($PathAtomsSetsRef)  && defined($PathBondsSetsRef) && defined($PathsProcessingStatusRef))) {
3040     return ($PathAtomsRef, $PathBondsRef);
3041   }
3042 
3043   $MaxNumOfPathBondsProcessed = -999;
3044   $AvailablePathIndex = undef;
3045 
3046   PATHINDEX: for $PathIndex (0 .. $#{$PathsProcessingStatusRef}) {
3047     if ($PathsProcessingStatusRef->[$PathIndex] =~ /^Processed$/i) {
3048       next PATHINDEX;
3049     }
3050 
3051     # Count of already processed bonds in an unprocessed path bonds through
3052     # their participation in any fused bonds sets...
3053     #
3054     $NumOfPathBondsProcessed = 0;
3055     PATHBOND: for $PathBond (@{$PathBondsSetsRef->[$PathIndex]}) {
3056       $PathBondID = $PathBond->GetID();
3057       if ($BondProcessingStatusMapRef->{$PathBondID} =~ /^NotProcessed$/i) {
3058         next PATHBOND;
3059       }
3060       $NumOfPathBondsProcessed++;
3061     }
3062 
3063     if ($NumOfPathBondsProcessed > $MaxNumOfPathBondsProcessed) {
3064       $AvailablePathIndex = $PathIndex;
3065       $MaxNumOfPathBondsProcessed = $NumOfPathBondsProcessed;
3066     }
3067 
3068   }
3069 
3070   # Is any path available?
3071   if (!$AvailablePathIndex) {
3072     return ($PathAtomsRef, $PathBondsRef);
3073   }
3074 
3075   $PathsProcessingStatusRef->[$AvailablePathIndex] = 'Processed';
3076 
3077   $PathAtomsRef = $PathAtomsSetsRef->[$AvailablePathIndex];
3078   $PathBondsRef = $PathBondsSetsRef->[$AvailablePathIndex];
3079 
3080   return ($PathAtomsRef, $PathBondsRef);
3081 }
3082 
3083 # Check for kekulization in a specific set of path bonds. For successful kekulization, all
3084 # all path atoms marked with DoubleBondPossible must be involved in a path double bond...
3085 #
3086 sub _DidKekulizationSucceedForPathBonds {
3087   my($This, $PathAtomsRef, $PathBondsRef, $AtomProcessingStatusMapRef, $BondProcessingStatusMapRef) = @_;
3088   my($PathAtom, $PathAtomID);
3089 
3090   for $PathAtom (@{$PathAtomsRef}) {
3091     $PathAtomID = $PathAtom->GetID();
3092     if (exists $AtomProcessingStatusMapRef->{$PathAtomID} && $AtomProcessingStatusMapRef->{$PathAtomID} =~ /^DoubleBondPossible$/i) {
3093       return 0;
3094     }
3095   }
3096   return 1;
3097 }
3098 
3099 # Assign bond orders to the bonds in a molecule which have been successfully
3100 # kekulized along with optional clearing of aromaticty property...
3101 #
3102 sub _ProcessBondOrdersAssignedDuringSuccessfulKekulization {
3103   my($This, $BondsRef, $BondsProcessingStatusMapRef, $DeleteBondsAromaticity) = @_;
3104   my($Bond, $BondID, $BondOrder);
3105 
3106   $DeleteBondsAromaticity = defined $DeleteBondsAromaticity ? $DeleteBondsAromaticity : 0;
3107 
3108   BOND: for $Bond (@{$BondsRef}) {
3109     $BondID = $Bond->GetID();
3110 
3111     if (!exists $BondsProcessingStatusMapRef->{$BondID}) {
3112       carp "Warning: ${ClassName}->_ProcessBondOrdersAssignedDuringSuccessfulKekulization: Couldn't process bond with bond ID, $BondID: It's not available in the list of bonds processed for kekulization...";
3113       next BOND;
3114     }
3115 
3116     $BondOrder = ($BondsProcessingStatusMapRef->{$BondID} =~ /^DoubleBondAssigned$/i) ? 2 : 1;
3117     $Bond->SetBondOrder($BondOrder);
3118 
3119     if ($DeleteBondsAromaticity) {
3120       $Bond->DeleteAromatic();
3121     }
3122   }
3123   return $This;
3124 }
3125 
3126 # Does molecule contains aromatic rings?
3127 #
3128 sub HasAromaticRings {
3129   my($This) = @_;
3130 
3131   return $This->GetNumOfAromaticRings() ? 1 : 0;
3132 }
3133 
3134 # Does molecule contains any aromatic atom in a ring?
3135 #
3136 sub HasAromaticAtomsInRings {
3137   my($This) = @_;
3138   my($Atom);
3139 
3140   ATOM: for $Atom ($This->GetAtoms()) {
3141     if (!$Atom->IsAromatic()) {
3142       next ATOM;
3143     }
3144     if ($Atom->IsInRing()) {
3145       return 1;
3146     }
3147   }
3148   return 0;
3149 }
3150 
3151 # Does molecule contains any aromatic atom not in a ring?
3152 #
3153 sub HasAromaticAtomsNotInRings {
3154   my($This) = @_;
3155   my($Atom);
3156 
3157   ATOM: for $Atom ($This->GetAtoms()) {
3158     if (!$Atom->IsAromatic()) {
3159       next ATOM;
3160     }
3161     if ($Atom->IsNotInRing()) {
3162       return 1;
3163     }
3164   }
3165   return 0;
3166 }
3167 
3168 # Does molecule contains rings?
3169 #
3170 sub HasRings {
3171   my($This) = @_;
3172 
3173   return $This->IsCyclic();
3174 }
3175 
3176 # Does molecule contains only one ring?
3177 #
3178 sub HasOnlyOneRing {
3179   my($This) = @_;
3180 
3181   return $This->IsUnicyclic();
3182 }
3183 
3184 # Does molecule contains any rings?
3185 #
3186 sub HasNoRings {
3187   my($This) = @_;
3188 
3189   return $This->IsAcyclic();
3190 }
3191 
3192 # Get size of smallest ring...
3193 #
3194 sub GetSizeOfSmallestRing {
3195   my($This) = @_;
3196 
3197   return $This->GetSizeOfSmallestCycle();
3198 }
3199 
3200 # Get size of largest ring...
3201 #
3202 sub GetSizeOfLargestRing {
3203   my($This) = @_;
3204 
3205   return $This->GetSizeOfLargestCycle();
3206 }
3207 
3208 # Get number of rings...
3209 #
3210 sub GetNumOfRings {
3211   my($This) = @_;
3212 
3213   return $This->GetNumOfCycles();
3214 }
3215 
3216 # Get number of aromatic rings...
3217 #
3218 sub GetNumOfAromaticRings {
3219   my($This) = @_;
3220   my($NumOfRings);
3221 
3222   $NumOfRings = scalar $This->GetAromaticRings();
3223 
3224   return $NumOfRings;
3225 }
3226 
3227 # Get num of rings with odd size...
3228 #
3229 sub GetNumOfRingsWithOddSize {
3230   my($This) = @_;
3231 
3232   return $This->GetNumOfCyclesWithOddSize();
3233 }
3234 
3235 # Get num of rings with even size...
3236 #
3237 sub GetNumOfRingsWithEvenSize {
3238   my($This) = @_;
3239 
3240   return $This->GetNumOfCyclesWithEvenSize();
3241 }
3242 
3243 # Get num of rings with specified size...
3244 #
3245 sub GetNumOfRingsWithSize {
3246   my($This, $RingSize) = @_;
3247 
3248   return $This->GetNumOfCyclesWithSize($RingSize);
3249 }
3250 
3251 # Get num of rings with size less than a specified size...
3252 #
3253 sub GetNumOfRingsWithSizeLessThan {
3254   my($This, $RingSize) = @_;
3255 
3256   return $This->GetNumOfCyclesWithSizeLessThan($RingSize);
3257 }
3258 
3259 # Get num of rings with size greater than a specified size...
3260 #
3261 sub GetNumOfRingsWithSizeGreaterThan {
3262   my($This, $RingSize) = @_;
3263 
3264   return $This->GetNumOfCyclesWithSizeGreaterThan($RingSize);
3265 }
3266 
3267 # Get largest ring as an array containing ring atoms...
3268 #
3269 sub GetLargestRing {
3270   my($This) = @_;
3271 
3272   return $This->_GetRing($This->GetLargestCycle());
3273 }
3274 
3275 # Get smallest ring as an array containing ring atoms...
3276 #
3277 sub GetSmallestRing {
3278   my($This) = @_;
3279 
3280   return $This->_GetRing($This->GetSmallestCycle());
3281 }
3282 
3283 # Get rings as an array containing references to arrays with ring atoms...
3284 #
3285 sub GetRings {
3286   my($This) = @_;
3287 
3288   return $This->_GetRings($This->GetCycles());
3289 }
3290 
3291 # Get aromatic rings as an array containing references to arrays with ring atoms...
3292 #
3293 sub GetAromaticRings {
3294   my($This) = @_;
3295 
3296   return $This->_GetAromaticRings($This->GetCycles());
3297 }
3298 
3299 # Get odd size rings as an array containing references to arrays with ring atoms...
3300 #
3301 sub GetRingsWithOddSize {
3302   my($This) = @_;
3303 
3304   return $This->_GetRings($This->GetCyclesWithOddSize());
3305 }
3306 
3307 # Get even size rings as an array containing references to arrays with ring atoms...
3308 #
3309 sub GetRingsWithEvenSize {
3310   my($This) = @_;
3311 
3312   return $This->_GetRings($This->GetCyclesWithEvenSize());
3313 }
3314 
3315 # Get rings with a specific size as an array containing references to arrays with ring atoms...
3316 #
3317 sub GetRingsWithSize {
3318   my($This, $RingSize) = @_;
3319 
3320   return $This->_GetRings($This->GetCyclesWithSize($RingSize));
3321 }
3322 
3323 # Get rings with size less than a specific size as an array containing references to arrays with ring atoms...
3324 #
3325 sub GetRingsWithSizeLessThan {
3326   my($This, $RingSize) = @_;
3327 
3328   return $This->_GetRings($This->GetCyclesWithSizeLessThan($RingSize));
3329 }
3330 
3331 # Get rings with size greater than a specific size as an array containing references to arrays with ring atoms...
3332 #
3333 sub GetRingsWithSizeGreaterThan {
3334   my($This, $RingSize) = @_;
3335 
3336   return $This->_GetRings($This->GetCyclesWithSizeGreaterThan($RingSize));
3337 }
3338 
3339 # Generate an array of bond objects for an array of ring atoms and return an array
3340 # of bond objects...
3341 #
3342 sub GetRingBonds {
3343   my($This, @RingAtoms) = @_;
3344   my(@Bonds);
3345 
3346   @Bonds = ();
3347   if (!@RingAtoms) {
3348     # Return an empty ring bonds list...
3349     return @Bonds;
3350   }
3351 
3352   my(@RingAtomIDs);
3353 
3354   @RingAtomIDs = ();
3355   @RingAtomIDs = $This->_GetAtomsIDsFromAtoms(@RingAtoms);
3356   if (!@RingAtomIDs) {
3357     carp "Warning: ${ClassName}->GetRingBonds: No ring bonds retrieved: Atom IDs couldn't be retrieved for specified atoms...";
3358     return @Bonds;
3359   }
3360 
3361   # Add start atom to the end to make it a cyclic path for ring: It's taken out during conversion
3362   # of cyclic path to a ring...
3363   push @RingAtomIDs, $RingAtomIDs[0];
3364 
3365   return $This->_GetPathBonds(@RingAtomIDs);
3366 }
3367 
3368 # Generate an array containing references to arrays of ring bond objects for rings specified
3369 # in an array of references to ring atoms...
3370 #
3371 sub GetRingBondsFromRings {
3372   my($This, @RingAtomsSets) = @_;
3373   my($RingAtomsRef, @RingBondsSets);
3374 
3375   @RingBondsSets = ();
3376   for $RingAtomsRef  (@RingAtomsSets) {
3377     my(@RingBonds);
3378     @RingBonds = $This->GetRingBonds(@{$RingAtomsRef});
3379 
3380     push @RingBondsSets, \@RingBonds;
3381   }
3382 
3383   return @RingBondsSets;
3384 }
3385 
3386 # Does molecule has any fused rings?
3387 #
3388 sub HasFusedRings {
3389   my($This) = @_;
3390 
3391   return $This->HasFusedCycles();
3392 }
3393 
3394 # Get references to array of fused ring sets and non-fused rings. Fused ring sets array reference
3395 # contains refernces to arrays of rings; Non-fused rings array reference contains references to
3396 # arrays of ring atoms...
3397 # rings.
3398 #
3399 sub GetFusedAndNonFusedRings {
3400   my($This) = @_;
3401   my($FusedCyclesSetsRef, $NonFusedCyclesRef, @FusedRingSets, @NonFusedRings);
3402 
3403   @FusedRingSets = (); @NonFusedRings = ();
3404   ($FusedCyclesSetsRef, $NonFusedCyclesRef) = $This->GetFusedAndNonFusedCycles();
3405   if (!(defined($FusedCyclesSetsRef) && defined($NonFusedCyclesRef))) {
3406     return (\@FusedRingSets, \@NonFusedRings);
3407   }
3408   my($FusedCyclesSetRef);
3409 
3410   for $FusedCyclesSetRef (@{$FusedCyclesSetsRef}) {
3411     my(@FusedRingSet);
3412     @FusedRingSet = ();
3413     @FusedRingSet = $This->_GetRings(@{$FusedCyclesSetRef});
3414     push @FusedRingSets, \@FusedRingSet;
3415   }
3416 
3417   @NonFusedRings = $This->_GetRings(@{$NonFusedCyclesRef});
3418 
3419   return (\@FusedRingSets, \@NonFusedRings);
3420 }
3421 
3422 # Get rings as an array containing references to arrays with ring atoms...
3423 #
3424 sub _GetRings {
3425   my($This, @CyclicPaths) = @_;
3426   my($CyclicPath, @Rings);
3427 
3428   @Rings = ();
3429   if (!@CyclicPaths) {
3430     return @Rings;
3431   }
3432   if (!@CyclicPaths) {
3433     # Return an empty ring list...
3434     return @Rings;
3435   }
3436 
3437   for $CyclicPath (@CyclicPaths) {
3438     my(@RingAtoms);
3439     @RingAtoms = ();
3440     push @RingAtoms, $This->_GetRing($CyclicPath);
3441 
3442     push @Rings, \@RingAtoms;
3443   }
3444   return @Rings;
3445 }
3446 
3447 # Get aromatic rings as an array containing references to arrays with ring atoms...
3448 #
3449 sub _GetAromaticRings {
3450   my($This, @CyclicPaths) = @_;
3451   my($RingAtomsRef, @Rings, @AromaticRings);
3452 
3453   @AromaticRings = ();
3454   @Rings = $This->_GetRings(@CyclicPaths);
3455 
3456   if (!@Rings) {
3457     return @AromaticRings;
3458   }
3459   RING: for $RingAtomsRef (@Rings) {
3460     if (!$This->IsRingAromatic(@{$RingAtomsRef})) {
3461       next RING;
3462     }
3463     my(@RingAtoms);
3464     @RingAtoms = ();
3465     push @RingAtoms, @{$RingAtomsRef};
3466 
3467     push @AromaticRings, \@RingAtoms;
3468   }
3469   return @AromaticRings;
3470 }
3471 
3472 # Map atom IDs in cyclic path to atoms and return a reference to an array containing ring atoms...
3473 #
3474 # Note:
3475 #   . Start and end vertex is same for cyclic paths. So end atom is removed before
3476 #     returning atoms array as ring atoms...
3477 #
3478 sub _GetRing {
3479   my($This, $CyclicPath) = @_;
3480   my(@RingAtoms);
3481 
3482   @RingAtoms = ();
3483   if (!defined $CyclicPath) {
3484     # Return an empty atoms list...
3485     return @RingAtoms;
3486   }
3487 
3488   @RingAtoms = $This->_GetPathAtoms($CyclicPath);
3489   if (@RingAtoms) {
3490     pop @RingAtoms;
3491   }
3492   return @RingAtoms;
3493 }
3494 
3495 # Map atom IDs to atoms and return a reference to an array containing these atoms...
3496 #
3497 sub _GetPathAtoms {
3498   my($This, $Path) = @_;
3499   my(@PathAtoms);
3500 
3501   @PathAtoms = ();
3502   if (!defined $Path) {
3503     carp "Warning: ${ClassName}->_GetPathAtoms: No path atoms retrieved: Path must be defined...";
3504     return @PathAtoms;
3505   }
3506   my(@AtomIDs);
3507 
3508   @AtomIDs = ();
3509   @AtomIDs = $Path->GetVertices();
3510 
3511   @PathAtoms = $This->_GetAtomsFromAtomIDs(@AtomIDs);
3512 
3513   return @PathAtoms;
3514 }
3515 
3516 # Get bonds for a path specified by atom IDs...
3517 #
3518 sub _GetPathBonds {
3519   my($This, @AtomIDs) = @_;
3520   my($Index, $AtomID1, $AtomID2, @Bonds, @EdgesAtomIDs);
3521 
3522   @Bonds = (); @EdgesAtomIDs = ();
3523 
3524   if (!@AtomIDs || @AtomIDs == 1) {
3525     return @Bonds;
3526   }
3527 
3528   # Setup edges...
3529   for $Index (0 .. ($#AtomIDs - 1) ) {
3530     $AtomID1 = $AtomIDs[$Index];
3531     $AtomID2 = $AtomIDs[$Index + 1];
3532     push @EdgesAtomIDs, ($AtomID1, $AtomID2);
3533   }
3534   @Bonds =  $This->GetEdgesProperty('Bond', @EdgesAtomIDs);
3535 
3536   return @Bonds;
3537 }
3538 
3539 # Map atom ID to an atom...
3540 #
3541 sub _GetAtomFromAtomID {
3542   my($This, $AtomID) = @_;
3543 
3544   return $This->GetVertexProperty('Atom', $AtomID);
3545 }
3546 
3547 # Map atom IDs to atoms and return an array containing these atoms...
3548 #
3549 sub _GetAtomsFromAtomIDs {
3550   my($This, @AtomIDs) = @_;
3551 
3552   return $This->GetVerticesProperty('Atom', @AtomIDs);
3553 }
3554 
3555 # Map atoms to atom IDs and return an array containing these atoms...
3556 #
3557 sub _GetAtomsIDsFromAtoms {
3558   my($This, @Atoms) = @_;
3559 
3560   return map { $_->GetID() } @Atoms;
3561 }
3562 
3563 # Get bonded atom pair atom IDs for specified list of atom IDs...
3564 #
3565 sub _GetBondedAtomPairAtomIDsFromAtomIDs {
3566   my($This, @AtomIDs) = @_;
3567   my($AtomIndex1, $AtomID1, $Atom1, $AtomIndex2, $AtomID2, $Atom2, @Atoms, @BondedAtomPairIDs);
3568 
3569   @BondedAtomPairIDs = ();
3570   @Atoms = $This->_GetAtomsFromAtomIDs(@AtomIDs);
3571 
3572   for $AtomIndex1 ( 0 .. $#Atoms) {
3573     $Atom1 = $Atoms[$AtomIndex1];
3574     $AtomID1 = $Atom1->GetID();
3575 
3576     ATOMINDEX2: for $AtomIndex2 ( ($AtomIndex1 + 1) .. $#Atoms) {
3577       $Atom2 = $Atoms[$AtomIndex2];
3578       if (!$Atom1->IsBondedToAtom($Atom2)) {
3579         next ATOMINDEX2;
3580        }
3581       $AtomID2 = $Atom2->GetID();
3582 
3583       push @BondedAtomPairIDs, ($AtomID1, $AtomID2);
3584     }
3585   }
3586 
3587   return @BondedAtomPairIDs;
3588 }
3589 
3590 # Get bonded atom pair atoms for specified list of atoms...
3591 #
3592 sub _GetBondedAtomPairAtomsFromAtoms {
3593   my($This, @Atoms) = @_;
3594   my($AtomIndex1, $Atom1, $AtomIndex2, $Atom2, @BondedAtomPairAtoms);
3595 
3596   @BondedAtomPairAtoms = ();
3597 
3598   for $AtomIndex1 ( 0 .. $#Atoms) {
3599     $Atom1 = $Atoms[$AtomIndex1];
3600 
3601     ATOMINDEX2: for $AtomIndex2 ( ($AtomIndex1 + 1) .. $#Atoms) {
3602       $Atom2 = $Atoms[$AtomIndex2];
3603       if ($Atom1->IsBondedToAtom($Atom2)) {
3604         next ATOMINDEX2;
3605        }
3606 
3607       push @BondedAtomPairAtoms, ($Atom1, $Atom2);
3608     }
3609   }
3610 
3611   return @BondedAtomPairAtoms;
3612 }
3613 
3614 # Is atom in a ring?
3615 #
3616 sub _IsAtomInRing {
3617   my($This, $Atom) = @_;
3618 
3619   return $This->IsCyclicVertex($Atom->GetID());
3620 }
3621 
3622 # Is atom not in a ring?
3623 #
3624 sub _IsAtomNotInRing {
3625   my($This, $Atom) = @_;
3626 
3627   return $This->IsAcyclicVertex($Atom->GetID());
3628 }
3629 
3630 # Is atom only in one ring?
3631 #
3632 sub _IsAtomInOnlyOneRing {
3633   my($This, $Atom) = @_;
3634 
3635   return $This->IsUnicyclicVertex($Atom->GetID());
3636 }
3637 
3638 # Is atom in a ring of specified size?
3639 #
3640 sub _IsAtomInRingOfSize {
3641   my($This, $Atom, $RingSize) = @_;
3642 
3643   return $This->GetNumOfVertexCyclesWithSize($Atom->GetID(), $RingSize) ? 1 : 0;
3644 }
3645 
3646 # Get size of smallest ring containing specified atom...
3647 #
3648 sub _GetSizeOfSmallestAtomRing {
3649   my($This, $Atom) = @_;
3650 
3651   return $This->GetSizeOfSmallestVertexCycle($Atom->GetID());
3652 }
3653 
3654 # Get size of largest ring containing specified atom...
3655 #
3656 sub _GetSizeOfLargestAtomRing {
3657   my($This, $Atom) = @_;
3658 
3659   return $This->GetSizeOfLargestVertexCycle($Atom->GetID());
3660 }
3661 
3662 # Get number of  rings containing specified atom...
3663 #
3664 sub _GetNumOfAtomRings {
3665   my($This, $Atom) = @_;
3666 
3667   return $This->GetNumOfVertexCycles($Atom->GetID());
3668 }
3669 
3670 # Get number of  rings with odd size containing specified atom...
3671 #
3672 sub _GetNumOfAtomRingsWithOddSize {
3673   my($This, $Atom) = @_;
3674 
3675   return $This->GetNumOfVertexCyclesWithOddSize($Atom->GetID());
3676 }
3677 
3678 # Get number of  rings with even size containing specified atom...
3679 #
3680 sub _GetNumOfAtomRingsWithEvenSize {
3681   my($This, $Atom) = @_;
3682 
3683   return $This->GetNumOfVertexCyclesWithEvenSize($Atom->GetID());
3684 }
3685 
3686 # Get number of  rings with specified size containing specified atom...
3687 #
3688 sub _GetNumOfAtomRingsWithSize {
3689   my($This, $Atom, $RingSize) = @_;
3690 
3691   return $This->GetNumOfVertexCyclesWithSize($Atom->GetID(), $RingSize);
3692 }
3693 
3694 # Get number of  rings with size less than specified containing specified atom...
3695 #
3696 sub _GetNumOfAtomRingsWithSizeLessThan {
3697   my($This, $Atom, $RingSize) = @_;
3698 
3699   return $This->GetNumOfVertexCyclesWithSizeLessThan($Atom->GetID(), $RingSize);
3700 }
3701 
3702 # Get number of  rings with size greater than specified containing specified atom...
3703 #
3704 sub _GetNumOfAtomRingsWithSizeGreaterThan {
3705   my($This, $Atom, $RingSize) = @_;
3706 
3707   return $This->GetNumOfVertexCyclesWithSizeGreaterThan($Atom->GetID(), $RingSize);
3708 }
3709 
3710 # Get smallest ring as an array containing ring atoms...
3711 #
3712 sub _GetSmallestAtomRing {
3713   my($This, $Atom) = @_;
3714 
3715   return $This->_GetRing($This->GetSmallestVertexCycle($Atom->GetID()));
3716 }
3717 
3718 # Get odd size rings an array of references to arrays containing ring atoms...
3719 #
3720 sub _GetLargestAtomRing {
3721   my($This, $Atom) = @_;
3722 
3723   return $This->_GetRing($This->GetLargestVertexCycle($Atom->GetID()));
3724 }
3725 
3726 # Get all rings an array of references to arrays containing ring atoms...
3727 #
3728 sub _GetAtomRings {
3729   my($This, $Atom) = @_;
3730 
3731   return $This->_GetRings($This->GetVertexCycles($Atom->GetID()));
3732 }
3733 
3734 # Get odd size rings an array of references to arrays containing ring atoms...
3735 #
3736 sub _GetAtomRingsWithOddSize {
3737   my($This, $Atom) = @_;
3738 
3739   return $This->_GetRings($This->GetVertexCyclesWithOddSize($Atom->GetID()));
3740 }
3741 
3742 # Get even size rings an array of references to arrays containing ring atoms...
3743 #
3744 sub _GetAtomRingsWithEvenSize {
3745   my($This, $Atom) = @_;
3746 
3747   return $This->_GetRings($This->GetVertexCyclesWithEvenSize($Atom->GetID()));
3748 }
3749 
3750 # Get rings with specified size  an array of references to arrays containing ring atoms...
3751 #
3752 sub _GetAtomRingsWithSize {
3753   my($This, $Atom, $RingSize) = @_;
3754 
3755   return $This->_GetRings($This->GetVertexCyclesWithSize($Atom->GetID(), $RingSize));
3756 }
3757 
3758 # Get rings with size less than specfied size as an array of references to arrays containing ring atoms...
3759 #
3760 sub _GetAtomRingsWithSizeLessThan {
3761   my($This, $Atom, $RingSize) = @_;
3762 
3763   return $This->_GetRings($This->GetVertexCyclesWithSizeLessThan($Atom->GetID(), $RingSize));
3764 }
3765 
3766 # Get rings with size less than specfied size as an array of references to arrays containing ring atoms...
3767 #
3768 sub _GetAtomRingsWithSizeGreaterThan {
3769   my($This, $Atom, $RingSize) = @_;
3770 
3771   return $This->_GetRings($This->GetVertexCyclesWithSizeGreaterThan($Atom->GetID(), $RingSize));
3772 }
3773 
3774 # Is bond in a ring?
3775 #
3776 sub _IsBondInRing {
3777   my($This, $Bond) = @_;
3778   my($Atom1, $Atom2);
3779 
3780   ($Atom1, $Atom2) = $Bond->GetAtoms();
3781 
3782   return $This->IsCyclicEdge($Atom1->GetID(), $Atom2->GetID());
3783 }
3784 
3785 # Is bond not in a ring?
3786 #
3787 sub _IsBondNotInRing {
3788   my($This, $Bond) = @_;
3789   my($Atom1, $Atom2);
3790 
3791   ($Atom1, $Atom2) = $Bond->GetAtoms();
3792 
3793   return $This->IsAcyclicEdge($Atom1->GetID(), $Atom2->GetID());
3794 }
3795 
3796 # Is bond only in one ring?
3797 #
3798 sub _IsBondInOnlyOneRing {
3799   my($This, $Bond) = @_;
3800   my($Atom1, $Atom2);
3801 
3802   ($Atom1, $Atom2) = $Bond->GetAtoms();
3803 
3804   return $This->IsUnicyclicEdge($Atom1->GetID(), $Atom2->GetID());
3805 }
3806 
3807 # Is bond in a ring of specified size?
3808 #
3809 sub _IsBondInRingOfSize {
3810   my($This, $Bond, $RingSize) = @_;
3811   my($Atom1, $Atom2);
3812 
3813   ($Atom1, $Atom2) = $Bond->GetAtoms();
3814 
3815   return $This->GetNumOfEdgeCyclesWithSize($Atom1->GetID(), $Atom2->GetID(), $RingSize) ? 1 : 0;
3816 }
3817 
3818 # Get size of smallest ring containing specified bond...
3819 #
3820 sub _GetSizeOfSmallestBondRing {
3821   my($This, $Bond) = @_;
3822   my($Atom1, $Atom2);
3823 
3824   ($Atom1, $Atom2) = $Bond->GetAtoms();
3825 
3826   return $This->GetSizeOfSmallestEdgeCycle($Atom1->GetID(), $Atom2->GetID());
3827 }
3828 
3829 # Get size of largest ring containing specified bond...
3830 #
3831 sub _GetSizeOfLargestBondRing {
3832   my($This, $Bond) = @_;
3833   my($Atom1, $Atom2);
3834 
3835   ($Atom1, $Atom2) = $Bond->GetAtoms();
3836 
3837   return $This->GetSizeOfLargestEdgeCycle($Atom1->GetID(), $Atom2->GetID());
3838 }
3839 
3840 # Get number of  rings containing specified bond...
3841 #
3842 sub _GetNumOfBondRings {
3843   my($This, $Bond) = @_;
3844   my($Atom1, $Atom2);
3845 
3846   ($Atom1, $Atom2) = $Bond->GetAtoms();
3847 
3848   return $This->GetNumOfEdgeCycles($Atom1->GetID(), $Atom2->GetID());
3849 }
3850 
3851 # Get number of  rings with odd size containing specified bond...
3852 #
3853 sub _GetNumOfBondRingsWithOddSize {
3854   my($This, $Bond) = @_;
3855   my($Atom1, $Atom2);
3856 
3857   ($Atom1, $Atom2) = $Bond->GetAtoms();
3858 
3859   return $This->GetNumOfEdgeCyclesWithOddSize($Atom1->GetID(), $Atom2->GetID());
3860 }
3861 
3862 # Get number of  rings with even size containing specified bond...
3863 #
3864 sub _GetNumOfBondRingsWithEvenSize {
3865   my($This, $Bond) = @_;
3866   my($Atom1, $Atom2);
3867 
3868   ($Atom1, $Atom2) = $Bond->GetAtoms();
3869 
3870   return $This->GetNumOfEdgeCyclesWithEvenSize($Atom1->GetID(), $Atom2->GetID());
3871 }
3872 
3873 # Get number of  rings with specified size containing specified bond...
3874 #
3875 sub _GetNumOfBondRingsWithSize {
3876   my($This, $Bond, $RingSize) = @_;
3877   my($Atom1, $Atom2);
3878 
3879   ($Atom1, $Atom2) = $Bond->GetAtoms();
3880 
3881   return $This->GetNumOfEdgeCyclesWithSize($Atom1->GetID(), $Atom2->GetID(), $RingSize);
3882 }
3883 
3884 # Get number of  rings with size less than specified containing specified bond...
3885 #
3886 sub _GetNumOfBondRingsWithSizeLessThan {
3887   my($This, $Bond, $RingSize) = @_;
3888   my($Atom1, $Atom2);
3889 
3890   ($Atom1, $Atom2) = $Bond->GetAtoms();
3891 
3892   return $This->GetNumOfEdgeCyclesWithSizeLessThan($Atom1->GetID(), $Atom2->GetID(), $RingSize);
3893 }
3894 
3895 # Get number of  rings with size greater than specified containing specified bond...
3896 #
3897 sub _GetNumOfBondRingsWithSizeGreaterThan {
3898   my($This, $Bond, $RingSize) = @_;
3899   my($Atom1, $Atom2);
3900 
3901   ($Atom1, $Atom2) = $Bond->GetAtoms();
3902 
3903   return $This->GetNumOfEdgeCyclesWithSizeGreaterThan($Atom1->GetID(), $Atom2->GetID(), $RingSize);
3904 }
3905 
3906 # Get smallest ring as an array containing ring atoms...
3907 #
3908 sub _GetSmallestBondRing {
3909   my($This, $Bond) = @_;
3910   my($Atom1, $Atom2);
3911 
3912   ($Atom1, $Atom2) = $Bond->GetAtoms();
3913 
3914   return $This->_GetRing($This->GetSmallestEdgeCycle($Atom1->GetID(), $Atom2->GetID()));
3915 }
3916 
3917 # Get odd size rings an array of references to arrays containing ring atoms...
3918 #
3919 sub _GetLargestBondRing {
3920   my($This, $Bond) = @_;
3921   my($Atom1, $Atom2);
3922 
3923   ($Atom1, $Atom2) = $Bond->GetAtoms();
3924 
3925   return $This->_GetRing($This->GetLargestEdgeCycle($Atom1->GetID(), $Atom2->GetID()));
3926 }
3927 
3928 # Get all rings an array of references to arrays containing ring atoms...
3929 #
3930 sub _GetBondRings {
3931   my($This, $Bond) = @_;
3932   my($Atom1, $Atom2);
3933 
3934   ($Atom1, $Atom2) = $Bond->GetAtoms();
3935 
3936   return $This->_GetRings($This->GetEdgeCycles($Atom1->GetID(), $Atom2->GetID()));
3937 }
3938 
3939 # Get odd size rings an array of references to arrays containing ring atoms...
3940 #
3941 sub _GetBondRingsWithOddSize {
3942   my($This, $Bond) = @_;
3943   my($Atom1, $Atom2);
3944 
3945   ($Atom1, $Atom2) = $Bond->GetAtoms();
3946 
3947   return $This->_GetRings($This->GetEdgeCyclesWithOddSize($Atom1->GetID(), $Atom2->GetID()));
3948 }
3949 
3950 # Get even size rings an array of references to arrays containing ring atoms...
3951 #
3952 sub _GetBondRingsWithEvenSize {
3953   my($This, $Bond) = @_;
3954   my($Atom1, $Atom2);
3955 
3956   ($Atom1, $Atom2) = $Bond->GetAtoms();
3957 
3958   return $This->_GetRings($This->GetEdgeCyclesWithEvenSize($Atom1->GetID(), $Atom2->GetID()));
3959 }
3960 
3961 # Get rings with specified size  an array of references to arrays containing ring atoms...
3962 #
3963 sub _GetBondRingsWithSize {
3964   my($This, $Bond, $RingSize) = @_;
3965   my($Atom1, $Atom2);
3966 
3967   ($Atom1, $Atom2) = $Bond->GetAtoms();
3968 
3969   return $This->_GetRings($This->GetEdgeCyclesWithSize($Atom1->GetID(), $Atom2->GetID(), $RingSize));
3970 }
3971 
3972 # Get rings with size less than specfied size as an array of references to arrays containing ring atoms...
3973 #
3974 sub _GetBondRingsWithSizeLessThan {
3975   my($This, $Bond, $RingSize) = @_;
3976   my($Atom1, $Atom2);
3977 
3978   ($Atom1, $Atom2) = $Bond->GetAtoms();
3979 
3980   return $This->_GetRings($This->GetEdgeCyclesWithSizeLessThan($Atom1->GetID(), $Atom2->GetID(), $RingSize));
3981 }
3982 
3983 # Get rings with size less than specfied size as an array of references to arrays containing ring atoms...
3984 #
3985 sub _GetBondRingsWithSizeGreaterThan {
3986   my($This, $Bond, $RingSize) = @_;
3987   my($Atom1, $Atom2);
3988 
3989   ($Atom1, $Atom2) = $Bond->GetAtoms();
3990 
3991   return $This->_GetRings($This->GetEdgeCyclesWithSizeGreaterThan($Atom1->GetID(), $Atom2->GetID(), $RingSize));
3992 }
3993 
3994 
3995 # Get atom paths starting from a specified atom as a reference to an array containing references
3996 # to arrays with path atoms.
3997 #
3998 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with length
3999 # upto a specified length and sharing of bonds in paths traversed. By default, rings are
4000 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4001 #
4002 # Note:
4003 #    . For molecule without any rings, this method returns the same set of atom paths
4004 #      as GetAtomPathsStartingAtWithLengthUpto method.
4005 #
4006 sub GetAllAtomPathsStartingAtWithLengthUpto {
4007   my($This, $StartAtom, $Length, $AllowCycles) = @_;
4008 
4009   return $This->_GetAtomPathsStartingAt('AllAtomPathsWithLengthUpto', $StartAtom, $Length, $AllowCycles);
4010 }
4011 
4012 # Get atom paths starting from a specified atom as a reference to an array containing references
4013 # to arrays with path atoms.
4014 #
4015 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with
4016 # specified length and sharing of bonds in paths traversed. By default, rings are
4017 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4018 #
4019 # Note:
4020 #    . For molecule without any rings, this method returns the same set of atom paths
4021 #      as GetAtomPathsStartingAtWithLengthUpto method.
4022 #
4023 sub GetAllAtomPathsStartingAtWithLength {
4024   my($This, $StartAtom, $Length, $AllowCycles) = @_;
4025 
4026   return $This->_GetAtomPathsStartingAt('AllAtomPathsWithLength', $StartAtom, $Length, $AllowCycles);
4027 }
4028 
4029 # Get atom paths starting from a specified atom as a reference to an array containing references
4030 # to arrays with path atoms.
4031 #
4032 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with all
4033 # possible lengths and sharing of bonds in paths traversed. By default, rings are
4034 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4035 #
4036 # Note:
4037 #    . For molecule without any rings, this method returns the same set of atom paths
4038 #      as GetAtomPathsStartingAt method.
4039 #
4040 sub GetAllAtomPathsStartingAt {
4041   my($This, $StartAtom, $AllowCycles) = @_;
4042 
4043   return $This->_GetAtomPathsStartingAt('AllAtomPathsWithAllLengths', $StartAtom, undef, $AllowCycles);
4044 }
4045 
4046 # Get atom paths starting from a specified atom as a reference to an array containing references
4047 # to arrays with path atoms.
4048 #
4049 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with length
4050 # upto a specified length and no sharing of bonds in paths traversed. By default, rings are
4051 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4052 #
4053 sub GetAtomPathsStartingAtWithLengthUpto {
4054   my($This, $StartAtom, $Length, $AllowCycles) = @_;
4055 
4056   return $This->_GetAtomPathsStartingAt('AtomPathsWithLengthUpto', $StartAtom, $Length, $AllowCycles);
4057 }
4058 
4059 # Get atom paths starting from a specified atom as a reference to an array containing references
4060 # to arrays with path atoms.
4061 #
4062 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with
4063 # specified length and no sharing of bonds in paths traversed. By default, rings are
4064 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4065 #
4066 sub GetAtomPathsStartingAtWithLength {
4067   my($This, $StartAtom, $Length, $AllowCycles) = @_;
4068 
4069   return $This->_GetAtomPathsStartingAt('AtomPathsWithLength', $StartAtom, $Length, $AllowCycles);
4070 }
4071 
4072 # Get atom paths starting from a specified atom as a reference to an array containing references
4073 # to arrays with path atoms.
4074 #
4075 # Path atoms atoms correspond to to all possible paths for specified atom in molecule with all
4076 # possible lengths and no sharing of bonds in paths traversed. By default, rings are
4077 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4078 #
4079 #
4080 sub GetAtomPathsStartingAt {
4081   my($This, $StartAtom, $AllowCycles) = @_;
4082 
4083   return $This->_GetAtomPathsStartingAt('AtomPathsWithAllLengths', $StartAtom, undef, $AllowCycles);
4084 }
4085 
4086 # Get atom paths as an array containing references to arrays with path atoms...
4087 #
4088 sub _GetAtomPathsStartingAt {
4089   my($This, $Mode, $StartAtom, $Length, $AllowCycles) = @_;
4090   my(@AtomPaths);
4091 
4092   @AtomPaths = ();
4093   if (!defined $StartAtom) {
4094     carp "Warning: ${ClassName}->_GetAtomPathsStartingAt: No atom paths retrieved: Start atom is not defined...";
4095     return @AtomPaths;
4096   }
4097   if (!$This->HasAtom($StartAtom)) {
4098     carp "Warning: ${ClassName}->_GetAtomPathsStartingAt: No atom paths retrieved: Start atom doesn't exist...";
4099     return @AtomPaths;
4100   }
4101   my($StartAtomID, @Paths);
4102 
4103   $StartAtomID = $StartAtom->GetID();
4104   @Paths = ();
4105 
4106   # Collect appropriate atom paths...
4107   MODE: {
4108     if ($Mode =~ /^AtomPathsWithLengthUpto$/i) { @Paths = $This->GetPathsStartingAtWithLengthUpto($StartAtomID, $Length, $AllowCycles); last MODE; }
4109     if ($Mode =~ /^AtomPathsWithLength$/i) { @Paths = $This->GetPathsStartingAtWithLength($StartAtomID, $Length, $AllowCycles); last MODE; }
4110     if ($Mode =~ /^AtomPathsWithAllLengths$/i) { @Paths = $This->GetPathsStartingAt($StartAtomID, $AllowCycles); last MODE; }
4111 
4112     if ($Mode =~ /^AllAtomPathsWithLengthUpto$/i) { @Paths = $This->GetAllPathsStartingAtWithLengthUpto($StartAtomID, $Length, $AllowCycles); last MODE; }
4113     if ($Mode =~ /^AllAtomPathsWithLength$/i) { @Paths = $This->GetAllPathsStartingAtWithLength($StartAtomID, $Length, $AllowCycles); last MODE; }
4114     if ($Mode =~ /^AllAtomPathsWithAllLengths$/i) { @Paths = $This->GetAllPathsStartingAt($StartAtomID, $AllowCycles); last MODE; }
4115 
4116     print "Warn: ${ClassName}->_GetAtomPathsStartingAt: No atom paths retrieved: Mode, $Mode, is not supported...";
4117     return @AtomPaths;
4118   }
4119   return $This->_GetAtomPathsFromPaths(\@Paths);
4120 }
4121 
4122 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4123 # path atoms.
4124 #
4125 # Path atoms correspond to to all possible paths for each atom in molecule with length
4126 # upto a specified length and sharing of bonds in paths traversed. By default, rings are
4127 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4128 #
4129 # Notes:
4130 #    . For molecule without any rings, this method returns the same set of atom paths
4131 #      as GetAtomPathsWithLengthUpto method.
4132 #
4133 sub GetAllAtomPathsWithLengthUpto {
4134   my($This, $Length, $AllowCycles) = @_;
4135 
4136   return $This->_GetAtomPaths('AllAtomPathsWithLengthUpto', $Length, $AllowCycles);
4137 }
4138 
4139 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4140 # path atoms.
4141 #
4142 # Path atoms correspond to to all possible paths for each atom in molecule with
4143 # a specified length and sharing of bonds in paths traversed. By default, rings are
4144 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4145 #
4146 # Notes:
4147 #    . For molecule without any rings, this method returns the same set of atom paths
4148 #      as GetAtomPathsWithLengthUpto method.
4149 #
4150 sub GetAllAtomPathsWithLength {
4151   my($This, $Length, $AllowCycles) = @_;
4152 
4153   return $This->_GetAtomPaths('AllAtomPathsWithLength', $Length, $AllowCycles);
4154 }
4155 
4156 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4157 # path atoms.
4158 #
4159 # Path atoms correspond to to all possible paths for each atom in molecule with all
4160 # possible lengths and sharing of bonds in paths traversed. By default, rings are
4161 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4162 #
4163 # Notes:
4164 #    . For molecule without any rings, this method returns the same set of atom paths
4165 #      as GetAtomPaths method.
4166 #
4167 sub GetAllAtomPaths {
4168   my($This, $AllowCycles) = @_;
4169 
4170   return $This->_GetAtomPaths('AllAtomPathsWithAllLengths', undef, $AllowCycles);
4171 }
4172 
4173 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4174 # path atoms.
4175 #
4176 # Path atoms correspond to to all possible paths for each atom in molecule with length
4177 # upto a specified length and no sharing of bonds in paths traversed. By default, rings are
4178 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4179 #
4180 sub GetAtomPathsWithLengthUpto {
4181   my($This, $Length, $AllowCycles) = @_;
4182 
4183   return $This->_GetAtomPaths('AtomPathsWithLengthUpto', $Length, $AllowCycles);
4184 }
4185 
4186 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4187 # path atoms.
4188 #
4189 # Path atoms correspond to to all possible paths for each atom in molecule with
4190 # a specified length and no sharing of bonds in paths traversed. By default, rings are
4191 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4192 #
4193 sub GetAtomPathsWithLength {
4194   my($This, $Length, $AllowCycles) = @_;
4195 
4196   return $This->_GetAtomPaths('AtomPathsWithLength', $Length, $AllowCycles);
4197 }
4198 
4199 
4200 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4201 # path atoms.
4202 #
4203 # Path atoms correspond to to all possible paths for each atom in molecule with all
4204 # possible lengths and no sharing of bonds in paths traversed. By default, rings are
4205 # included in paths. A path containing a ring is terminated at an atom completing the ring.
4206 #
4207 sub GetAtomPaths {
4208   my($This, $AllowCycles) = @_;
4209 
4210   return $This->_GetAtomPaths('AtomPathsWithAllLengths', undef, $AllowCycles);
4211 }
4212 
4213 # Get atom paths for all atoms as a reference to an array containing references to arrays with
4214 # path atoms.
4215 #
4216 sub _GetAtomPaths {
4217   my($This, $Mode, $Length, $AllowCycles) = @_;
4218   my($PathsRef, @AtomPaths);
4219 
4220   @AtomPaths = ();
4221   # Collect appropriate atom paths...
4222   MODE: {
4223     if ($Mode =~ /^AtomPathsWithLengthUpto$/i) { $PathsRef = $This->GetPathsWithLengthUpto($Length, $AllowCycles); last MODE; }
4224     if ($Mode =~ /^AtomPathsWithLength$/i) { $PathsRef = $This->GetPathsWithLength($Length, $AllowCycles); last MODE; }
4225     if ($Mode =~ /^AtomPathsWithAllLengths$/i) { $PathsRef = $This->GetPaths($AllowCycles); last MODE; }
4226 
4227     if ($Mode =~ /^AllAtomPathsWithLengthUpto$/i) { $PathsRef = $This->GetAllPathsWithLengthUpto($Length, $AllowCycles); last MODE; }
4228     if ($Mode =~ /^AllAtomPathsWithLength$/i) { $PathsRef = $This->GetAllPathsWithLength($Length, $AllowCycles); last MODE; }
4229     if ($Mode =~ /^AllAtomPathsWithAllLengths$/i) { $PathsRef = $This->GetAllPaths($AllowCycles); last MODE; }
4230 
4231     print "Warn: ${ClassName}->_GetAtomPaths: No atom paths retrieved: Mode, $Mode, is not supported...";
4232     return \@AtomPaths;
4233   }
4234   return $This->_GetAtomPathsFromPaths($PathsRef);
4235 }
4236 
4237 # Get atom paths as an array reference containing references to arrays with path atoms...
4238 #
4239 sub _GetAtomPathsFromPaths {
4240   my($This, $PathsRef) = @_;
4241   my($Path, @AtomPaths);
4242 
4243   @AtomPaths = ();
4244   if (!defined $PathsRef) {
4245     return \@AtomPaths;
4246   }
4247   if (!@{$PathsRef}) {
4248     # Return an empty atom paths list...
4249     return \@AtomPaths;
4250   }
4251   for $Path (@{$PathsRef}) {
4252     my(@PathAtoms);
4253     @PathAtoms = ();
4254     @PathAtoms = $This->_GetAtomPathFromPath($Path);
4255 
4256     push @AtomPaths, \@PathAtoms;
4257   }
4258   return \@AtomPaths;
4259 }
4260 
4261 # Generate an array of bond objects for an array of path atoms and return an array
4262 # of bond objects...
4263 #
4264 sub GetAtomPathBonds {
4265   my($This, @PathAtoms) = @_;
4266   my(@Bonds);
4267 
4268   if (!@PathAtoms) {
4269     # Return an empty ring bonds list...
4270     return @Bonds;
4271   }
4272   my(@PathAtomIDs);
4273 
4274   @PathAtomIDs = ();
4275   @PathAtomIDs = $This->_GetAtomsIDsFromAtoms(@PathAtoms);
4276 
4277   return $This->_GetPathBonds(@PathAtomIDs);
4278 }
4279 
4280 # Map atom IDs in path to atoms and return a reference to an array containing ring atoms...
4281 #
4282 sub _GetAtomPathFromPath {
4283   my($This, $Path) = @_;
4284   my(@PathAtoms);
4285 
4286   @PathAtoms = ();
4287   if (!defined $Path) {
4288     # Return an empty atoms list...
4289     return @PathAtoms;
4290   }
4291 
4292   return $This->_GetPathAtoms($Path);
4293 }
4294 
4295 # Get atom paths between two specified atoms as a reference to an array containing references
4296 # to arrays with path atoms. For molecules with rings, atom paths array contains may contain
4297 # two paths.
4298 #
4299 sub GetAtomPathsBetween {
4300   my($This, $StartAtom, $EndAtom) = @_;
4301   my(@AtomPaths);
4302 
4303   @AtomPaths = ();
4304   if (!(defined($StartAtom) && $This->HasAtom($StartAtom))) {
4305     carp "Warning: ${ClassName}->_GetAtomPathsBetween: No atom paths retrieved: Start atom is not defined  or it doesn't exist...";
4306     return @AtomPaths;
4307   }
4308   if (!(defined($EndAtom) && $This->HasAtom($EndAtom))) {
4309     carp "Warning: ${ClassName}->_GetAtomPathsBetween: No atom paths retrieved: End atom is not defined  or it doesn't exist...";
4310     return @AtomPaths;
4311   }
4312   return $This->_GetAtomPathsBetween($StartAtom, $EndAtom);
4313 }
4314 
4315 # Get atom paths between two specified atoms as a reference to an array containing references
4316 # to arrays with path atoms.
4317 #
4318 sub _GetAtomPathsBetween {
4319   my($This, $StartAtom, $EndAtom) = @_;
4320   my($StartAtomID, $EndAtomID, @Paths);
4321 
4322   $StartAtomID = $StartAtom->GetID();
4323   $EndAtomID = $EndAtom->GetID();
4324 
4325   @Paths = ();
4326   @Paths = $This->GetPathsBetween($StartAtomID, $EndAtomID);
4327 
4328   return $This->_GetAtomPathsFromPaths(\@Paths);
4329 }
4330 
4331 # Get atom neighborhoods around a specified atom as an array containing references
4332 # to arrays with neighborhood atoms at different radii upto specified radius...
4333 #
4334 sub GetAtomNeighborhoodsWithRadiusUpto {
4335   my($This, $StartAtom, $Radius) = @_;
4336 
4337   return $This->_GetAtomNeighborhoods('RadiusUpto', $StartAtom, $Radius);
4338 }
4339 
4340 # Get atom neighborhoods around a specified atom as an array containing references
4341 # to arrays with neighborhood atoms at possible radii...
4342 #
4343 sub GetAtomNeighborhoods {
4344   my($This, $StartAtom) = @_;
4345 
4346   return $This->_GetAtomNeighborhoods('AllRadii', $StartAtom, undef);
4347 }
4348 
4349 # Get atom neighborhood around a specified atom, along with their successor connected atoms, collected
4350 # with in a specified radius as a list containing references to lists with first value corresponding to neighborhood
4351 # atom at a specific radius and second value as reference to a list containing its successor connected atoms.
4352 #
4353 # For a neighborhood atom at each radius level, the successor connected atoms correspond to the
4354 # neighborhood atoms at the next radius level. Consequently, the neighborhood atoms at the last
4355 # radius level don't contain any successor atoms which fall outside the range of specified radius.
4356 #
4357 sub GetAtomNeighborhoodsWithSuccessorAtomsAndRadiusUpto {
4358   my($This, $StartAtom, $Radius) = @_;
4359 
4360   return $This->_GetAtomNeighborhoods('WithSuccessorsAndRadiusUpto', $StartAtom, $Radius);
4361 }
4362 
4363 # Get atom neighborhood around a specified atom, along with their successor connected atoms, collected
4364 # at all radii as a list containing references to lists with first value corresponding to neighborhood
4365 # atom at a specific radius and second value as reference to a list containing its successor connected atoms.
4366 #
4367 # For a neighborhood atom at each radius level, the successor connected atoms correspond to the
4368 # neighborhood atoms at the next radius level. Consequently, the neighborhood atoms at the last
4369 # radius level don't contain any successor atoms which fall outside the range of specified radius.
4370 #
4371 #
4372 sub GetAtomNeighborhoodsWithSuccessorAtoms {
4373   my($This, $StartAtom) = @_;
4374 
4375   return $This->_GetAtomNeighborhoods('WithSuccessorsAndAllRadii', $StartAtom, undef);
4376 }
4377 
4378 # Get atom neighborhoods...
4379 #
4380 sub _GetAtomNeighborhoods {
4381   my($This, $Mode, $StartAtom, $Radius) = @_;
4382   my(@AtomNeighborhoods);
4383 
4384   @AtomNeighborhoods = ();
4385 
4386   if (!(defined($StartAtom) && $This->HasAtom($StartAtom))) {
4387     carp "Warning: ${ClassName}->_GetAtomNeighborhoods: No atom neighborhoods retrieved: Start atom is not defined  or it doesn't exist...";
4388     return @AtomNeighborhoods;
4389   }
4390   if ($Mode =~ /^(RadiusUpto|WithSuccessorsAndRadiusUpto)$/i) {
4391     if (!(defined($Radius) && $Radius > 0)) {
4392       carp "Warning: ${ClassName}->_GetAtomNeighborhoods: No atom neighborhoods retrieved: Radius is not defined or it's <= 0 ...";
4393       return @AtomNeighborhoods;
4394     }
4395   }
4396 
4397   # Collect neighborhood atom IDs...
4398   my($StartAtomID, @NeighborhoodAtomIDs, @NeighborhoodAtomIDsWithSuccessors);
4399 
4400   @NeighborhoodAtomIDs = (); @NeighborhoodAtomIDsWithSuccessors = ();
4401   $StartAtomID = $StartAtom->GetID();
4402 
4403   MODE: {
4404     if ($Mode =~ /^RadiusUpto$/i) { @NeighborhoodAtomIDs = $This->GetNeighborhoodVerticesWithRadiusUpto($StartAtomID, $Radius); last MODE; }
4405     if ($Mode =~ /^AllRadii$/i) { @NeighborhoodAtomIDs = $This->GetNeighborhoodVertices($StartAtomID); last MODE; }
4406 
4407     if ($Mode =~ /^WithSuccessorsAndRadiusUpto$/i) { @NeighborhoodAtomIDsWithSuccessors = $This->GetNeighborhoodVerticesWithSuccessorsAndRadiusUpto($StartAtomID, $Radius); last MODE; }
4408     if ($Mode =~ /^WithSuccessorsAndAllRadii$/i) { @NeighborhoodAtomIDsWithSuccessors = $This->GetNeighborhoodVerticesWithSuccessors($StartAtomID); last MODE; }
4409 
4410     print "Warn: ${ClassName}->_GetAtomNeighborhood: No atom neighborhoods retrieved: Mode, $Mode, is not supported...";
4411     return @AtomNeighborhoods;
4412   }
4413   if ($Mode =~ /^(RadiusUpto|AllRadii)$/i) {
4414     return $This->_GetNeighborhoodAtomsFromAtomIDs(\@NeighborhoodAtomIDs);
4415   }
4416   elsif ($Mode =~ /^(WithSuccessorsAndRadiusUpto|WithSuccessorsAndAllRadii)$/i) {
4417     return $This->_GetNeighborhoodAtomsWithSuccessorsFromAtomIDs(\@NeighborhoodAtomIDsWithSuccessors);
4418   }
4419 
4420   return @AtomNeighborhoods;
4421 }
4422 
4423 # Map neighborhood atom IDs to atoms...
4424 #
4425 sub _GetNeighborhoodAtomsFromAtomIDs {
4426   my($This, $NeighborhoodsAtomIDsRef) = @_;
4427   my($NeighborhoodAtomIDsRef, @AtomNeighborhoods);
4428 
4429   @AtomNeighborhoods = ();
4430   for $NeighborhoodAtomIDsRef (@{$NeighborhoodsAtomIDsRef}) {
4431     my(@AtomNeighborhood);
4432 
4433     @AtomNeighborhood = ();
4434     @AtomNeighborhood = $This->_GetAtomsFromAtomIDs(@{$NeighborhoodAtomIDsRef});
4435     push @AtomNeighborhoods, \@AtomNeighborhood;
4436   }
4437   return @AtomNeighborhoods;
4438 }
4439 
4440 # Map neighborhood atom IDs with successors to atoms...
4441 #
4442 sub _GetNeighborhoodAtomsWithSuccessorsFromAtomIDs {
4443   my($This, $NeighborhoodsAtomIDsWithSuccessorsRef) = @_;
4444   my($Depth, $NeighborhoodAtomIDsWithSuccessorsRef, $NeighborhoodAtomIDWithSuccessorsRef, $NeighborhoodAtomID, $NeighborhoodAtomSuccessorsIDsRef, @AtomNeighborhoods);
4445 
4446   $Depth = 0;
4447   @AtomNeighborhoods = ();
4448 
4449   # Go over neighborhoods at each level...
4450   for $NeighborhoodAtomIDsWithSuccessorsRef (@{$NeighborhoodsAtomIDsWithSuccessorsRef}) {
4451     @{$AtomNeighborhoods[$Depth]} = ();
4452 
4453     # Go over the neighborhood atoms and their successors at a specific level..
4454     for $NeighborhoodAtomIDWithSuccessorsRef (@{$NeighborhoodAtomIDsWithSuccessorsRef}) {
4455       my($NeighborhoodAtom, @NeighborhoodAtomWithSuccessors, @NeighborhoodAtomSuccessorAtoms);
4456 
4457       @NeighborhoodAtomWithSuccessors = (); @NeighborhoodAtomSuccessorAtoms = ();
4458       ($NeighborhoodAtomID, $NeighborhoodAtomSuccessorsIDsRef) = @{$NeighborhoodAtomIDWithSuccessorsRef};
4459 
4460       # Map atom IDs to atoms...
4461       $NeighborhoodAtom = $This->_GetAtomFromAtomID($NeighborhoodAtomID);
4462       if (@{$NeighborhoodAtomSuccessorsIDsRef}) {
4463         @NeighborhoodAtomSuccessorAtoms = $This->_GetAtomsFromAtomIDs(@{$NeighborhoodAtomSuccessorsIDsRef});
4464       }
4465 
4466       # Store an atom and its successors at each level in an array...
4467       push @NeighborhoodAtomWithSuccessors, ($NeighborhoodAtom, \@NeighborhoodAtomSuccessorAtoms);
4468 
4469       push @{$AtomNeighborhoods[$Depth]} , \@NeighborhoodAtomWithSuccessors;
4470     }
4471       $Depth++;
4472   }
4473   return @AtomNeighborhoods;
4474 }
4475 
4476 # Get next object ID...
4477 sub _GetNewObjectID {
4478   $ObjectID++;
4479   return $ObjectID;
4480 }
4481 
4482 # Is aromatic property set for the molecule?
4483 sub IsAromatic {
4484   my($This) = @_;
4485   my($Aromatic);
4486 
4487   $Aromatic = $This->GetAromatic();
4488 
4489   return (defined($Aromatic) && $Aromatic) ? 1 : 0;
4490 }
4491 
4492 # Does molecule contains any atoms with non-zero Z coordiantes?
4493 sub IsThreeDimensional {
4494   my($This) = @_;
4495   my($Atom, @Atoms);
4496 
4497   @Atoms = $This->GetAtoms();
4498   ATOM: for $Atom (@Atoms) {
4499       if ($Atom->GetZ() != 0) {
4500         return 1;
4501       }
4502   }
4503   return 0;
4504 }
4505 
4506 # Does molecule contains any atoms with non-zero X or Y coordinates
4507 # and only zero Z-coordinates?
4508 sub IsTwoDimensional {
4509   my($This) = @_;
4510   my($Atom, @Atoms);
4511 
4512   @Atoms = $This->GetAtoms();
4513   ATOM: for $Atom (@Atoms) {
4514       if ($Atom->GetZ() != 0) {
4515         return 0;
4516       }
4517       if ($Atom->GetX() != 0 || $Atom->GetY() != 0) {
4518         return 1;
4519       }
4520   }
4521   return 0;
4522 }
4523 
4524 # Get dimensionality of the molecule using one of the following two methods:
4525 #   . Using explicitly set Dimensionality
4526 #   . Going over atomic coordinates
4527 #
4528 # The valid dimensionality values are:
4529 #   . 3D - Three dimensional: One of X, Y or Z coordinate is non-zero
4530 #   . 2D - Two dimensional: One of X or Y coordinate is non-zero; All Z coordinates are zero
4531 #   . 0D - Zero dimensional: All atomic coordinates are zero
4532 #
4533 sub GetDimensionality {
4534   my($This) = @_;
4535 
4536   # Is Dimensionality property explicitly set?
4537   if ($This->HasProperty('Dimensionality')) {
4538     return $This->GetProperty('Dimensionality');
4539   }
4540   my($Atom, @Atoms);
4541 
4542   @Atoms = $This->GetAtoms();
4543   ATOM: for $Atom (@Atoms) {
4544       if ($Atom->GetZ() != 0) {
4545         return '3D';
4546       }
4547       if ($Atom->GetX() != 0 || $Atom->GetY() != 0) {
4548         return '2D';
4549       }
4550   }
4551   return '0D';
4552 }
4553 
4554 # Is it a molecule object?
4555 sub IsMolecule ($) {
4556   my($Object) = @_;
4557 
4558   return _IsMolecule($Object);
4559 }
4560 
4561 # Return a string containing vertices, edges and other properties...
4562 sub StringifyMolecule {
4563   my($This) = @_;
4564   my($MoleculeString, $ID, $Name, $NumOfAtoms, $NumOfBonds, $MolecularFormula, $NumOfRings, $MolecularWeight, $ExactMass, $FormalCharge, $SpinMultiplicity, $FreeRadicalElectrons, $Charge, $ElementsRef, $ElementsCompositionRef, $ElementalComposition);
4565 
4566   $ID = $This->GetID();
4567   $Name = $This->GetName();
4568   $NumOfAtoms = $This->GetNumOfAtoms();
4569   $NumOfBonds = $This->GetNumOfBonds();
4570 
4571   $NumOfRings = $This->GetNumOfRings();
4572   if (!defined $NumOfRings) {
4573     $NumOfRings = 'undefined';
4574   }
4575 
4576   $MolecularFormula = $This->GetMolecularFormula();
4577 
4578   $MolecularWeight = $This->GetMolecularWeight();
4579   $MolecularWeight = round($MolecularWeight, 4) + 0;
4580 
4581   $ExactMass = $This->GetExactMass();
4582   $ExactMass = round($ExactMass, 4) + 0;
4583 
4584   $FormalCharge = $This->GetFormalCharge();
4585   $Charge = $This->GetCharge();
4586 
4587   $SpinMultiplicity = $This->GetSpinMultiplicity();
4588   $FreeRadicalElectrons = $This->GetFreeRadicalElectrons();
4589 
4590   ($ElementsRef, $ElementsCompositionRef) = $This->GetElementalComposition();
4591   $ElementalComposition = 'None';
4592   if (defined($ElementsRef) && @{$ElementsRef}) {
4593     $ElementalComposition = "[ " . FormatElementalCompositionInformation($ElementsRef, $ElementsCompositionRef) . " ]";
4594   }
4595 
4596   $MoleculeString = "Molecule: ID: $ID; Name: \"$Name\"; NumOfAtoms: $NumOfAtoms; NumOfBonds: $NumOfBonds; NumOfRings: $NumOfRings; MolecularFormula: $MolecularFormula; MolecularWeight: $MolecularWeight; ExactMass: $ExactMass; FormalCharge: $FormalCharge; Charge: $Charge; SpinMultiplicity: $SpinMultiplicity; FreeRadicalElectrons: $FreeRadicalElectrons; ElementalComposition: $ElementalComposition";
4597 
4598   return $MoleculeString;
4599 }
4600 
4601 # Load appropriate atom data files from <MayaChemTools>/lib directory used by various
4602 # object methods in the current class...
4603 #
4604 sub _LoadMoleculeClassData {
4605   my($MayaChemToolsLibDir);
4606 
4607   $MayaChemToolsLibDir = GetMayaChemToolsLibDirName();
4608 
4609   # Load and process data for aromaticity models...
4610   _LoadAromaticityModelsData($MayaChemToolsLibDir);
4611   _ProcessAromaticityModelsData();
4612 }
4613 
4614 #
4615 # Load data for supported aromaticity models...
4616 #
4617 sub _LoadAromaticityModelsData {
4618   my($MayaChemToolsLibDir) = @_;
4619   my($DataFile, $Index, $InDelim, $Line, $NumOfCols, $ParameterName, $ParameterValue, $ModelName, @ColLabels, @LineWords, %ParameterNames, %ColIndexToModelName, %SupportedParameterNames);
4620 
4621   %AromaticityModelsDataMap = ();
4622   %CanonicalAromaticityModelNamesMap = ();
4623 
4624   # File format:
4625   #
4626   # "ParameterName","MDLAromaticityModel","TriposAromaticityModel","MMFFAromaticityModel","ChemAxonBasicAromaticityModel","ChemAxonGeneralAromaticityModel","DaylightAromaticityModel","MayaChemToolsAromaticityModel"
4627   # "AllowHeteroRingAtoms","No","No","Yes","Yes","Yes","Yes","Yes"
4628   #
4629   $DataFile = $MayaChemToolsLibDir . "/data/AromaticityModelsData.csv";
4630   if (! -e "$DataFile") {
4631     croak "Error: ${ClassName}::_LoadAromaticityModelsData: MayaChemTools package file, $DataFile, is missing: Possible installation problems...";
4632   }
4633 
4634   # Setup a list of currently supported aromaticity parameters...
4635   #
4636   my(@KnownNames);
4637   @KnownNames = qw(AllowHeteroRingAtoms HeteroRingAtomsList AllowExocyclicDoubleBonds AllowHomoNuclearExocyclicDoubleBonds AllowElectronegativeRingAtomExocyclicDoubleBonds AllowRingAtomFormalCharge AllowHeteroRingAtomFormalCharge MinimumRingSize);
4638 
4639   %SupportedParameterNames = ();
4640   for $ParameterName (@KnownNames) {
4641     $SupportedParameterNames{$ParameterName} = $ParameterName;
4642   }
4643 
4644   $InDelim = "\,";
4645   open DATAFILE, "$DataFile" or croak "Couldn't open $DataFile: $! ...";
4646 
4647   # Skip lines up to column labels...
4648   LINE: while ($Line = GetTextLine(\*DATAFILE)) {
4649     if ($Line !~ /^#/) {
4650       last LINE;
4651     }
4652   }
4653   @ColLabels= quotewords($InDelim, 0, $Line);
4654   $NumOfCols = @ColLabels;
4655 
4656   %ColIndexToModelName = ();
4657 
4658   # Process names of aromaticity models...
4659   for $Index (1 .. $#ColLabels) {
4660     $ModelName = $ColLabels[$Index];
4661     $ModelName =~ s/ //g;
4662 
4663     if (exists $AromaticityModelsDataMap{$ModelName}) {
4664       croak "Error: ${ClassName}::_LoadAromaticityModelsData: The aromaticity model name, $ModelName, in $DataFile has already exists.\nLine: $Line...";
4665     }
4666     %{$AromaticityModelsDataMap{$ModelName}} = ();
4667 
4668     # Cannonicalize aromatic model name by converting into all lowercase...
4669     $CanonicalAromaticityModelNamesMap{lc($ModelName)} = $ModelName;
4670 
4671     $ColIndexToModelName{$Index} = $ModelName;
4672   }
4673 
4674   # Process paramater name and their values for specified aromaticity models...
4675   #
4676   %ParameterNames = ();
4677   LINE: while ($Line = GetTextLine(\*DATAFILE)) {
4678     if ($Line =~ /^#/) {
4679       next LINE;
4680     }
4681     @LineWords = ();
4682     @LineWords = quotewords($InDelim, 0, $Line);
4683     if (@LineWords != $NumOfCols) {
4684       croak "Error: ${ClassName}::_LoadAromaticityModelsData: The number of data fields, @LineWords, in $DataFile must be $NumOfCols.\nLine: $Line...";
4685     }
4686 
4687     # Process parameter name and values for aromaticity models...
4688     #
4689     $ParameterName = $LineWords[0];
4690 
4691     if (!exists $SupportedParameterNames{$ParameterName}) {
4692       carp "Warning: ${ClassName}::_LoadAromaticityModelsData: The current release of MayaChemTools doesn't support aromaticity model parameter name, $ParameterName, specified in $DataFile. It would be ignore during aromaticity detection.\nLine: $Line...";
4693     }
4694 
4695     if (exists $ParameterNames{$ParameterName}) {
4696       carp "Warning: ${ClassName}::_LoadAromaticityModelsData: Ignoring aromaticity model data for parameter name, $ParameterName, in $DataFile. It has already been loaded.\nLine: $Line...";
4697       next LINE;
4698     }
4699     $ParameterNames{$ParameterName} = $ParameterName;
4700 
4701     for $Index (1 .. $#LineWords) {
4702       $ModelName = $ColIndexToModelName{$Index};
4703       $ParameterValue = $LineWords[$Index];
4704       $AromaticityModelsDataMap{$ModelName}{$ParameterName} = $ParameterValue;
4705     }
4706   }
4707   close DATAFILE;
4708 }
4709 
4710 # Process already loaded aromaticity model data...
4711 #
4712 sub _ProcessAromaticityModelsData {
4713   my($ParameterName, $ParameterValue, $ModelName, $NewParameterValue);
4714 
4715   for $ModelName (keys %AromaticityModelsDataMap) {
4716     for $ParameterName (keys %{$AromaticityModelsDataMap{$ModelName}}) {
4717       $ParameterValue = $AromaticityModelsDataMap{$ModelName}{$ParameterName};
4718       $ParameterValue =~ s/ //g;
4719 
4720       VALUE: {
4721         if ($ParameterValue =~ /^Yes$/i) {
4722           $NewParameterValue = 1;
4723           last VALUE;
4724         }
4725         if ($ParameterValue =~ /^(NA|No)$/i) {
4726           $NewParameterValue = 0;
4727           last VALUE;
4728         }
4729         if ($ParameterValue =~ /^None$/i) {
4730           $NewParameterValue = '';
4731           last VALUE;
4732         }
4733         $NewParameterValue = $ParameterValue;
4734       }
4735       $AromaticityModelsDataMap{$ModelName}{$ParameterName} = $NewParameterValue;
4736 
4737       if ($ParameterName =~ /List/i) {
4738         # Setup a new parameter conatining a reference to a hash for the specified values...
4739         my($DataMapRefName, $DataValue, %DataMap);
4740 
4741         $DataMapRefName = "${ParameterName}MapRef";
4742 
4743         %DataMap = ();
4744         for $DataValue (split /\,/, $NewParameterValue) {
4745           $DataMap{$DataValue} = $DataValue;
4746         }
4747         $AromaticityModelsDataMap{$ModelName}{$DataMapRefName} = \%DataMap;
4748       }
4749     }
4750   }
4751 }
4752 
4753 # Is it a molecule object?
4754 sub _IsMolecule {
4755   my($Object) = @_;
4756 
4757   return (Scalar::Util::blessed($Object) && $Object->isa($ClassName)) ? 1 : 0;
4758 }
4759