       1 package AtomicDescriptors::EStateValuesDescriptors;
       2 #
       3 # File: EStateValuesDescriptors.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 Scalar::Util ();
      30 use Matrix;
      31 use Constants;
      32 use TextUtil ();
      33 use MathUtil ();
      34 use StatisticsUtil ();
      35 use Atom;
      36 use Molecule;
      37 use AtomicDescriptors::AtomicDescriptors;
      38
      39 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
      40
      41 @ISA = qw(AtomicDescriptors::AtomicDescriptors Exporter);
      42 @EXPORT = qw();
      43 @EXPORT_OK = qw();
      44
      45 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
      46
      47 # Setup class variables...
      48 my($ClassName);
      49 _InitializeClass();
      50
      51 # Overload Perl functions...
      52 use overload '""' => 'StringifyEStateValuesDescriptors';
      53
      54 # Class constructor...
      55 sub new {
      56   my($Class, %NamesAndValues) = @_;
      57
      58   # Initialize object...
      59   my $This = $Class->SUPER::new();
      60   bless $This, ref($Class) || $Class;
      61   $This->_InitializeEStateValuesDescriptors();
      62
      63   $This->_InitializeEStateValuesDescriptorsProperties(%NamesAndValues);
      64
      65   return $This;
      66 }
      67
      68 # Initialize class ...
      69 sub _InitializeClass {
      70   #Class name...
      71   $ClassName = __PACKAGE__;
      72 }
      73
      74
      75 # Initialize object data...
      76 #
      77 sub _InitializeEStateValuesDescriptors {
      78   my($This) = @_;
      79
      80   # Type of AtomicDescriptor...
      81   $This->{Type} = 'EStateValue';
      82
      83   # Intrinsic state values calculated for non-hydrogen atom...
      84   #
      85   %{$This->{IStateValues}} = ();
      86
      87   # Calculatetion of E-state values for types for hydrogens is not supported...
      88   $This->{IgnoreHydrogens} = 1;
      89
      90   # Perturbation to intrinsic state values of atoms corresponding to all non-hydrogen
      91   # atom pairs...
      92   #
      93   %{$This->{DeltaIStateMatrix}} = ();
      94
      95   # E-state values calculated for non-hydrogen atoms...
      96   #
      97   %{$This->{EStateValues}} = ();
      98
      99   return $This;
     100 }
     101
     102 # Initialize object properties...
     103 #
     104 sub _InitializeEStateValuesDescriptorsProperties {
     105   my($This, %NamesAndValues) = @_;
     106
     107   my($Name, $Value, $MethodName);
     108   while (($Name, $Value) = each  %NamesAndValues) {
     109     $MethodName = "Set${Name}";
     110     $This->$MethodName($Value);
     111   }
     112
     113   # Make sure molecule object was specified...
     114   if (!exists $NamesAndValues{Molecule}) {
     115     croak "Error: ${ClassName}->New: Object can't be instantiated without specifying molecule...";
     116   }
     117
     118   # Intialize atomic descriptor values...
     119   $This->_InitializeDescriptorValues();
     120
     121   return $This;
     122 }
     123
     124 # Disable change of ignore hydrogens...
     125 #
     126 sub SetIgnoreHydrogens {
     127   my($This, $IgnoreHydrogens) = @_;
     128
     129   carp "Warning: ${ClassName}->SetIgnoreHydrogens: Ignore hydrogens value can't be changed: It's not supported...";
     130
     131   return $This;
     132 }
     133
     134 # Generate electrotopological state (E-state) values [ Ref 75-78 ] for all atoms
     135 # in the molecule...
     136 #
     137 # Calculation of E-state values for non-hydrogen atoms:
     138 #
     139 # Let:
     140 #
     141 #   N = Principal quantum number or period number corresponding to element symbol
     142 #
     143 #   Sigma = Number of sigma electrons involves in bonds to hydrogen and non-hydrogen atoms
     144 #           attached to atom
     145 #         = Number of sigma bonds to hydrogen and non-hydrogen atoms attached to atom
     146 #   PI = Number of PI electrons involved in bonds to non-hydrogen atoms attached to atom
     147 #       = Number of PI bonds to non-hydrogen atoms attached to atom
     148 #
     149 #   LP = Number of lone pair electrons on atom
     150 #
     151 #   Zv = Number of electrons in valence shell of atom
     152 #
     153 #   X = Number of non-hydrogen atom neighbors or heavy atoms attached to atom
     154 #   H = Number of implicit and explicit hydrogens for atom
     155 #
     156 #   Delta = Number of sigma electrons involved to bonds to non-hydrogen atoms
     157 #   DeltaV = ValenceDelta = Number of valence shell electrons not involved in bonding to hydrogen atoms
     158 #
     159 #   Ii = Intrinsic state value for atom i
     160 #
     161 #   DeltaIi = Sum of perturbations to intrinsic state value Ii of atom i by all other atoms besides atom i
     162 #
     163 #   DeltaIij = Perturbation to intrinsic state value Ii of atom i by atom j
     164 #
     165 #   Dij = Graph/bond distance between atom i and j
     166 #   Rij = Dij + 1
     167 #
     168 #   Si = E-state value for atom i
     169 #
     170 #
     171 # Then:
     172 #
     173 #   Delta = Sigma - H = X
     174 #
     175 #   DeltaV = Zv - H
     176 #      = Sigma + PI + LP - H
     177 #
     178 #   Ii = ( ( ( 2 / N ) ** 2 ) * DeltaV + 1 ) / Delta
     179 #
     180 #   Si = Ii + DeltaIi
     181 #
     182 #   DeltaIij = (Ii - Ij) / (Rij ** 2) for j not equal to i
     183 #
     184 #   DeltaIji = - DeltaIij
     185 #
     186 #   DeltaIi = SUM ( (Ii - Ij) / (Rij ** 2) ) for j = 1 to num of atoms skipping atom i
     187 #
     188 # Methodology:
     189 #   . Calculate intrinsic state values for atoms.
     190 #   . Generate a distance matrix.
     191 #   . Use distance matrix to calculate DeltaIij matrix with each row i containing
     192 #     DeltaIij values corresponding to perturbation to intrinsic state value of atom
     193 #     i by atom j.
     194 #   . Calculate E-state values using DeltaIij matrix.
     195 #   . Assign E-state values to atoms.
     196 #
     197 # Notes:
     198 #   . The current release of MayaChemTools doesn't support calculation of Hydrogen
     199 #     E-state values.
     200 #
     201 sub GenerateDescriptors {
     202   my($This) = @_;
     203
     204   # Cache appropriate molecule data...
     205   $This->_SetupMoleculeDataCache();
     206
     207   # Generate distance matrix...
     208   if (!$This->_SetupDistanceMatrix()) {
     209     carp "Warning: ${ClassName}->GenerateDescriptors: E-state values description generation didn't succeed: Couldn't generate a distance matrix...";
     210     return $This;
     211   }
     212
     213   # Calculate EState values..
     214   if (!$This->_CalculateEStateValuesDescriptors()) {
     215     carp "Warning: ${ClassName}->GenerateDescriptors: E-state values description generation didn't succeed: Couldn't calculate IState values for all atoms...";
     216     return $This;
     217   }
     218
     219   # Set final descriptor values...
     220   $This->_SetFinalValues();
     221
     222   # Clear cached molecule data...
     223   $This->_ClearMoleculeDataCache();
     224
     225   return $This;
     226 }
     227
     228 # Calculate E-state values for non-hydrogen atoms...
     229 #
     230 sub _CalculateEStateValuesDescriptors {
     231   my($This) = @_;
     232   my($DeltaIStateMatrix, $NumOfRows, $NumOfCols, $RowIndex, $AtomID, $EStateValue, $IStateValue, $DeltaIStateValue, @DeltaIStateMatrixRowValues);
     233
     234   %{$This->{EStateValues}} = ();
     235
     236   # Calculate intrinsic state values for non-hydrogen atoms...
     237   if (!$This->_CalculateIStateValues()) {
     238     return undef;
     239   }
     240
     241   # Calculate delta intrinsic state matrix for non-hydrogen atoms...
     242   $This->_CalculateDeltaIStateMatrix();
     243
     244   # Get DeltaIState matrix information...
     245   $DeltaIStateMatrix = $This->{DeltaIStateMatrix};
     246   ($NumOfRows, $NumOfCols) = $DeltaIStateMatrix->GetSize();
     247
     248   # Calculate EState values...
     249   ROWINDEX: for $RowIndex (0 .. ($NumOfRows - 1) ) {
     250     $AtomID = $This->{AtomIndexToID}{$RowIndex};
     251     if ( !(exists($This->{IStateValues}{$AtomID})) ) {
     252       next ROWINDEX;
     253     }
     254     $IStateValue = $This->{IStateValues}{$AtomID};
     255
     256     @DeltaIStateMatrixRowValues = $DeltaIStateMatrix->GetRowValues($RowIndex);
     257     $DeltaIStateValue = StatisticsUtil::Sum(\@DeltaIStateMatrixRowValues);
     258
     259     $EStateValue = $IStateValue + $DeltaIStateValue;
     260
     261     $This->{EStateValues}{$AtomID} = $EStateValue;
     262   }
     263   return $This;
     264 }
     265
     266 # Calculate intrinsic state values for non-hydrogen atoms...
     267 #
     268 sub _CalculateIStateValues {
     269   my($This) = @_;
     270   my($Atom, $AtomID);
     271
     272   %{$This->{IStateValues}} = ();
     273
     274   ATOM: for $Atom (@{$This->{Atoms}}) {
     275     # Irrespective of IgoreHydrogens value, just ignore hydrogens...
     276     if ($Atom->IsHydrogen()) {
     277       next ATOM;
     278     }
     279     $AtomID = $Atom->GetID();
     280     $This->{IStateValues}{$AtomID} = $This->_CalculateIStateValue($Atom);
     281     if (!defined($This->{IStateValues}{$AtomID})) {
     282       return undef;
     283     }
     284   }
     285   return $This;
     286 }
     287
     288 # Calculation intrinsic state value for non-hydrogen...
     289 #
     290 sub _CalculateIStateValue {
     291   my($This, $Atom) = @_;
     292   my($IStateValue, $Delta, $DeltaV, $PeriodNumber);
     293
     294   $PeriodNumber = $Atom->GetPeriodNumber();
     295   ($Delta, $DeltaV) = $This->_GetDeltaValues($Atom);
     296
     297   if (!(defined($Delta) && defined($PeriodNumber) && defined($DeltaV))) {
     298     return undef;
     299   }
     300
     301   $IStateValue = ($PeriodNumber && $Delta) ? (((2/$PeriodNumber)**2)*$DeltaV + 1)/$Delta : 0;
     302
     303   return $IStateValue;
     304 }
     305
     306 # Get Delta and DeltaV values for atom...
     307 #
     308 sub _GetDeltaValues {
     309   my($This, $Atom) = @_;
     310   my($Delta, $DeltaV, $ValenceElectrons, $NumOfHydrogens);
     311
     312   ($Delta, $DeltaV) = (undef, undef);
     313
     314   $ValenceElectrons = $Atom->GetValenceElectrons();
     315   $NumOfHydrogens = $Atom->GetAtomicInvariantValue('H');
     316
     317   $Delta = $Atom->GetAtomicInvariantValue('X');
     318   if (defined($ValenceElectrons) && defined($NumOfHydrogens)) {
     319     $DeltaV = $ValenceElectrons - $NumOfHydrogens;
     320   }
     321
     322   return ($Delta, $DeltaV);
     323 }
     324
     325 # Calculate DeltaIState matrix for atoms with each row i containing DeltaIij values
     326 # corresponding atom atoms i and j.
     327 #
     328 # Notes:
     329 #   . Matrix elements corresponding to hydrogen atoms or unconnected
     330 #     are assigned zero value.
     331 #
     332 sub _CalculateDeltaIStateMatrix {
     333   my($This) = @_;
     334   my($DistanceMatrix, $NumOfRows, $NumOfCols, $RowIndex, $ColIndex, $AtomID1, $AtomID2, $DeltaIStateMatrix, $IStateValue1, $IStateValue2, $GraphDistance, $DeltaIState12, $DeltaIState21, $SkipIndexCheck);
     335
     336   # Get distance matrix information...
     337   $DistanceMatrix = $This->{DistanceMatrix};
     338   ($NumOfRows, $NumOfCols) = $DistanceMatrix->GetSize();
     339
     340   # Initialize DeltaIState matrix...
     341   $This->{DeltaIStateMatrix} = new Matrix($NumOfRows, $NumOfCols);
     342   $DeltaIStateMatrix = $This->{DeltaIStateMatrix};
     343
     344   $SkipIndexCheck = 1;
     345
     346   # Calculate DeltaIState matrix values...
     347   ROWINDEX: for $RowIndex (0 .. ($NumOfRows - 1) ) {
     348     $AtomID1 = $This->{AtomIndexToID}{$RowIndex};
     349     if (!exists($This->{IStateValues}{$AtomID1})) {
     350       next ROWINDEX;
     351     }
     352     $IStateValue1 = $This->{IStateValues}{$AtomID1};
     353
     354     COLINDEX: for $ColIndex (($RowIndex + 1) .. ($NumOfCols - 1) ) {
     355       $AtomID2 = $This->{AtomIndexToID}{$ColIndex};
     356       if (!exists($This->{IStateValues}{$AtomID2})) {
     357         next COLINDEX;
     358       }
     359       $IStateValue2 = $This->{IStateValues}{$AtomID2};
     360
     361       # Make sure it's a connected atom...
     362       $GraphDistance = $DistanceMatrix->GetValue($RowIndex, $ColIndex, $SkipIndexCheck);
     363       if ($GraphDistance >= BigNumber) {
     364         next COLINDEX;
     365       }
     366
     367       $DeltaIState12 = ($IStateValue1 - $IStateValue2)/(($GraphDistance + 1)**2);
     368       $DeltaIState21 = -$DeltaIState12;
     369
     370       # Set DeltaIState values...
     371       $DeltaIStateMatrix->SetValue($RowIndex, $ColIndex, $DeltaIState12, $SkipIndexCheck);
     372       $DeltaIStateMatrix->SetValue($ColIndex, $RowIndex, $DeltaIState21, $SkipIndexCheck);
     373     }
     374   }
     375 }
     376
     377 # Setup distance matrix...
     378 #
     379 sub _SetupDistanceMatrix {
     380   my($This) = @_;
     381
     382   $This->{DistanceMatrix} = $This->GetMolecule()->GetDistanceMatrix();
     383
     384   if (!defined($This->{DistanceMatrix})) {
     385     return undef;
     386   }
     387
     388   return $This;
     389 }
     390
     391 # Setup final descriptor values...
     392 #
     393 sub _SetFinalValues {
     394   my($This) = @_;
     395   my($Atom, $AtomID, $EStateValue);
     396
     397   ATOM: for $Atom (@{$This->{Atoms}}) {
     398     $AtomID = $Atom->GetID();
     399     if (!exists $This->{EStateValues}{$AtomID}) {
     400       next ATOM;
     401     }
     402     $EStateValue = $This->{EStateValues}{$AtomID};
     403     $This->SetDescriptorValue($Atom, $EStateValue);
     404   }
     405
     406   return $This;
     407 }
     408
     409 # Cache  appropriate molecule data...
     410 #
     411 sub _SetupMoleculeDataCache {
     412   my($This) = @_;
     413
     414   # Get all atoms including hydrogens to correctly map atom indicies to atom IDs for
     415   # usage of distance matrix.
     416   #
     417   @{$This->{Atoms}} = $This->GetMolecule()->GetAtoms();
     418
     419   # Get all atom IDs...
     420   my(@AtomIDs);
     421   @AtomIDs = ();
     422   @AtomIDs =  map { $_->GetID() } @{$This->{Atoms}};
     423
     424   # Set AtomIndex to AtomID hash...
     425   %{$This->{AtomIndexToID}} = ();
     426   @{$This->{AtomIndexToID}}{ (0 .. $#AtomIDs) } = @AtomIDs;
     427
     428   return $This;
     429 }
     430
     431 # Clear cached molecule data...
     432 #
     433 sub _ClearMoleculeDataCache {
     434   my($This) = @_;
     435
     436   @{$This->{Atoms}} = ();
     437
     438   return $This;
     439 }
     440
     441 # Return a string containg data for EStateValuesDescriptors object...
     442 #
     443 sub StringifyEStateValuesDescriptors {
     444   my($This) = @_;
     445   my($EStateValuesDescriptorsString);
     446
     447   # Type of AtomicValues...
     448   $EStateValuesDescriptorsString = "AtomicDescriptorType: $This->{Type}; IgnoreHydrogens: " . ($This->{IgnoreHydrogens} ? "Yes" : "No");
     449
     450   # Setup atomic descriptor information...
     451   my($AtomID, $DescriptorValue, @DescriptorValuesInfo, %DescriptorValues);
     452
     453   @DescriptorValuesInfo = ();
     454   %DescriptorValues = $This->GetDescriptorValues();
     455
     456   for $AtomID (sort { $a <=> $b } keys %DescriptorValues) {
     457     $DescriptorValue = $DescriptorValues{$AtomID};
     458     $DescriptorValue = (TextUtil::IsEmpty($DescriptorValue) || $DescriptorValue =~ /^None$/i) ? 'None' : MathUtil::round($DescriptorValue, 3) + 0;
     459     push @DescriptorValuesInfo, "$AtomID:$DescriptorValue";
     460   }
     461   $EStateValuesDescriptorsString .= "; AtomIDs:EStateValuesDescriptors: <" . TextUtil::JoinWords(\@DescriptorValuesInfo, ", ", 0) . ">";
     462
     463   return $EStateValuesDescriptorsString;
     464 }
     465
     466 # Is it a EStateValuesDescriptors object?
     467 sub _IsEStateValuesDescriptors {
     468   my($Object) = @_;
     469
     470   return (Scalar::Util::blessed($Object) && $Object->isa($ClassName)) ? 1 : 0;
     471 }
     472