       1 package Atom;
       2 #
       3 # File: Atom.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 PeriodicTable;
      33 use Vector;
      34 use MathUtil;
      35 use Text::ParseWords;
      36 use TextUtil;
      37 use FileUtil;
      38
      39 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
      40
      41 @ISA = qw(ObjectProperty Exporter);
      42 @EXPORT = qw();
      43 @EXPORT_OK = qw();
      44
      45 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
      46
      47 # Setup class variables...
      48 my($ClassName, $ObjectID, %MDLValenceModelDataMap, %DaylightValenceModelDataMap);
      49 _InitializeClass();
      50
      51 # Overload Perl functions...
      52 use overload '""' => 'StringifyAtom';
      53
      54 # Class constructor...
      55 sub new {
      56   my($Class, %NamesAndValues) = @_;
      57
      58   # Initialize object...
      59   my $This = {};
      60   bless $This, ref($Class) || $Class;
      61   $This->_InitializeAtom();
      62
      63   $This->_InitializeAtomProperties(%NamesAndValues);
      64
      65   return $This;
      66 }
      67
      68 # Initialize object data...
      69 sub _InitializeAtom {
      70   my($This) = @_;
      71   my($ObjectID) = _GetNewObjectID();
      72
      73   # All other property names and values along with all Set/Get<PropertyName> methods
      74   # are implemented on-demand using ObjectProperty class.
      75   $This->{ID} = $ObjectID;
      76   $This->{Name} = "Atom ${ObjectID}";
      77   $This->{AtomSymbol} = '';
      78   $This->{AtomicNumber} = 0;
      79   $This->{XYZ} = Vector::ZeroVector;
      80 }
      81
      82 # Initialize atom properties...
      83 sub _InitializeAtomProperties {
      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   if (!exists $NamesAndValues{'AtomSymbol'}) {
      92     carp "Warning: ${ClassName}->new: Atom object instantiated without setting atom symbol...";
      93   }
      94
      95   return $This;
      96 }
      97
      98 # Initialize class ...
      99 sub _InitializeClass {
     100   #Class name...
     101   $ClassName = __PACKAGE__;
     102
     103   # ID to keep track of objects...
     104   $ObjectID = 0;
     105
     106   # Load atom class data...
     107   _LoadAtomClassData();
     108 }
     109
     110 # Setup an explicit SetID method to block setting of ID by AUTOLOAD function...
     111 sub SetID {
     112   my($This, $Value) = @_;
     113
     114   carp "Warning: ${ClassName}->SetID: Object ID can't be changed: it's used for internal tracking...";
     115
     116   return $This;
     117 }
     118
     119 # Setup an explicit SetMolecule method to block setting of ID by AUTOLOAD function...
     120 sub SetMolecule {
     121   my($This, $Value) = @_;
     122
     123   carp "Warning: ${ClassName}->SetMolecule: Molecule property can't be changed: it's used for internal tracking...";
     124
     125   return $This;
     126 }
     127
     128 # Assign atom to  molecule...
     129 sub _SetMolecule {
     130   my($This, $Molecule) = @_;
     131
     132   $This->{Molecule} = $Molecule;
     133
     134   # Weaken the reference to disable increment of reference count; otherwise,
     135   # it it becomes a circular reference and destruction of Molecule object doesn't
     136   # get initiated which in turn disables destruction of atom object.
     137   #
     138   Scalar::Util::weaken($This->{Molecule});
     139
     140   return $This;
     141 }
     142
     143 # Setup atom symbol and atomic number for the element...
     144 #
     145 # Possible atom symbol values:
     146 #    . An element symbol or some other type of atom: L - Atom list; LP - Lone pair; R# - R group;
     147 #       A, Q, * - unknown atom; or something else?
     148 #
     149 # Default mass number corresponds to the most abundant natural isotope unless it's explicity
     150 # set using "MassNumber" property.
     151 #
     152 sub SetAtomSymbol {
     153   my($This, $AtomSymbol) = @_;
     154   my($AtomicNumber);
     155
     156   $This->{AtomSymbol} = $AtomSymbol;
     157
     158   $AtomicNumber = PeriodicTable::GetElementAtomicNumber($AtomSymbol);
     159   $This->{AtomicNumber} = (defined $AtomicNumber) ? $AtomicNumber : 0;
     160
     161   return $This;
     162 }
     163
     164 # Setup atom symbol and atomic number for the element...
     165 sub SetAtomicNumber {
     166   my($This, $AtomicNumber) = @_;
     167   my($AtomSymbol);
     168
     169   $AtomSymbol = PeriodicTable::GetElementAtomSymbol($AtomicNumber);
     170   if (!defined $AtomSymbol) {
     171     carp "Warning: ${ClassName}->SetAtomicNumber: Didn't set atomic number: Invalid atomic number, $AtomicNumber, specified...";
     172     return;
     173   }
     174   $This->{AtomicNumber} = $AtomicNumber;
     175   $This->{AtomSymbol} = $AtomSymbol;
     176
     177   return $This;
     178 }
     179
     180 # Set atom as stereo center...
     181 #
     182 sub SetStereoCenter {
     183   my($This, $StereoCenter) = @_;
     184
     185   $This->SetProperty('StereoCenter', $StereoCenter);
     186
     187   return $This;
     188 }
     189
     190 # Is it a stereo center?
     191 #
     192 sub IsStereoCenter {
     193   my($This) = @_;
     194   my($StereoCenter);
     195
     196   $StereoCenter = $This->GetProperty('StereoCenter');
     197
     198   return (defined($StereoCenter) && $StereoCenter) ? 1 : 0;
     199 }
     200
     201 # Set atom stereochemistry.
     202 #
     203 # Supported values are: R, S.
     204 #
     205 # Notes:
     206 #
     207 # . After the ligands around a central stereocenter has been ranked using CIP priority scheme and
     208 # the lowest ranked ligand lies behind the center atom, then R and S values correspond to:
     209 #
     210 # R: Clockwise arrangement of remaining ligands around the central atom going from highest to lowest ranked ligand
     211 # S: CounterClockwise arrangement of remaining ligands around the central atom going from highest to lowest ranked ligand
     212 #
     213 # . Assignment of any other arbitray values besides R and S is also allowed; however, a warning is printed.
     214 #
     215 sub SetStereochemistry {
     216   my($This, $Stereochemistry) = @_;
     217
     218   if ($Stereochemistry !~ /^(R|S)$/i) {
     219     carp "Warning: ${ClassName}->SetStereochemistry: Assigning non-supported Stereochemistry value of $Stereochemistry. Supported values: R, S...";
     220   }
     221
     222   $This->SetProperty('StereoCenter', 1);
     223   $This->SetProperty('Stereochemistry', $Stereochemistry);
     224
     225   return $This;
     226 }
     227
     228 # Setup mass number for atom...
     229 sub SetMassNumber {
     230   my($This, $MassNumber) = @_;
     231   my($AtomicNumber, $AtomSymbol);
     232
     233   $AtomicNumber = $This->{AtomicNumber};
     234   $AtomSymbol = $This->{AtomSymbol};
     235   if (!$AtomicNumber) {
     236     carp "Warning: ${ClassName}->SetMassNumber: Didn't set mass number: Non standard atom with atomic number, $AtomicNumber, and atomic symbol, $AtomSymbol...";
     237     return;
     238   }
     239   if (!PeriodicTable::IsElementNaturalIsotopeMassNumber($AtomicNumber, $MassNumber)) {
     240     carp "Warning: ${ClassName}->SetMassNumber: Unknown mass number, $MassNumber, specified for atom with atomic number, $AtomicNumber, and atomic symbol, $AtomSymbol. Don't forget to Set ExactMass property explicitly; otherwise, GetExactMass method would return mass of most abundant isotope...";
     241   }
     242   $This->SetProperty('MassNumber', $MassNumber);
     243
     244   return $This;
     245 }
     246
     247 # Get mass number...
     248 #
     249 sub GetMassNumber {
     250   my($This) = @_;
     251
     252   # Is mass number explicity set?
     253   if ($This->HasProperty('MassNumber')) {
     254     return $This->GetProperty('MassNumber');
     255   }
     256
     257   # Is it an element symbol?
     258   my($AtomicNumber) = $This->{AtomicNumber};
     259   if (!$AtomicNumber) {
     260     return 0;
     261   }
     262
     263   # Return most abundant mass number...
     264   return PeriodicTable::GetElementMostAbundantNaturalIsotopeMassNumber($AtomicNumber);
     265 }
     266
     267 # Get atomic weight:
     268 #   . Explicitly set by the caller
     269 #   . Using atomic number
     270 #
     271 sub GetAtomicWeight {
     272   my($This) = @_;
     273
     274   # Is atomic weight explicity set?
     275   if ($This->HasProperty('AtomicWeight')) {
     276     return $This->GetProperty('AtomicWeight');
     277   }
     278
     279   # Is it an element symbol?
     280   my($AtomicNumber) = $This->{AtomicNumber};
     281   if (!$AtomicNumber) {
     282     return 0;
     283   }
     284
     285   # Return its atomic weight...
     286   return PeriodicTable::GetElementAtomicWeight($AtomicNumber);
     287 }
     288
     289 # Get exact mass weight:
     290 #   . Explicitly set by the caller
     291 #   . Using atomic number and mass number explicity set by the caller
     292 #   . Using atomic number and most abundant isotope
     293 #
     294 sub GetExactMass {
     295   my($This) = @_;
     296
     297   # Is exact mass explicity set?
     298   if ($This->HasProperty('ExactMass')) {
     299     return $This->GetProperty('ExactMass');
     300   }
     301
     302   # Is it an element symbol?
     303   my($AtomicNumber) = $This->{AtomicNumber};
     304   if (!$AtomicNumber) {
     305     return 0;
     306   }
     307
     308   # Is mass number explicitly set?
     309   if ($This->HasProperty('MassNumber')) {
     310     my($MassNumber) = $This->GetProperty('MassNumber');
     311     if (PeriodicTable::IsElementNaturalIsotopeMassNumber($AtomicNumber, $MassNumber)) {
     312       return PeriodicTable::GetElementNaturalIsotopeMass($AtomicNumber, $MassNumber);
     313     }
     314   }
     315
     316   # Return most abundant isotope mass...
     317   return PeriodicTable::GetElementMostAbundantNaturalIsotopeMass($AtomicNumber);
     318 }
     319
     320 # Get formal charge:
     321 #   . Explicitly set by the caller
     322 #   . Or return zero insetad of undef
     323 #
     324 sub GetFormalCharge {
     325   my($This) = @_;
     326   my($FormalCharge);
     327
     328   $FormalCharge = 0;
     329   if ($This->HasProperty('FormalCharge')) {
     330     $FormalCharge = $This->GetProperty('FormalCharge');
     331   }
     332
     333   return defined($FormalCharge) ? $FormalCharge : 0;
     334 }
     335
     336 # Get spin multiplicity:
     337 #   . Explicitly set by the caller
     338 #   . From FreeRadicalElectrons value explicitly set by the caller
     339 #   . Or return zero insetad of undef
     340 #
     341 sub GetSpinMultiplicity {
     342   my($This) = @_;
     343   my($SpinMultiplicity);
     344
     345   $SpinMultiplicity = 0;
     346   if ($This->HasProperty('SpinMultiplicity')) {
     347     $SpinMultiplicity = $This->GetProperty('SpinMultiplicity');
     348     return defined($SpinMultiplicity) ? $SpinMultiplicity : 0;
     349   }
     350
     351   if ($This->HasProperty('FreeRadicalElectrons')) {
     352     my($FreeRadicalElectrons);
     353     $FreeRadicalElectrons = $This->GetProperty('FreeRadicalElectrons');
     354
     355     SPINMULTIPLICITY: {
     356       if ($FreeRadicalElectrons == 1) { $SpinMultiplicity = 2; last SPINMULTIPLICITY;}
     357       if ($FreeRadicalElectrons == 2) { $SpinMultiplicity = 1; last SPINMULTIPLICITY;}
     358       carp "Warning: ${ClassName}->GetSpinMultiplicity: It's not possible to determine spin multiplicity from the specified free radical electrons value, $FreeRadicalElectrons. It has been set to 0...";
     359       $SpinMultiplicity = 0;
     360     }
     361   }
     362
     363   return $SpinMultiplicity;
     364 }
     365
     366 # Get number of free radical electrons:
     367 #   . Explicitly set by the caller
     368 #   . From SpinMultiplicity value explicitly set by the caller
     369 #   . Or return zero insetad of undef
     370 #
     371 # Notes:
     372 #  . For atoms with explicit assignment of SpinMultiplicity property values corresponding to
     373 #    Singlet (two unpaired electrons corresponding to one spin state), Doublet (free radical; an unpaired
     374 #    electron corresponding to two spin states), and Triplet (two unparied electrons corresponding to
     375 #    three spin states; divalent carbon atoms (carbenes)), FreeRadicalElectrons are calculated as follows:
     376 #
     377 #       SpinMultiplicity: Doublet(2); FreeRadicalElectrons: 1
     378 #       SpinMultiplicity: Singlet(1)/Triplet(3); FreeRadicalElectrons: 2
     379 #
     380 sub GetFreeRadicalElectrons {
     381   my($This) = @_;
     382   my($FreeRadicalElectrons);
     383
     384   $FreeRadicalElectrons = 0;
     385
     386   if ($This->HasProperty('FreeRadicalElectrons')) {
     387     $FreeRadicalElectrons = $This->GetProperty('FreeRadicalElectrons');
     388     return defined($FreeRadicalElectrons) ? $FreeRadicalElectrons : 0;
     389   }
     390
     391   if ($This->HasProperty('SpinMultiplicity')) {
     392     my($SpinMultiplicity);
     393     $SpinMultiplicity = $This->GetProperty('SpinMultiplicity');
     394
     395     SPINMULTIPLICITY: {
     396       if ($SpinMultiplicity == 1) { $FreeRadicalElectrons = 2; last SPINMULTIPLICITY;}
     397       if ($SpinMultiplicity == 2) { $FreeRadicalElectrons = 1; last SPINMULTIPLICITY;}
     398       if ($SpinMultiplicity == 3) { $FreeRadicalElectrons = 2; last SPINMULTIPLICITY;}
     399       carp "Warning: ${ClassName}->GetFreeRadicalElectrons: It's not possible to determine free radical electrons from the specified spin multiplicity value, $FreeRadicalElectrons. It has been set to 0...";
     400       $FreeRadicalElectrons = 0;
     401     }
     402   }
     403
     404   return $FreeRadicalElectrons;
     405 }
     406
     407 # Set atom coordinates using:
     408 # . An array reference with three values
     409 # . An array containg three values
     410 # . A 3D vector
     411 #
     412 sub SetXYZ {
     413   my($This, @Values) = @_;
     414
     415   if (!@Values) {
     416     carp "Warning: ${ClassName}->SetXYZ: No values specified...";
     417     return;
     418   }
     419
     420   $This->{XYZ}->SetXYZ(@Values);
     421   return $This;
     422 }
     423
     424 # Set X value...
     425 sub SetX {
     426   my($This, $Value) = @_;
     427
     428   if (!defined $Value) {
     429     carp "Warning: ${ClassName}->SetX: Undefined X value...";
     430     return;
     431   }
     432   $This->{XYZ}->SetX($Value);
     433   return $This;
     434 }
     435
     436 # Set Y value...
     437 sub SetY {
     438   my($This, $Value) = @_;
     439
     440   if (!defined $Value) {
     441     carp "Warning: ${ClassName}->SetY: Undefined Y value...";
     442     return;
     443   }
     444   $This->{XYZ}->SetY($Value);
     445   return $This;
     446 }
     447
     448 # Set Z value...
     449 sub SetZ {
     450   my($This, $Value) = @_;
     451
     452   if (!defined $Value) {
     453     carp "Warning: ${ClassName}->SetZ: Undefined Z value...";
     454     return;
     455   }
     456   $This->{XYZ}->SetZ($Value);
     457   return $This;
     458 }
     459
     460 # Return XYZ as:
     461 # . Reference to an array
     462 # . An array
     463 #
     464 sub GetXYZ {
     465   my($This) = @_;
     466
     467   return $This->{XYZ}->GetXYZ();
     468 }
     469
     470 # Return XYZ as a vector object...
     471 #
     472 sub GetXYZVector {
     473   my($This) = @_;
     474
     475   return $This->{XYZ};
     476 }
     477
     478 # Get X value...
     479 sub GetX {
     480   my($This) = @_;
     481
     482   return $This->{XYZ}->GetX();
     483 }
     484
     485 # Get Y value...
     486 sub GetY {
     487   my($This) = @_;
     488
     489   return $This->{XYZ}->GetY();
     490 }
     491
     492 # Get Z value...
     493 sub GetZ {
     494   my($This) = @_;
     495
     496   return $This->{XYZ}->GetZ();
     497 }
     498
     499 # Delete atom...
     500 sub DeleteAtom {
     501   my($This) = @_;
     502
     503   # Is this atom in a molecule?
     504   if (!$This->HasProperty('Molecule')) {
     505     # Nothing to do...
     506     return $This;
     507   }
     508   my($Molecule) = $This->GetProperty('Molecule');
     509
     510   return $Molecule->_DeleteAtom($This);
     511 }
     512
     513 # Get atom neighbor objects as array. In scalar conetxt, return number of neighbors...
     514 sub GetNeighbors {
     515   my($This, @ExcludeNeighbors) = @_;
     516
     517   # Is this atom in a molecule?
     518   if (!$This->HasProperty('Molecule')) {
     519     return undef;
     520   }
     521   my($Molecule) = $This->GetProperty('Molecule');
     522
     523   if (@ExcludeNeighbors) {
     524     return $This->_GetAtomNeighbors(@ExcludeNeighbors);
     525   }
     526   else {
     527     return $This->_GetAtomNeighbors();
     528   }
     529 }
     530
     531 # Get atom neighbor objects as array. In scalar conetxt, return number of neighbors...
     532 sub _GetAtomNeighbors {
     533   my($This, @ExcludeNeighbors) = @_;
     534   my($Molecule) = $This->GetProperty('Molecule');
     535
     536   if (!@ExcludeNeighbors) {
     537     return $Molecule->_GetAtomNeighbors($This);
     538   }
     539
     540   # Setup a map for neigbhors to exclude...
     541   my($ExcludeNeighbor, $ExcludeNeighborID, %ExcludeNeighborsIDsMap);
     542
     543   %ExcludeNeighborsIDsMap = ();
     544   for $ExcludeNeighbor (@ExcludeNeighbors) {
     545     $ExcludeNeighborID = $ExcludeNeighbor->GetID();
     546     $ExcludeNeighborsIDsMap{$ExcludeNeighborID} = $ExcludeNeighborID;
     547   }
     548
     549   # Generate a filtered neighbors list...
     550   my($Neighbor, $NeighborID, @FilteredAtomNeighbors);
     551   @FilteredAtomNeighbors = ();
     552   NEIGHBOR: for $Neighbor ($Molecule->_GetAtomNeighbors($This)) {
     553       $NeighborID = $Neighbor->GetID();
     554       if (exists $ExcludeNeighborsIDsMap{$NeighborID}) {
     555         next NEIGHBOR;
     556       }
     557     push @FilteredAtomNeighbors, $Neighbor;
     558   }
     559
     560   return wantarray ? @FilteredAtomNeighbors : scalar @FilteredAtomNeighbors;
     561 }
     562
     563 # Get specific atom neighbor objects as array. In scalar conetxt, return number of neighbors.
     564 #
     565 # Notes:
     566 #   . AtomSpecification correspond to any valid AtomicInvariant based atomic specifications
     567 #     as implemented in DoesAtomNeighborhoodMatch method.
     568 #   . Multiple atom specifications can be used in a string delimited by comma.
     569 #
     570 sub GetNeighborsUsingAtomSpecification {
     571   my($This, $AtomSpecification, @ExcludeNeighbors) = @_;
     572   my(@AtomNeighbors);
     573
     574   @AtomNeighbors = ();
     575   @AtomNeighbors = $This->GetNeighbors(@ExcludeNeighbors);
     576
     577   # Does atom has any neighbors and do they need to be filtered?
     578   if (!(@AtomNeighbors && defined($AtomSpecification) && $AtomSpecification)) {
     579     return wantarray ? @AtomNeighbors : scalar @AtomNeighbors;
     580   }
     581
     582   # Filter neighbors using atom specification...
     583   my($AtomNeighbor, @FilteredAtomNeighbors);
     584
     585   @FilteredAtomNeighbors = ();
     586   NEIGHBOR: for $AtomNeighbor (@AtomNeighbors) {
     587     if (!$AtomNeighbor->_DoesAtomSpecificationMatch($AtomSpecification)) {
     588       next NEIGHBOR;
     589     }
     590     push @FilteredAtomNeighbors, $AtomNeighbor;
     591   }
     592
     593   return wantarray ? @FilteredAtomNeighbors : scalar @FilteredAtomNeighbors;
     594 }
     595
     596
     597 # Get non-hydrogen atom neighbor objects as array. In scalar context, return number of neighbors...
     598 sub GetHeavyAtomNeighbors {
     599   my($This) = @_;
     600
     601   return $This->GetNonHydrogenAtomNeighbors();
     602 }
     603
     604 # Get non-hydrogen atom neighbor objects as array. In scalar context, return number of neighbors...
     605 sub GetNonHydrogenAtomNeighbors {
     606   my($This) = @_;
     607
     608   # Is this atom in a molecule?
     609   if (!$This->HasProperty('Molecule')) {
     610     return undef;
     611   }
     612   my($NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = (1, 0);
     613
     614   return $This->_GetFilteredAtomNeighbors($NonHydrogenAtomsOnly, $HydrogenAtomsOnly);
     615 }
     616
     617 # Get hydrogen atom neighbor objects as array. In scalar context, return numbe of neighbors...
     618 sub GetHydrogenAtomNeighbors {
     619   my($This) = @_;
     620
     621   # Is this atom in a molecule?
     622   if (!$This->HasProperty('Molecule')) {
     623     return undef;
     624   }
     625   my($NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = (0, 1);
     626
     627   return $This->_GetFilteredAtomNeighbors($NonHydrogenAtomsOnly, $HydrogenAtomsOnly);
     628 }
     629
     630 # Get non-hydrogen neighbor of hydrogen atom...
     631 #
     632 sub GetNonHydrogenNeighborOfHydrogenAtom {
     633   my($This) = @_;
     634
     635   # Is it Hydrogen?
     636   if (!$This->IsHydrogen()) {
     637     return undef;
     638   }
     639   my(@Neighbors);
     640
     641   @Neighbors = $This->GetNonHydrogenAtomNeighbors();
     642
     643   return (@Neighbors == 1) ? $Neighbors[0] : undef;
     644 }
     645
     646 # Get filtered atom atom neighbors
     647 sub _GetFilteredAtomNeighbors {
     648   my($This, $NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = @_;
     649
     650   # Check flags...
     651   if (!defined $NonHydrogenAtomsOnly) {
     652     $NonHydrogenAtomsOnly = 0;
     653   }
     654   if (!defined $HydrogenAtomsOnly) {
     655     $HydrogenAtomsOnly = 0;
     656   }
     657   my($Neighbor, @FilteredAtomNeighbors);
     658
     659   @FilteredAtomNeighbors = ();
     660   NEIGHBOR: for $Neighbor ($This->GetNeighbors()) {
     661     if ($NonHydrogenAtomsOnly && $Neighbor->IsHydrogen()) {
     662       next NEIGHBOR;
     663     }
     664     if ($HydrogenAtomsOnly && (!$Neighbor->IsHydrogen())) {
     665       next NEIGHBOR;
     666     }
     667     push @FilteredAtomNeighbors, $Neighbor;
     668   }
     669
     670   return wantarray ? @FilteredAtomNeighbors : scalar @FilteredAtomNeighbors;
     671 }
     672
     673 # Get number of neighbors...
     674 #
     675 sub GetNumOfNeighbors {
     676   my($This) = @_;
     677   my($NumOfNeighbors);
     678
     679   $NumOfNeighbors = $This->GetNeighbors();
     680
     681   return (defined $NumOfNeighbors) ? $NumOfNeighbors : undef;
     682 }
     683
     684 # Get number of neighbors which are non-hydrogen atoms...
     685 sub GetNumOfHeavyAtomNeighbors {
     686   my($This) = @_;
     687
     688   return $This->GetNumOfNonHydrogenAtomNeighbors();
     689 }
     690
     691 # Get number of neighbors which are non-hydrogen atoms...
     692 sub GetNumOfNonHydrogenAtomNeighbors {
     693   my($This) = @_;
     694   my($NumOfNeighbors);
     695
     696   $NumOfNeighbors = $This->GetNonHydrogenAtomNeighbors();
     697
     698   return (defined $NumOfNeighbors) ? $NumOfNeighbors : undef;
     699 }
     700
     701 # Get number of neighbors which are hydrogen atoms...
     702 sub GetNumOfHydrogenAtomNeighbors {
     703   my($This) = @_;
     704   my($NumOfNeighbors);
     705
     706   $NumOfNeighbors = $This->GetHydrogenAtomNeighbors();
     707
     708   return (defined $NumOfNeighbors) ? $NumOfNeighbors : undef;
     709 }
     710
     711 # Get bond objects as array. In scalar context, return number of bonds...
     712 sub GetBonds {
     713   my($This) = @_;
     714
     715   # Is this atom in a molecule?
     716   if (!$This->HasProperty('Molecule')) {
     717     return undef;
     718   }
     719   my($Molecule) = $This->GetProperty('Molecule');
     720
     721   return $Molecule->_GetAtomBonds($This);
     722 }
     723
     724 # Get bond to specified atom...
     725 sub GetBondToAtom {
     726   my($This, $Other) = @_;
     727
     728   # Is this atom in a molecule?
     729   if (!$This->HasProperty('Molecule')) {
     730     return undef;
     731   }
     732   my($Molecule) = $This->GetProperty('Molecule');
     733
     734   return $Molecule->_GetBondToAtom($This, $Other);
     735 }
     736
     737 # It it bonded to a specified atom?
     738 sub IsBondedToAtom {
     739   my($This, $Other) = @_;
     740
     741   # Is this atom in a molecule?
     742   if (!$This->HasProperty('Molecule')) {
     743     return undef;
     744   }
     745   my($Molecule) = $This->GetProperty('Molecule');
     746
     747   return $Molecule->_IsBondedToAtom($This, $Other);
     748 }
     749
     750 # Get bond objects to non-hydrogen atoms as array. In scalar context, return number of bonds...
     751 sub GetBondsToHeavyAtoms {
     752   my($This) = @_;
     753
     754   return $This->GetBondsToNonHydrogenAtoms();
     755 }
     756
     757 # Get bond objects to non-hydrogen atoms as array. In scalar context, return number of bonds...
     758 sub GetBondsToNonHydrogenAtoms {
     759   my($This) = @_;
     760
     761   # Is this atom in a molecule?
     762   if (!$This->HasProperty('Molecule')) {
     763     return undef;
     764   }
     765   my($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = (1, 0);
     766
     767   return $This->_GetFilteredBonds($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly);
     768 }
     769
     770 # Get bond objects to hydrogen atoms as array. In scalar context, return number of bonds...
     771 sub GetBondsToHydrogenAtoms {
     772   my($This) = @_;
     773
     774   # Is this atom in a molecule?
     775   if (!$This->HasProperty('Molecule')) {
     776     return undef;
     777   }
     778   my($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = (0, 1);
     779
     780   return $This->_GetFilteredBonds($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly);
     781 }
     782
     783 # Get filtered bonds...
     784 sub _GetFilteredBonds {
     785   my($This, $BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = @_;
     786
     787   # Check flags...
     788   if (!defined $BondsToNonHydrogenAtomsOnly) {
     789     $BondsToNonHydrogenAtomsOnly = 0;
     790   }
     791   if (!defined $BondsToHydrogenAtomsOnly) {
     792     $BondsToHydrogenAtomsOnly = 0;
     793   }
     794
     795   my($Bond, $BondedAtom, @FilteredBonds);
     796
     797   @FilteredBonds = ();
     798   BOND: for $Bond ($This->GetBonds()) {
     799     $BondedAtom = $Bond->GetBondedAtom($This);
     800     if ($BondsToNonHydrogenAtomsOnly && $BondedAtom->IsHydrogen()) {
     801       next BOND;
     802     }
     803     if ($BondsToHydrogenAtomsOnly && (!$BondedAtom->IsHydrogen())) {
     804       next BOND;
     805     }
     806     push @FilteredBonds, $Bond;
     807   }
     808
     809   return wantarray ? @FilteredBonds : (scalar @FilteredBonds);
     810 }
     811
     812 # Get number of bonds...
     813 #
     814 sub GetNumOfBonds {
     815   my($This) = @_;
     816   my($NumOfBonds);
     817
     818   $NumOfBonds = $This->GetBonds();
     819
     820   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
     821 }
     822
     823 # Get number of bonds to non-hydrogen atoms...
     824 sub GetNumOfBondsToHeavyAtoms {
     825   my($This) = @_;
     826
     827   return $This->GetNumOfBondsToNonHydrogenAtoms();
     828 }
     829
     830 # Get number of bonds to non-hydrogen atoms...
     831 sub GetNumOfBondsToNonHydrogenAtoms {
     832   my($This) = @_;
     833   my($NumOfBonds);
     834
     835   $NumOfBonds = $This->GetBondsToNonHydrogenAtoms();
     836
     837   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
     838 }
     839
     840 # Get number of single bonds to heavy atoms...
     841 sub GetNumOfSingleBondsToHeavyAtoms {
     842   my($This) = @_;
     843
     844   return $This->GetNumOfSingleBondsToNonHydrogenAtoms();
     845 }
     846
     847 # Get number of single bonds to non-hydrogen atoms...
     848 sub GetNumOfSingleBondsToNonHydrogenAtoms {
     849   my($This) = @_;
     850
     851   # Is this atom in a molecule?
     852   if (!$This->HasProperty('Molecule')) {
     853     return undef;
     854   }
     855   return $This->_GetNumOfBondsWithSpecifiedBondOrderToNonHydrogenAtoms(1);
     856 }
     857
     858 # Get number of double bonds to heavy atoms...
     859 sub GetNumOfDoubleBondsToHeavyAtoms {
     860   my($This) = @_;
     861
     862   return $This->GetNumOfDoubleBondsToNonHydrogenAtoms();
     863 }
     864
     865 # Get number of double bonds to non-hydrogen atoms...
     866 sub GetNumOfDoubleBondsToNonHydrogenAtoms {
     867   my($This) = @_;
     868
     869   # Is this atom in a molecule?
     870   if (!$This->HasProperty('Molecule')) {
     871     return undef;
     872   }
     873   return $This->_GetNumOfBondsWithSpecifiedBondOrderToNonHydrogenAtoms(2);
     874 }
     875
     876 # Get number of triple bonds to heavy atoms...
     877 sub GetNumOfTripleBondsToHeavyAtoms {
     878   my($This) = @_;
     879
     880   return $This->GetNumOfTripleBondsToNonHydrogenAtoms();
     881 }
     882
     883 # Get number of triple bonds to non-hydrogen atoms...
     884 sub GetNumOfTripleBondsToNonHydrogenAtoms {
     885   my($This) = @_;
     886
     887   # Is this atom in a molecule?
     888   if (!$This->HasProperty('Molecule')) {
     889     return undef;
     890   }
     891   return $This->_GetNumOfBondsWithSpecifiedBondOrderToNonHydrogenAtoms(3);
     892 }
     893
     894 # Get number of bonds of specified bond order to non-hydrogen atoms...
     895 sub _GetNumOfBondsWithSpecifiedBondOrderToNonHydrogenAtoms {
     896   my($This, $SpecifiedBondOrder) = @_;
     897   my($NumOfBonds, $Bond, $BondOrder, @Bonds);
     898
     899   $NumOfBonds = 0;
     900   @Bonds = $This->GetBondsToNonHydrogenAtoms();
     901   for $Bond (@Bonds) {
     902     $BondOrder = $Bond->GetBondOrder();
     903     if ($SpecifiedBondOrder == $BondOrder) {
     904       $NumOfBonds++;
     905     }
     906   }
     907   return $NumOfBonds;
     908 }
     909
     910 # Get number of aromatic bonds to heavy atoms...
     911 sub GetNumOfAromaticBondsToHeavyAtoms {
     912   my($This) = @_;
     913
     914   return $This->GetNumOfAromaticBondsToNonHydrogenAtoms();
     915 }
     916
     917 # Get number of aromatic bonds to non-hydrogen atoms...
     918 sub GetNumOfAromaticBondsToNonHydrogenAtoms {
     919   my($This) = @_;
     920   my($NumOfBonds, $Bond, @Bonds);
     921
     922   # Is this atom in a molecule?
     923   if (!$This->HasProperty('Molecule')) {
     924     return undef;
     925   }
     926
     927   $NumOfBonds = 0;
     928   @Bonds = $This->GetBondsToNonHydrogenAtoms();
     929   for $Bond (@Bonds) {
     930     if ($Bond->IsAromatic()) { $NumOfBonds++; }
     931   }
     932   return $NumOfBonds;
     933 }
     934
     935 # Get number of different bond types to non-hydrogen atoms...
     936 #
     937 sub GetNumOfBondTypesToHeavyAtoms {
     938   my($This, $CountAromaticBonds) = @_;
     939
     940   return $This->GetNumOfBondTypesToNonHydrogenAtoms($CountAromaticBonds);
     941 }
     942
     943 # Get number of single, double, triple, and aromatic bonds from an atom to all other
     944 # non-hydrogen atoms. Value of CountAtomaticBonds parameter controls whether
     945 # number of aromatic bonds is returned; default is not to count aromatic bonds. During
     946 # counting of aromatic bonds, the bond marked aromatic is not included in the count
     947 # of other bond types.
     948 #
     949 sub GetNumOfBondTypesToNonHydrogenAtoms {
     950   my($This, $CountAromaticBonds) = @_;
     951   my($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfAromaticBonds, $None, $Bond, @Bonds);
     952
     953   $CountAromaticBonds = defined($CountAromaticBonds) ? $CountAromaticBonds : 0;
     954
     955   ($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds) = ('0') x 3;
     956   $NumOfAromaticBonds = $CountAromaticBonds ? 0 : undef;
     957
     958   # Is this atom in a molecule?
     959   if (!$This->HasProperty('Molecule')) {
     960     return ($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfAromaticBonds);
     961   }
     962
     963   @Bonds = $This->GetBondsToNonHydrogenAtoms();
     964
     965   for $Bond (@Bonds) {
     966     BONDTYPE: {
     967       if ($CountAromaticBonds) {
     968         if ($Bond->IsAromatic()) { $NumOfAromaticBonds++; last BONDTYPE; }
     969       }
     970       if ($Bond->IsSingle()) { $NumOfSingleBonds++; last BONDTYPE; }
     971       if ($Bond->IsDouble()) { $NumOfDoubleBonds++; last BONDTYPE; }
     972       if ($Bond->IsTriple()) { $NumOfTripleBonds++; last BONDTYPE; }
     973       $None = 1;
     974     }
     975   }
     976   return ($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfAromaticBonds);
     977 }
     978
     979 # Get number of sigma and pi bonds to heavy atoms...
     980 #
     981 sub GetNumOfSigmaAndPiBondsToHeavyAtoms {
     982   my($This) = @_;
     983
     984   return $This->GetNumOfSigmaAndPiBondsToNonHydrogenAtoms();
     985 }
     986
     987 # Get number of sigma and pi bonds from an atom to all other non-hydrogen atoms.
     988 # Sigma and pi bonds are counted using the following methodology: a single bond
     989 # correspond to one sigma bond; a double bond contributes one to sigma bond count
     990 # and one to pi bond count; a triple bond contributes one to sigma bond count and
     991 # two to pi bond count.
     992 #
     993 sub GetNumOfSigmaAndPiBondsToNonHydrogenAtoms {
     994   my($This) = @_;
     995   my($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfSigmaBonds, $NumOfPiBonds);
     996
     997   ($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds) = $This->GetNumOfBondTypesToNonHydrogenAtoms();
     998
     999   $NumOfSigmaBonds = $NumOfSingleBonds + $NumOfDoubleBonds + $NumOfTripleBonds;
    1000   $NumOfPiBonds = $NumOfDoubleBonds + 2*$NumOfTripleBonds;
    1001
    1002   return ($NumOfSigmaBonds, $NumOfPiBonds);
    1003 }
    1004
    1005 # Get information related to atoms for all heavy atoms attached to an atom..
    1006 #
    1007 sub GetHeavyAtomNeighborsAtomInformation {
    1008   my($This) = @_;
    1009
    1010   return $This->GetNonHydrogenAtomNeighborsAtomInformation();
    1011 }
    1012
    1013 # Get information related to atoms for all non-hydrogen atoms attached to an atom..
    1014 #
    1015 # The following values are returned:
    1016 #   . Number of non-hydrogen atom neighbors
    1017 #   . A reference to an array containing atom objects correpsonding to non-hydrogen
    1018 #     atom neighbors
    1019 #   . Number of different types of non-hydrogen atom neighbors
    1020 #   . A reference to a hash containing atom symbol as key with value corresponding
    1021 #     to its count for non-hydrogen atom neighbors
    1022 #
    1023 sub GetNonHydrogenAtomNeighborsAtomInformation {
    1024   my($This) = @_;
    1025
    1026   # Is this atom in a molecule?
    1027   if (!$This->HasProperty('Molecule')) {
    1028     return (undef, undef, undef, undef);
    1029   }
    1030   my($AtomSymbol, $AtomNeighbor, $NumOfAtomNeighbors, $NumOfAtomNeighborsType, @AtomNeighbors, %AtomNeighborsTypeMap);
    1031
    1032   $NumOfAtomNeighbors = 0; @AtomNeighbors = ();
    1033   $NumOfAtomNeighborsType = 0; %AtomNeighborsTypeMap = ();
    1034
    1035   @AtomNeighbors = $This->GetNonHydrogenAtomNeighbors();
    1036   $NumOfAtomNeighbors = scalar @AtomNeighbors;
    1037
    1038   for $AtomNeighbor (@AtomNeighbors) {
    1039     $AtomSymbol = $AtomNeighbor->{AtomSymbol};
    1040     if (exists $AtomNeighborsTypeMap{$AtomSymbol}) {
    1041       $AtomNeighborsTypeMap{$AtomSymbol} += 1;
    1042     }
    1043     else {
    1044       $AtomNeighborsTypeMap{$AtomSymbol} = 1;
    1045       $NumOfAtomNeighborsType++;
    1046     }
    1047   }
    1048
    1049   return ($NumOfAtomNeighbors, \@AtomNeighbors, $NumOfAtomNeighborsType, \%AtomNeighborsTypeMap);
    1050 }
    1051
    1052 # Get information related to bonds for all heavy atoms attached to an atom..
    1053 #
    1054 sub GetHeavyAtomNeighborsBondformation {
    1055   my($This) = @_;
    1056
    1057   return $This->GetNonHydrogenAtomNeighborsBondInformation();
    1058 }
    1059
    1060 # Get information related to bonds for all non-hydrogen atoms attached to an atom..
    1061 #
    1062 # The following values are returned:
    1063 #   . Number of bonds to non-hydrogen atom neighbors
    1064 #   . A reference to an array containing bond objects correpsonding to non-hydrogen
    1065 #     atom neighbors
    1066 #   . A reference to a hash containing bond type as key with value corresponding
    1067 #     to its count for non-hydrogen atom neighbors. Bond types are: Single, Double or Triple
    1068 #   . A reference to a hash containing atom symbol as key pointing to bond type as second
    1069 #     key with values correponding to count of bond types for atom symbol for non-hydrogen
    1070 #     atom neighbors
    1071 #   . A reference to a hash containing atom symbol as key pointing to bond type as second
    1072 #     key with values correponding to atom objects array involved in corresponding bond type for
    1073 #     atom symbol for non-hydrogen atom neighbors
    1074 #
    1075 sub GetNonHydrogenAtomNeighborsBondInformation {
    1076   my($This) = @_;
    1077
    1078   # Is this atom in a molecule?
    1079   if (!$This->HasProperty('Molecule')) {
    1080     return (undef, undef, undef, undef, undef);
    1081   }
    1082   my($BondedAtom, $BondedAtomSymbol, $BondType, $None, $Bond, $NumOfBonds, @Bonds, %BondTypeCountMap, %AtomsBondTypesCountMap, %AtomsBondTypeAtomsMap);
    1083
    1084   $NumOfBonds = 0; @Bonds = ();
    1085   %BondTypeCountMap = ();
    1086   %AtomsBondTypesCountMap = (); %AtomsBondTypeAtomsMap = ();
    1087
    1088   $BondTypeCountMap{Single} = 0;
    1089   $BondTypeCountMap{Double} = 0;
    1090   $BondTypeCountMap{Triple} = 0;
    1091
    1092   @Bonds = $This->GetBondsToNonHydrogenAtoms();
    1093   $NumOfBonds = scalar @Bonds;
    1094
    1095   BOND: for $Bond (@Bonds) {
    1096     $BondType = $Bond->IsSingle() ? "Single" : ($Bond->IsDouble() ? "Double" : ($Bond->IsTriple() ? "Triple" : ""));
    1097     if (!$BondType) {
    1098       next BOND;
    1099     }
    1100
    1101     # Track bond types...
    1102     if (exists $BondTypeCountMap{$BondType}) {
    1103       $BondTypeCountMap{$BondType} += 1;
    1104     }
    1105     else {
    1106       $BondTypeCountMap{$BondType} = 1;
    1107     }
    1108
    1109     $BondedAtom = $Bond->GetBondedAtom($This);
    1110     $BondedAtomSymbol = $BondedAtom->{AtomSymbol};
    1111
    1112     # Track bond types count for atom types involved in specific bond types...
    1113     if (!exists $AtomsBondTypesCountMap{$BondedAtomSymbol}) {
    1114       %{$AtomsBondTypesCountMap{$BondedAtomSymbol}} = ();
    1115     }
    1116     if (exists $AtomsBondTypesCountMap{$BondedAtomSymbol}{$BondType}) {
    1117       $AtomsBondTypesCountMap{$BondedAtomSymbol}{$BondType} += 1;
    1118     }
    1119     else {
    1120       $AtomsBondTypesCountMap{$BondedAtomSymbol}{$BondType} = 1;
    1121     }
    1122
    1123     # Track atoms involved in specific bond types for specific atom types...
    1124     if (!exists $AtomsBondTypeAtomsMap{$BondedAtomSymbol}) {
    1125       %{$AtomsBondTypeAtomsMap{$BondedAtomSymbol}} = ();
    1126     }
    1127     if (!exists $AtomsBondTypeAtomsMap{$BondedAtomSymbol}{$BondType}) {
    1128       @{$AtomsBondTypeAtomsMap{$BondedAtomSymbol}{$BondType}} = ();
    1129     }
    1130     push @{$AtomsBondTypeAtomsMap{$BondedAtomSymbol}{$BondType}}, $BondedAtom;
    1131   }
    1132
    1133   return ($NumOfBonds, \@Bonds, \%BondTypeCountMap, \%AtomsBondTypesCountMap, \%AtomsBondTypeAtomsMap);
    1134 }
    1135
    1136 # Get number of bonds to hydrogen atoms...
    1137 sub GetNumOfBondsToHydrogenAtoms {
    1138   my($This) = @_;
    1139   my($NumOfBonds);
    1140
    1141   $NumOfBonds = $This->GetBondsToHydrogenAtoms();
    1142
    1143   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
    1144 }
    1145
    1146 # Get sum of bond orders to all bonded atoms...
    1147 #
    1148 sub GetSumOfBondOrders {
    1149   my($This) = @_;
    1150
    1151   # Is this atom in a molecule?
    1152   if (!$This->HasProperty('Molecule')) {
    1153     return undef;
    1154   }
    1155
    1156   return $This->_GetSumOfBondOrders();
    1157 }
    1158
    1159 # Get sum of bond orders to non-hydrogen atoms only...
    1160 #
    1161 sub GetSumOfBondOrdersToHeavyAtoms {
    1162   my($This) = @_;
    1163
    1164   return $This->GetSumOfBondOrdersToNonHydrogenAtoms();
    1165 }
    1166
    1167 # Get sum of bond orders to non-hydrogen atoms only...
    1168 #
    1169 sub GetSumOfBondOrdersToNonHydrogenAtoms {
    1170   my($This) = @_;
    1171
    1172   # Is this atom in a molecule?
    1173   if (!$This->HasProperty('Molecule')) {
    1174     return undef;
    1175   }
    1176   my($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = (1, 0);
    1177
    1178   return $This->_GetSumOfBondOrders($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly);
    1179 }
    1180
    1181 # Get sum of bond orders to hydrogen atoms only...
    1182 #
    1183 sub GetSumOfBondOrdersToHydrogenAtoms {
    1184   my($This) = @_;
    1185
    1186   # Is this atom in a molecule?
    1187   if (!$This->HasProperty('Molecule')) {
    1188     return undef;
    1189   }
    1190   my($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = (0, 1);
    1191
    1192   return $This->_GetSumOfBondOrders($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly);
    1193 }
    1194
    1195 # Get sum of bond orders to all bonded atoms,  non-hydrogen or hydrogen bonded atoms...
    1196 #
    1197 sub _GetSumOfBondOrders {
    1198   my($This, $ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = @_;
    1199
    1200   # Check flags...
    1201   if (!defined $ToNonHydrogenAtomsOnly) {
    1202     $ToNonHydrogenAtomsOnly = 0;
    1203   }
    1204   if (!defined $ToHydrogenAtomsOnly) {
    1205     $ToHydrogenAtomsOnly = 0;
    1206   }
    1207   my($Bond, $SumOfBondOrders, @Bonds);
    1208   @Bonds = ();
    1209
    1210   if ($ToNonHydrogenAtomsOnly) {
    1211     @Bonds = $This->GetBondsToNonHydrogenAtoms();
    1212   }
    1213   elsif ($ToHydrogenAtomsOnly) {
    1214     @Bonds = $This->GetBondsToHydrogenAtoms();
    1215   }
    1216   else {
    1217     # All bonds...
    1218     @Bonds = $This->GetBonds();
    1219   }
    1220
    1221   $SumOfBondOrders = 0;
    1222   for $Bond (@Bonds) {
    1223     $SumOfBondOrders += $Bond->GetBondOrder();
    1224   }
    1225
    1226   if ($SumOfBondOrders =~ /\./) {
    1227     #
    1228     # Change any fractional bond order to next largest integer...
    1229     #
    1230     # As long as aromatic bond orders in a ring are correctly using using 4n + 2 Huckel rule
    1231     # (BondOrder: 1.5) or explicity set as Kekule  bonds (alternate single/double),
    1232     # SumOfBondOrders should add up to an integer.
    1233     #
    1234     $SumOfBondOrders = ceil($SumOfBondOrders);
    1235   }
    1236
    1237   return $SumOfBondOrders;
    1238 }
    1239
    1240 # Get largest bond order to any bonded atoms...
    1241 #
    1242 sub GetLargestBondOrder {
    1243   my($This) = @_;
    1244
    1245   # Is this atom in a molecule?
    1246   if (!$This->HasProperty('Molecule')) {
    1247     return undef;
    1248   }
    1249
    1250   return $This->_GetLargestBondOrder();
    1251 }
    1252
    1253 # Get largest bond order to bonded non-hydrogen atoms...
    1254 #
    1255 sub GetLargestBondOrderToHeavyAtoms {
    1256   my($This) = @_;
    1257
    1258   return $This->GetLargestBondOrderToNonHydrogenAtoms();
    1259 }
    1260
    1261 # Get largest bond order to bonded non-hydrogen atoms...
    1262 #
    1263 sub GetLargestBondOrderToNonHydrogenAtoms {
    1264   my($This) = @_;
    1265
    1266   # Is this atom in a molecule?
    1267   if (!$This->HasProperty('Molecule')) {
    1268     return undef;
    1269   }
    1270
    1271   my($ToNonHydrogenAtomsOnly) = (1);
    1272
    1273   return $This->_GetLargestBondOrder($ToNonHydrogenAtomsOnly);
    1274 }
    1275
    1276 # Get largest bond order to all bonded atoms, non-hydrogen or hydrogen bonded atoms...
    1277 #
    1278 sub _GetLargestBondOrder {
    1279   my($This, $ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = @_;
    1280
    1281   # Check flags...
    1282   if (!defined $ToNonHydrogenAtomsOnly) {
    1283     $ToNonHydrogenAtomsOnly = 0;
    1284   }
    1285   if (!defined $ToHydrogenAtomsOnly) {
    1286     $ToHydrogenAtomsOnly = 0;
    1287   }
    1288   my($Bond, $LargestBondOrder, $BondOrder, @Bonds);
    1289   @Bonds = ();
    1290
    1291   if ($ToNonHydrogenAtomsOnly) {
    1292     @Bonds = $This->GetBondsToNonHydrogenAtoms();
    1293   }
    1294   elsif ($ToHydrogenAtomsOnly) {
    1295     @Bonds = $This->GetBondsToHydrogenAtoms();
    1296   }
    1297   else {
    1298     # All bonds...
    1299     @Bonds = $This->GetBonds();
    1300   }
    1301
    1302   $LargestBondOrder = 0;
    1303   for $Bond (@Bonds) {
    1304     $BondOrder = $Bond->GetBondOrder();
    1305     if ($BondOrder > $LargestBondOrder) {
    1306       $LargestBondOrder = $BondOrder;
    1307     }
    1308   }
    1309
    1310   return $LargestBondOrder;
    1311 }
    1312
    1313 # Get number of implicit hydrogen for atom...
    1314 #
    1315 sub GetImplicitHydrogens {
    1316   my($This) = @_;
    1317
    1318   return $This->GetNumOfImplicitHydrogens();
    1319 }
    1320
    1321 # Get number of implicit hydrogen for atom...
    1322 #
    1323 sub GetNumOfImplicitHydrogens {
    1324   my($This) = @_;
    1325
    1326   # Is this atom in a molecule?
    1327   if (!$This->HasProperty('Molecule')) {
    1328     return undef;
    1329   }
    1330
    1331   # Is ImplicitHydrogens property explicitly set?
    1332   if ($This->HasProperty('ImplicitHydrogens')) {
    1333     return $This->GetProperty('ImplicitHydrogens');
    1334   }
    1335
    1336   # Is it an element symbol?
    1337   if (!$This->{AtomicNumber}) {
    1338     return 0;
    1339   }
    1340
    1341   my($ImplicitHydrogens, $PotentialTotalValence, $SumOfBondOrders);
    1342
    1343   $ImplicitHydrogens = 0;
    1344   $SumOfBondOrders = $This->GetSumOfBondOrders();
    1345   $PotentialTotalValence = $This->GetPotentialTotalCommonValence();
    1346
    1347   if (defined($PotentialTotalValence) && defined($SumOfBondOrders)) {
    1348     # Subtract sum of bond orders to non-hydrogen and hydrogen atom neighbors...
    1349     $ImplicitHydrogens = $PotentialTotalValence - $SumOfBondOrders;
    1350   }
    1351
    1352   return $ImplicitHydrogens > 0 ? $ImplicitHydrogens : 0;
    1353 }
    1354
    1355 # Get number of bonds available to form additional bonds with heavy atoms, excluding
    1356 # any implicit bonds to hydrogens set using ImplicitHydrogens property.
    1357 #
    1358 # It's different from number of implicit or missing hydrogens, both of which are equivalent.
    1359 #
    1360 # For example, in a SMILES string, [nH] ring atom corresponds to an aromatic nitrogen.
    1361 # Although the hydrogen specified for n is treated internally as implicit hydrogen and shows
    1362 # up in missing hydrogen count, it's not available to participate in double bonds to additional
    1363 # heavy atoms.
    1364 #
    1365 sub GetNumOfBondsAvailableForHeavyAtoms {
    1366   my($This) = @_;
    1367
    1368   return $This->GetNumOfBondsAvailableForNonHydrogenAtoms();
    1369 }
    1370
    1371 # It's another name for GetNumOfBondsAvailableForHeavyAtoms
    1372 #
    1373 sub GetNumOfBondsAvailableForNonHydrogenAtoms {
    1374   my($This) = @_;
    1375   my($NumOfAvailableBonds, $PotentialTotalValence, $SumOfBondOrders);
    1376
    1377   $NumOfAvailableBonds = 0;
    1378
    1379   $SumOfBondOrders = $This->GetSumOfBondOrders();
    1380   $PotentialTotalValence = $This->GetPotentialTotalCommonValence();
    1381
    1382   if (defined($PotentialTotalValence) && defined($SumOfBondOrders)) {
    1383     # Subtract sum of bond orders to non-hydrogen and hydrogen atom neighbors...
    1384     $NumOfAvailableBonds = $PotentialTotalValence - $SumOfBondOrders;
    1385   }
    1386
    1387   if ($This->HasProperty('ImplicitHydrogens')) {
    1388     $NumOfAvailableBonds -= $This->GetProperty('ImplicitHydrogens');
    1389   }
    1390
    1391   return $NumOfAvailableBonds > 0 ? $NumOfAvailableBonds : 0;
    1392 }
    1393
    1394 # Disable setting of explicit hydrogens property...
    1395 sub SetExplicitHydrogens {
    1396   my($This, $Value) = @_;
    1397
    1398   carp "Warning: ${ClassName}->SetExplicitHydrogens: Setting of explicit hydrogens is not supported...";
    1399
    1400   return $This;
    1401 }
    1402
    1403 # Get number of explicit hydrogens for atom...
    1404 #
    1405 sub GetExplicitHydrogens {
    1406   my($This) = @_;
    1407
    1408   return $This->GetNumOfExplicitHydrogens();
    1409 }
    1410
    1411 # Get number of explicit hydrogens for atom...
    1412 #
    1413 sub GetNumOfExplicitHydrogens {
    1414   my($This) = @_;
    1415   my($HydrogenAtomNbrs);
    1416
    1417   # Is this atom in a molecule?
    1418   if (!$This->HasProperty('Molecule')) {
    1419     return undef;
    1420   }
    1421
    1422   $HydrogenAtomNbrs = $This->GetNumOfHydrogenAtomNeighbors();
    1423
    1424   return defined $HydrogenAtomNbrs ? $HydrogenAtomNbrs : 0;
    1425 }
    1426
    1427 # Get num of missing hydrogens...
    1428 #
    1429 sub GetMissingHydrogens {
    1430   my($This) = @_;
    1431
    1432   return $This->GetNumOfMissingHydrogens();
    1433 }
    1434
    1435 # Get num of missing hydrogens...
    1436 #
    1437 sub GetNumOfMissingHydrogens {
    1438   my($This) = @_;
    1439
    1440   # Is this atom in a molecule?
    1441   if (!$This->HasProperty('Molecule')) {
    1442     return undef;
    1443   }
    1444
    1445   return $This->GetNumOfImplicitHydrogens();
    1446 }
    1447
    1448 # Get total number of hydrogens...
    1449 #
    1450 sub GetHydrogens {
    1451   my($This) = @_;
    1452
    1453   return $This->GetNumOfHydrogens();
    1454 }
    1455
    1456 # Get total number of hydrogens...
    1457 #
    1458 sub GetNumOfHydrogens {
    1459   my($This) = @_;
    1460
    1461   # Is this atom in a molecule?
    1462   if (!$This->HasProperty('Molecule')) {
    1463     return undef;
    1464   }
    1465
    1466   return $This->GetNumOfImplicitHydrogens() + $This->GetNumOfExplicitHydrogens();
    1467 }
    1468
    1469 # Valence corresponds to the number of electrons used by an atom in bonding:
    1470 #
    1471 #   Valence = ValenceElectrons - ValenceFreeElectrons = BondingElectrons
    1472 #
    1473 # Single, double, triple bonds with bond orders of 1, 2 and 3 correspond to contribution of
    1474 # 1, 2, and 3 bonding electrons. So Valence can be computed using:
    1475 #
    1476 #   Valence = SumOfBondOrders + NumOfMissingHydrogens + FormalCharge
    1477 #
    1478 # where positive and negative values of FormalCharge increase and decrease the number
    1479 # of bonding electrons respectively.
    1480 #
    1481 # The current release of MayaChemTools supports the following three valence models, which
    1482 # are used during calculation of implicit hydrogens: MDLValenceModel, DaylightValenceModel,
    1483 # InternalValenceModel or MayaChemToolsValenceModel.
    1484 #
    1485 # Notes:
    1486 #  . This doesn't always corresponds to explicit valence.
    1487 #  . Missing hydrogens are included in the valence.
    1488 #  . For neutral molecules, valence and sum of bond order are equal.
    1489 #  . For molecules containing only single bonds, SumOfBondOrders and NumOfBonds are equal.
    1490 #  . Free radical electrons lead to the decrease in valence. For atoms with explicit assignment
    1491 #    of SpinMultiplicity property values corresponding to Singlet (two unparied electrons
    1492 #    corresponding to one spin state), Doublet (free radical; an unpaired electron corresponding
    1493 #    to two spin states), and Triplet (two unparied electrons corresponding to three spin states;
    1494 #    divalent carbon atoms (carbenes)), FreeRadicalElectrons are calculated as follows:
    1495 #
    1496 #       SpinMultiplicity: Doublet(2); FreeRadicalElectrons: 1 (one valence electron not available for bonding)
    1497 #       SpinMultiplicity: Singlet(1)/Triplet(3); FreeRadicalElectrons: 2 (two valence electrons not available for bonding)
    1498 #
    1499 sub GetValence {
    1500   my($This) = @_;
    1501
    1502   # Is this atom in a molecule?
    1503   if (!$This->HasProperty('Molecule')) {
    1504     return undef;
    1505   }
    1506
    1507   # Is Valence property explicitly set?
    1508   if ($This->HasProperty('Valence')) {
    1509     return $This->GetProperty('Valence');
    1510   }
    1511   my($Valence);
    1512
    1513   $Valence = $This->GetSumOfBondOrders() + $This->GetNumOfMissingHydrogens() + $This->GetFormalCharge();
    1514
    1515   return $Valence > 0 ? $Valence : 0;
    1516 }
    1517
    1518 # Get free non-bodning valence electrons left on atom after taking into account
    1519 # sum of bond orders, missing hydrogens and formal charged on the atom. Free
    1520 # radical electrons are included in the valence free electrons count by default.
    1521 #
    1522 # Valence corresponds to number of electrons used by atom in bonding:
    1523 #
    1524 #   Valence = ValenceElectrons - ValenceFreeElectrons
    1525 #
    1526 # Additionally, valence can also be calculated by:
    1527 #
    1528 #   Valence = SumOfBondOrders + NumOfMissingHydrogens + FormalCharge
    1529 #
    1530 # Valence and SumOfBondOrders are equal for neutral molecules.
    1531 #
    1532 # From two formulas for Valence described above, non-bonding free electrons
    1533 # left can be computed by:
    1534 #
    1535 #  ValenceFreeElectrons = ValenceElectrons - Valence
    1536 #                       = ValenceElectrons - SumOfBondOrders -
    1537 #                         NumOfMissingHydrogens - FormalCharge
    1538 #
    1539 # . Notes:
    1540 #    . Missing hydrogens are excluded from the valence free electrons.
    1541 #    . Any free radical electrons are considered part of the valence free electrons
    1542 #      by default.
    1543 #
    1544 # Examples:
    1545 #
    1546 # o NH3: ValenceFreeElectrons = 5 - 3 = 5 - 3 - 0 - 0 = 2
    1547 # o NH2: ValenceFreeElectrons = 5 - 3 = 5 - 2 - 1 - 0 = 2
    1548 # o NH4+; ValenceFreeElectrons = 5 - 5 = 5 - 4 - 0 - 1 = 0
    1549 # o NH3+; ValenceFreeElectrons = 5 - 5 = 5 - 3 - 1 - 1 = 0
    1550 # o C(=O)O- : ValenceFreeElectrons on O- = 6 - 0 = 6 - 1 - 0 - (-1) = 6
    1551 # o C(=O)O- : ValenceFreeElectrons on =O = 6 - 2 = 6 - 2 - 0 - 0 = 4
    1552 #
    1553 #
    1554 sub GetValenceFreeElectrons {
    1555   my($This, $ExcludeFreeRadicalElectrons) = @_;
    1556
    1557   # Is this atom in a molecule?
    1558   if (!$This->HasProperty('Molecule')) {
    1559     return undef;
    1560   }
    1561
    1562   # Is ValenceFreeElectrons property explicitly set?
    1563   if ($This->HasProperty('ValenceFreeElectrons')) {
    1564     return $This->GetProperty('ValenceFreeElectrons');
    1565   }
    1566
    1567   if (!$This->{AtomicNumber}) {
    1568     return 0;
    1569   }
    1570
    1571   my($ValenceFreeElectrons);
    1572
    1573   $ValenceFreeElectrons = $This->GetValenceElectrons() - $This->GetValence();
    1574   if ($ExcludeFreeRadicalElectrons) {
    1575     $ValenceFreeElectrons -= $This->GetFreeRadicalElectrons();
    1576   }
    1577
    1578   return $ValenceFreeElectrons > 0 ? $ValenceFreeElectrons : 0;
    1579 }
    1580
    1581 # Get potential total common valence for calculating the number of implicit hydrogens
    1582 # using the specified common valence model or default internal model for a molecule...
    1583 #
    1584 sub GetPotentialTotalCommonValence {
    1585   my($This) = @_;
    1586
    1587   # Is this atom in a molecule?
    1588   if (!$This->HasProperty('Molecule')) {
    1589     return undef;
    1590   }
    1591   my($PotentialTotalValence, $ValenceModel);
    1592
    1593   $PotentialTotalValence = 0;
    1594   $ValenceModel = $This->GetProperty('Molecule')->GetValenceModel();
    1595
    1596   VALENCEMODEL: {
    1597     if ($ValenceModel =~ /^MDLValenceModel$/i) {
    1598       $PotentialTotalValence = $This->_GetPotentialTotalCommonValenceUsingMDLValenceModel();
    1599       last VALENCEMODEL;
    1600     }
    1601     if ($ValenceModel =~ /^DaylightValenceModel$/i) {
    1602        $PotentialTotalValence = $This->_GetPotentialTotalCommonValenceUsingDaylightValenceModel();
    1603       last VALENCEMODEL;
    1604     }
    1605     if ($ValenceModel !~ /^(InternalValenceModel|MayaChemToolsValenceModel)$/i) {
    1606       carp "Warning: ${ClassName}->GetPotentialTotalCommonValence: The current release of MayaChemTools doesn't support the specified valence model $ValenceModel. Supported valence models: MDLValenceModel, DaylightValenceModel, InternalValenceModel. Using internal valence model...";
    1607     }
    1608     # Use internal valence model as the default valence model...
    1609     $PotentialTotalValence = $This->_GetPotentialTotalCommonValenceUsingInternalValenceModel();
    1610   }
    1611
    1612   return $PotentialTotalValence;
    1613 }
    1614
    1615 # Get potential total common valence using data for MDL valence model available in file,
    1616 # lib/data/MDLValenceModelData.csv, distributed with the package...
    1617 #
    1618 sub _GetPotentialTotalCommonValenceUsingMDLValenceModel {
    1619   my($This) = @_;
    1620
    1621   return $This->_GetPotentialTotalCommonValenceUsingValenceModelData(\%MDLValenceModelDataMap);
    1622
    1623 }
    1624
    1625 # Get potential total common valence using data for Daylight valence model available in file,
    1626 # lib/data/DaylightValenceModelData.csv, distributed with the release...
    1627 #
    1628 sub _GetPotentialTotalCommonValenceUsingDaylightValenceModel {
    1629   my($This) = @_;
    1630
    1631   return $This->_GetPotentialTotalCommonValenceUsingValenceModelData(\%DaylightValenceModelDataMap);
    1632 }
    1633
    1634 # Get potential total common valence using data for a specific valence model...
    1635 #
    1636 sub _GetPotentialTotalCommonValenceUsingValenceModelData {
    1637   my($This, $ValenceModelDataRef) = @_;
    1638   my($AtomicNumber, $FormalCharge);
    1639
    1640   $AtomicNumber = $This->{AtomicNumber};
    1641   if (!$AtomicNumber) {
    1642     return 0;
    1643   }
    1644
    1645   $FormalCharge = $This->GetFormalCharge();
    1646
    1647   # Is any valence model data available for atomic number and formal charge?
    1648   if (!exists $ValenceModelDataRef->{$AtomicNumber}) {
    1649     return 0;
    1650   }
    1651   if (!exists $ValenceModelDataRef->{$AtomicNumber}{$FormalCharge}) {
    1652     return 0;
    1653   }
    1654
    1655   my($PotentialTotalValence, $SumOfBondOrders, $CurrentEffectiveValence, $AvailableCommonValence);
    1656
    1657   $SumOfBondOrders = $This->GetSumOfBondOrders();
    1658   if (!defined $SumOfBondOrders) {
    1659     $SumOfBondOrders = 0;
    1660   }
    1661   $CurrentEffectiveValence = $SumOfBondOrders + $This->GetFreeRadicalElectrons();
    1662
    1663   $PotentialTotalValence = 0;
    1664   VALENCE: for $AvailableCommonValence (@{$ValenceModelDataRef->{$AtomicNumber}{$FormalCharge}{CommonValences}}) {
    1665       if ($CurrentEffectiveValence <= $AvailableCommonValence) {
    1666         $PotentialTotalValence = $AvailableCommonValence;
    1667         last VALENCE;
    1668       }
    1669   }
    1670
    1671   return $PotentialTotalValence;
    1672 }
    1673
    1674 #
    1675 # For elements with one one common valence, potential total common valence used
    1676 # during the calculation for number of implicit hydrogens during InternalValenceMode
    1677 # corresponds to:
    1678 #
    1679 #   CommonValence + FormalCharge - FreeRadicalElectrons
    1680 #
    1681 # For elements with multiple common valences, each common valence is used to
    1682 # calculate total potential common valence as shown above, and the first total potential
    1683 # common valence gerater than the sum of bond orderes is selected as the final total
    1684 # common valence.
    1685 #
    1686 # Group numbers > 14 - Group numbers 15 (N), 16 (O), 17 (F), 18 (He)
    1687 #
    1688 # Formal charge sign is not adjusted. Positive and negative values result in the
    1689 # increase and decrease of valence.
    1690 #
    1691 # Group 14 containing C, Si, Ge, Sn, Pb...
    1692 #
    1693 # Formal charge sign is reversed for positive values. Both positive and negative
    1694 # values result in the decrease of valence.
    1695 #
    1696 # Group 13 containing B, Al, Ga, In, Tl...
    1697 #
    1698 # Formal charge sign is always reversed. Positive and negative values result in the
    1699 # decrease and increase of valence.
    1700 #
    1701 # Groups 1 (H) through 12 (Zn)...
    1702 #
    1703 # Formal charge sign is reversed for positive values. Both positive and negative
    1704 # values result in the decrease of valence.
    1705 #
    1706 # Lanthanides and actinides...
    1707 #
    1708 # Formal charge sign is reversed for positive values. Both positive and negative
    1709 # values result in the decrease of valence.
    1710 #
    1711 # Notes:
    1712 #  . CommonValence and HighestCommonValence available from PeriodicTable module
    1713 #    are equivalent to most common and highest sum of bond orders for an element. For
    1714 #    neutral atoms involved only in single bonds, it corresponds to highest number of
    1715 #    allowed bonds for the atom.
    1716 #  . FormalCharge sign is reversed for electropositive elements with positive formal charge
    1717 #    during common valence calculations. Electropositive elements, metals and transition elements,
    1718 #    have usually plus formal charge and it leads to decrease in common valence; the negative
    1719 #    formal charge should result in the decrease of common valence.
    1720 #  . For carbon, both plus/minus formal charge cause decrease in common valence
    1721 #  . For elements on the right of carbon in periodic table, electronegative elements, plus formal
    1722 #    charge causes common valence to increase and minus formal charge cause it to decrease.
    1723 #
    1724 sub _GetPotentialTotalCommonValenceUsingInternalValenceModel {
    1725   my($This) = @_;
    1726   my($AtomicNumber, $CommonValences);
    1727
    1728   $AtomicNumber = $This->{AtomicNumber};
    1729   if (!$AtomicNumber) {
    1730     return 0;
    1731   }
    1732
    1733   $CommonValences = PeriodicTable::GetElementCommonValences($AtomicNumber);
    1734   if (!$CommonValences) {
    1735     return 0;
    1736   }
    1737
    1738   my($PotentialTotalValence, $AdjustedFormalCharge, $FreeRadicalElectrons, $SumOfBondOrders, $AvailableCommonValence, @AvailableCommonValences);
    1739
    1740   $AdjustedFormalCharge = $This->_GetFormalChargeAdjustedForInternalValenceModel();
    1741   $FreeRadicalElectrons = $This->GetFreeRadicalElectrons();
    1742
    1743   $SumOfBondOrders = $This->GetSumOfBondOrders();
    1744   if (!defined $SumOfBondOrders) {
    1745     $SumOfBondOrders = 0;
    1746   }
    1747
    1748   @AvailableCommonValences = split /\,/, $CommonValences;
    1749
    1750   if (@AvailableCommonValences == 1) {
    1751     # Calculate potential total valence using the only available common valence...
    1752     $PotentialTotalValence = $AvailableCommonValences[0] + $AdjustedFormalCharge - $FreeRadicalElectrons;
    1753   }
    1754   else {
    1755     # Calculate potential total valence using common valence from a list of available valences
    1756     # that makes it higher than sum of bond orders or using the highest common valence...
    1757     VALENCE: for $AvailableCommonValence (@AvailableCommonValences) {
    1758       $PotentialTotalValence = $AvailableCommonValence + $AdjustedFormalCharge - $FreeRadicalElectrons;
    1759
    1760       if ($PotentialTotalValence < 0 || $PotentialTotalValence >= $SumOfBondOrders) {
    1761         last VALENCE;
    1762       }
    1763     }
    1764   }
    1765
    1766   return $PotentialTotalValence > 0 ? $PotentialTotalValence : 0;
    1767 }
    1768
    1769 # Adjust sign of the formal charge for potential total common valence calculation
    1770 # used during internal valence model to figure out number of implicit hydrogens.
    1771 #
    1772 sub _GetFormalChargeAdjustedForInternalValenceModel {
    1773   my($This) = @_;
    1774   my($FormalCharge, $GroupNumber, $SwitchSign);
    1775
    1776   $FormalCharge = $This->GetFormalCharge();
    1777   if ($FormalCharge == 0) {
    1778     return 0;
    1779   }
    1780
    1781   $GroupNumber = $This->GetGroupNumber();
    1782   if (!defined $GroupNumber) {
    1783     return $FormalCharge;
    1784   }
    1785
    1786   # Group numbers > 14 - Group numbers 15 (N), 16 (O), 17 (F), 18 (He)
    1787   #
    1788   # Formal charge sign is not adjusted. Positive and negative values result in the
    1789   # increase and decrease of valence.
    1790   #
    1791   # Group 14 containing C, Si, Ge, Sn, Pb...
    1792   #
    1793   # Formal charge sign is reversed for positive values. Both positive and negative
    1794   # values result in the decrease of valence.
    1795   #
    1796   # Group 13 containing B, Al, Ga, In, Tl...
    1797   #
    1798   # Formal charge sign is always reversed. Positive and negative values result in the
    1799   # decrease and increase of valence.
    1800   #
    1801   # Groups 1 (H) through 12 (Zn)...
    1802   #
    1803   # Formal charge sign is reversed for positive values. Both positive and negative
    1804   # values result in the decrease of valence.
    1805   #
    1806   # Lanthanides and actinides...
    1807   #
    1808   # Formal charge sign is reversed for positive values. Both positive and negative
    1809   # values result in the decrease of valence.
    1810   #
    1811
    1812   $SwitchSign = 0;
    1813   if (length $GroupNumber) {
    1814     GROUPNUMBER: {
    1815       if ($GroupNumber > 14) {
    1816         # Groups on the right side of C group in the periodic table...
    1817         $SwitchSign = 0;
    1818         last GROUPNUMBER;
    1819       }
    1820       if ($GroupNumber == 14) {
    1821         # Group containing C, Si, Ge, Sn, Pb...
    1822         $SwitchSign = ($FormalCharge > 0) ? 1 : 0;
    1823         last GROUPNUMBER;
    1824       }
    1825       if ($GroupNumber == 13) {
    1826         # Group containing B, Al, Ga, In, Tl...
    1827         $SwitchSign = 1;
    1828         last GROUPNUMBER;
    1829       }
    1830       # Groups 1 (H) through 12 (Zn)...
    1831       if ($GroupNumber >=1 && $GroupNumber <= 12) {
    1832         # Groups 1 (H) through 12 (Zn)...
    1833         $SwitchSign = ($FormalCharge > 0) ? 1 : 0;
    1834         last GROUPNUMBER;
    1835       }
    1836     }
    1837   }
    1838   else {
    1839     # Lanthanides and actinides...
    1840     $SwitchSign = ($FormalCharge > 0) ? 1 : 0;
    1841   }
    1842
    1843   if ($SwitchSign) {
    1844     $FormalCharge *= -1.0;
    1845   }
    1846
    1847   return $FormalCharge;
    1848 }
    1849
    1850 # Get lowest common valence...
    1851 sub GetLowestCommonValence {
    1852   my($This) = @_;
    1853
    1854   # Is LowestCommonValence property explicitly set?
    1855   if ($This->HasProperty('LowestCommonValence')) {
    1856     return $This->GetProperty('LowestCommonValence');
    1857   }
    1858   my($AtomicNumber, $LowestCommonValence);
    1859
    1860   $AtomicNumber = $This->{AtomicNumber};
    1861   if (!$AtomicNumber) {
    1862     return 0;
    1863   }
    1864   # Any need to differentiate between internal and other valence models...
    1865
    1866   # LowestCommonValence is not set for all elements...
    1867   $LowestCommonValence = PeriodicTable::GetElementLowestCommonValence($AtomicNumber);
    1868   if (!$LowestCommonValence) {
    1869     $LowestCommonValence = undef;
    1870   }
    1871
    1872   return $LowestCommonValence;
    1873 }
    1874
    1875 # Get highest common valence...
    1876 sub GetHighestCommonValence {
    1877   my($This) = @_;
    1878
    1879   # Is HighestCommonValence property explicitly set?
    1880   if ($This->HasProperty('HighestCommonValence')) {
    1881     return $This->GetProperty('HighestCommonValence');
    1882   }
    1883   my($AtomicNumber, $HighestCommonValence);
    1884
    1885   $AtomicNumber = $This->{AtomicNumber};
    1886   if (!$AtomicNumber) {
    1887     return 0;
    1888   }
    1889
    1890   # Any need to differentiate between internal and other valence models...
    1891
    1892   # HighestCommonValence is not set for all elements...
    1893   $HighestCommonValence = PeriodicTable::GetElementHighestCommonValence($AtomicNumber);
    1894   if (!$HighestCommonValence) {
    1895     $HighestCommonValence = undef;
    1896   }
    1897
    1898   return $HighestCommonValence;
    1899 }
    1900
    1901 # Get valence electrons...
    1902 sub GetValenceElectrons {
    1903   my($This) = @_;
    1904
    1905   # Is ValenceElectrons property explicitly set?
    1906   if ($This->HasProperty('ValenceElectrons')) {
    1907     return $This->GetProperty('ValenceElectrons');
    1908   }
    1909   my($AtomicNumber, $ValenceElectrons);
    1910
    1911   $AtomicNumber = $This->{AtomicNumber};
    1912   if (!$AtomicNumber) {
    1913     return 0;
    1914   }
    1915
    1916   $ValenceElectrons = PeriodicTable::GetElementValenceElectrons($AtomicNumber);
    1917
    1918   return $ValenceElectrons;
    1919 }
    1920
    1921 # Add hydrogens to specified atom in molecule and return number of hydrogens added:
    1922 #
    1923 #   o HydrogensToAdd = ImplicitHydrogenCount - ExplicitHydrogenCount
    1924 #
    1925 #   o XYZ are set to ZeroVector
    1926 #
    1927 sub AddHydrogens {
    1928   my($This, $HydrogenPositionsWarning) = @_;
    1929
    1930   # Is this atom in a molecule?
    1931   if (!$This->HasProperty('Molecule')) {
    1932     return undef;
    1933   }
    1934   if (!defined $HydrogenPositionsWarning) {
    1935     $HydrogenPositionsWarning = 1;
    1936   }
    1937   if ($HydrogenPositionsWarning) {
    1938     carp "Warning: ${ClassName}->AddHydrogens: The current release of MayaChemTools doesn't assign any hydrogen positions...";
    1939   }
    1940
    1941   # Is it an element symbol?
    1942   if (!$This->{AtomicNumber}) {
    1943     return 0;
    1944   }
    1945
    1946   my($Molecule, $HydrogensAdded, $HydrogensToAdd);
    1947
    1948   $Molecule = $This->GetProperty('Molecule');
    1949   $HydrogensAdded = 0;
    1950   $HydrogensToAdd = $This->GetNumOfMissingHydrogens();
    1951   if ($HydrogensToAdd <= 0) {
    1952     return $HydrogensAdded;
    1953   }
    1954
    1955   my($Count, $Hydrogen);
    1956
    1957   for $Count (1 .. $HydrogensToAdd) {
    1958     $HydrogensAdded++;
    1959
    1960     $Hydrogen = $Molecule->NewAtom('AtomSymbol' => 'H', 'XYZ' => [0, 0, 0]);
    1961     $Molecule->NewBond('Atoms' => [$This, $Hydrogen], 'BondOrder' => 1);
    1962   }
    1963
    1964   return $HydrogensAdded;
    1965 }
    1966
    1967 # Delete hydrogens attached to atom in molecule and return total number of hydrogens deleted...
    1968 sub DeleteHydrogens {
    1969   my($This) = @_;
    1970
    1971   # Is this atom in a molecule?
    1972   if (!$This->HasProperty('Molecule')) {
    1973     return undef;
    1974   }
    1975
    1976   # Is it an element symbol?
    1977   if (!$This->{AtomicNumber}) {
    1978     return 0;
    1979   }
    1980
    1981   my($Molecule, $Neighbor, $HydrogensDeleted, @Neighbors);
    1982
    1983   $Molecule = $This->GetProperty('Molecule');
    1984   $HydrogensDeleted = 0;
    1985   @Neighbors = $This->GetNeighbors();
    1986
    1987   NEIGHBOR: for $Neighbor (@Neighbors) {
    1988     if (!$Neighbor->IsHydrogen()) {
    1989       next NEIGHBOR;
    1990     }
    1991     $Molecule->_DeleteAtom($Neighbor);
    1992     $HydrogensDeleted++;
    1993   }
    1994
    1995   return $HydrogensDeleted;
    1996 }
    1997
    1998 # Copy atom and all its associated data...
    1999 sub Copy {
    2000   my($This) = @_;
    2001   my($Atom);
    2002
    2003   $Atom = Storable::dclone($This);
    2004
    2005   return $Atom;
    2006 }
    2007
    2008 # Get atomic invariant value...
    2009 #
    2010 sub GetAtomicInvariantValue {
    2011   my($This, $AtomicInvariant) = @_;
    2012   my($Value);
    2013
    2014   $Value = "";
    2015
    2016   ATOMICVARIANT: {
    2017     if ($AtomicInvariant =~ /^(AS|AtomSymbol|ElementSymbol)$/i) {
    2018       $Value = $This->GetAtomSymbol();
    2019       last ATOMICVARIANT;
    2020     }
    2021     if ($AtomicInvariant =~ /^(X|NumOfNonHydrogenAtomNeighbors|NumOfHeavyAtomNeighbors)$/i) {
    2022       $Value = $This->GetNumOfNonHydrogenAtomNeighbors();
    2023       last ATOMICVARIANT;
    2024     }
    2025     if ($AtomicInvariant =~ /^(BO|SumOfBondOrdersToNonHydrogenAtoms|SumOfBondOrdersToHeavyAtoms)$/i) {
    2026       $Value = $This->GetSumOfBondOrdersToNonHydrogenAtoms();
    2027       last ATOMICVARIANT;
    2028     }
    2029     if ($AtomicInvariant =~ /^(LBO|LargestBondOrderToNonHydrogenAtoms|LargestBondOrderToHeavyAtoms)$/i) {
    2030       $Value = $This->GetLargestBondOrderToNonHydrogenAtoms();
    2031       last ATOMICVARIANT;
    2032     }
    2033     if ($AtomicInvariant =~ /^(H|NumOfImplicitAndExplicitHydrogens)$/i) {
    2034       $Value = $This->GetNumOfHydrogens();
    2035       last ATOMICVARIANT;
    2036     }
    2037     if ($AtomicInvariant =~ /^(SB|NumOfSingleBondsToNonHydrogenAtoms|NumOfSingleBondsToHeavyAtoms)$/i) {
    2038       $Value = $This->GetNumOfSingleBondsToNonHydrogenAtoms();
    2039       last ATOMICVARIANT;
    2040     }
    2041     if ($AtomicInvariant =~ /^(DB|NumOfDoubleBondsToNonHydrogenAtoms|NumOfDoubleBondsToHeavyAtoms)$/i) {
    2042       $Value = $This->GetNumOfDoubleBondsToNonHydrogenAtoms();
    2043       last ATOMICVARIANT;
    2044     }
    2045     if ($AtomicInvariant =~ /^(TB|NumOfTripleBondsToNonHydrogenAtoms|NumOfTripleBondsToHeavyAtoms)$/i) {
    2046       $Value = $This->GetNumOfTripleBondsToNonHydrogenAtoms();
    2047       last ATOMICVARIANT;
    2048     }
    2049     if ($AtomicInvariant =~ /^(AB|NumOfAromaticBondsToNonHydrogenAtoms|NumOfAromaticBondsToHeavyAtoms)$/i) {
    2050       $Value = $This->GetNumOfAromaticBondsToNonHydrogenAtoms();
    2051       last ATOMICVARIANT;
    2052     }
    2053     if ($AtomicInvariant =~ /^(FC|FormalCharge)$/i) {
    2054       $Value = $This->GetFormalCharge();
    2055       $Value = defined $Value ? $Value : 0;
    2056       last ATOMICVARIANT;
    2057     }
    2058     if ($AtomicInvariant =~ /^(T|TotalNumOfAtomNeighbors)$/i) {
    2059       $Value = $This->GetNumOfNonHydrogenAtomNeighbors() + $This->GetNumOfHydrogens();
    2060       last ATOMICVARIANT;
    2061     }
    2062     if ($AtomicInvariant =~ /^(TSB|TotalNumOfSingleBonds)$/i) {
    2063       $Value = $This->GetNumOfSingleBondsToNonHydrogenAtoms() + $This->GetNumOfHydrogens();
    2064       last ATOMICVARIANT;
    2065     }
    2066     if ($AtomicInvariant =~ /^(Ar|Aromatic)$/i) {
    2067       $Value = $This->IsAromatic() ? 1 : 0;
    2068       last ATOMICVARIANT;
    2069     }
    2070     if ($AtomicInvariant =~ /^(RA|RingAtom)$/i) {
    2071       $Value = $This->IsInRing() ? 1 : 0;
    2072       last ATOMICVARIANT;
    2073     }
    2074     if ($AtomicInvariant =~ /^(Str|Stereochemistry)$/i) {
    2075       $Value = $This->GetStereochemistry();
    2076       $Value= (defined($Value) && ($Value =~ /^(R|S)$/i)) ? $Value : '';
    2077       last ATOMICVARIANT;
    2078     }
    2079     if ($AtomicInvariant =~ /^(AN|AtomicNumber)$/i) {
    2080       $Value = $This->GetAtomicNumber();
    2081       last ATOMICVARIANT;
    2082     }
    2083     if ($AtomicInvariant =~ /^(AM|AtomicMass)$/i) {
    2084       $Value = round($This->GetExactMass(), 4) + 0;
    2085       last ATOMICVARIANT;
    2086     }
    2087     if ($AtomicInvariant =~ /^(MN|MassNumber)$/i) {
    2088       $Value = $This->GetMassNumber();
    2089       last ATOMICVARIANT;
    2090     }
    2091     if ($AtomicInvariant =~ /^(SM|SpinMultiplicity)$/i) {
    2092       $Value = $This->GetSpinMultiplicity();
    2093       $Value = defined $Value ? $Value : '';
    2094       last ATOMICVARIANT;
    2095     }
    2096     $Value = "";
    2097     carp "Warning: ${ClassName}->GetAtomicInvariantValue: Unknown atomic invariant $AtomicInvariant...";
    2098   }
    2099
    2100   return $Value;
    2101 }
    2102
    2103 # Get period number of the atom..
    2104 #
    2105 sub GetPeriodNumber {
    2106   my($This) = @_;
    2107
    2108   # Is PeriodNumber property explicitly set?
    2109   if ($This->HasProperty('PeriodNumber')) {
    2110     return $This->GetProperty('PeriodNumber');
    2111   }
    2112   my($AtomicNumber, $PeriodNumber);
    2113
    2114   $AtomicNumber = $This->{AtomicNumber};
    2115   if (!$AtomicNumber) {
    2116     return 0;
    2117   }
    2118
    2119   $PeriodNumber = PeriodicTable::GetElementPeriodNumber($AtomicNumber);
    2120
    2121   return $PeriodNumber;
    2122 }
    2123
    2124 # Get group number of the atom..
    2125 #
    2126 sub GetGroupNumber {
    2127   my($This) = @_;
    2128
    2129   # Is GroupNumber property explicitly set?
    2130   if ($This->HasProperty('GroupNumber')) {
    2131     return $This->GetProperty('GroupNumber');
    2132   }
    2133   my($AtomicNumber, $GroupNumber);
    2134
    2135   $AtomicNumber = $This->{AtomicNumber};
    2136   if (!$AtomicNumber) {
    2137     return 0;
    2138   }
    2139
    2140   $GroupNumber = PeriodicTable::GetElementGroupNumber($AtomicNumber);
    2141
    2142   return $GroupNumber;
    2143 }
    2144
    2145 # Is it a specified topological pharmacophore atom type?
    2146 #
    2147 sub IsTopologicalPharmacophoreType {
    2148   my($This, $Type) = @_;
    2149
    2150   return $This->_IsFunctionalClassType($Type);
    2151 }
    2152
    2153 # Is it a specified functional class atom type?
    2154 #
    2155 sub IsFunctionalClassType {
    2156   my($This, $Type) = @_;
    2157
    2158   return $This->_IsFunctionalClassType($Type);
    2159 }
    2160
    2161 # Is it a specified functional/topological pharmacophore atom type?
    2162 #
    2163 sub _IsFunctionalClassType {
    2164   my($This, $Type) = @_;
    2165   my($Value);
    2166
    2167   $Value = 0;
    2168
    2169   TYPE: {
    2170     if ($Type =~ /^(HBD|HydrogenBondDonor)$/i) {
    2171       $Value = $This->IsHydrogenBondDonor();
    2172       last TYPE;
    2173     }
    2174     if ($Type =~ /^(HBA|HydrogenBondAcceptor)$/i) {
    2175       $Value = $This->IsHydrogenBondAcceptor();
    2176       last TYPE;
    2177     }
    2178     if ($Type =~ /^(PI|PositivelyIonizable)$/i) {
    2179       $Value = $This->IsPositivelyIonizable();
    2180       last TYPE;
    2181     }
    2182     if ($Type =~ /^(NI|NegativelyIonizable)$/i) {
    2183       $Value = $This->IsNegativelyIonizable();
    2184       last TYPE;
    2185     }
    2186     if ($Type =~ /^(H|Hydrophobic)$/i) {
    2187       $Value = $This->IsHydrophobic();
    2188       last TYPE;
    2189     }
    2190     if ($Type =~ /^(Ar|Aromatic)$/i) {
    2191       $Value = $This->IsAromatic();
    2192       last TYPE;
    2193     }
    2194     if ($Type =~ /^(Hal|Halogen)$/i) {
    2195       $Value = $This->IsHalogen();
    2196       last TYPE;
    2197     }
    2198     if ($Type =~ /^(RA|RingAtom)$/i) {
    2199       $Value = $This->IsInRing();
    2200       last TYPE;
    2201     }
    2202     if ($Type =~ /^(CA|ChainAtom)$/i) {
    2203       $Value = $This->IsNotInRing();
    2204       last TYPE;
    2205     }
    2206     $Value = 0;
    2207     carp "Warning: ${ClassName}->_IsType: Unknown functional/pharmacohore type $Type...";
    2208   }
    2209   return $Value;
    2210 }
    2211
    2212 # Is it a Hydrogen atom?
    2213 sub IsHydrogen {
    2214   my($This) = @_;
    2215
    2216   return ($This->{AtomicNumber} == 1) ? 1 : 0;
    2217 }
    2218
    2219 # Is it a Carbon atom?
    2220 sub IsCarbon {
    2221   my($This) = @_;
    2222
    2223   return ($This->{AtomicNumber} == 6) ? 1 : 0;
    2224 }
    2225
    2226 # Is it a Nitrogen atom?
    2227 sub IsNitrogen {
    2228   my($This) = @_;
    2229
    2230   return ($This->{AtomicNumber} == 7) ? 1 : 0;
    2231 }
    2232
    2233 # Is it a Oxygen atom?
    2234 sub IsOxygen {
    2235   my($This) = @_;
    2236
    2237   return ($This->{AtomicNumber} == 8) ? 1 : 0;
    2238 }
    2239
    2240 # Is it a Fluorine atom?
    2241 sub IsFluorine {
    2242   my($This) = @_;
    2243
    2244   return ($This->{AtomicNumber} == 9) ? 1 : 0;
    2245 }
    2246
    2247 # Is it a Silicon atom?
    2248 sub IsSilicon {
    2249   my($This) = @_;
    2250
    2251   return ($This->{AtomicNumber} == 14) ? 1 : 0;
    2252 }
    2253
    2254 # Is it a Phosphorus atom?
    2255 sub IsPhosphorus {
    2256   my($This) = @_;
    2257
    2258   return ($This->{AtomicNumber} == 15) ? 1 : 0;
    2259 }
    2260
    2261 # Is it a Sulphur atom?
    2262 sub IsSulphur {
    2263   my($This) = @_;
    2264
    2265   return $This->IsSulfur();
    2266 }
    2267
    2268 # Is it a Sulfur atom?
    2269 sub IsSulfur {
    2270   my($This) = @_;
    2271
    2272   return ($This->{AtomicNumber} == 16) ? 1 : 0;
    2273 }
    2274
    2275 # Is it a Chlorine atom?
    2276 sub IsChlorine {
    2277   my($This) = @_;
    2278
    2279   return ($This->{AtomicNumber} == 17) ? 1 : 0;
    2280 }
    2281
    2282 # Is it a Arsenic atom?
    2283 sub IsArsenic {
    2284   my($This) = @_;
    2285
    2286   return ($This->{AtomicNumber} == 33) ? 1 : 0;
    2287 }
    2288
    2289 # Is it a Selenium atom?
    2290 sub IsSelenium {
    2291   my($This) = @_;
    2292
    2293   return ($This->{AtomicNumber} == 34) ? 1 : 0;
    2294 }
    2295
    2296 # Is it a Bromine atom?
    2297 sub IsBromine {
    2298   my($This) = @_;
    2299
    2300   return ($This->{AtomicNumber} == 35) ? 1 : 0;
    2301 }
    2302
    2303 # Is it a Tellurium atom?
    2304 sub IsTellurium {
    2305   my($This) = @_;
    2306
    2307   return ($This->{AtomicNumber} == 52) ? 1 : 0;
    2308 }
    2309
    2310 # Is it a Iodine atom?
    2311 sub IsIodine {
    2312   my($This) = @_;
    2313
    2314   return ($This->{AtomicNumber} == 53) ? 1 : 0;
    2315 }
    2316
    2317 # Is it a hetro atom? (N, O, F, P, S, Cl, Br, I)
    2318 sub IsHeteroAtom {
    2319   my($This) = @_;
    2320
    2321   return ($This->{AtomicNumber} =~ /^(7|8|9|15|16|17|35|53)$/) ? 1 : 0;
    2322 }
    2323
    2324 # Is it a halogen atom? (F, Cl, Br, I)
    2325 sub IsHalogen {
    2326   my($This) = @_;
    2327
    2328   return ($This->{AtomicNumber} =~ /^(9|17|35|53)$/) ? 1 : 0;
    2329 }
    2330
    2331 # Is it classified as metallic?
    2332 sub IsMetallic {
    2333   my($This) = @_;
    2334   my($Classification);
    2335
    2336   $Classification = PeriodicTable::GetElementClassification($This->{AtomicNumber});
    2337
    2338   return ($Classification =~ /^Metallic$/i) ? 1 : 0;
    2339 }
    2340
    2341 # Is it a non carbon or hydrogen atom? (C, H)
    2342 sub IsNonCarbonOrHydrogen {
    2343   my($This) = @_;
    2344
    2345   return ($This->{AtomicNumber} =~ /^(1|6)$/) ? 0 : 1;
    2346 }
    2347
    2348 # Is it a polar atom? ( N, O,  P, S)
    2349 sub IsPolarAtom {
    2350   my($This) = @_;
    2351
    2352   return ($This->{AtomicNumber} =~ /^(7|8|15|16)$/) ? 1 : 0;
    2353 }
    2354
    2355 # Is it an isotope?
    2356 sub IsIsotope {
    2357   my($This) = @_;
    2358
    2359   my($AtomicNumber) = $This->{AtomicNumber};
    2360   if (!$AtomicNumber) {
    2361     return 0;
    2362   }
    2363
    2364   if (!$This->HasProperty('MassNumber')) {
    2365     return 0;
    2366   }
    2367   my($MassNumber, $MostAbundantMassNumber);
    2368
    2369   $MassNumber = $This->GetProperty('MassNumber');
    2370   $MostAbundantMassNumber = PeriodicTable::GetElementMostAbundantNaturalIsotopeMassNumber($AtomicNumber);
    2371
    2372   return ($MassNumber == $MostAbundantMassNumber) ? 0 : 1;
    2373 }
    2374
    2375 # Is it a terminal atom?
    2376 sub IsTerminal {
    2377   my($This) = @_;
    2378
    2379   # Is this atom in a molecule?
    2380   if (!$This->HasProperty('Molecule')) {
    2381     return undef;
    2382   }
    2383
    2384   return ($This->GetNumOfNonHydrogenAtomNeighbors() <= 1) ? 1 : 0
    2385
    2386 }
    2387
    2388 # Is aromatic property set for the atom?
    2389 sub IsAromatic {
    2390   my($This) = @_;
    2391   my($Aromatic);
    2392
    2393   $Aromatic = $This->GetAromatic();
    2394
    2395   return (defined($Aromatic) && $Aromatic) ? 1 : 0;
    2396 }
    2397
    2398 # Is this a hydrogen atom and attached to one of these atoms: N, O, P, S
    2399 sub IsPolarHydrogen {
    2400   my($This) = @_;
    2401
    2402   if (!$This->IsHydrogen()) {
    2403     return 0;
    2404   }
    2405
    2406   my(@Bonds);
    2407   @Bonds = $This->GetBonds();
    2408   if (@Bonds > 1) {
    2409     return 0;
    2410   }
    2411
    2412   my($Bond, $BondedAtom);
    2413   ($Bond) = @Bonds;
    2414   $BondedAtom = $Bond->GetBondedAtom($This);
    2415
    2416   return $BondedAtom->IsPolarAtom() ? 1 : 0;
    2417 }
    2418
    2419 # Is it a hydrogen bond donor atom?
    2420 #
    2421 sub IsHBondDonor {
    2422   my($This, $HydrogenBondsType) = @_;
    2423
    2424   return $This->IsHydrogenBondDonor($HydrogenBondsType);
    2425 }
    2426
    2427 # The currrent release of MayaChemTools supports identification of two types of
    2428 # hydrogen bond donor and acceptor atoms with these names:
    2429 #
    2430 # HBondsType1 or HydrogenBondsType1
    2431 # HBondsType2 or HydrogenBondsType2
    2432 #
    2433 # The names of these hydrogen bond types are rather arbitrary. However, their
    2434 # definitions have specific meaning and are as follows:
    2435 #
    2436 # HydrogenBondsType1 [ Ref 60-61, Ref 65-66 ]:
    2437 #   . Donor: NH, NH2, NH3, OH - Any N and O with available H
    2438 #   . Acceptor: N[!H], O - Any N without available H and any O
    2439 #
    2440 # HydrogenBondsType2 [ Ref 91 ]:
    2441 #   . Donor: NH, NH2, NH3, OH - Any N and O with availabe H
    2442 #   . Acceptor: N, O - Any N and O
    2443 #
    2444 # Note:
    2445 #   . HydrogenBondsType2 definition corresponds to Rule of 5.
    2446 #
    2447
    2448 # Is it a hydrogen bond donor atom?
    2449 #
    2450 # The currrent release of MayaChemTools supports identification of two types of
    2451 sub IsHydrogenBondDonor {
    2452   my($This, $HydrogenBondsType) = @_;
    2453   my($Status);
    2454
    2455   $HydrogenBondsType = defined $HydrogenBondsType ? $HydrogenBondsType : 'HBondsType1';
    2456   $Status = 0;
    2457
    2458   HYDROGENBONDSTYPE: {
    2459
    2460       if ($HydrogenBondsType =~ /^(HBondsType1|HydrogenBondsType1)$/i) {
    2461         $Status = $This->_IsHydrogenBondDonorOfType1();
    2462         last HYDROGENBONDSTYPE;
    2463       }
    2464
    2465       if ($HydrogenBondsType =~ /^(HBondsType2|HydrogenBondsType2)$/i) {
    2466         $Status = $This->_IsHydrogenBondDonorOfType2();
    2467         last HYDROGENBONDSTYPE;
    2468       }
    2469
    2470       $Status = 0;
    2471       carp "Warning: ${ClassName}->IsHydrogenBondDonor: The current release of MayaChemTools doesn't support specified value, $HydrogenBondsType, for HydrogenBondsType. Valid values: HBondsType1, HydrogenBondsType1, HBondsType2 HydrogenBondsType2 ...";
    2472   }
    2473
    2474   return $Status;
    2475 }
    2476
    2477 # Is it a MayaChemTools HBondType1 hydrogen bond donor atom?
    2478 #
    2479 sub _IsHydrogenBondDonorOfType1 {
    2480   my($This) = @_;
    2481
    2482   return $This->_IsHydrogenBondDonorOfType1OrType2();
    2483 }
    2484
    2485 # Is it a MayaChemTools HBondType2 hydrogen bond donor atom?
    2486 #
    2487 sub _IsHydrogenBondDonorOfType2 {
    2488   my($This) = @_;
    2489
    2490   return $This->_IsHydrogenBondDonorOfType1OrType2();
    2491 }
    2492
    2493 # Is it a hydrogen bond donor atom of MayaChemTools Type1 or Type2?
    2494 #
    2495 # HydrogenBondDonor definition [ Ref 60-61, Ref 65-66, Ref 91 ]: NH, NH2, OH
    2496 #
    2497 # In other words:
    2498 #   . NH, NH2 - Nitrogen atom with available hydrogen
    2499 #   . OH - Oxygen atom with avilable hydrogen
    2500 #
    2501 sub _IsHydrogenBondDonorOfType1OrType2 {
    2502   my($This) = @_;
    2503
    2504   # Is this atom in a molecule?
    2505   if (!$This->HasProperty('Molecule')) {
    2506     return 0;
    2507   }
    2508
    2509   # Is it N or O?
    2510   if ($This->{AtomicNumber} !~ /^(7|8)$/) {
    2511     return 0;
    2512   }
    2513
    2514   # Any explicitly attached hydrogens?
    2515   if ($This->GetExplicitHydrogens()) {
    2516     return 1;
    2517   }
    2518
    2519   # Any missing hydrogens?
    2520   return $This->GetNumOfMissingHydrogens() ? 1 : 0;
    2521 }
    2522
    2523 # Is it a hydrogen bond acceptor atom?
    2524 #
    2525 sub IsHBondAcceptor {
    2526   my($This, $HydrogenBondsType) = @_;
    2527
    2528   return $This->IsHydrogenBondAcceptor($HydrogenBondsType);
    2529 }
    2530
    2531 # Is it a hydrogen bond acceptor atom?
    2532 #
    2533 sub IsHydrogenBondAcceptor {
    2534   my($This, $HydrogenBondsType) = @_;
    2535   my($Status);
    2536
    2537   $HydrogenBondsType = defined $HydrogenBondsType ? $HydrogenBondsType : 'HBondsType1';
    2538   $Status = 0;
    2539
    2540   HYDROGENBONDSTYPE: {
    2541
    2542       if ($HydrogenBondsType =~ /^(HBondsType1|HydrogenBondsType1)$/i) {
    2543         $Status = $This->_IsHydrogenBondAcceptorOfType1();
    2544         last HYDROGENBONDSTYPE;
    2545       }
    2546
    2547       if ($HydrogenBondsType =~ /^(HBondsType2|HydrogenBondsType2)$/i) {
    2548         $Status = $This->_IsHydrogenBondAcceptorOfType2();
    2549         last HYDROGENBONDSTYPE;
    2550       }
    2551
    2552       $Status = 0;
    2553       carp "Warning: ${ClassName}->IsHydrogenBondAcceptor: The current release of MayaChemTools doesn't support specified value, $HydrogenBondsType, for HydrogenBondsType. Valid values: HBondsType1, HydrogenBondsType1, HBondsType2 HydrogenBondsType2 ...";
    2554   }
    2555
    2556   return $Status;
    2557 }
    2558
    2559 # Is it a MayaChemTools HBondType1 hydrogen bond acceptor atom?
    2560 #
    2561 # HydrogenBondAcceptor definition [ Ref 60-61, Ref 65-66 ]: N[!H], O
    2562 #
    2563 # In other words:
    2564 #   . N[!H] - Nitrogen atom with no hydrogen
    2565 #   . O - Oxygen atom
    2566 #
    2567 sub _IsHydrogenBondAcceptorOfType1 {
    2568   my($This) = @_;
    2569
    2570   # Is this atom in a molecule?
    2571   if (!$This->HasProperty('Molecule')) {
    2572     return 0;
    2573   }
    2574
    2575   # Is it N or O?
    2576   if ($This->{AtomicNumber} !~ /^(7|8)$/) {
    2577     return 0;
    2578   }
    2579
    2580   # Is it O?
    2581   if ($This->{AtomicNumber} == 8 ) {
    2582     return 1;
    2583   }
    2584
    2585   # Any explicitly attached hydrogens?
    2586   if ($This->GetExplicitHydrogens()) {
    2587     return 0;
    2588   }
    2589
    2590   # Any missing hydrogens?
    2591   return $This->GetNumOfMissingHydrogens() ? 0 : 1;
    2592 }
    2593
    2594 # Is it a MayaChemTools HBondType2 hydrogen bond acceptor atom?
    2595 #
    2596 # HydrogenBondAcceptor definition [ Ref 91 ]: N, O
    2597 #
    2598 # In other words:
    2599 #   . Any Nitrogen or Oxygen atom
    2600 #
    2601 # Note:
    2602 #   . HydrogenBondsType2 definition corresponds to Rule of 5.
    2603 #
    2604 sub _IsHydrogenBondAcceptorOfType2 {
    2605   my($This) = @_;
    2606
    2607   # Is this atom in a molecule?
    2608   if (!$This->HasProperty('Molecule')) {
    2609     return 0;
    2610   }
    2611
    2612   return ($This->{AtomicNumber} =~ /^(7|8)$/) ? 1 : 0;
    2613 }
    2614
    2615 # Is it a positively ionizable atom?
    2616 #
    2617 # PositivelyIonizable defintion [ Ref 60-61, Ref 65-66 ]: +, NH2
    2618 #
    2619 # In other words:
    2620 #   . Any atom with positve formal charge
    2621 #   . NH2 - Nitogen atom in amino group
    2622 #
    2623 sub IsPositivelyIonizable {
    2624   my($This) = @_;
    2625   my($FormalCharge);
    2626
    2627   # Is this atom in a molecule?
    2628   if (!$This->HasProperty('Molecule')) {
    2629     return 0;
    2630   }
    2631
    2632   # Any explicit positive formal charge?
    2633   $FormalCharge = $This->GetFormalCharge();
    2634   if (defined($FormalCharge) && $FormalCharge > 0) {
    2635     return 1;
    2636   }
    2637
    2638   # Is it  N?
    2639   if ($This->{AtomicNumber} != 7 ) {
    2640     return 0;
    2641   }
    2642
    2643   return ($This->GetNumOfHydrogens() == 2) ? 1 : 0;
    2644 }
    2645
    2646 # Is it a negatively ionizable atom?
    2647 #
    2648 # NegativelyIonizable definition [ Ref 60-61, Ref 65-66 ]: -, C(=O)OH, S(=O)OH, P(=O)OH
    2649 #
    2650 # In other words:
    2651 #   . Any atom with negative formal charge
    2652 #   . Carbon atom in C(=O)OH group
    2653 #   . Phosphorous in P(=O)OH group
    2654 #   . Sulfur atom in S(=O)OH group
    2655 #
    2656 sub IsNegativelyIonizable {
    2657   my($This) = @_;
    2658   my($FormalCharge);
    2659
    2660   # Is this atom in a molecule?
    2661   if (!$This->HasProperty('Molecule')) {
    2662     return 0;
    2663   }
    2664
    2665   # Any explicit negative formal charge?
    2666   $FormalCharge = $This->GetFormalCharge();
    2667   if (defined($FormalCharge) && $FormalCharge < 0) {
    2668     return 1;
    2669   }
    2670
    2671   # Is it C, P or S?
    2672   if ($This->{AtomicNumber} !~ /^(6|15|16)$/ ) {
    2673     return 0;
    2674   }
    2675
    2676   # Collect oxygens connected to C, P or S with single or double bonds and not connected to
    2677   # any other heavy atom...
    2678   my($Neighbor, $NeighborOxygenBondOrder, $NumOfNeighborOxygensWithSingleBonds, $NumOfNeighborOxygensWithDoubleBonds);
    2679
    2680   $NumOfNeighborOxygensWithSingleBonds = 0; $NumOfNeighborOxygensWithDoubleBonds = 0;
    2681
    2682   NEIGHBOR: for $Neighbor ($This->GetNeighbors()) {
    2683     # Is it an oxygen?
    2684     if ($Neighbor->{AtomicNumber} != 8) {
    2685       next NEIGHBOR;
    2686     }
    2687     # Is oxygent connected to only heavy atom?
    2688     if ($Neighbor->GetNumOfHeavyAtomNeighbors() != 1) {
    2689       next NEIGHBOR;
    2690     }
    2691     $NeighborOxygenBondOrder = $This->GetBondToAtom($Neighbor)->GetBondOrder();
    2692
    2693     if ($NeighborOxygenBondOrder == 2) {
    2694       $NumOfNeighborOxygensWithDoubleBonds++;
    2695     }
    2696     elsif ($NeighborOxygenBondOrder == 1) {
    2697       $NumOfNeighborOxygensWithSingleBonds++;
    2698     }
    2699   }
    2700   return ($NumOfNeighborOxygensWithDoubleBonds >= 1 && $NumOfNeighborOxygensWithSingleBonds >= 1) ? 1 : 0;
    2701 }
    2702
    2703 # Is it a liphophilic atom?
    2704 #
    2705 # Lipophilic definition [ Ref 60-61, Ref 65-66 ]: C(C)(C)(C)(C), Cl, Br, I, S(C)(C)
    2706 #
    2707 # In other words:
    2708 #   . C(C)(C)(C)(C) - Carbon atom connected to only other carbons
    2709 #   . Chlorine, Bromine or Iodine atom
    2710 #   . S(C)(C) - Sulfur connected to two carbons
    2711 #
    2712 sub IsLipophilic {
    2713   my($This) = @_;
    2714
    2715   # Is this atom in a molecule?
    2716   if (!$This->HasProperty('Molecule')) {
    2717     return 0;
    2718   }
    2719
    2720   # Is it Cl, Br, I?
    2721   if ($This->{AtomicNumber} =~ /^(17|35|53)$/) {
    2722     return 1;
    2723   }
    2724
    2725   # Is it C, S?
    2726   if ($This->{AtomicNumber} !~ /^(6|16)$/) {
    2727     return 0;
    2728   }
    2729
    2730   # Are all heavy atom neighbors Carbons?
    2731   my($HeavyAtomNeighbor, @HeavyAtomNeighbors);
    2732   @HeavyAtomNeighbors = ();
    2733   @HeavyAtomNeighbors = $This->GetHeavyAtomNeighbors();
    2734
    2735   for $HeavyAtomNeighbor (@HeavyAtomNeighbors) {
    2736     if ($HeavyAtomNeighbor->{AtomicNumber} != 6) {
    2737       return 0;
    2738     }
    2739   }
    2740
    2741   # Does sulfur has two carbon neighbors?
    2742   if ($This->{AtomicNumber} == 16) {
    2743     if (@HeavyAtomNeighbors != 2) {
    2744       return 0;
    2745     }
    2746   }
    2747   return 1;
    2748 }
    2749
    2750 # Is it hydrophobic?
    2751 #
    2752 sub IsHydrophobic {
    2753   my($This) = @_;
    2754
    2755   return $This->IsLipophilic();
    2756 }
    2757
    2758 # Is it a Nitrogen atom in Guadinium group?
    2759 #
    2760 sub IsGuadiniumNitrogen {
    2761   my($This) = @_;
    2762
    2763   # Is it Nitrogen?
    2764   if (!$This->IsNitrogen()) {
    2765     return 0;
    2766   }
    2767
    2768   # Is it connected to a Guadinium Carbon?
    2769   my($AtomNeighbor);
    2770
    2771   for $AtomNeighbor ($This->GetNonHydrogenAtomNeighbors()) {
    2772     if ($AtomNeighbor->IsGuadiniumCarbon()) {
    2773       return 1;
    2774     }
    2775   }
    2776
    2777   return 0;
    2778 }
    2779
    2780 # Is it a Carbon atom in Guadinium group?
    2781 #
    2782 # Guadinium group definition:
    2783 #
    2784 #   R2N-C(=NR)-(NR2) or R2N-C(=NR2+)-(NR2)
    2785 #
    2786 #   where:
    2787 #      . R = Hydrogens or group of atoms attached through Carbon
    2788 #      . Only one of the three Nitrogens has a double bond to Carbon and has optional
    2789 #        formal charge allowing it to be neutral or charged state
    2790 #
    2791 sub IsGuadiniumCarbon {
    2792   my($This) = @_;
    2793
    2794   # Is it Carbon?
    2795   if (!$This->IsCarbon()) {
    2796     return 0;
    2797   }
    2798
    2799   # Match atom neighborhood...
    2800   my($CentralAtomSpec, @NbrAtomSpecsRef, @NbrBondSpecsRef, @NbrOfNbrAtomSpecsRef);
    2801
    2802   $CentralAtomSpec = 'C.X3.BO4';
    2803   @NbrAtomSpecsRef = ('N.FC0', 'N.FC0', 'N.FC0,N.FC+1');
    2804   @NbrBondSpecsRef = ('-', '-', '=');
    2805   @NbrOfNbrAtomSpecsRef = ('C,H', 'C,H', 'C,H');
    2806
    2807   if ($This->DoesAtomNeighborhoodMatch($CentralAtomSpec, \@NbrAtomSpecsRef, \@NbrBondSpecsRef, \@NbrOfNbrAtomSpecsRef)) {
    2808     return 1;
    2809   }
    2810
    2811   return 0;
    2812 }
    2813
    2814 # Is it a Nitrogen atom in Amide group?
    2815 #
    2816 sub IsAmideNitrogen {
    2817   my($This) = @_;
    2818
    2819   # Is it Nitrogen?
    2820   if (!$This->IsNitrogen()) {
    2821     return 0;
    2822   }
    2823
    2824   # Is it connected to a Amide Carbon?
    2825   my($AtomNeighbor);
    2826
    2827   for $AtomNeighbor ($This->GetNonHydrogenAtomNeighbors()) {
    2828     if ($AtomNeighbor->IsAmideCarbon()) {
    2829       return 1;
    2830     }
    2831   }
    2832
    2833   return 0;
    2834 }
    2835
    2836 # Is it a Carbon atom in Amide group?
    2837 #
    2838 # Amide group definition: R-C(=O)-N(-R')-R''
    2839 #
    2840 #   where:
    2841 #      . R = Hydrogen or groups of atoms attached through Carbon
    2842 #      . R' = Hydrogens or groups of atoms attached through Carbon or hetro atoms
    2843 #      . R'' = Hydrogens or groups of atoms attached through Carbon or hetro atoms
    2844 #
    2845 sub IsAmideCarbon {
    2846   my($This) = @_;
    2847
    2848   # Is this atom in a molecule?
    2849   if (!$This->HasProperty('Molecule')) {
    2850     return 0;
    2851   }
    2852
    2853   # Is it Carbon?
    2854   if (!$This->IsCarbon()) {
    2855     return 0;
    2856   }
    2857
    2858   # Match atom neighborhood...
    2859   my($CentralAtomSpec, @NbrAtomSpecsRef, @NbrBondSpecsRef, @NbrOfNbrAtomSpecsRef);
    2860
    2861   $CentralAtomSpec = 'C.X3.BO4,C.X2.BO3';
    2862   @NbrAtomSpecsRef = ('C,H', 'O', 'N');
    2863   @NbrBondSpecsRef = ('-', '=', '-');
    2864   @NbrOfNbrAtomSpecsRef = ('C,H', 'C', 'C,H,N,O,S,P,F,Cl,Br,I');
    2865
    2866   if ($This->DoesAtomNeighborhoodMatch($CentralAtomSpec, \@NbrAtomSpecsRef, \@NbrBondSpecsRef, \@NbrOfNbrAtomSpecsRef)) {
    2867     return 1;
    2868   }
    2869
    2870   return 0;
    2871 }
    2872
    2873 # Is it a Oxygen atom in Carboxylate group?
    2874 #
    2875 sub IsCarboxylateOxygen {
    2876   my($This) = @_;
    2877
    2878   return $This->_MatchCarboxylateAndOrCarboxylOxygen('Carboxylate');
    2879 }
    2880
    2881 # Is it a Carbon atom in Carboxylate group?
    2882 #
    2883 # Carboxyl group definition: R-C(=O)-O-
    2884 #
    2885 sub IsCarboxylateCarbon {
    2886   my($This) = @_;
    2887
    2888   return $This->_MatchCarboxylateAndOrCarboxylCarbon('Carboxylate');
    2889 }
    2890
    2891 # Is it a Oxygen atom in Carboxyl group?
    2892 #
    2893 sub IsCarboxylOxygen {
    2894   my($This) = @_;
    2895
    2896   return $This->_MatchCarboxylateAndOrCarboxylOxygen('Carboxyl');
    2897 }
    2898
    2899 # Is it a Carbon atom in Carboxyl group?
    2900 #
    2901 # Carboxyl group definition: R-C(=O)-OH
    2902 #
    2903 sub IsCarboxylCarbon {
    2904   my($This) = @_;
    2905
    2906   return $This->_MatchCarboxylateAndOrCarboxylCarbon('Carboxyl');
    2907 }
    2908
    2909 # Match Carboxylate and/or Carboxyl oxygen...
    2910 #
    2911 sub _MatchCarboxylateAndOrCarboxylOxygen {
    2912   my($This, $Mode) = @_;
    2913
    2914   # Is it Oxygen?
    2915   if (!$This->IsOxygen()) {
    2916     return 0;
    2917   }
    2918
    2919   # Is it connected to a Carboxylate Carbon?
    2920   my($AtomNeighbor);
    2921
    2922   for $AtomNeighbor ($This->GetNonHydrogenAtomNeighbors()) {
    2923     if ($AtomNeighbor->_MatchCarboxylateAndOrCarboxylCarbon($Mode)) {
    2924       return 1;
    2925     }
    2926   }
    2927
    2928   return 0;
    2929 }
    2930
    2931 # Match Carboxylate and Carboxyl Carbon
    2932 #
    2933 # Carboxylate group definition: R-C(=O)-O-
    2934 # Carboxyl group definition: R-C(=O)-OH
    2935 #
    2936 #   where:
    2937 #      . R = Hydrogens or groups of atoms attached through Carbon
    2938 #
    2939 sub _MatchCarboxylateAndOrCarboxylCarbon {
    2940   my($This, $Mode) = @_;
    2941
    2942   # Is this atom in a molecule?
    2943   if (!$This->HasProperty('Molecule')) {
    2944     return 0;
    2945   }
    2946
    2947   # Is it Carbon?
    2948   if (!$This->IsCarbon()) {
    2949     return 0;
    2950   }
    2951
    2952   # Match atom neighborhood...
    2953   my($CentralAtomSpec, @NbrAtomSpecsRef, @NbrBondSpecsRef, @NbrOfNbrAtomSpecsRef);
    2954
    2955   $CentralAtomSpec = 'C.X3.BO4,C.X2.BO3';
    2956   MODE: {
    2957     if ($Mode =~ /^Carboxylate$/i) {
    2958       @NbrAtomSpecsRef = ('C,H', 'O', 'O.X1.FC-1');
    2959       last MODE;
    2960     }
    2961     if ($Mode =~ /^Carboxyl$/i) {
    2962       @NbrAtomSpecsRef = ('C,H', 'O', 'O.X1.FC0');
    2963       last MODE;
    2964     }
    2965     if ($Mode =~ /^CarboxylateOrCarboxyl$/i) {
    2966       @NbrAtomSpecsRef = ('C,H', 'O', 'O.X1.FC-1,O.X1.FC0');
    2967       last MODE;
    2968     }
    2969     carp "Warning: ${ClassName}->_MatchCarboxylateAndCarboxylCarbon.: Unknown mode $Mode...";
    2970     return 0;
    2971   }
    2972   @NbrBondSpecsRef = ('-', '=', '-');
    2973   @NbrOfNbrAtomSpecsRef = ('C,H', 'C', 'C');
    2974
    2975   if ($This->DoesAtomNeighborhoodMatch($CentralAtomSpec, \@NbrAtomSpecsRef, \@NbrBondSpecsRef, \@NbrOfNbrAtomSpecsRef)) {
    2976     return 1;
    2977   }
    2978
    2979   return 0;
    2980 }
    2981
    2982 # Is it a Oxygen atom in Phosphate group?
    2983 #
    2984 sub IsPhosphateOxygen {
    2985   my($This) = @_;
    2986
    2987   # Is it Oxygen?
    2988   if (!$This->IsOxygen()) {
    2989     return 0;
    2990   }
    2991
    2992   # Is it connected to a Phosphate Phosphorus?
    2993   my($AtomNeighbor);
    2994
    2995   for $AtomNeighbor ($This->GetNonHydrogenAtomNeighbors()) {
    2996     if ($AtomNeighbor->IsPhosphatePhosphorus()) {
    2997       return 1;
    2998     }
    2999   }
    3000
    3001   return 0;
    3002 }
    3003
    3004 # Is it a Phosphorus atom in Phosphate group?
    3005 #
    3006 # Phosphate group definition: AO-(O=)P(-OA)-OA
    3007 #
    3008 #   where:
    3009 #      . A = Any Groups of atoms including hydrogens
    3010 #
    3011 sub IsPhosphatePhosphorus {
    3012   my($This) = @_;
    3013
    3014   # Is this atom in a molecule?
    3015   if (!$This->HasProperty('Molecule')) {
    3016     return 0;
    3017   }
    3018
    3019   # Is it Phosphorus?
    3020   if (!$This->IsPhosphorus()) {
    3021     return 0;
    3022   }
    3023
    3024   # Match atom neighborhood...
    3025   my($CentralAtomSpec, @NbrAtomSpecsRef, @NbrBondSpecsRef, @NbrOfNbrAtomSpecsRef);
    3026
    3027   $CentralAtomSpec = 'P.X4.BO5';
    3028   @NbrAtomSpecsRef = ('O', 'O', 'O', 'O');
    3029   @NbrBondSpecsRef = ('-', '=', '-', '-');
    3030   @NbrOfNbrAtomSpecsRef = (undef, undef, undef, undef);
    3031
    3032   if ($This->DoesAtomNeighborhoodMatch($CentralAtomSpec, \@NbrAtomSpecsRef, \@NbrBondSpecsRef, \@NbrOfNbrAtomSpecsRef)) {
    3033     return 1;
    3034   }
    3035
    3036   return 0;
    3037 }
    3038
    3039
    3040 # Match central atom and its neighborhood using specified atom and bonds specifications...
    3041 #
    3042 # Let:
    3043 #   AS = Atom symbol corresponding to element symbol, atomic number (#n) or any
    3044 #        atom (A)
    3045 #
    3046 #   X<n>   = Number of non-hydrogen atom neighbors or heavy atoms attached to atom
    3047 #   T<n>   = Total number of atom neighbors including implcit and explicit hydrogens
    3048 #   BO<n> = Sum of bond orders to non-hydrogen atom neighbors or heavy atoms attached to atom
    3049 #   LBO<n> = Largest bond order of non-hydrogen atom neighbors or heavy atoms attached to atom
    3050 #   SB<n> = Number of single bonds to non-hydrogen atom neighbors or heavy atoms attached to atom
    3051 #   TSB<n> = Total number of single bonds to atom neighbors including implcit and explicit hydrogens
    3052 #   DB<n> = Number of double bonds to non-hydrogen atom neighbors or heavy atoms attached to atom
    3053 #   TB<n> = Number of triple bonds to non-hydrogen atom neighbors or heavy atoms attached to atom
    3054 #   H<n>   = Number of implicit and explicit hydrogens for atom
    3055 #   Ar     = Aromatic annotation indicating whether atom is aromatic
    3056 #   RA or RA<n>  = Ring atom annotation indicating whether atom is a ring
    3057 #   TR<n>  = Total number of rings containing atom
    3058 #   FC<+n/-n> = Formal charge assigned to atom
    3059 #   MN<n> = Mass number indicating isotope other than most abundant isotope
    3060 #   SM<n> = Spin multiplicity of atom. Possible values: 1 (singlet), 2 (doublet) or 3 (triplet)
    3061 #
    3062 # Then:
    3063 #
    3064 #   Atom specification corresponds to:
    3065 #
    3066 #     AS.X<n>.T<n>.BO<n>.LBO<n>.<SB><n>.TSB<n>.<DB><n>.<TB><n>.H<n>.Ar.RA<n>.TR<n>FC<+n/-n>.MN<n>.SM<n>
    3067 #
    3068 # Except for AS which is a required atomic invariant in atom specification, all other atomic invariants are
    3069 # optional. For an atom specification to match an atom, the values of all specified atomic invariants must
    3070 # match. Exclamation in from of atomic invariant can be used to negate its effect during the match.
    3071 #
    3072 # A comma delimited atom specification string is used to match any one of the specifed atom specification.
    3073 #
    3074 # Notes:
    3075 #   . During atom specification match to an atom, the first atomic invariant is always assumed to
    3076 #     atom symbol.
    3077 #   . Atom match specfication is based on AtomicInvariantAtomTypes implemented in
    3078 #     AotmTypes::AtomicInvariantAtomType.pm module
    3079 #
    3080 # Examples:
    3081 #     . ('N', 'N', 'N')
    3082 #     . ('N.FC0', 'N.FC0', 'N,N.FC+1.H1')
    3083 #     . ('N.H2', 'N.H2', 'N.H1')
    3084 #     . ('C,N', '!N', '!H')
    3085 #     . ('C,N', 'N.Ar', 'N.R5')
    3086 #
    3087 # Let:
    3088 #   -|1|s|Single = Single bond
    3089 #   =|2|d|Double = Double bond
    3090 #   #|3|t|Triple  = Triple bond
    3091 #   :|1.5|a|Ar|Aromatic = Aromatic bond
    3092 #
    3093 #   @|RB|Ring = Ring bond
    3094 #   ~|*|Any = Any bond
    3095 #
    3096 # Then:
    3097 #
    3098 #   Bond specification corresponds to:
    3099 #
    3100 #     -.:
    3101 #     =.@
    3102 #     Double.Aromatic
    3103 #
    3104 # For a bond specification to match bond between two atoms, the values of all specified bond symbols must
    3105 # match. Exclamation in from of bond symbol can be used to negate its effect during the match.
    3106 #
    3107 # A comma delimited bond specification string is used to match any one of the specifed atom specification.
    3108 #
    3109 sub DoesAtomNeighborhoodMatch {
    3110   my($CentralAtom, $CentralAtomSpec, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef) = @_;
    3111   my($NumOfNbrAtomSpecs, $NumOfNbrBondSpecs, $NumOfNbrOfNbrAtomSpecs);
    3112
    3113   # Is this atom in a molecule?
    3114   if (!$CentralAtom->HasProperty('Molecule')) {
    3115     return 0;
    3116   }
    3117
    3118   $NumOfNbrAtomSpecs = defined $NbrAtomSpecsRef ? scalar @{$NbrAtomSpecsRef} : 0;
    3119   $NumOfNbrBondSpecs = defined $NbrBondSpecsRef ? scalar @{$NbrBondSpecsRef} : 0;
    3120   $NumOfNbrOfNbrAtomSpecs = defined $NbrOfNbrAtomSpecsRef ? scalar @{$NbrOfNbrAtomSpecsRef} : 0;
    3121
    3122   # Validate number of specifications...
    3123   if ($NumOfNbrBondSpecs && ($NumOfNbrAtomSpecs != $NumOfNbrBondSpecs)) {
    3124     carp "Warning: ${ClassName}->DoesAtomNeighborhoodMatch: Number of specified central atom, $NumOfNbrAtomSpecs, and bond, $NumOfNbrBondSpecs, specifications must be same; No neighborhood match performed ...";
    3125     return 0;
    3126   }
    3127
    3128   if ($NumOfNbrOfNbrAtomSpecs && ($NumOfNbrOfNbrAtomSpecs != $NumOfNbrAtomSpecs)) {
    3129     carp "Warning: ${ClassName}->DoesAtomNeighborhoodMatch: Number of specified central atom, $NumOfNbrAtomSpecs, and neighbor of neighbor atoms specifications, $NumOfNbrOfNbrAtomSpecs, must be same; No neighborhood match performed ...";
    3130     return 0;
    3131   }
    3132
    3133   # Sort atom and bond specifications in terms of being most specific to least specific..
    3134   ($NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef) = $CentralAtom->_SortSpecificationsForAtomNeighborhoodMatch($NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef);
    3135
    3136   # Does central atom specification match?
    3137   if (!$CentralAtom->_DoesAtomSpecificationMatch($CentralAtomSpec)) {
    3138     return 0;
    3139   }
    3140
    3141   # No neighbors to match...
    3142   if (!$NumOfNbrAtomSpecs) {
    3143     return 1;
    3144   }
    3145
    3146   # Match neighbors...
    3147   my($NbrSpecsMatched, $NbrSpecCount, $NbrSpecMatchCount, %NbrSpecAlreadyMatchedMap);
    3148
    3149   $NbrSpecCount = $NumOfNbrAtomSpecs;
    3150   $NbrSpecMatchCount = 0;
    3151
    3152   %NbrSpecAlreadyMatchedMap = ();
    3153   ($NbrSpecsMatched, $NbrSpecMatchCount) = $CentralAtom->_MatchAtomNeighborhoodUsingAtomBondSpecs($NbrSpecCount, $NbrSpecMatchCount, \%NbrSpecAlreadyMatchedMap, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef);
    3154
    3155   if ($NbrSpecsMatched) {
    3156     # It's match...
    3157     return 1;
    3158   }
    3159
    3160   # Match central atom's missing hydrogens with any unmatched atom
    3161   # and bond specifications...
    3162   #
    3163   ($NbrSpecsMatched, $NbrSpecMatchCount) = $CentralAtom->_MatchAtomNeighborhoodUsingMissingHydrogens($NbrSpecCount, $NbrSpecMatchCount, \%NbrSpecAlreadyMatchedMap, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef);
    3164
    3165   if ($NbrSpecsMatched) {
    3166     # It's match...
    3167     return 1;
    3168   }
    3169
    3170   # No match...
    3171   return 0;
    3172 }
    3173
    3174 # Match central atom neighborhood atom and bond specifications...
    3175 #
    3176 sub _MatchAtomNeighborhoodUsingAtomBondSpecs {
    3177   my($CentralAtom, $NbrSpecCount, $NbrSpecMatchCount, $NbrSpecAlreadyMatchedRef, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef) = @_;
    3178   my($Index, $NbrAtom, $NbrAtomSpec, $NbrBondSpec, $NbrOfNbrAtom, $NbrOfNbrAtomSpec, $MatchNbrOfNbrAtomSpecs, $NbrSpecsMatched);
    3179
    3180   $MatchNbrOfNbrAtomSpecs = (defined $NbrOfNbrAtomSpecsRef && scalar @{$NbrOfNbrAtomSpecsRef}) ? 1 : 0;
    3181
    3182   $NbrSpecsMatched = 0;
    3183
    3184   # Match central atom's  immediate neighbors atom and bond specifications...
    3185   NBRATOM:  for $NbrAtom ($CentralAtom->GetNeighbors()) {
    3186     NBRATOMSPEC: for $Index (0 .. ($NbrSpecCount - 1)) {
    3187       if (exists $NbrSpecAlreadyMatchedRef->{$Index}) {
    3188         next NBRATOMSPEC;
    3189       }
    3190       $NbrAtomSpec = $NbrAtomSpecsRef->[$Index];
    3191       $NbrBondSpec = $NbrBondSpecsRef->[$Index];
    3192
    3193       $NbrOfNbrAtomSpec = $MatchNbrOfNbrAtomSpecs ? $NbrOfNbrAtomSpecsRef->[$Index] : undef;
    3194
    3195       # Match neighbor atom specification...
    3196       if (!$NbrAtom->_DoesAtomSpecificationMatch($NbrAtomSpec)) {
    3197         next NBRATOMSPEC;
    3198       }
    3199
    3200       # Match central atom to neighbor atom bond specification...
    3201       if (!$CentralAtom->_DoesBondSpecificationMatch($NbrAtom, $NbrBondSpec)) {
    3202         next NBRATOMSPEC;
    3203       }
    3204
    3205       # Match any neighbor of neighbor atom specifications...
    3206       if (defined $NbrOfNbrAtomSpec) {
    3207         # Go over the neighbors of central atom skipping the central atom...
    3208         for $NbrOfNbrAtom ($NbrAtom->GetNeighbors($CentralAtom)) {
    3209           if (!$NbrOfNbrAtom->_DoesAtomSpecificationMatch($NbrOfNbrAtomSpec)) {
    3210             next NBRATOMSPEC;
    3211           }
    3212         }
    3213       }
    3214
    3215       # It's a match for a neighbor atom specification...
    3216       $NbrSpecAlreadyMatchedRef->{$Index} = $Index;
    3217       $NbrSpecMatchCount++;
    3218
    3219       if ($NbrSpecMatchCount == $NbrSpecCount) {
    3220         # It's match...
    3221         $NbrSpecsMatched = 1;
    3222         last NBRATOM;
    3223       }
    3224       # Match next neighbor atom...
    3225       next NBRATOM;
    3226     }
    3227   }
    3228   return ($NbrSpecsMatched, $NbrSpecMatchCount);
    3229 }
    3230
    3231 # Match central atom's missing hydrogens with any unmatched atom and bond
    3232 # specifications...
    3233 #
    3234 sub _MatchAtomNeighborhoodUsingMissingHydrogens {
    3235   my($CentralAtom, $NbrSpecCount, $NbrSpecMatchCount, $NbrSpecAlreadyMatchedRef, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef) = @_;
    3236   my($Index, $NbrAtom, $NbrAtomSpec, $NbrBondSpec, $NumOfMissingHydrogens, $MissingHydrogensIndex, $NbrSpecsMatched, $AtomSpecMatched, $AtomSpec, $AtomSymbol);
    3237
    3238   $NbrSpecsMatched = 0;
    3239
    3240   $NumOfMissingHydrogens = $CentralAtom->GetNumOfMissingHydrogens();
    3241   if (($NbrSpecCount - $NbrSpecMatchCount) > $NumOfMissingHydrogens) {
    3242     # It won't match...
    3243     return ($NbrSpecsMatched, $NbrSpecMatchCount);
    3244   }
    3245
    3246   MISSINGHYDROGENNBR: for $MissingHydrogensIndex (0 .. ($NumOfMissingHydrogens - 1)) {
    3247     NBRATOMSPEC: for $Index (0 .. ($NbrSpecCount - 1)) {
    3248       if (exists $NbrSpecAlreadyMatchedRef->{$Index}) {
    3249         next NBRATOMSPEC;
    3250       }
    3251       $NbrAtomSpec = $NbrAtomSpecsRef->[$Index];
    3252       $NbrBondSpec = $NbrBondSpecsRef->[$Index];
    3253
    3254       $NbrAtomSpec =~ s/ //g;
    3255
    3256       # Match neighbor atom specification hydrogen atom symbol...
    3257       $AtomSpecMatched = 0;
    3258       ATOMSPEC: for $AtomSpec (split /\,/, $NbrAtomSpec) {
    3259         ($AtomSymbol) = split /\./, $AtomSpec;
    3260         if ($AtomSymbol =~ /^(H|A|\*)$/i) {
    3261           $AtomSpecMatched = 1;
    3262           last ATOMSPEC;
    3263         }
    3264       }
    3265       if (!$AtomSpecMatched) {
    3266         next NBRATOMSPEC;
    3267       }
    3268
    3269       # Match neighbor atom bond specification to singal bond...
    3270       if (defined $NbrBondSpec) {
    3271         $NbrBondSpec =~ s/ //g;
    3272         if ($NbrBondSpec !~ /^(-|1|s|Single|\~|\*|Any)/i) {
    3273           next NBRATOMSPEC;
    3274         }
    3275       }
    3276
    3277       # It's a match for a neighbor atom specification...
    3278       $NbrSpecAlreadyMatchedRef->{$Index} = $Index;
    3279       $NbrSpecMatchCount++;
    3280
    3281       if ($NbrSpecMatchCount == $NbrSpecCount) {
    3282         # It's match...
    3283         $NbrSpecsMatched = 1;
    3284         last MISSINGHYDROGENNBR;
    3285       }
    3286       # Match next missing hydrogen neighbor...
    3287       next MISSINGHYDROGENNBR;
    3288     }
    3289   }
    3290
    3291   return ($NbrSpecsMatched, $NbrSpecMatchCount);
    3292 }
    3293
    3294 # Sort atom and bond specifications base on neighbor atom specifications going
    3295 # from most to least specific atom specifications.
    3296 #
    3297 # Atom specifications are sorted at the following two levels:
    3298 #
    3299 #   o By atom specification count with in each specification going from most specific
    3300 #      to least specific, where count is determined by the number of "," in each
    3301 #      specification. Wild card containing specifications are considered least specific
    3302 #      and end up at the end of the sorted list.
    3303 #   o By atomic invariant count with in each sorted list going from most specific to
    3304 #      least specific, where count is determined by the number of "." in each atom
    3305 #      specification.
    3306 #
    3307 #A single atom specification,
    3308 # without any commas in atom specification, is is considered most specific...
    3309 #
    3310 sub _SortSpecificationsForAtomNeighborhoodMatch {
    3311   my($This, $NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef) = @_;
    3312   my($Index, $NeedToSort, $NumOfNbrAtomSpecs, $NbrAtomSpecCount,  $NbrAtomSpecAtomicInvarintCount, $NbrAtomSpecToMatch, $NbrAtomSpec, $FirstAtomicInvariant, $WildCardInNbrAtomSpec, @NbrAtomSpecs, @NbrAtomSpecAtomicInvariants, @SortedNbrAtomSpecs, @SortedNbrBondSpecs, @SortedNbrOfNbrAtomSpecs, %NbrAtomSpecDataMap);
    3313
    3314   $NumOfNbrAtomSpecs = defined $NbrAtomSpecsRef ? scalar @{$NbrAtomSpecsRef} : 0;
    3315
    3316   # Figure out whether sorting is necessary...
    3317   $NeedToSort = 0;
    3318   if ($NumOfNbrAtomSpecs > 1) {
    3319     ATOMSPEC: for $NbrAtomSpecToMatch (@{$NbrAtomSpecsRef}) {
    3320       if ($NbrAtomSpecToMatch =~ /(,|\.|A|\*)/i) {
    3321         $NeedToSort = 1;
    3322         last ATOMSPEC;
    3323       }
    3324     }
    3325   }
    3326   if (!$NeedToSort) {
    3327     # Nothing to do...
    3328     return ($NbrAtomSpecsRef, $NbrBondSpecsRef, $NbrOfNbrAtomSpecsRef);
    3329   }
    3330
    3331   %NbrAtomSpecDataMap = ();
    3332
    3333   for $Index (0 .. ($NumOfNbrAtomSpecs - 1)) {
    3334     $NbrAtomSpecToMatch = $NbrAtomSpecsRef->[$Index];
    3335     $NbrAtomSpecToMatch =~ s/ //g;
    3336
    3337     @NbrAtomSpecs = split /\,/, $NbrAtomSpecToMatch;
    3338     $NbrAtomSpecCount = scalar @NbrAtomSpecs;
    3339
    3340     # Does neighbor specification contains a wild card in atom symbol specification?
    3341     #
    3342     if ($NbrAtomSpecToMatch =~ /(A|\*)/i) {
    3343       $WildCardInNbrAtomSpec = 0;
    3344       NBRATOMSPEC: for $NbrAtomSpec (@NbrAtomSpecs) {
    3345         ($FirstAtomicInvariant) = split /\./, $NbrAtomSpec;
    3346         if ($FirstAtomicInvariant =~ /^!/) {
    3347           $FirstAtomicInvariant =~ s/^!//;
    3348         }
    3349         $WildCardInNbrAtomSpec = $This->_IsWildCardAtomSymbolAtomicInvariant($FirstAtomicInvariant);
    3350         if ($WildCardInNbrAtomSpec) {
    3351           last NBRATOMSPEC;
    3352         }
    3353       }
    3354       if ($WildCardInNbrAtomSpec) {
    3355         # Set NbrAtomSpecCount arbitrarily high to make the spec containing wild
    3356         # card last on the sorted list while maintaining its original order in the list...
    3357         $NbrAtomSpecCount = 999;
    3358       }
    3359     }
    3360
    3361     if (!exists $NbrAtomSpecDataMap{$NbrAtomSpecCount}) {
    3362       %{$NbrAtomSpecDataMap{$NbrAtomSpecCount}} = ();
    3363     }
    3364
    3365     # Use first NbrAtomSpec available in @NbrAtomSpecs to determine atomic invariant count
    3366     # with in each NbrAtomSpecToMatch, as @NbrAtomSpecs derived from $NbrAtomSpecToMatch
    3367     # simply corresponds to a list of possible matches...
    3368     #
    3369     ($NbrAtomSpec) = @NbrAtomSpecs;
    3370     @NbrAtomSpecAtomicInvariants = split /\./, $NbrAtomSpec;
    3371     $NbrAtomSpecAtomicInvarintCount = scalar @NbrAtomSpecAtomicInvariants;
    3372
    3373     if (!exists $NbrAtomSpecDataMap{$NbrAtomSpecCount}{$NbrAtomSpecAtomicInvarintCount}) {
    3374       @{$NbrAtomSpecDataMap{$NbrAtomSpecCount}{$NbrAtomSpecAtomicInvarintCount}} = ();
    3375     }
    3376     push @{$NbrAtomSpecDataMap{$NbrAtomSpecCount}{$NbrAtomSpecAtomicInvarintCount}}, $Index;
    3377
    3378   }
    3379
    3380   @SortedNbrAtomSpecs = (); @SortedNbrBondSpecs = ();
    3381   @SortedNbrOfNbrAtomSpecs = ();
    3382
    3383   for $NbrAtomSpecCount ( sort { $a <=> $b } keys %NbrAtomSpecDataMap) {
    3384     for $NbrAtomSpecAtomicInvarintCount ( sort { $b <=> $a } keys %{$NbrAtomSpecDataMap{$NbrAtomSpecCount}}) {
    3385       for $Index (@{$NbrAtomSpecDataMap{$NbrAtomSpecCount}{$NbrAtomSpecAtomicInvarintCount}}) {
    3386         push @SortedNbrAtomSpecs, $NbrAtomSpecsRef->[$Index];
    3387         if (defined $NbrBondSpecsRef) {
    3388           push @SortedNbrBondSpecs, $NbrBondSpecsRef->[$Index];
    3389         }
    3390         if (defined $NbrOfNbrAtomSpecsRef) {
    3391           push @SortedNbrOfNbrAtomSpecs, $NbrOfNbrAtomSpecsRef->[$Index];
    3392         }
    3393       }
    3394     }
    3395   }
    3396
    3397   return (\@SortedNbrAtomSpecs, defined $NbrBondSpecsRef ? \@SortedNbrBondSpecs : undef, defined $NbrOfNbrAtomSpecsRef ? \@SortedNbrOfNbrAtomSpecs : undef);
    3398 }
    3399
    3400 # Check whether atom matches supported atom specification...
    3401 #
    3402 sub _DoesAtomSpecificationMatch {
    3403   my($This, $AtomSpecificationToMatch) = @_;
    3404   my($AtomSpecification, $AtomicInvariant, $AtomSpecificationMatched, $AtomicInvariantMatched, $FirstMatch);
    3405
    3406   # Anything to match...
    3407   if (!(defined($AtomSpecificationToMatch) && $AtomSpecificationToMatch)) {
    3408     return 1;
    3409   }
    3410
    3411   # Take out any spaces...
    3412   $AtomSpecificationToMatch =~ s/ //g;
    3413
    3414   # Match specified atom specifications. For multiple atom specifications in a comma delimited string,
    3415   # only one atom specification needs to match for a successful match. It's up to the caller to make
    3416   # sure that the specificaton list is ordered from least to most specific atom specification...
    3417   #
    3418   for $AtomSpecification (split /\,/, $AtomSpecificationToMatch) {
    3419     $AtomSpecificationMatched = 1;
    3420     $FirstMatch = 1;
    3421
    3422     # Match all atom symbol atomic invariants...
    3423     ATOMICINVARIANT: for $AtomicInvariant (split /\./, $AtomSpecification) {
    3424       if ($FirstMatch) {
    3425         # Match atom symbol atomic invariant...
    3426         $FirstMatch = 0;
    3427         $AtomicInvariantMatched = $This->_MatchAtomSymbolAtomicInvariant($AtomicInvariant);
    3428       }
    3429       else {
    3430         # Match non atom symbol atomic invariant...
    3431         $AtomicInvariantMatched = $This->_MatchNonAtomSymbolAtomicInvariant($AtomicInvariant);
    3432       }
    3433
    3434       if (!$AtomicInvariantMatched) {
    3435         # No need to match other atomic invariants...
    3436         $AtomSpecificationMatched = 0;
    3437         last ATOMICINVARIANT;
    3438       }
    3439     }
    3440
    3441     if ($AtomSpecificationMatched) {
    3442       # No need to match other atom specifications...
    3443       return 1;
    3444     }
    3445   }
    3446
    3447   # Nothing matched...
    3448   return 0;
    3449 }
    3450
    3451 # Check whether atom matches atom symbol atomic invariant...
    3452 #
    3453 sub _MatchAtomSymbolAtomicInvariant {
    3454   my($This, $AtomicInvariant) = @_;
    3455   my($NegateMatch, $Status, $AtomicNumber);
    3456
    3457   $Status = 0;
    3458   $NegateMatch = 0;
    3459
    3460   # Does match needs to be negated?
    3461   if ($AtomicInvariant =~ /^!/) {
    3462     $NegateMatch = 1;
    3463     $AtomicInvariant =~ s/^!//;
    3464   }
    3465
    3466   ATOMICINVARIANT: {
    3467     # Any atom match...
    3468     if ($This->_IsWildCardAtomSymbolAtomicInvariant($AtomicInvariant)) {
    3469       $Status = 1;
    3470       last ATOMICINVARIANT;
    3471     }
    3472
    3473     # Atomic number match...
    3474     if ($AtomicInvariant =~ /^#/) {
    3475       $AtomicNumber = $AtomicInvariant; $AtomicNumber =~ s/^#//;
    3476       $Status = ($This->{AtomicNumber} == $AtomicNumber) ? 1 : 0;
    3477       last ATOMICINVARIANT;
    3478     }
    3479
    3480     # Atom symbol match...
    3481     $Status = ($This->{AtomSymbol} =~ /^$AtomicInvariant$/i) ? 1 : 0;
    3482   }
    3483
    3484   if ($NegateMatch) {
    3485     $Status = $Status ? 0 : 1;
    3486   }
    3487
    3488   return $Status;
    3489 }
    3490
    3491 # Is it a wild card atom symbol atomic invariant?
    3492 #
    3493 sub _IsWildCardAtomSymbolAtomicInvariant {
    3494   my($This, $AtomicInvariant) = @_;
    3495
    3496   return ($AtomicInvariant =~ /^(A|\*)$/i) ? 1 : 0;
    3497 }
    3498
    3499 # Check whether atom matches non atom symbol atomic invariants...
    3500 #
    3501 sub _MatchNonAtomSymbolAtomicInvariant {
    3502   my($This, $AtomicInvariant) = @_;
    3503   my($NegateMatch, $Status, $Name, $Value, $UnknownName);
    3504
    3505   ($Status, $NegateMatch, $UnknownName) = ('0') x 3;
    3506
    3507   # Does match needs to be negated?
    3508   if ($AtomicInvariant =~ /^!/) {
    3509     $NegateMatch = 1;
    3510     $AtomicInvariant =~ s/^!//;
    3511   }
    3512
    3513   # Extract atomic invariant name and any value...
    3514   if ($AtomicInvariant =~ /[0-9\*]+/) {
    3515     ($Name, $Value) = $AtomicInvariant =~ /^([a-zA-Z]+)([0-9\-\+\*\>\<\=]+)$/;
    3516   }
    3517   else {
    3518     ($Name, $Value) = ($AtomicInvariant, undef);
    3519   }
    3520
    3521   NAME: {
    3522     # Match number of non-hydrogen atom neighbors
    3523     if ($Name =~ /^X$/i) {
    3524       $Status = (defined($Value) && $This->GetNumOfNonHydrogenAtomNeighbors() == $Value) ? 1 : 0;
    3525       last NAME;
    3526     }
    3527
    3528     # Match total number of atom neighbors including missing hydrogens...
    3529     if ($Name =~ /^T$/i) {
    3530       $Status = (defined($Value) && ($This->GetNumOfNonHydrogenAtomNeighbors() + $This->GetNumOfHydrogens()) == $Value) ? 1 : 0;
    3531       last NAME;
    3532     }
    3533
    3534     # Match formal charge...
    3535     if ($Name =~ /^FC$/i) {
    3536       my $FormalCharge = $This->GetFormalCharge();
    3537       $Status = $This->_MatchNonAtomSymbolAtomicInvariantValue($FormalCharge, $Value);
    3538       last NAME;
    3539     }
    3540
    3541     # Match aromatic annotation indicating whether atom is aromatic...
    3542     if ($Name =~ /^Ar$/i) {
    3543       $Status = $This->IsAromatic() ? 1 : 0;
    3544       last NAME;
    3545     }
    3546
    3547     # Match number of implicit and explicit hydrogens...
    3548     if ($Name =~ /^H$/i) {
    3549       $Status = (defined($Value) && ($This->GetNumOfHydrogens() == $Value)) ? 1 : 0;
    3550       last NAME;
    3551     }
    3552
    3553     # Match ring atom annotation indicating whether atom is in ring...
    3554     if ($Name =~ /^RA$/i) {
    3555       $Status = defined($Value) ? $This->IsInRingOfSize($Value) : ($This->IsInRing() ? 1 : 0);
    3556       last NAME;
    3557     }
    3558
    3559     # Match number of rings for atom..
    3560     if ($Name =~ /^TR$/i) {
    3561       $Status = (defined($Value) && ($Value == $This->GetNumOfRings())) ? 1 : 0;
    3562       last NAME;
    3563     }
    3564
    3565     # Match sum of bond orders to non-hydrogen atom neighbors...
    3566     if ($Name =~ /^BO$/i) {
    3567       $Status = (defined($Value) && $This->GetSumOfBondOrdersToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3568       last NAME;
    3569     }
    3570
    3571     # Match largest bond order of non-hydrogen atom neighbors...
    3572     if ($Name =~ /^LBO$/i) {
    3573       $Status = (defined($Value) && $This->GetLargestBondOrderToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3574       last NAME;
    3575     }
    3576
    3577     # Match number of single bonds to non-hydrogen atom neighbors...
    3578     if ($Name =~ /^SB$/i) {
    3579       $Status = (defined($Value) && $This->GetNumOfSingleBondsToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3580       last NAME;
    3581     }
    3582
    3583     # Match total number of single bonds to atom neighbors including missing and explicit hydrogens...
    3584     if ($Name =~ /^TSB$/i) {
    3585       $Status = (defined($Value) && ($This->GetNumOfSingleBondsToNonHydrogenAtoms() + $This->GetNumOfHydrogens()) == $Value) ? 1 : 0;
    3586       last NAME;
    3587     }
    3588
    3589     # Match number of double bonds to non-hydrogen atom neighbors...
    3590     if ($Name =~ /^DB$/i) {
    3591       $Status = (defined($Value) && $This->GetNumOfDoubleBondsToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3592       last NAME;
    3593     }
    3594
    3595     # Match number of triple bonds to non-hydrogen atom neighbors...
    3596     if ($Name =~ /^TB$/i) {
    3597       $Status = (defined($Value) && $This->GetNumOfTripleBondsToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3598       last NAME;
    3599     }
    3600
    3601     # Match number of aromatic bonds to non-hydrogen atom neighbors...
    3602     if ($Name =~ /^AB$/i) {
    3603       $Status = (defined($Value) && $This->GetNumOfAromaticBondsToNonHydrogenAtoms() == $Value) ? 1 : 0;
    3604       last NAME;
    3605     }
    3606
    3607
    3608     # Match mass number indicating isotope other than most abundant isotope...
    3609     if ($Name =~ /^MN$/i) {
    3610       $Status = (defined($Value) && $This->GetMassNumber() == $Value) ? 1 : 0;
    3611       last NAME;
    3612     }
    3613
    3614     # Match spin multiplicity...
    3615     if ($Name =~ /^SM$/i) {
    3616       my $SpinMultiplicity = $This->GetSpinMultiplicity();
    3617       if (!defined $SpinMultiplicity) { $SpinMultiplicity = 0; }
    3618       $Status = (defined($Value) && defined($SpinMultiplicity) && $Value == $SpinMultiplicity) ? 1 : 0;
    3619       last NAME;
    3620     }
    3621
    3622     $UnknownName = 1;
    3623     carp "Warning: ${ClassName}->_MatchNonAtomSymbolAtomicInvariant: Unknown atomic invariant $AtomicInvariant...";
    3624   }
    3625
    3626   if (!$UnknownName) {
    3627     if ($NegateMatch) {
    3628       $Status = $Status ? 0 : 1;
    3629     }
    3630   }
    3631
    3632   return $Status;
    3633 }
    3634
    3635 # Match atomic invariant value...
    3636 #
    3637 # Specified value format:
    3638 #   . +* : Any positive value
    3639 #   . -* : Any negative value
    3640 #   . >ValidNumber or >=ValidNumber
    3641 #   . <ValidNumber or <=ValidNumber
    3642 #   . Any valid number
    3643 #
    3644 sub _MatchNonAtomSymbolAtomicInvariantValue {
    3645   my($This, $TargetValue, $SpecifiedValue) = @_;
    3646   my($Status);
    3647
    3648   $Status = 0;
    3649
    3650   if (!(defined($TargetValue) && defined($SpecifiedValue))) {
    3651     return $Status;
    3652   }
    3653
    3654   VALUE: {
    3655     if ($SpecifiedValue =~ /^\+\*/) {
    3656       $Status = ($TargetValue > 0) ? 1 : 0;
    3657       last VALUE;
    3658     }
    3659     if ($SpecifiedValue =~ /^\-\*/) {
    3660       $Status = ($TargetValue < 0) ? 1 : 0;
    3661       last VALUE;
    3662     }
    3663     if ($SpecifiedValue =~ /^>/) {
    3664       if ($SpecifiedValue =~ /^>=/) {
    3665         $SpecifiedValue =~ s/^>=//;
    3666         $Status = ($SpecifiedValue >= $TargetValue) ? 1 : 0;
    3667       }
    3668       else {
    3669         $SpecifiedValue =~ s/^>//;
    3670         $Status = ($SpecifiedValue > $TargetValue) ? 1 : 0;
    3671       }
    3672       last VALUE;
    3673     }
    3674     if ($SpecifiedValue =~ /^</) {
    3675       if ($SpecifiedValue =~ /^<=/) {
    3676         $SpecifiedValue =~ s/^<=//;
    3677         $Status = ($SpecifiedValue <= $TargetValue) ? 1 : 0;
    3678       }
    3679       else {
    3680         $SpecifiedValue =~ s/^<//;
    3681         $Status = ($SpecifiedValue < $TargetValue) ? 1 : 0;
    3682       }
    3683       last VALUE;
    3684     }
    3685     # Default is do perform an equality match...
    3686     $Status = ($SpecifiedValue == $TargetValue) ? 1 : 0;
    3687   }
    3688
    3689   return $Status;
    3690 }
    3691
    3692 # Check whether atoms match  bond specifications...
    3693 #
    3694 sub _DoesBondSpecificationMatch {
    3695   my($This, $BondedAtom, $BondSpecificationToMatch) = @_;
    3696   my($BondSpecification, $BondSymbolSpecification, $BondSpecificationMatched);
    3697
    3698   # Anything to match...
    3699   if (!(defined($BondSpecificationToMatch) && $BondSpecificationToMatch)) {
    3700     return 1;
    3701   }
    3702
    3703   # Take out any spaces...
    3704   $BondSpecificationToMatch =~ s/ //g;
    3705
    3706   # Match specified bond specifications. For multiple bond specifications in a comma delimited string,
    3707   # only one bond specification needs to match for a successful match...
    3708   #
    3709   for $BondSpecification (split /\,/, $BondSpecificationToMatch) {
    3710     $BondSpecificationMatched = 1;
    3711
    3712     # Match all specified bond symbol specifications...
    3713     BONDSYMBOL: for $BondSymbolSpecification (split /\./, $BondSpecification) {
    3714       if (!$This->_MatchBondSymbolSpecification($BondedAtom, $BondSymbolSpecification)) {
    3715         # No need to match other bond symbol specifications...
    3716         $BondSpecificationMatched = 0;
    3717         last BONDSYMBOL;
    3718       }
    3719     }
    3720     if ($BondSpecificationMatched) {
    3721       # No need to try matching other bond specifications...
    3722       return 1;
    3723     }
    3724   }
    3725
    3726   # Nothing matched...
    3727   return 0;
    3728 }
    3729
    3730 # Check whether atoms match  bond symbol specification...
    3731 #
    3732 sub _MatchBondSymbolSpecification {
    3733   my($This, $BondedAtom, $BondSymbolSpecification) = @_;
    3734   my($NegateMatch, $Status, $Bond, $BondSymbol, $UnknownBondSymbol);
    3735
    3736   ($Status, $NegateMatch, $UnknownBondSymbol) = ('0') x 3;
    3737
    3738   # Does match needs to be negated?
    3739   if ($BondSymbolSpecification =~ /^!/) {
    3740     $NegateMatch = 1;
    3741     $BondSymbolSpecification =~ s/^!//;
    3742   }
    3743   $BondSymbol = $BondSymbolSpecification;
    3744   $Bond = $This->GetBondToAtom($BondedAtom);
    3745
    3746   BONDSYMBOL: {
    3747     if ($BondSymbol =~ /^(-|1|s|Single)$/i) { $Status = $Bond->IsSingle() ? 1 : 0; last BONDSYMBOL; }
    3748     if ($BondSymbol =~ /^(=|2|d|Double)$/i) { $Status = $Bond->IsDouble() ? 1 : 0; last BONDSYMBOL; }
    3749     if ($BondSymbol =~ /^(#|3|t|Triple)$/i) { $Status = $Bond->IsTriple() ? 1 : 0; last BONDSYMBOL; }
    3750     if ($BondSymbol =~ /^(:|a|Ar|Aromatic)$/i) { $Status = $Bond->IsAromatic() ? 1 : 0; last BONDSYMBOL; }
    3751
    3752     if ($BondSymbol =~ /^(\@|RB|Ring)$/i) { $Status = $Bond->IsInRing() ? 1 : 0; last BONDSYMBOL; }
    3753
    3754     if ($BondSymbol =~ /^(\~|\*|Any)$/i) { $Status = 1; last BONDSYMBOL; }
    3755
    3756     $UnknownBondSymbol = 1;
    3757     carp "Warning: ${ClassName}->_MatchBondSpecification: Unknown bond specification $BondSymbolSpecification...";
    3758   }
    3759
    3760   if (!$UnknownBondSymbol) {
    3761     if ($NegateMatch) {
    3762       $Status = $Status ? 0 : 1;
    3763     }
    3764   }
    3765
    3766   return $Status;
    3767 }
    3768
    3769 # Is it a saturated atom?
    3770 #
    3771 sub IsSaturated {
    3772   my($This) = @_;
    3773
    3774   return !$This->IsUnsaturated();
    3775 }
    3776
    3777 # Is it an unsaturated atom containing at least one non-single bond?
    3778 #
    3779 sub IsUnsaturated {
    3780   my($This) = @_;
    3781   my($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfAromaticBonds);
    3782
    3783   ($NumOfSingleBonds, $NumOfDoubleBonds, $NumOfTripleBonds, $NumOfAromaticBonds) = $This->GetNumOfBondTypesToNonHydrogenAtoms();
    3784
    3785   return ($NumOfDoubleBonds || $NumOfTripleBonds || $NumOfAromaticBonds) ? 1 : 0;
    3786 }
    3787
    3788 # Is atom in a ring?
    3789 #
    3790 sub IsInRing {
    3791   my($This) = @_;
    3792
    3793   # Is this atom in a molecule?
    3794   if (!$This->HasProperty('Molecule')) {
    3795     return undef;
    3796   }
    3797   my($Molecule);
    3798   $Molecule = $This->GetProperty('Molecule');
    3799
    3800   return $Molecule->_IsAtomInRing($This);
    3801 }
    3802
    3803 # Is atom not in a ring?
    3804 #
    3805 sub IsNotInRing {
    3806   my($This) = @_;
    3807
    3808   # Is this atom in a molecule?
    3809   if (!$This->HasProperty('Molecule')) {
    3810     return undef;
    3811   }
    3812   my($Molecule);
    3813   $Molecule = $This->GetProperty('Molecule');
    3814
    3815   return $Molecule->_IsAtomNotInRing($This);
    3816 }
    3817
    3818 # Is atom only in one ring?
    3819 #
    3820 sub IsOnlyInOneRing {
    3821   my($This) = @_;
    3822
    3823   # Is this atom in a molecule?
    3824   if (!$This->HasProperty('Molecule')) {
    3825     return undef;
    3826   }
    3827   my($Molecule);
    3828   $Molecule = $This->GetProperty('Molecule');
    3829
    3830   return $Molecule->_IsAtomInOnlyOneRing($This);
    3831 }
    3832
    3833 # Is atom in a ring of specific size?
    3834 #
    3835 sub IsInRingOfSize {
    3836   my($This, $RingSize) = @_;
    3837
    3838   # Is this atom in a molecule?
    3839   if (!$This->HasProperty('Molecule')) {
    3840     return undef;
    3841   }
    3842   my($Molecule);
    3843   $Molecule = $This->GetProperty('Molecule');
    3844
    3845   return $Molecule->_IsAtomInRingOfSize($This, $RingSize);
    3846 }
    3847
    3848 # Get size of smallest ring containing the atom...
    3849 #
    3850 sub GetSizeOfSmallestRing {
    3851   my($This) = @_;
    3852
    3853   # Is this atom in a molecule?
    3854   if (!$This->HasProperty('Molecule')) {
    3855     return undef;
    3856   }
    3857   my($Molecule);
    3858   $Molecule = $This->GetProperty('Molecule');
    3859
    3860   return $Molecule->_GetSizeOfSmallestAtomRing($This);
    3861 }
    3862
    3863 # Get size of largest ring containing the atom...
    3864 #
    3865 sub GetSizeOfLargestRing {
    3866   my($This) = @_;
    3867
    3868   # Is this atom in a molecule?
    3869   if (!$This->HasProperty('Molecule')) {
    3870     return undef;
    3871   }
    3872   my($Molecule);
    3873   $Molecule = $This->GetProperty('Molecule');
    3874
    3875   return $Molecule->_GetSizeOfLargestAtomRing($This);
    3876 }
    3877
    3878 # Get number of  rings containing the atom...
    3879 #
    3880 sub GetNumOfRings {
    3881   my($This) = @_;
    3882
    3883   # Is this atom in a molecule?
    3884   if (!$This->HasProperty('Molecule')) {
    3885     return undef;
    3886   }
    3887   my($Molecule);
    3888   $Molecule = $This->GetProperty('Molecule');
    3889
    3890   return $Molecule->_GetNumOfAtomRings($This);
    3891 }
    3892
    3893 # Get number of  rings with odd size containing the atom...
    3894 #
    3895 sub GetNumOfRingsWithOddSize {
    3896   my($This) = @_;
    3897
    3898   # Is this atom in a molecule?
    3899   if (!$This->HasProperty('Molecule')) {
    3900     return undef;
    3901   }
    3902   my($Molecule);
    3903   $Molecule = $This->GetProperty('Molecule');
    3904
    3905   return $Molecule->_GetNumOfAtomRingsWithOddSize($This);
    3906 }
    3907
    3908 # Get number of  rings with even size containing the atom...
    3909 #
    3910 sub GetNumOfRingsWithEvenSize {
    3911   my($This) = @_;
    3912
    3913   # Is this atom in a molecule?
    3914   if (!$This->HasProperty('Molecule')) {
    3915     return undef;
    3916   }
    3917   my($Molecule);
    3918   $Molecule = $This->GetProperty('Molecule');
    3919
    3920   return $Molecule->_GetNumOfAtomRingsWithEvenSize($This);
    3921 }
    3922
    3923 # Get number of  rings with specified size containing the atom...
    3924 #
    3925 sub GetNumOfRingsWithSize {
    3926   my($This, $RingSize) = @_;
    3927
    3928   # Is this atom in a molecule?
    3929   if (!$This->HasProperty('Molecule')) {
    3930     return undef;
    3931   }
    3932   my($Molecule);
    3933   $Molecule = $This->GetProperty('Molecule');
    3934
    3935   return $Molecule->_GetNumOfAtomRingsWithSize($This, $RingSize);
    3936
    3937 }
    3938
    3939 # Get number of  rings with size less than specified containing the atom...
    3940 #
    3941 sub GetNumOfRingsWithSizeLessThan {
    3942   my($This, $RingSize) = @_;
    3943
    3944   # Is this atom in a molecule?
    3945   if (!$This->HasProperty('Molecule')) {
    3946     return undef;
    3947   }
    3948   my($Molecule);
    3949   $Molecule = $This->GetProperty('Molecule');
    3950
    3951   return $Molecule->_GetNumOfAtomRingsWithSizeLessThan($This, $RingSize);
    3952 }
    3953
    3954 # Get number of  rings with size greater than specified size containing the atom...
    3955 #
    3956 sub GetNumOfRingsWithSizeGreaterThan {
    3957   my($This, $RingSize) = @_;
    3958
    3959   # Is this atom in a molecule?
    3960   if (!$This->HasProperty('Molecule')) {
    3961     return undef;
    3962   }
    3963   my($Molecule);
    3964   $Molecule = $This->GetProperty('Molecule');
    3965
    3966   return $Molecule->_GetNumOfAtomRingsWithSizeGreaterThan($This, $RingSize);
    3967 }
    3968
    3969 # Get all rings an array of references to arrays containing ring atoms...
    3970 #
    3971 sub GetRings {
    3972   my($This) = @_;
    3973
    3974   # Is this atom in a molecule?
    3975   if (!$This->HasProperty('Molecule')) {
    3976     return undef;
    3977   }
    3978   my($Molecule);
    3979   $Molecule = $This->GetProperty('Molecule');
    3980
    3981   return $Molecule->_GetAtomRings($This);
    3982 }
    3983
    3984 # Get smallest ring as an array containing ring atoms...
    3985 #
    3986 sub GetSmallestRing {
    3987   my($This) = @_;
    3988
    3989   # Is this atom in a molecule?
    3990   if (!$This->HasProperty('Molecule')) {
    3991     return undef;
    3992   }
    3993   my($Molecule);
    3994   $Molecule = $This->GetProperty('Molecule');
    3995
    3996   return $Molecule->_GetSmallestAtomRing($This);
    3997 }
    3998
    3999 # Get largest ring as an array containing ring atoms...
    4000 #
    4001 sub GetLargestRing {
    4002   my($This) = @_;
    4003
    4004   # Is this atom in a molecule?
    4005   if (!$This->HasProperty('Molecule')) {
    4006     return undef;
    4007   }
    4008   my($Molecule);
    4009   $Molecule = $This->GetProperty('Molecule');
    4010
    4011   return $Molecule->_GetLargestAtomRing($This);
    4012 }
    4013
    4014 # Get odd size rings an array of references to arrays containing ring atoms...
    4015 #
    4016 sub GetRingsWithOddSize {
    4017   my($This) = @_;
    4018
    4019   # Is this atom in a molecule?
    4020   if (!$This->HasProperty('Molecule')) {
    4021     return undef;
    4022   }
    4023   my($Molecule);
    4024   $Molecule = $This->GetProperty('Molecule');
    4025
    4026   return $Molecule->_GetAtomRingsWithOddSize($This);
    4027 }
    4028
    4029 # Get even size rings an array of references to arrays containing ring atoms...
    4030 #
    4031 sub GetRingsWithEvenSize {
    4032   my($This) = @_;
    4033
    4034   # Is this atom in a molecule?
    4035   if (!$This->HasProperty('Molecule')) {
    4036     return undef;
    4037   }
    4038   my($Molecule);
    4039   $Molecule = $This->GetProperty('Molecule');
    4040
    4041   return $Molecule->_GetAtomRingsWithEvenSize($This);
    4042 }
    4043
    4044 # Get rings with specified size as an array of references to arrays containing ring atoms...
    4045 #
    4046 sub GetRingsWithSize {
    4047   my($This, $RingSize) = @_;
    4048
    4049   # Is this atom in a molecule?
    4050   if (!$This->HasProperty('Molecule')) {
    4051     return undef;
    4052   }
    4053   my($Molecule);
    4054   $Molecule = $This->GetProperty('Molecule');
    4055
    4056   return $Molecule->_GetAtomRingsWithSize($This, $RingSize);
    4057 }
    4058
    4059 # Get rings with size less than specfied size as an array of references to arrays containing ring atoms...
    4060 #
    4061 sub GetRingsWithSizeLessThan {
    4062   my($This, $RingSize) = @_;
    4063
    4064   # Is this atom in a molecule?
    4065   if (!$This->HasProperty('Molecule')) {
    4066     return undef;
    4067   }
    4068   my($Molecule);
    4069   $Molecule = $This->GetProperty('Molecule');
    4070
    4071   return $Molecule->_GetAtomRingsWithSizeLessThan($This, $RingSize);
    4072 }
    4073
    4074 # Get rings with size greater than specfied size as an array of references to arrays containing ring atoms...
    4075 #
    4076 sub GetRingsWithSizeGreaterThan {
    4077   my($This, $RingSize) = @_;
    4078
    4079   # Is this atom in a molecule?
    4080   if (!$This->HasProperty('Molecule')) {
    4081     return undef;
    4082   }
    4083   my($Molecule);
    4084   $Molecule = $This->GetProperty('Molecule');
    4085
    4086   return $Molecule->_GetAtomRingsWithSizeGreaterThan($This, $RingSize);
    4087 }
    4088
    4089 # Get next object ID...
    4090 sub _GetNewObjectID {
    4091   $ObjectID++;
    4092   return $ObjectID;
    4093 }
    4094
    4095 # Return a string containing vertices, edges and other properties...
    4096 sub StringifyAtom {
    4097   my($This) = @_;
    4098   my($AtomString, $ID, $Name, $AtomSymbol, $AtomicNumber, $XYZVector, $AtomicWeight, $ExactMass, $NumOfNeighbors, $NumOfBonds, $Valence, $MissingHydrogens, $TotalHydrogens, $ImplicitHydrogens, $ExplicitHydrogens, $FormalCharge, $Charge, $SpinMultiplicity, $FreeRadicalElectrons, $StereoCenter, $StereoCenterStatus, $StereoChemistry, $StereochemistryString, $RingAtom, $NumOfRings, $AromaticAtom);
    4099
    4100   $ID = $This->GetID();
    4101   $Name = $This->GetName();
    4102   $AtomSymbol = $This->GetAtomSymbol();
    4103   $AtomicNumber = $This->GetAtomicNumber();
    4104   $XYZVector = $This->GetXYZVector();
    4105
    4106   $AtomicWeight = $This->GetAtomicWeight();
    4107   if (!defined $AtomicWeight) {
    4108     $AtomicWeight = 'undefined';
    4109   }
    4110   $ExactMass = $This->GetExactMass();
    4111   if (!defined $ExactMass) {
    4112     $ExactMass = 'undefined';
    4113   }
    4114   $NumOfNeighbors = $This->GetNumOfNeighbors();
    4115   if (!defined $NumOfNeighbors) {
    4116     $NumOfNeighbors = 'undefined';
    4117   }
    4118   $NumOfBonds = $This->GetNumOfBonds();
    4119   if (!defined $NumOfBonds) {
    4120     $NumOfBonds = 'undefined';
    4121   }
    4122   $Valence = $This->GetValence();
    4123   if (!defined $Valence) {
    4124     $Valence = 'undefined';
    4125   }
    4126
    4127   $MissingHydrogens = $This->GetNumOfMissingHydrogens();
    4128   if (!defined $MissingHydrogens) {
    4129     $MissingHydrogens = 'undefined';
    4130   }
    4131   $TotalHydrogens = $This->GetNumOfHydrogens();
    4132   if (!defined $TotalHydrogens) {
    4133     $TotalHydrogens = 'undefined';
    4134   }
    4135   $ImplicitHydrogens = $This->GetNumOfImplicitHydrogens();
    4136   if (!defined $ImplicitHydrogens) {
    4137     $ImplicitHydrogens = 'undefined';
    4138   }
    4139   $ExplicitHydrogens = $This->GetNumOfExplicitHydrogens();
    4140   if (!defined $ExplicitHydrogens) {
    4141     $ExplicitHydrogens = 'undefined';
    4142   }
    4143
    4144   $FormalCharge = $This->GetFormalCharge();
    4145   if (!defined $FormalCharge) {
    4146     $FormalCharge = 'undefined';
    4147   }
    4148   $Charge = $This->GetCharge();
    4149   if (!defined $Charge) {
    4150     $Charge = 'undefined';
    4151   }
    4152
    4153   $SpinMultiplicity = $This->GetSpinMultiplicity();
    4154   if (!defined $SpinMultiplicity) {
    4155     $SpinMultiplicity = 'undefined';
    4156   }
    4157   $FreeRadicalElectrons = $This->GetFreeRadicalElectrons();
    4158   if (!defined $FreeRadicalElectrons) {
    4159     $FreeRadicalElectrons = 'undefined';
    4160   }
    4161
    4162   $RingAtom = $This->IsInRing();
    4163   if (defined $RingAtom) {
    4164     $RingAtom = $RingAtom  ? 'Yes' : 'No';
    4165     $NumOfRings = $This->GetNumOfRings();
    4166   }
    4167   else {
    4168     $RingAtom = 'undefined';
    4169     $NumOfRings = 'undefined';
    4170   }
    4171
    4172   $AromaticAtom = $This->GetAromatic();
    4173   if (defined $AromaticAtom) {
    4174     $AromaticAtom = $AromaticAtom  ? 'Yes' : 'No';
    4175   }
    4176   else {
    4177     $AromaticAtom = 'undefined';
    4178   }
    4179
    4180   $StereochemistryString = '';
    4181   $StereoCenter = $This->GetStereoCenter();
    4182   if (defined $StereoCenter) {
    4183     $StereoCenterStatus = $This->IsStereoCenter() ? 'Yes' : 'No';
    4184     $StereoChemistry = $This->GetStereochemistry();
    4185     if (!defined $StereoChemistry) {
    4186       $StereoChemistry = 'undefined';
    4187     }
    4188     $StereochemistryString = "StereoCenter: $StereoCenterStatus; Stereochemistry: $StereoChemistry";
    4189   }
    4190
    4191   $AtomString = "Atom: ID: $ID; Name: \"$Name\"; AtomSymbol: \"$AtomSymbol\"; AtomicNumber: $AtomicNumber; XYZ: $XYZVector; AtomicWeight: $AtomicWeight; ExactMass: $ExactMass; NumOfNeighbors: $NumOfNeighbors;  NumOfBonds: $NumOfBonds; Valence: $Valence; MissingHydrogens: $MissingHydrogens; TotalHydrogens: $TotalHydrogens; ImplicitHydrogens: $ImplicitHydrogens; ExplicitHydrogens: $ExplicitHydrogens; FormalCharge: $FormalCharge; Charge: $Charge; SpinMultiplicity: $SpinMultiplicity; FreeRadicalElectrons: $FreeRadicalElectrons; RingAtom: $RingAtom; NumOfAtomRings: $NumOfRings; AromaticAtom: $AromaticAtom";
    4192
    4193   if ($StereochemistryString) {
    4194     $AtomString .= "; $StereochemistryString";
    4195   }
    4196
    4197   return $AtomString;
    4198 }
    4199
    4200 # Load appropriate atom data files from <MayaChemTools>/lib directory used by various
    4201 # object methods in the current class...
    4202 #
    4203 sub _LoadAtomClassData {
    4204   my($MayaChemToolsLibDir);
    4205
    4206   $MayaChemToolsLibDir = GetMayaChemToolsLibDirName();
    4207
    4208   _LoadAtomValenceModelData($MayaChemToolsLibDir);
    4209
    4210 }
    4211
    4212 #
    4213 # Load data for supported valence models...
    4214 #
    4215 sub _LoadAtomValenceModelData {
    4216   my($MayaChemToolsLibDir) = @_;
    4217   my($MDLValenceModelDataFile, $DaylightValenceModelDataFile);
    4218
    4219   %MDLValenceModelDataMap = ();
    4220   %DaylightValenceModelDataMap = ();
    4221
    4222   $MDLValenceModelDataFile = $MayaChemToolsLibDir . "/data/MDLValenceModelData.csv";
    4223   $DaylightValenceModelDataFile = $MayaChemToolsLibDir . "/data/DaylightValenceModelData.csv";
    4224
    4225   if (! -e "$MDLValenceModelDataFile") {
    4226     croak "Error: ${ClassName}::_LoadAtomValenceModelData: MayaChemTools package file, $MDLValenceModelDataFile, is missing: Possible installation problems...";
    4227   }
    4228
    4229   if (! -e "$DaylightValenceModelDataFile") {
    4230     croak "Error: ${ClassName}::_LoadAtomValenceModelData: MayaChemTools package file, $DaylightValenceModelDataFile, is missing: Possible installation problems...";
    4231   }
    4232
    4233   _LoadValenceModelDataFile($MDLValenceModelDataFile, \%MDLValenceModelDataMap);
    4234   _LoadValenceModelDataFile($DaylightValenceModelDataFile, \%DaylightValenceModelDataMap);
    4235
    4236 }
    4237
    4238 #
    4239 # Load valence model data file...
    4240 #
    4241 sub _LoadValenceModelDataFile {
    4242   my($DataFile, $DataMapRef) = @_;
    4243
    4244   # File format:
    4245   #
    4246   # "AtomicNumber","ElementSymbol","FormalCharge","CommomValences"
    4247   #
    4248   my($InDelim, $Line, $NumOfCols, @ColLabels, @LineWords);
    4249
    4250   $InDelim = "\,";
    4251   open DATAFILE, "$DataFile" or croak "Couldn't open $DataFile: $! ...";
    4252
    4253   # Skip lines up to column labels...
    4254   LINE: while ($Line = GetTextLine(\*DATAFILE)) {
    4255     if ($Line !~ /^#/) {
    4256       last LINE;
    4257     }
    4258   }
    4259   @ColLabels= quotewords($InDelim, 0, $Line);
    4260   $NumOfCols = @ColLabels;
    4261
    4262   my($AtomicNumber, $ElementSymbol, $FormalCharge, $CommonValences);
    4263
    4264   # Process element data...
    4265   LINE: while ($Line = GetTextLine(\*DATAFILE)) {
    4266     if ($Line =~ /^#/) {
    4267       next LINE;
    4268     }
    4269     @LineWords = ();
    4270     @LineWords = quotewords($InDelim, 0, $Line);
    4271     if (@LineWords != $NumOfCols) {
    4272       croak "Error: ${ClassName}::_LoadValenceModelDataFile: The number of data fields, @LineWords, in $DataFile must be $NumOfCols.\nLine: $Line...";
    4273     }
    4274
    4275     ($AtomicNumber, $ElementSymbol, $FormalCharge, $CommonValences) = @LineWords;
    4276
    4277     if (exists $DataMapRef->{$AtomicNumber}) {
    4278       # Additional data for an element...
    4279       if (exists $DataMapRef->{$AtomicNumber}{$FormalCharge}) {
    4280         # Duplicate data entries for an element...
    4281         carp "Warning: ${ClassName}::_LoadValenceModelDataFile: Ignoring valence data for element with atomic number $AtomicNumber and formal charge $FormalCharge in data file $DataFile: It has already been loaded.\nLine: $Line...";
    4282         next LINE;
    4283       }
    4284     }
    4285     else {
    4286       # Data for a new element...
    4287       %{$DataMapRef->{$AtomicNumber}} = ();
    4288     }
    4289
    4290     %{$DataMapRef->{$AtomicNumber}{$FormalCharge}} = ();
    4291     $DataMapRef->{$AtomicNumber}{$FormalCharge}{ElementSymbol} = $ElementSymbol;
    4292
    4293     @{$DataMapRef->{$AtomicNumber}{$FormalCharge}{CommonValences}} = ();
    4294     @{$DataMapRef->{$AtomicNumber}{$FormalCharge}{CommonValences}} = sort { $a <=> $b } split /\,/, $CommonValences;
    4295   }
    4296   close DATAFILE;
    4297 }
    4298