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