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