1 package Atom;
   2 #
   3 # $RCSfile: Atom.pm,v $
   4 # $Date: 2008/04/25 00:00:44 $
   5 # $Revision: 1.28 $
   6 #
   7 # Author: Manish Sud <msud@san.rr.com>
   8 #
   9 # Copyright (C) 2004-2008 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 use 5.006;
  29 use strict;
  30 use Carp;
  31 use Exporter;
  32 use Storable ();
  33 use ObjectProperty;
  34 use PeriodicTable;
  35 use Vector;
  36 
  37 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  38 
  39 $VERSION = '1.00';
  40 @ISA = qw(ObjectProperty Exporter);
  41 @EXPORT = qw();
  42 @EXPORT_OK = qw();
  43 
  44 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
  45 
  46 # Setup class variables...
  47 my($ClassName, $ObjectID);
  48 _InitializeClass();
  49 
  50 # Overload Perl functions...
  51 use overload '""' => 'StringifyAtom';
  52 
  53 # Class constructor...
  54 sub new {
  55   my($Class, %NamesAndValues) = @_;
  56 
  57   # Initialize object...
  58   my $This = {};
  59   bless $This, ref($Class) || $Class;
  60   $This->_InitializeAtom();
  61 
  62   $This->_InitializeAtomProperties(%NamesAndValues);
  63 
  64   return $This;
  65 }
  66 
  67 # Initialize object data...
  68 sub _InitializeAtom {
  69   my($This) = @_;
  70   my($ObjectID) = _GetNewObjectID();
  71 
  72   # All other property names and values along with all Set/Get<PropertyName> methods
  73   # are implemented on-demand using ObjectProperty class.
  74   $This->{ID} = $ObjectID;
  75   $This->{Name} = "Atom ${ObjectID}";
  76   $This->{AtomSymbol} = '';
  77   $This->{AtomicNumber} = 0;
  78   $This->{XYZ} = Vector::ZeroVector;
  79 }
  80 
  81 # Initialize atom properties...
  82 sub _InitializeAtomProperties {
  83   my($This, %NamesAndValues) = @_;
  84 
  85   my($Name, $Value, $MethodName);
  86   while (($Name, $Value) = each  %NamesAndValues) {
  87     $MethodName = "Set${Name}";
  88     $This->$MethodName($Value);
  89   }
  90   if (!exists $NamesAndValues{'AtomSymbol'}) {
  91     carp "Warning: ${ClassName}->new: Atom object instantiated without setting atom symbol...";
  92   }
  93 
  94   return $This;
  95 }
  96 
  97 # Initialize class ...
  98 sub _InitializeClass {
  99   #Class name...
 100   $ClassName = __PACKAGE__;
 101 
 102   # ID to keep track of objects...
 103   $ObjectID = 0;
 104 }
 105 
 106 # Setup an explicit SetID method to block setting of ID by AUTOLOAD function...
 107 sub SetID {
 108   my($This, $Value) = @_;
 109 
 110   carp "Warning: ${ClassName}->SetID: Object ID can't be changed: it's used for internal tracking...";
 111 
 112   return $This;
 113 }
 114 
 115 # Setup atom symbol and atomic number for the element...
 116 #
 117 # Possible atom symbol values:
 118 #    . An element symbol or some other type of atom: L - Atom list; LP - Lone pair; R# - R group;
 119 #       A, Q, * - unknown atom; or something else?
 120 #
 121 # Default mass number corresponds to the most abundant natural isotope unless it's explicity
 122 # set using "MassNumber" property.
 123 #
 124 sub SetAtomSymbol {
 125   my($This, $AtomSymbol) = @_;
 126   my($AtomicNumber);
 127 
 128   $This->{AtomSymbol} = $AtomSymbol;
 129 
 130   $AtomicNumber = PeriodicTable::GetElementAtomicNumber($AtomSymbol);
 131   $This->{AtomicNumber} = (defined $AtomicNumber) ? $AtomicNumber : 0;
 132 
 133   return $This;
 134 }
 135 
 136 # Setup atom symbol and atomic number for the element...
 137 sub SetAtomicNumber {
 138   my($This, $AtomicNumber) = @_;
 139   my($AtomSymbol);
 140 
 141   $AtomSymbol = PeriodicTable::GetElementAtomSymbol($AtomicNumber);
 142   if (!defined $AtomSymbol) {
 143     carp "Warning: ${ClassName}->SetAtomicNumber: Didn't set atomic number: Invalid atomic number, $AtomicNumber, specified...";
 144     return;
 145   }
 146   $This->{AtomicNumber} = $AtomicNumber;
 147   $This->{AtomSymbol} = $AtomSymbol;
 148 
 149   return $This;
 150 }
 151 
 152 # Set atom as stereo center...
 153 #
 154 sub SetStereoCenter {
 155   my($This, $StereoCenter) = @_;
 156 
 157   $This->SetProperty('StereoCenter', $StereoCenter);
 158 
 159   return $This;
 160 }
 161 
 162 # Is it a stereo center?
 163 #
 164 sub IsStereoCenter {
 165   my($This) = @_;
 166   my($StereoCenter);
 167 
 168   $StereoCenter = $This->GetProperty('StereoCenter');
 169 
 170   return (defined($StereoCenter) && $StereoCenter) ? 1 : 0;
 171 }
 172 
 173 # Set atom stereochemistry.
 174 #
 175 # Supported values are: R, S.
 176 #
 177 # Notes:
 178 #
 179 # . After the ligands around a central stereocenter has been ranked using CIP priority scheme and
 180 # the lowest ranked ligand lies behind the center atom, then R and S values correspond to:
 181 #
 182 # R: Clockwise arrangement of remaining ligands around the central atom going from highest to lowest ranked ligand
 183 # S: CounterClockwise arrangement of remaining ligands around the central atom going from highest to lowest ranked ligand
 184 #
 185 # . Assignment of any other arbitray values besides R and S is also allowed; however, a warning is printed.
 186 #
 187 sub SetStereochemistry {
 188   my($This, $Stereochemistry) = @_;
 189 
 190   if ($Stereochemistry !~ /^(R|S)$/i) {
 191     carp "Warning: ${ClassName}->SetStereochemistry: Assigning non-supported Stereochemistry value of $Stereochemistry. Supported values: R, S...";
 192   }
 193 
 194   $This->SetProperty('StereoCenter', 1);
 195   $This->SetProperty('Stereochemistry', $Stereochemistry);
 196 
 197   return $This;
 198 }
 199 
 200 # Setup mass number for atom...
 201 sub SetMassNumber {
 202   my($This, $MassNumber) = @_;
 203   my($AtomicNumber, $AtomSymbol);
 204 
 205   $AtomicNumber = $This->{AtomicNumber};
 206   $AtomSymbol = $This->{AtomSymbol};
 207   if (!$AtomicNumber) {
 208     carp "Warning: ${ClassName}->SetMassNumber: Didn't set mass number: Non standard atom with atomic number, $AtomicNumber, and atomic symbol, $AtomSymbol...";
 209     return;
 210   }
 211   if (!PeriodicTable::IsElementNaturalIsotopeMassNumber($AtomicNumber, $MassNumber)) {
 212     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...";
 213   }
 214   $This->SetProperty('MassNumber', $MassNumber);
 215 
 216   return $This;
 217 }
 218 
 219 # Get mass number...
 220 #
 221 sub GetMassNumber {
 222   my($This) = @_;
 223 
 224   # Is mass number explicity set?
 225   if ($This->HasProperty('MassNumber')) {
 226     return $This->GetProperty('MassNumber');
 227   }
 228 
 229   # Is it an element symbol?
 230   my($AtomicNumber) = $This->{AtomicNumber};
 231   if (!$AtomicNumber) {
 232     return undef;
 233   }
 234 
 235   # Return most abundant mass number...
 236   return PeriodicTable::GetElementMostAbundantNaturalIsotopeMassNumber($AtomicNumber);
 237 }
 238 
 239 # Get atomic weight:
 240 #   . Explicitly set by the caller
 241 #   . Using atomic number
 242 #
 243 sub GetAtomicWeight {
 244   my($This) = @_;
 245 
 246   # Is atomic weight explicity set?
 247   if ($This->HasProperty('AtomicWeight')) {
 248     return $This->GetProperty('AtomicWeight');
 249   }
 250 
 251   # Is it an element symbol?
 252   my($AtomicNumber) = $This->{AtomicNumber};
 253   if (!$AtomicNumber) {
 254     return undef;
 255   }
 256 
 257   # Return its atomic weight...
 258   return PeriodicTable::GetElementAtomicWeight($AtomicNumber);
 259 }
 260 
 261 # Get exact mass weight:
 262 #   . Explicitly set by the caller
 263 #   . Using atomic number and mass number explicity set by the caller
 264 #   . Using atomic number and most abundant isotope
 265 #
 266 sub GetExactMass {
 267   my($This) = @_;
 268 
 269   # Is exact mass explicity set?
 270   if ($This->HasProperty('ExactMass')) {
 271     return $This->GetProperty('ExactMass');
 272   }
 273 
 274   # Is it an element symbol?
 275   my($AtomicNumber) = $This->{AtomicNumber};
 276   if (!$AtomicNumber) {
 277     return undef;
 278   }
 279 
 280   # Is mass number explicitly set?
 281   if ($This->HasProperty('MassNumber')) {
 282     my($MassNumber) = $This->GetProperty('MassNumber');
 283     if (PeriodicTable::IsElementNaturalIsotopeMassNumber($AtomicNumber, $MassNumber)) {
 284       return PeriodicTable::GetElementNaturalIsotopeMass($AtomicNumber, $MassNumber);
 285     }
 286   }
 287 
 288   # Return most abundant isotope mass...
 289   return PeriodicTable::GetElementMostAbundantNaturalIsotopeMass($AtomicNumber);
 290 }
 291 
 292 
 293 # Set atom coordinates using:
 294 # . An array reference with three values
 295 # . An array containg three values
 296 # . A 3D vector
 297 #
 298 sub SetXYZ {
 299   my($This, @Values) = @_;
 300 
 301   if (!@Values) {
 302     carp "Warning: ${ClassName}->SetXYZ: No values specified...";
 303     return;
 304   }
 305 
 306   $This->{XYZ}->SetXYZ(@Values);
 307   return $This;
 308 }
 309 
 310 # Set X value...
 311 sub SetX {
 312   my($This, $Value) = @_;
 313 
 314   if (!defined $Value) {
 315     carp "Warning: ${ClassName}->SetX: Undefined X value...";
 316     return;
 317   }
 318   $This->{XYZ}->SetX($Value);
 319   return $This;
 320 }
 321 
 322 # Set Y value...
 323 sub SetY {
 324   my($This, $Value) = @_;
 325 
 326   if (!defined $Value) {
 327     carp "Warning: ${ClassName}->SetY: Undefined Y value...";
 328     return;
 329   }
 330   $This->{XYZ}->SetY($Value);
 331   return $This;
 332 }
 333 
 334 # Set Z value...
 335 sub SetZ {
 336   my($This, $Value) = @_;
 337 
 338   if (!defined $Value) {
 339     carp "Warning: ${ClassName}->SetZ: Undefined Z value...";
 340     return;
 341   }
 342   $This->{XYZ}->SetZ($Value);
 343   return $This;
 344 }
 345 
 346 # Return XYZ as:
 347 # . Reference to an array
 348 # . An array
 349 #
 350 sub GetXYZ {
 351   my($This) = @_;
 352 
 353   return $This->{XYZ}->GetXYZ();
 354 }
 355 
 356 # Return XYZ as a vector object...
 357 #
 358 sub GetXYZVector {
 359   my($This) = @_;
 360 
 361   return $This->{XYZ};
 362 }
 363 
 364 # Get X value...
 365 sub GetX {
 366   my($This) = @_;
 367 
 368   return $This->{XYZ}->GetX();
 369 }
 370 
 371 # Get Y value...
 372 sub GetY {
 373   my($This) = @_;
 374 
 375   return $This->{XYZ}->GetY();
 376 }
 377 
 378 # Get Z value...
 379 sub GetZ {
 380   my($This) = @_;
 381 
 382   return $This->{XYZ}->GetZ();
 383 }
 384 
 385 # Delete atom...
 386 sub DeleteAtom {
 387   my($This) = @_;
 388 
 389   # Is this atom in a molecule?
 390   if (!$This->HasProperty('Molecule')) {
 391     # Nothing to do...
 392     return $This;
 393   }
 394   my($Molecule) = $This->GetProperty('Molecule');
 395 
 396   return $Molecule->_DeleteAtom($This);
 397 }
 398 
 399 # Get atom neighbor objects as array. In scalar conetxt, return number of neighbors...
 400 sub GetNeighbors {
 401   my($This) = @_;
 402 
 403   # Is this atom in a molecule?
 404   if (!$This->HasProperty('Molecule')) {
 405     return undef;
 406   }
 407   my($Molecule) = $This->GetProperty('Molecule');
 408 
 409   return $Molecule->_GetAtomNeighbors($This);
 410 }
 411 
 412 # Get non-hydrogen atom neighbor objects as array. In scalar context, return number of neighbors...
 413 sub GetHeavyAtomNeighbors {
 414   my($This) = @_;
 415 
 416   return $This->GetNonHydrogenAtomNeighbors();
 417 }
 418 
 419 # Get non-hydrogen atom neighbor objects as array. In scalar context, return number of neighbors...
 420 sub GetNonHydrogenAtomNeighbors {
 421   my($This) = @_;
 422 
 423   # Is this atom in a molecule?
 424   if (!$This->HasProperty('Molecule')) {
 425     return undef;
 426   }
 427   my($NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = (1, 0);
 428 
 429   return $This->_GetFilteredAtomNeighbors($NonHydrogenAtomsOnly, $HydrogenAtomsOnly);
 430 }
 431 
 432 # Get hydrogen atom neighbor objects as array. In scalar context, return numbe of neighbors...
 433 sub GetHydrogenAtomNeighbors {
 434   my($This) = @_;
 435 
 436   # Is this atom in a molecule?
 437   if (!$This->HasProperty('Molecule')) {
 438     return undef;
 439   }
 440   my($NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = (0, 1);
 441 
 442   return $This->_GetFilteredAtomNeighbors($NonHydrogenAtomsOnly, $HydrogenAtomsOnly);
 443 }
 444 
 445 # Get filtered atom atom neighbors
 446 sub _GetFilteredAtomNeighbors {
 447   my($This, $NonHydrogenAtomsOnly, $HydrogenAtomsOnly) = @_;
 448 
 449   # Check flags...
 450   if (!defined $NonHydrogenAtomsOnly) {
 451     $NonHydrogenAtomsOnly = 0;
 452   }
 453   if (!defined $HydrogenAtomsOnly) {
 454     $HydrogenAtomsOnly = 0;
 455   }
 456   my($Neighbor, @FilteredAtomNeighbors);
 457 
 458   @FilteredAtomNeighbors = ();
 459   NEIGHBOR: for $Neighbor ($This->GetNeighbors()) {
 460     if ($NonHydrogenAtomsOnly && $Neighbor->IsHydrogen()) {
 461       next NEIGHBOR;
 462     }
 463     if ($HydrogenAtomsOnly && (!$Neighbor->IsHydrogen())) {
 464       next NEIGHBOR;
 465     }
 466     push @FilteredAtomNeighbors, $Neighbor;
 467   }
 468 
 469   return wantarray ? @FilteredAtomNeighbors : scalar @FilteredAtomNeighbors;
 470 }
 471 
 472 # Get number of neighbors...
 473 #
 474 sub GetNumOfNeighbors {
 475   my($This) = @_;
 476   my($NumOfNeighbors);
 477 
 478   $NumOfNeighbors = $This->GetNeighbors();
 479 
 480   return (defined $NumOfNeighbors) ? ($NumOfNeighbors) : undef;
 481 }
 482 
 483 # Get number of neighbors which are non-hydrogen atoms...
 484 sub GetNumOfHeavyAtomNeighbors {
 485   my($This) = @_;
 486 
 487   return $This->GetNumOfNonHydrogenAtomNeighbors();
 488 }
 489 
 490 # Get number of neighbors which are non-hydrogen atoms...
 491 sub GetNumOfNonHydrogenAtomNeighbors {
 492   my($This) = @_;
 493   my($NumOfNeighbors);
 494 
 495   $NumOfNeighbors = $This->GetNonHydrogenAtomNeighbors();
 496 
 497   return (defined $NumOfNeighbors) ? ($NumOfNeighbors) : undef;
 498 }
 499 
 500 # Get number of neighbors which are hydrogen atoms...
 501 sub GetNumOfHydrogenAtomNeighbors {
 502   my($This) = @_;
 503   my($NumOfNeighbors);
 504 
 505   $NumOfNeighbors = $This->GetHydrogenAtomNeighbors();
 506 
 507   return (defined $NumOfNeighbors) ? ($NumOfNeighbors) : undef;
 508 }
 509 
 510 # Get bond objects as array. In scalar context, return number of bonds...
 511 sub GetBonds {
 512   my($This) = @_;
 513 
 514   # Is this atom in a molecule?
 515   if (!$This->HasProperty('Molecule')) {
 516     return undef;
 517   }
 518   my($Molecule) = $This->GetProperty('Molecule');
 519 
 520   return $Molecule->_GetAtomBonds($This);
 521 }
 522 
 523 # Get bond to specified atom...
 524 sub GetBondToAtom {
 525   my($This, $Other) = @_;
 526 
 527   # Is this atom in a molecule?
 528   if (!$This->HasProperty('Molecule')) {
 529     return undef;
 530   }
 531   my($Molecule) = $This->GetProperty('Molecule');
 532 
 533   return $Molecule->_GetBondToAtom($This, $Other);
 534 }
 535 
 536 # Get bond objects to non-hydrogen atoms as array. In scalar context, return number of bonds...
 537 sub GetBondsToHeavyAtoms {
 538   my($This) = @_;
 539 
 540   return $This->GetBondsToNonHydrogenAtoms();
 541 }
 542 
 543 # Get bond objects to non-hydrogen atoms as array. In scalar context, return number of bonds...
 544 sub GetBondsToNonHydrogenAtoms {
 545   my($This) = @_;
 546 
 547   # Is this atom in a molecule?
 548   if (!$This->HasProperty('Molecule')) {
 549     return undef;
 550   }
 551   my($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = (1, 0);
 552 
 553   return $This->_GetFilteredBonds($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly);
 554 }
 555 
 556 # Get bond objects to hydrogen atoms as array. In scalar context, return number of bonds...
 557 sub GetBondsToHydrogenAtoms {
 558   my($This) = @_;
 559 
 560   # Is this atom in a molecule?
 561   if (!$This->HasProperty('Molecule')) {
 562     return undef;
 563   }
 564   my($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = (0, 1);
 565 
 566   return $This->_GetFilteredBonds($BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly);
 567 }
 568 
 569 # Get filtered bonds...
 570 sub _GetFilteredBonds {
 571   my($This, $BondsToNonHydrogenAtomsOnly, $BondsToHydrogenAtomsOnly) = @_;
 572 
 573   # Check flags...
 574   if (!defined $BondsToNonHydrogenAtomsOnly) {
 575     $BondsToNonHydrogenAtomsOnly = 0;
 576   }
 577   if (!defined $BondsToHydrogenAtomsOnly) {
 578     $BondsToHydrogenAtomsOnly = 0;
 579   }
 580 
 581   my($Bond, $BondedAtom, @FilteredBonds);
 582 
 583   @FilteredBonds = ();
 584   BOND: for $Bond ($This->GetBonds()) {
 585     $BondedAtom = $Bond->GetBondedAtom($This);
 586     if ($BondsToNonHydrogenAtomsOnly && $BondedAtom->IsHydrogen()) {
 587       next BOND;
 588     }
 589     if ($BondsToHydrogenAtomsOnly && (!$BondedAtom->IsHydrogen())) {
 590       next BOND;
 591     }
 592     push @FilteredBonds, $Bond;
 593   }
 594 
 595   return wantarray ? @FilteredBonds : (scalar @FilteredBonds);
 596 }
 597 
 598 # Get number of bonds...
 599 #
 600 sub GetNumOfBonds {
 601   my($This) = @_;
 602   my($NumOfBonds);
 603 
 604   $NumOfBonds = $This->GetBonds();
 605 
 606   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
 607 }
 608 
 609 # Get number of bonds to non-hydrogen atoms...
 610 sub GetNumOfBondsToHeavyAtoms {
 611   my($This) = @_;
 612 
 613   return $This->GetNumOfBondsToNonHydrogenAtoms();
 614 }
 615 
 616 # Get number of bonds to non-hydrogen atoms...
 617 sub GetNumOfBondsToNonHydrogenAtoms {
 618   my($This) = @_;
 619   my($NumOfBonds);
 620 
 621   $NumOfBonds = $This->GetBondsToNonHydrogenAtoms();
 622 
 623   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
 624 }
 625 
 626 # Get number of bonds to hydrogen atoms...
 627 sub GetNumOfBondsToHydrogenAtoms {
 628   my($This) = @_;
 629   my($NumOfBonds);
 630 
 631   $NumOfBonds = $This->GetBondsToHydrogenAtoms();
 632 
 633   return (defined $NumOfBonds) ? ($NumOfBonds) : undef;
 634 }
 635 
 636 # Get sum of bond orders to all bonded atoms...
 637 #
 638 sub GetSumOfBondOrders {
 639   my($This) = @_;
 640 
 641   # Is this atom in a molecule?
 642   if (!$This->HasProperty('Molecule')) {
 643     return undef;
 644   }
 645 
 646   return $This->_GetSumOfBondOrders();
 647 }
 648 
 649 # Get sum of bond orders to non-hydrogen atoms only...
 650 #
 651 sub GetSumOfBondOrdersToHeavyAtoms {
 652   my($This) = @_;
 653 
 654   return $This->GetSumOfBondOrdersToNonHydrogenAtoms();
 655 }
 656 
 657 # Get sum of bond orders to non-hydrogen atoms only...
 658 #
 659 sub GetSumOfBondOrdersToNonHydrogenAtoms {
 660   my($This) = @_;
 661 
 662   # Is this atom in a molecule?
 663   if (!$This->HasProperty('Molecule')) {
 664     return undef;
 665   }
 666   my($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = (1, 0);
 667 
 668   return $This->_GetSumOfBondOrders($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly);
 669 }
 670 
 671 # Get sum of bond orders to hydrogen atoms only...
 672 #
 673 sub GetSumOfBondOrdersToHydrogenAtoms {
 674   my($This) = @_;
 675 
 676   # Is this atom in a molecule?
 677   if (!$This->HasProperty('Molecule')) {
 678     return undef;
 679   }
 680   my($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = (0, 1);
 681 
 682   return $This->_GetSumOfBondOrders($ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly);
 683 }
 684 
 685 # Get sum of bond orders to all bonded atoms,  non-hydrogen or hydrogen bonded atoms...
 686 #
 687 sub _GetSumOfBondOrders {
 688   my($This, $ToNonHydrogenAtomsOnly, $ToHydrogenAtomsOnly) = @_;
 689 
 690   # Check flags...
 691   if (!defined $ToNonHydrogenAtomsOnly) {
 692     $ToNonHydrogenAtomsOnly = 0;
 693   }
 694   if (!defined $ToHydrogenAtomsOnly) {
 695     $ToHydrogenAtomsOnly = 0;
 696   }
 697   my($Bond, $SumOfBondOrders, $BondOrder, $NumOfAromaticBonds, @Bonds);
 698   @Bonds = ();
 699 
 700   if ($ToNonHydrogenAtomsOnly) {
 701     @Bonds = $This->GetBondsToNonHydrogenAtoms();
 702   }
 703   elsif ($ToHydrogenAtomsOnly) {
 704     @Bonds = $This->GetBondsToHydrogenAtoms();
 705   }
 706   else {
 707     # All bonds...
 708     @Bonds = $This->GetBonds();
 709   }
 710 
 711   $SumOfBondOrders = 0;
 712   $NumOfAromaticBonds = 0;
 713   for $Bond (@Bonds) {
 714     $BondOrder = $Bond->GetBondOrder();
 715     $SumOfBondOrders += $BondOrder;
 716     if ($BondOrder == 1.5) {
 717       $NumOfAromaticBonds++;
 718     }
 719   }
 720 
 721   if ($NumOfAromaticBonds) {
 722     # As long as aromatic bond orders in a ring are correctly assigned in a ring using
 723     # using 4n + 2 Huckel rule (BondOrder: 1.5) or explicity set as Kekule  bonds (alternate single/double),
 724     # SumOfBondOrders should add up to an integer. Just in case it's not, turn it into an integer...
 725     $SumOfBondOrders = int $SumOfBondOrders;
 726   }
 727 
 728   return $SumOfBondOrders;
 729 }
 730 
 731 # Valence corresponds to number of electrons used by an atom in bonding:
 732 #
 733 #   Valence = ValenceElectrons - ValenceFreeElectrons = BondingElectrons
 734 #
 735 # Single, double, triple bonds with bond orders of 1, 2 and 3 correspond to contribution of
 736 # 1, 2, and 3 bonding electrons. So Valence can be computed using:
 737 #
 738 #   Valence = SumOfBondOrders + FormalCharge
 739 #
 740 # where positive and negative values of FormalCharge increase and decrease the number
 741 # of bonding electrons respectively.
 742 #
 743 # Notes:
 744 #    . For neutral molecules, valence and sum of bond order are equal.
 745 #    . For molecues containing only single bonds, SumOfBondOrders and NumOfBonds are equal.
 746 #
 747 sub GetValence {
 748   my($This) = @_;
 749 
 750   # Is this atom in a molecule?
 751   if (!$This->HasProperty('Molecule')) {
 752     return undef;
 753   }
 754 
 755   # Is Valence property explicitly set?
 756   if ($This->HasProperty('Valence')) {
 757     return $This->GetProperty('Valence');
 758   }
 759   my($FormalCharge, $Valence);
 760 
 761   $Valence = $This->GetSumOfBondOrders();
 762   if ($This->HasProperty('FormalCharge')) {
 763     $Valence += $This->GetProperty('FormalCharge');
 764   }
 765   return $Valence;
 766 }
 767 
 768 # For elements with one one common valence, ImplictValence corresponds to:
 769 #
 770 #   . ImplicitValence = CommonValence + FormalCharge + SpinMultiplicityCorrection
 771 #
 772 # For elements with multiple common valences and no explicit FormalCharge assignment,
 773 # ImplicitValence is determined based on SumOfBondOrdersToNonHydrogenAtoms and next
 774 # available common valence.
 775 #
 776 # For elements with multiple common valences and explicit FormalCharge assignment,
 777 # ImplicitValence corresponds to:
 778 #
 779 #   . ImplicitValence = HighestCommonValence + FormalCharge + SpinMultiplicityCorrection
 780 #
 781 # Notes:
 782 #  . For atoms with explicit assignment of SpinMultiplicity property values corresponding to
 783 #    Singlet (two unparied electrons corresponding to one spin state), Doublet (free radical; an unpaired
 784 #    electron corresponding to two spin states), and Triplet (two unparied electrons corresponding to
 785 #    three spin states; divalent carbon atoms (carbenes)), SpinMultiplicityCorrection is calculated as follows:
 786 #
 787 #       SpinMultiplicity: Doublet(2); SpinMultiplicityCorrection: -1 (one valence electron not available for bonding)
 788 #       SpinMultiplicity: Singlet(1)/Triplet(3); SpinMultiplicityCorrection: -2 (two valence electrons not available for bonding)
 789 #
 790 #
 791 sub GetImplicitValence {
 792   my($This) = @_;
 793 
 794   # Is this atom in a molecule?
 795   if (!$This->HasProperty('Molecule')) {
 796     return undef;
 797   }
 798 
 799   # Is ImplicitValence property explicitly set?
 800   if ($This->HasProperty('ImplicitValence')) {
 801     return $This->GetProperty('ImplicitValence');
 802   }
 803   # Assign implicit valence...
 804   my($AtomicNumber, $CommonValences);
 805 
 806   $AtomicNumber = $This->{AtomicNumber};
 807   if (!$AtomicNumber) {
 808     return undef;
 809   }
 810 
 811   $CommonValences = PeriodicTable::GetElementCommonValences($AtomicNumber);
 812   if (!$CommonValences) {
 813     return undef;
 814   }
 815   my($ImplicitValence, $FormalCharge, $SpinMultiplicityCorrection, @AvailableCommonValences);
 816 
 817   $FormalCharge = $This->GetFormalCharge();
 818   if (!defined $FormalCharge) {
 819     $FormalCharge = 0;
 820   }
 821 
 822   $SpinMultiplicityCorrection = 0;
 823   if ($This->HasProperty('SpinMultiplicity')) {
 824     my($SpinMultiplicity);
 825     $SpinMultiplicity = $This->GetSpinMultiplicity();
 826     SPINMULTIPLICITY: {
 827       if ($SpinMultiplicity =~ /^2$/i) { $SpinMultiplicityCorrection = -1; last SPINMULTIPLICITY; }
 828       if ($SpinMultiplicity =~ /^(1|3)$/i) { $SpinMultiplicityCorrection = -2; last SPINMULTIPLICITY; }
 829       $SpinMultiplicityCorrection = 0;
 830     }
 831   }
 832 
 833   @AvailableCommonValences = split /\,/, $CommonValences;
 834 
 835   # Only only available common valence...
 836   if (@AvailableCommonValences == 1) {
 837     my($CommonValence);
 838 
 839     $CommonValence = $CommonValences;
 840     $ImplicitValence = $CommonValence + $FormalCharge + $SpinMultiplicityCorrection;
 841 
 842     return $ImplicitValence;
 843   }
 844 
 845   # Non-zero formal charge...
 846   if ($FormalCharge != 0) {
 847     # Use HighestCommonValence to set ImplicitValence...
 848     my($HighestCommonValence);
 849 
 850     $HighestCommonValence = $This->GetHighestCommonValence();
 851     $ImplicitValence = $HighestCommonValence + $FormalCharge + $SpinMultiplicityCorrection;
 852 
 853     return $ImplicitValence;
 854   }
 855 
 856   # Set ImplicitValence using SumOfBondOrdersToNonHydrogenAtoms and appropriate CommonValence
 857   # value...
 858   my($SumOfBondOrders, $ValenceFound, $Valence, $AvailableCommonValence);
 859 
 860   $SumOfBondOrders = $This->GetSumOfBondOrdersToNonHydrogenAtoms();
 861   if (!defined $SumOfBondOrders) {
 862     $SumOfBondOrders = 0;
 863   }
 864 
 865   # Get first available valence higher than SumOfBondOrders...
 866   $ValenceFound = 0;
 867   $Valence = 0;
 868   VALENCE: for $AvailableCommonValence (@AvailableCommonValences) {
 869     if ($AvailableCommonValence >= $SumOfBondOrders) {
 870       $ValenceFound = 1;
 871       $Valence = $AvailableCommonValence;
 872       last VALENCE;
 873     }
 874   }
 875   if (!$ValenceFound) {
 876     # Use highest common valence...
 877     $Valence = $This->GetHighestCommonValence();
 878   }
 879   $ImplicitValence = $Valence + $SpinMultiplicityCorrection;
 880 
 881   return $ImplicitValence;
 882 }
 883 
 884 # Get lowest common valence...
 885 sub GetLowestCommonValence {
 886   my($This) = @_;
 887 
 888   # Is LowestCommonValence property explicitly set?
 889   if ($This->HasProperty('LowestCommonValence')) {
 890     return $This->GetProperty('LowestCommonValence');
 891   }
 892   my($AtomicNumber, $LowestCommonValence);
 893 
 894   $AtomicNumber = $This->{AtomicNumber};
 895   if (!$AtomicNumber) {
 896     return undef;
 897   }
 898 
 899   # LowestCommonValence is not set for all elements...
 900   $LowestCommonValence = PeriodicTable::GetElementLowestCommonValence($AtomicNumber);
 901   if (!$LowestCommonValence) {
 902     $LowestCommonValence = undef;
 903   }
 904 
 905   return $LowestCommonValence;
 906 }
 907 
 908 # Get highest common valence...
 909 sub GetHighestCommonValence {
 910   my($This) = @_;
 911 
 912   # Is HighestCommonValence property explicitly set?
 913   if ($This->HasProperty('HighestCommonValence')) {
 914     return $This->GetProperty('HighestCommonValence');
 915   }
 916   my($AtomicNumber, $HighestCommonValence);
 917 
 918   $AtomicNumber = $This->{AtomicNumber};
 919   if (!$AtomicNumber) {
 920     return undef;
 921   }
 922 
 923   # HighestCommonValence is not set for all elements...
 924   $HighestCommonValence = PeriodicTable::GetElementHighestCommonValence($AtomicNumber);
 925   if (!$HighestCommonValence) {
 926     $HighestCommonValence = undef;
 927   }
 928 
 929   return $HighestCommonValence;
 930 }
 931 
 932 # Get valence electrons...
 933 sub GetValenceElectrons {
 934   my($This) = @_;
 935 
 936   # Is ValenceElectrons property explicitly set?
 937   if ($This->HasProperty('ValenceElectrons')) {
 938     return $This->GetProperty('ValenceElectrons');
 939   }
 940   my($AtomicNumber, $ValenceElectrons);
 941 
 942   $AtomicNumber = $This->{AtomicNumber};
 943   if (!$AtomicNumber) {
 944     return undef;
 945   }
 946 
 947   $ValenceElectrons = PeriodicTable::GetElementValenceElectrons($AtomicNumber);
 948 
 949   return $ValenceElectrons;
 950 }
 951 
 952 # Get free non-bodning valence electrons left on atom after taking into account
 953 # sum of bond orders and formal charged on atom.
 954 #
 955 # Valence corresponds to number of electrons used by atom in bonding:
 956 #
 957 #   Valence = ValenceElectrons - ValenceFreeElectrons
 958 #
 959 # Additionally, valence can also be calculated by:
 960 #
 961 #   Valence = SumOfBondOrders + FormalCharge
 962 #
 963 # Implying for neutral molecules, Valence and SumOfBondOrders are equal.
 964 #
 965 # From two formulas for Valence described above, non-bonding free electrons
 966 # left can be computed by:
 967 #
 968 #  ValenceFreeElectrons = ValenceElectrons - SumOfBondOrders - FormalCharge
 969 #
 970 # Examples:
 971 #
 972 # o NH3: ValenceFreeElectrons = 5 - 3 - 0 = 2
 973 # o NH4+; ValenceFreeElectrons = 5 - 4 - 1 = 0
 974 # o C(=O)O- : ValenceFreeElectrons on O- = 6 - 1 + 1 = 6
 975 # o C(=O)O- : ValenceFreeElectrons on =O = 6 - 2 - 0 = 4
 976 #
 977 #
 978 sub GetValenceFreeElectrons {
 979   my($This) = @_;
 980 
 981   # Is this atom in a molecule?
 982   if (!$This->HasProperty('Molecule')) {
 983     return undef;
 984   }
 985 
 986   # Is ValenceFreeElectrons property explicitly set?
 987   if ($This->HasProperty('ValenceFreeElectrons')) {
 988     return $This->GetProperty('ValenceFreeElectrons');
 989   }
 990   my($ValenceElectrons, $SumOfBondOrders, $ValenceFreeElectrons);
 991 
 992   if (!$This->{AtomicNumber}) {
 993     return undef;
 994   }
 995   $ValenceElectrons = $This->GetValenceElectrons();
 996   $SumOfBondOrders = $This->GetSumOfBondOrders();
 997 
 998   $ValenceFreeElectrons = $ValenceElectrons - $SumOfBondOrders;
 999 
1000   if ($This->HasProperty('FormalCharge')) {
1001     $ValenceFreeElectrons -= $This->GetProperty('FormalCharge');
1002   }
1003   return $ValenceFreeElectrons;
1004 }
1005 
1006 # Get num of missing hydrogens...
1007 #
1008 sub GetNumOfMissingHydrogens {
1009   my($This) = @_;
1010   my($NumOfMissingHydrogens);
1011 
1012   $NumOfMissingHydrogens = $This->GetImplicitHydrogens() - $This->GetExplicitHydrogens();
1013 
1014   return $NumOfMissingHydrogens;
1015 }
1016 
1017 # Get number of implicit hydrogen for atom...
1018 #
1019 sub GetImplicitHydrogens {
1020   my($This) = @_;
1021 
1022   # Is this atom in a molecule?
1023   if (!$This->HasProperty('Molecule')) {
1024     return undef;
1025   }
1026 
1027   # Is ImplicitHydrogens property explicitly set?
1028   if ($This->HasProperty('ImplicitHydrogens')) {
1029     return $This->GetProperty('ImplicitHydrogens');
1030   }
1031 
1032   # Is it an element symbol?
1033   if (!$This->{AtomicNumber}) {
1034     return undef;
1035   }
1036   my($ImplicitHydrogens, $ImplicitValence, $SumOfBondOrders);
1037 
1038   $ImplicitHydrogens = 0;
1039   $ImplicitValence = $This->GetImplicitValence();
1040   $SumOfBondOrders = $This->GetSumOfBondOrdersToNonHydrogenAtoms();
1041 
1042   if (defined($ImplicitValence) && defined($SumOfBondOrders)) {
1043     $ImplicitHydrogens = ($SumOfBondOrders >= $ImplicitValence) ? 0 : ($ImplicitValence - $SumOfBondOrders);
1044   }
1045 
1046   return $ImplicitHydrogens;
1047 }
1048 
1049 # Get number of explicit hydrogens for atom...
1050 sub GetExplicitHydrogens {
1051   my($This) = @_;
1052 
1053   return $This->GetNumOfHydrogenAtomNeighbors();
1054 }
1055 
1056 # Add hydrogens to specified atom in molecule and return number of hydrogens added:
1057 #
1058 #   o HydrogensToAdd = ImplicitHydrogenCount - ExplicitHydrogenCount
1059 #
1060 #   o XYZ are set to ZeroVector
1061 #
1062 sub AddHydrogens {
1063   my($This, $HydrogenPositionsWarning) = @_;
1064 
1065   # Is this atom in a molecule?
1066   if (!$This->HasProperty('Molecule')) {
1067     return undef;
1068   }
1069   if (!defined $HydrogenPositionsWarning) {
1070     $HydrogenPositionsWarning = 1;
1071   }
1072   if ($HydrogenPositionsWarning) {
1073     carp "Warning: ${ClassName}->AddHydrogens: The current release of MayaChemTools doesn't assign any hydrogen positions...";
1074   }
1075 
1076   # Is it an element symbol?
1077   if (!$This->{AtomicNumber}) {
1078     return undef;
1079   }
1080 
1081   my($Molecule, $HydrogensAdded, $HydrogensToAdd);
1082 
1083   $Molecule = $This->GetProperty('Molecule');
1084   $HydrogensAdded = 0;
1085   $HydrogensToAdd = $This->GetNumOfMissingHydrogens();
1086   if ($HydrogensToAdd <= 0) {
1087     return $HydrogensAdded;
1088   }
1089 
1090   my($Count, $Hydrogen);
1091 
1092   for $Count (1 .. $HydrogensToAdd) {
1093     $HydrogensAdded++;
1094 
1095     $Hydrogen = $Molecule->NewAtom('AtomSymbol' => 'H', 'XYZ' => [0, 0, 0]);
1096     $Molecule->NewBond('Atoms' => [$This, $Hydrogen], 'BondOrder' => 1);
1097   }
1098 
1099   return $HydrogensAdded;
1100 }
1101 
1102 # Delete hydrogens attached to atom in molecule and return total number of hydrogens deleted...
1103 sub DeleteHydrogens {
1104   my($This) = @_;
1105 
1106   # Is this atom in a molecule?
1107   if (!$This->HasProperty('Molecule')) {
1108     return undef;
1109   }
1110 
1111   # Is it an element symbol?
1112   if (!$This->{AtomicNumber}) {
1113     return undef;
1114   }
1115 
1116   my($Molecule, $Neighbor, $HydrogensDeleted, @Neighbors);
1117 
1118   $Molecule = $This->GetProperty('Molecule');
1119   $HydrogensDeleted = 0;
1120   @Neighbors = $This->GetNeighbors();
1121 
1122   NEIGHBOR: for $Neighbor (@Neighbors) {
1123     if (!$Neighbor->IsHydrogen()) {
1124       next NEIGHBOR;
1125     }
1126     $Molecule->_DeleteAtom($Neighbor);
1127     $HydrogensDeleted++;
1128   }
1129 
1130   return $HydrogensDeleted;
1131 }
1132 
1133 # Copy atom and all its associated data...
1134 sub Copy {
1135   my($This) = @_;
1136   my($Atom);
1137 
1138   $Atom = Storable::dclone($This);
1139 
1140   return $Atom;
1141 }
1142 
1143 # Is it a Hydrogen atom?
1144 sub IsHydrogen {
1145   my($This) = @_;
1146 
1147   return ($This->{AtomicNumber} == 1) ? 1 : 0;
1148 }
1149 
1150 # Is it a Carbon atom?
1151 sub IsCarbon {
1152   my($This) = @_;
1153 
1154   return ($This->{AtomicNumber} == 6) ? 1 : 0;
1155 }
1156 
1157 # Is it a Nitrogen atom?
1158 sub IsNitrogen {
1159   my($This) = @_;
1160 
1161   return ($This->{AtomicNumber} == 7) ? 1 : 0;
1162 }
1163 
1164 # Is it a Oxygen atom?
1165 sub IsOxygen {
1166   my($This) = @_;
1167 
1168   return ($This->{AtomicNumber} == 8) ? 1 : 0;
1169 }
1170 
1171 # Is it a Fluorine atom?
1172 sub IsFluorine {
1173   my($This) = @_;
1174 
1175   return ($This->{AtomicNumber} == 9) ? 1 : 0;
1176 }
1177 
1178 # Is it a Phosphorus atom?
1179 sub IsPhosphorus {
1180   my($This) = @_;
1181 
1182   return ($This->{AtomicNumber} == 15) ? 1 : 0;
1183 }
1184 
1185 # Is it a Sulfur atom?
1186 sub IsSulfur {
1187   my($This) = @_;
1188 
1189   return ($This->{AtomicNumber} == 16) ? 1 : 0;
1190 }
1191 
1192 # Is it a Chlorine atom?
1193 sub IsChlorine {
1194   my($This) = @_;
1195 
1196   return ($This->{AtomicNumber} == 17) ? 1 : 0;
1197 }
1198 
1199 # Is it a Bromine atom?
1200 sub IsBromine {
1201   my($This) = @_;
1202 
1203   return ($This->{AtomicNumber} == 35) ? 1 : 0;
1204 }
1205 
1206 # Is it a Iodine atom?
1207 sub IsIodine {
1208   my($This) = @_;
1209 
1210   return ($This->{AtomicNumber} == 53) ? 1 : 0;
1211 }
1212 
1213 # Is it a hetro atom? (N, O, F, P, S, Cl, Br, I)
1214 sub IsHetroAtom {
1215   my($This) = @_;
1216 
1217   return ($This->{AtomicNumber} =~ /^(7|8|9|15|16|17|35|53)$/) ? 1 : 0;
1218 }
1219 
1220 # Is it a polar atom? ( N, O,  P, S)
1221 sub IsPolarAtom {
1222   my($This) = @_;
1223 
1224   return ($This->{AtomicNumber} =~ /^(7|8|15|16)$/) ? 1 : 0;
1225 }
1226 
1227 # Is aromatic property set for the atom?
1228 sub IsAromatic {
1229   my($This) = @_;
1230   my($Aromatic);
1231 
1232   $Aromatic = $This->GetAromatic();
1233 
1234   return (defined($Aromatic) && $Aromatic) ? 1 : 0;
1235 }
1236 
1237 # Is this a hydrogen atom and attached to one of these atoms: N, O, P, S
1238 sub IsPolarHydrogen {
1239   my($This) = @_;
1240 
1241   if (!$This->IsHydrogen()) {
1242     return 0;
1243   }
1244 
1245   my(@Bonds);
1246   @Bonds = $This->GetBonds();
1247   if (@Bonds > 1) {
1248     return 0;
1249   }
1250 
1251   my($Bond, $BondedAtom);
1252   ($Bond) = @Bonds;
1253   $BondedAtom = $Bond->GetBondedAtom($This);
1254 
1255   return $BondedAtom->IsPolarAtom() ? 1 : 0;
1256 }
1257 
1258 # Is atom in a ring?
1259 #
1260 sub IsInRing {
1261   my($This) = @_;
1262 
1263   # Is this atom in a molecule?
1264   if (!$This->HasProperty('Molecule')) {
1265     return undef;
1266   }
1267   my($Molecule);
1268   $Molecule = $This->GetProperty('Molecule');
1269 
1270   return $Molecule->_IsAtomInRing($This);
1271 }
1272 
1273 # Is atom not in a ring?
1274 #
1275 sub IsNotInRing {
1276   my($This) = @_;
1277 
1278   # Is this atom in a molecule?
1279   if (!$This->HasProperty('Molecule')) {
1280     return undef;
1281   }
1282   my($Molecule);
1283   $Molecule = $This->GetProperty('Molecule');
1284 
1285   return $Molecule->_IsAtomNotInRing($This);
1286 }
1287 
1288 # Is atom only in one ring?
1289 #
1290 sub IsOnlyInOneRing {
1291   my($This) = @_;
1292 
1293   # Is this atom in a molecule?
1294   if (!$This->HasProperty('Molecule')) {
1295     return undef;
1296   }
1297   my($Molecule);
1298   $Molecule = $This->GetProperty('Molecule');
1299 
1300   return $Molecule->_IsAtomInOnlyOneRing($This);
1301 }
1302 
1303 # Is atom in a ring of specific size?
1304 #
1305 sub IsInRingOfSize {
1306   my($This, $RingSize) = @_;
1307 
1308   # Is this atom in a molecule?
1309   if (!$This->HasProperty('Molecule')) {
1310     return undef;
1311   }
1312   my($Molecule);
1313   $Molecule = $This-&