       1 package MolecularFormula;
       2 #
       3 # File: MolecularFormula.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 Text::ParseWords;
      29 use TextUtil;
      30 use PeriodicTable;
      31
      32 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
      33
      34 @ISA = qw(Exporter);
      35 @EXPORT = qw();
      36 @EXPORT_OK = qw(CalculateMolecularWeight CalculateExactMass CalculateElementalComposition FormatCompositionInfomation GetElementsAndCount IsMolecularFormula);
      37
      38 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
      39
      40 #
      41 # Calculate molecular weight assuming its a valid molecular formula...
      42 #
      43 sub CalculateMolecularWeight {
      44   my($MolecularFormula) = @_;
      45   my($Index, $MolecularWeight, $ElementSymbol, $ElementCount, $AtomicWeight, $FormulaElementsRef, $FormulaElementCountRef);
      46
      47   ($FormulaElementsRef, $FormulaElementCountRef) = _ProcessMolecularFormula($MolecularFormula);
      48   if (!(defined($FormulaElementsRef) && defined($FormulaElementCountRef))) {
      49     return undef;
      50   }
      51
      52   $MolecularWeight = 0;
      53
      54   for $Index (0 .. $#{$FormulaElementsRef}) {
      55     $ElementSymbol = $FormulaElementsRef->[$Index];
      56     $ElementCount = $FormulaElementCountRef->[$Index];
      57     $AtomicWeight = PeriodicTable::GetElementAtomicWeight($ElementSymbol);
      58     $MolecularWeight += $AtomicWeight * $ElementCount;
      59   }
      60   return $MolecularWeight;
      61 }
      62
      63 #
      64 # Calculate exact mass assuming it's a valid formula...
      65 #
      66 sub CalculateExactMass {
      67   my($MolecularFormula) = @_;
      68   my($Index, $ElementSymbol, $ElementCount, $ExactMass, $RelativeAtomicMass, $FormulaElementsRef, $FormulaElementCountRef);
      69
      70   ($FormulaElementsRef, $FormulaElementCountRef) = _ProcessMolecularFormula($MolecularFormula);
      71   if (!(defined($FormulaElementsRef) && defined($FormulaElementCountRef))) {
      72     return undef;
      73   }
      74   $ExactMass = 0;
      75
      76   for $Index (0 .. $#{$FormulaElementsRef}) {
      77     $ElementSymbol = $FormulaElementsRef->[$Index];
      78     $ElementCount = $FormulaElementCountRef->[$Index];
      79     $RelativeAtomicMass = PeriodicTable::GetElementMostAbundantNaturalIsotopeMass($ElementSymbol);
      80     if (!defined($RelativeAtomicMass)) {
      81       next ELEMENT;
      82     }
      83     $ExactMass += $RelativeAtomicMass * $ElementCount;
      84   }
      85   return $ExactMass;
      86 }
      87
      88
      89 #
      90 # Calculate elemental composition and return reference to arrays
      91 # containing elements and their percent composition...
      92 #
      93 sub CalculateElementalComposition {
      94   my($MolecularFormula) = @_;
      95   my($Index, $MolecularWeight, $ElementSymbol, $ElementCount, $AtomicWeight, $Composition, $CompositionMultiplier, $FormulaElementsRef, $FormulaElementCountRef, @FormulaElements, @FormulaElementComposition);
      96
      97   $MolecularWeight = CalculateMolecularWeight($MolecularFormula);
      98   if (! defined $MolecularWeight) {
      99     return (undef, undef);
     100   }
     101   ($FormulaElementsRef, $FormulaElementCountRef) = _ProcessMolecularFormula($MolecularFormula);
     102
     103   @FormulaElements = ();
     104   @FormulaElementComposition = ();
     105
     106   if (!$MolecularWeight) {
     107     return ( \@FormulaElements, \@FormulaElementComposition);
     108   }
     109
     110   $CompositionMultiplier = 100 / $MolecularWeight;
     111
     112   for $Index (0 .. $#{$FormulaElementsRef}) {
     113     $ElementSymbol = $FormulaElementsRef->[$Index];
     114     $ElementCount = $FormulaElementCountRef->[$Index];
     115     $AtomicWeight = PeriodicTable::GetElementAtomicWeight($ElementSymbol);
     116     $Composition = ($AtomicWeight * $ElementCount) * $CompositionMultiplier;
     117
     118     push @FormulaElements, $ElementSymbol;
     119     push @FormulaElementComposition, $Composition;
     120   }
     121
     122   return ( \@FormulaElements, \@FormulaElementComposition);
     123 }
     124
     125 # Using refernece to element and its composition arrays, format composition information
     126 # as: Element: Composition;...
     127 #
     128 sub FormatCompositionInfomation {
     129   my($Index, $ElementSymbol, $ElementComposition, $ElementsRef, $ElementCompositionRef, $Precision, $Composition);
     130
     131   $Precision = 2;
     132   if (@_ == 3) {
     133     ($ElementsRef, $ElementCompositionRef, $Precision) = @_;
     134   }
     135   else {
     136     ($ElementsRef, $ElementCompositionRef) = @_;
     137   }
     138
     139   $Composition = '';
     140   for $Index (0 .. $#{$ElementsRef}) {
     141     $ElementSymbol = $ElementsRef->[$Index];
     142     $ElementComposition = $ElementCompositionRef->[$Index];
     143     $ElementComposition = sprintf("%.${Precision}f", $ElementComposition);
     144
     145     $Composition .= ($Composition) ? '; ' : '';
     146     $Composition .=  "${ElementSymbol}: ${ElementComposition}%";
     147   }
     148
     149   return $Composition;
     150 }
     151
     152 #
     153 # Get elements and their count...
     154 #
     155 sub GetElementsAndCount {
     156   my($MolecularFormula) = @_;
     157   my($FormulaElementsRef, $FormulaElementCountRef, $ErrorMsg);
     158
     159   ($FormulaElementsRef, $FormulaElementCountRef, $ErrorMsg) = _ProcessMolecularFormula($MolecularFormula);
     160
     161   return ($FormulaElementsRef, $FormulaElementCountRef);
     162 }
     163
     164 #
     165 # Is it a valid molecular formula?
     166 #
     167 sub IsMolecularFormula {
     168   my($MolecularFormula, $PrintErrorMsg, $Status, $FormulaElementsRef, $FormulaElementCountRef, $ErrorMsg);
     169
     170   ($MolecularFormula) = @_;
     171
     172   ($FormulaElementsRef, $FormulaElementCountRef, $ErrorMsg) = _ProcessMolecularFormula($MolecularFormula);
     173   $Status = (defined($FormulaElementsRef) && defined($FormulaElementCountRef)) ? 1 : 0;
     174
     175   return (wantarray ? ($Status, $ErrorMsg) : $Status);
     176 }
     177
     178 #
     179 # Process molecular formula. For a valid formula, return references to arrays conatining elements
     180 # and element count; otherwsie, return undef.
     181 #
     182 sub _ProcessMolecularFormula {
     183   my($MolecularFormula) = @_;
     184   my($ErrorMsg) = '';
     185
     186   $MolecularFormula = _CleanUpFormula($MolecularFormula);
     187
     188   # Make sure it only contains numbers and letters...
     189   if ($MolecularFormula =~ /[^a-zA-Z0-9\(\)\[\]]/) {
     190     $ErrorMsg = 'Molecular formula contains characters other than a-zA-Z0-9';
     191     return (undef, undef, $ErrorMsg);
     192   }
     193
     194   # Parse the formula...
     195   my($ElementSpec, $FormulaElementSpec, $Spec, $ElementSymbol, $ElementCount,  @FormulaElements, @ElementCount, %FormulaElementsToCountMap, @SubFormulaElements, %SubFormulaElementsToCountMap);
     196
     197   @FormulaElements = (); @ElementCount = ();
     198   %FormulaElementsToCountMap = ();
     199
     200 # Setup element symbol and count regular expression...
     201 # IUPAC: http://www.iupac.org/reports/provisional/abstract04/RB-prs310804/Chap4-3.04.pdf
     202 #
     203
     204   $FormulaElementSpec = qr/
     205                                       \G(                         # $1
     206                                                   (?:
     207                                                       ([A-Z][a-z]?)   # Two or one letter element symbol; $2
     208                                                       ([0-9]*)          # Optionally followed by element count; $3
     209                                                   )
     210                                                   | \( | \[
     211                                                   | \)[0-9]* | \][0-9]*
     212                                                   | .
     213                                             )
     214                                       /x;
     215
     216   my($ProcessingParenthesis);
     217   $ProcessingParenthesis = 0;
     218   # Go over the formula...
     219   FORMULA: while ($MolecularFormula =~ /$FormulaElementSpec/gx) {
     220     ($Spec, $ElementSymbol, $ElementCount) = ($1, $2, $3);
     221
     222     # Handle parenthesis in formula to indicate repeating units...
     223     if ($Spec =~ /^(\(|\[)/) {
     224       if ($ProcessingParenthesis) {
     225         $ErrorMsg = "Molecular formula contains multiple level of () or []";
     226         return (undef, undef, $ErrorMsg);
     227       }
     228       $ProcessingParenthesis = 1;
     229       @SubFormulaElements = ();
     230       %SubFormulaElementsToCountMap = ();
     231       next FORMULA;
     232     }
     233     elsif ($Spec =~ /^(\)|\])/) {
     234       $ProcessingParenthesis = 0;
     235
     236       # Retrieve repeat count and move data to @FormulaElements and %FormulaElementsToCountMap;
     237       my($RepeatCount, $Symbol, $Count);
     238       $RepeatCount = $Spec;
     239       $RepeatCount =~  s/(\)|\])//g;
     240       if (!$RepeatCount) {
     241         $RepeatCount = 1;
     242       }
     243       # Copy data...
     244       for $Symbol (@SubFormulaElements) {
     245         $Count = $SubFormulaElementsToCountMap{$Symbol} * $RepeatCount;
     246         _SetupFormulaElementData(\@FormulaElements, \%FormulaElementsToCountMap, $Symbol, $Count);
     247       }
     248
     249       # Get ready again...
     250       @SubFormulaElements = ();
     251       %SubFormulaElementsToCountMap = ();
     252
     253       next FORMULA;
     254     }
     255
     256     # Retrieve element symbol and count...
     257     $ElementSymbol = ($Spec && !$ElementSymbol) ? $Spec : ($ElementSymbol ? $ElementSymbol : '');
     258     $ElementCount = $ElementCount ? $ElementCount : 1;
     259     if (!PeriodicTable::IsElement($ElementSymbol)) {
     260       $ErrorMsg = "Molecular formula contains unknown elemental symbol $ElementSymbol";
     261       return (undef, undef, $ErrorMsg);
     262     }
     263
     264     if ($ProcessingParenthesis) {
     265       _SetupFormulaElementData(\@SubFormulaElements, \%SubFormulaElementsToCountMap, $ElementSymbol, $ElementCount);
     266     }
     267     else {
     268       _SetupFormulaElementData(\@FormulaElements, \%FormulaElementsToCountMap, $ElementSymbol, $ElementCount);
     269     }
     270   }
     271
     272   # Setup element count array...
     273   for $ElementSymbol (@FormulaElements) {
     274     $ElementCount = $FormulaElementsToCountMap{$ElementSymbol};
     275     push @ElementCount, $ElementCount;
     276   }
     277
     278   # Make sure it all adds up to 100%; otherwise, adjust the last value..
     279
     280   return (\@FormulaElements, \@ElementCount, $ErrorMsg);
     281 }
     282
     283 # Clean it up...
     284 sub _CleanUpFormula {
     285   my($MolecularFormula) = @_;
     286   #Take out any spaces...
     287   $MolecularFormula =~ s/ //g;
     288
     289   # Eliminate any charge specifications: +, - or [1-9]+[+-]
     290   # e.g NO+ [Al(H2O)6]3+ [H2NO3]+
     291   if ($MolecularFormula =~ /[\+\-]/) {
     292     if ($MolecularFormula =~ /\][0-9]+[\+\-]/) {
     293       # Bracket followed optionally by number and then, +/- ...
     294       # [Al(H2O)6]3+ ...
     295       $MolecularFormula =~ s/\][0-9]+[\+\-]/\]/g;
     296     }
     297     elsif ($MolecularFormula =~ /[\+\-][0-9]*/) {
     298       # +/- followed optionally by a number...
     299       # C37H42N2O6+2, Cu+
     300       $MolecularFormula =~ s/[\+\-][0-9]*//g;
     301     }
     302   }
     303
     304   # Eliminate any brackets - ] or ) - not followed by numbers:
     305   # e.g. Li[H2PO4]
     306   if ($MolecularFormula !~ /\][0-9]+/) {
     307     $MolecularFormula =~ s/[\[\]]//g;
     308   }
     309   if ($MolecularFormula !~ /\)[0-9]+/) {
     310     $MolecularFormula =~ s/[\(\)]//g;
     311   }
     312   # Change adducts to parenthesis format...
     313   # Na2CO3.10H2O -> Na2CO3(H2O)10
     314   # 3CdSO4.8H2O -> (CdSO4)3(H2O)8
     315   if ($MolecularFormula =~ /\./) {
     316     my($SubFormula, $Count, $Spec);
     317     my(@MolecularFormulaSplits) = split /\./, $MolecularFormula;
     318     $MolecularFormula = '';
     319     for $SubFormula (@MolecularFormulaSplits) {
     320       ($Count, $Spec) = $SubFormula =~ /^([0-9]*)(.*?)$/;
     321       if ($Count) {
     322         $MolecularFormula .= "(${Spec})${Count}";
     323       }
     324       else {
     325         $MolecularFormula .= $Spec;
     326       }
     327     }
     328   }
     329
     330   return $MolecularFormula;
     331 }
     332
     333 # Store the element and count...
     334 sub _SetupFormulaElementData {
     335   my($ElementsRef, $ElementsToCountMapRef, $Element, $Count) = @_;
     336
     337   if (exists $ElementsToCountMapRef->{$Element}) {
     338     $ElementsToCountMapRef->{$Element} += $Count;
     339   }
     340   else {
     341     push @{$ElementsRef}, $Element;
     342     $ElementsToCountMapRef->{$Element} = $Count;
     343   }
     344 }
     345