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