       1 package SDFileUtil;
       2 #
       3 # File: SDFileUtil.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 Exporter;
      28 use Carp;
      29 use PeriodicTable qw(IsElement);
      30 use TimeUtil ();
      31
      32 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
      33
      34 @ISA = qw(Exporter);
      35 @EXPORT = qw(GenerateCmpdAtomLine GenerateCmpdBondLine GenerateCmpdChargePropertyLines GenerateCmpdCommentsLine GenerateCmpdCountsLine GenerateCmpdAtomAliasPropertyLines GenerateCmpdIsotopePropertyLines GenerateCmpdDataHeaderLabelsAndValuesLines GenerateCmpdMiscInfoLine GenerateCmpdRadicalPropertyLines GenerateCmpdMolNameLine GenerateEmptyCtabBlockLines GenerateMiscLineDateStamp GetAllAndCommonCmpdDataHeaderLabels GetCmpdDataHeaderLabels GetCmpdDataHeaderLabelsAndValues GetCmpdFragments GetCtabLinesCount GetUnknownAtoms GetInvalidAtomNumbers MDLChargeToInternalCharge InternalChargeToMDLCharge MDLBondTypeToInternalBondOrder InternalBondOrderToMDLBondType MDLBondStereoToInternalBondStereochemistry InternalBondStereochemistryToMDLBondStereo InternalSpinMultiplicityToMDLRadical MDLRadicalToInternalSpinMultiplicity IsCmpd3D IsCmpd2D ParseCmpdAtomLine ParseCmpdBondLine ParseCmpdCommentsLine ParseCmpdCountsLine ParseCmpdMiscInfoLine ParseCmpdMolNameLine ParseCmpdAtomAliasPropertyLine ParseCmpdChargePropertyLine ParseCmpdIsotopePropertyLine ParseCmpdRadicalPropertyLine ReadCmpdString RemoveCmpdDataHeaderLabelAndValue WashCmpd);
      36 @EXPORT_OK = qw();
      37 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
      38
      39 # Format data for compounds count line...
      40 sub GenerateCmpdCountsLine {
      41   my($AtomCount, $BondCount, $ChiralFlag, $PropertyCount, $Version, $Line);
      42
      43   if (@_ == 5) {
      44     ($AtomCount, $BondCount, $ChiralFlag, $PropertyCount, $Version) = @_;
      45   }
      46   elsif (@_ == 3) {
      47     ($AtomCount, $BondCount, $ChiralFlag) = @_;
      48     $PropertyCount = 999;
      49     $Version = "V2000";
      50   }
      51   else {
      52     ($AtomCount, $BondCount) = @_;
      53     $ChiralFlag = 0;
      54     $PropertyCount = 999;
      55     $Version = "V2000";
      56   }
      57   if ($AtomCount > 999) {
      58     croak "Error: SDFileUtil::GenerateCmpdCountsLine: The atom count, $AtomCount, exceeds maximum of 999 allowed for CTAB version 2000. The Extended Connection Table (V3000) format in MDL MOL and SD files is not supported by the current release of MayaChemTools...";
      59   }
      60   $Line = sprintf "%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i%6s", $AtomCount, $BondCount, 0, 0, $ChiralFlag, 0, 0, 0, 0, 0, $PropertyCount, $Version;
      61
      62   return ($Line);
      63 }
      64
      65 # Generate comments line...
      66 sub GenerateCmpdCommentsLine {
      67   my($Comments) = @_;
      68   my($Line);
      69
      70   $Line = (length($Comments) > 80) ? substr($Comments, 0, 80) : $Comments;
      71
      72   return $Line;
      73 }
      74
      75 # Generate molname line...
      76 sub GenerateCmpdMolNameLine {
      77   my($MolName) = @_;
      78   my($Line);
      79
      80   $Line = (length($MolName) > 80) ? substr($MolName, 0, 80) : $MolName;
      81
      82   return $Line;
      83 }
      84
      85 # Generate data for compounds misc info line...
      86 sub GenerateCmpdMiscInfoLine {
      87   my($ProgramName, $UserInitial, $Code) = @_;
      88   my($Date, $Line);
      89
      90   if (!(defined($ProgramName) && $ProgramName)) {
      91     $ProgramName = "MayaChem";
      92   }
      93   if (!(defined($UserInitial) && $UserInitial)) {
      94     $UserInitial = "  ";
      95   }
      96   if (!(defined($Code) && $Code)) {
      97     $Code = "2D";
      98   }
      99
     100   if (length($ProgramName) > 8) {
     101     $ProgramName = substr($ProgramName, 0, 8);
     102   }
     103   if (length($UserInitial) > 2) {
     104     $UserInitial = substr($UserInitial, 0, 2);
     105   }
     106   if (length($Code) > 2) {
     107     $Code = substr($Code, 0, 2);
     108   }
     109   $Date = GenerateMiscLineDateStamp();
     110
     111   $Line = "${UserInitial}${ProgramName}${Date}${Code}";
     112
     113   return $Line;
     114 }
     115
     116 # Generate data for compounds misc info line...
     117 sub GenerateEmptyCtabBlockLines {
     118   my($Date, $Lines);
     119
     120   if (@_ == 1) {
     121     ($Date) = @_;
     122   }
     123   else {
     124     $Date = GenerateMiscLineDateStamp();
     125   }
     126   # First line: Blank molname line...
     127   # Second line: Misc info...
     128   # Third line: Blank comments line...
     129   # Fourth line: Counts line reflecting empty structure data block...
     130   $Lines = "\n";
     131   $Lines .= "  MayaChem${Date}2D\n";
     132   $Lines .= "\n";
     133   $Lines .= GenerateCmpdCountsLine(0, 0, 0) . "\n";
     134   $Lines .= "M  END";
     135
     136   return $Lines;
     137 }
     138
     139 # Generate SD file data stamp...
     140 sub GenerateMiscLineDateStamp {
     141   return TimeUtil::SDFileTimeStamp();
     142 }
     143
     144 # Generate data for compound atom line...
     145 #
     146 sub GenerateCmpdAtomLine {
     147   my($AtomSymbol, $AtomX, $AtomY, $AtomZ, $MassDifference, $Charge, $StereoParity) = @_;
     148   my($Line);
     149
     150   if (!defined $MassDifference) {
     151     $MassDifference = 0;
     152   }
     153   if (!defined $Charge) {
     154     $Charge = 0;
     155   }
     156   if (!defined $StereoParity) {
     157     $StereoParity = 0;
     158   }
     159   $Line = sprintf "%10.4f%10.4f%10.4f %-3s%2i%3i%3i  0  0  0  0  0  0  0  0  0", $AtomX, $AtomY, $AtomZ, $AtomSymbol, $MassDifference, $Charge, $StereoParity;
     160
     161   return $Line
     162 }
     163
     164 # Generate data for compound bond line...
     165 #
     166 sub GenerateCmpdBondLine {
     167   my($FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo) = @_;
     168   my($Line);
     169
     170   if (!defined $BondStereo) {
     171     $BondStereo = 0;
     172   }
     173   $Line = sprintf "%3i%3i%3i%3i  0  0  0", $FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo;
     174
     175   return $Line
     176 }
     177
     178 # Generate charge property lines for CTAB block...
     179 #
     180 sub GenerateCmpdChargePropertyLines {
     181   my($ChargeValuePairsRef) = @_;
     182
     183   return _GenerateCmpdGenericPropertyLines('Charge', $ChargeValuePairsRef);
     184 }
     185
     186 # Generate isotope property lines for CTAB block...
     187 #
     188 sub GenerateCmpdIsotopePropertyLines {
     189   my($IsotopeValuePairsRef) = @_;
     190
     191   return _GenerateCmpdGenericPropertyLines('Isotope', $IsotopeValuePairsRef);
     192 }
     193
     194 # Generate radical property line property lines for CTAB block...
     195 #
     196 sub GenerateCmpdRadicalPropertyLines {
     197   my($RadicalValuePairsRef) = @_;
     198
     199   return _GenerateCmpdGenericPropertyLines('Radical', $RadicalValuePairsRef);
     200 }
     201
     202 # Generate atom alias property line property lines for CTAB block...
     203 #
     204 # Atom alias property line format:
     205 #
     206 # A  aaa
     207 # x...
     208 #
     209 #    aaa: Atom number
     210 #    x: Atom alias in next line
     211 #
     212 sub GenerateCmpdAtomAliasPropertyLines {
     213   my($PropertyValuePairsRef) = @_;
     214   my($Index, $AtomNum, $AtomAlias, $Line, @PropertyLines);
     215
     216   @PropertyLines = ();
     217
     218   for ($Index = 0; $Index < $#{$PropertyValuePairsRef}; $Index += 2) {
     219     $AtomNum = $PropertyValuePairsRef->[$Index];
     220     $AtomAlias = $PropertyValuePairsRef->[$Index + 1];
     221
     222     $Line = "A  " . sprintf "%3i", $AtomNum;
     223
     224     push @PropertyLines, $Line;
     225     push @PropertyLines, $AtomAlias;
     226   }
     227
     228   return @PropertyLines;
     229 }
     230
     231 # Generate data header labels and values lines...
     232 #
     233 sub GenerateCmpdDataHeaderLabelsAndValuesLines {
     234   my($DataHeaderLabelsRef, $DataHeaderLabelsAndValuesRef, $SortDataLabels) = @_;
     235   my($DataLabel, $DataValue, @DataLabels, @DataLines);
     236
     237   if (!defined $SortDataLabels) {
     238     $SortDataLabels = 0;
     239   }
     240
     241   @DataLines = ();
     242   @DataLabels = ();
     243   if ($SortDataLabels) {
     244     push @DataLabels, sort @{$DataHeaderLabelsRef};
     245   }
     246   else {
     247     push @DataLabels,  @{$DataHeaderLabelsRef};
     248   }
     249   for $DataLabel (@DataLabels) {
     250     $DataValue = '';
     251     if (exists $DataHeaderLabelsAndValuesRef->{$DataLabel}) {
     252       $DataValue = $DataHeaderLabelsAndValuesRef->{$DataLabel};
     253     }
     254     push @DataLines, (">  <${DataLabel}>", "$DataValue", "");
     255   }
     256   return @DataLines;
     257 }
     258
     259 # Parse data field header in SD file and return lists of all and common data field
     260 # labels.
     261 sub GetAllAndCommonCmpdDataHeaderLabels {
     262   my($SDFileRef) = @_;
     263   my($CmpdCount, $CmpdString, $Label, @CmpdLines, @DataFieldLabels, @CommonDataFieldLabels, %DataFieldLabelsMap);
     264
     265   $CmpdCount = 0;
     266   @DataFieldLabels = ();
     267   @CommonDataFieldLabels = ();
     268   %DataFieldLabelsMap = ();
     269
     270   while ($CmpdString = ReadCmpdString($SDFileRef)) {
     271     $CmpdCount++;
     272     @CmpdLines = split "\n", $CmpdString;
     273     # Process compound data header labels and figure out which ones are present for
     274     # all the compounds...
     275     if (@DataFieldLabels) {
     276       my (@CmpdDataFieldLabels) = GetCmpdDataHeaderLabels(\@CmpdLines);
     277       my(%CmpdDataFieldLabelsMap) = ();
     278       # Setup a map for the current labels...
     279       for $Label (@CmpdDataFieldLabels) {
     280         $CmpdDataFieldLabelsMap{$Label} = "PresentInSome";
     281       }
     282       # Check the presence old labels for this compound; otherwise, mark 'em new...
     283       for $Label (@DataFieldLabels) {
     284         if (!$CmpdDataFieldLabelsMap{$Label}) {
     285           $DataFieldLabelsMap{$Label} = "PresentInSome";
     286         }
     287       }
     288       # Check the presence this compound in the old labels; otherwise, add 'em...
     289       for $Label (@CmpdDataFieldLabels ) {
     290         if (!$DataFieldLabelsMap{$Label}) {
     291           # It's a new label...
     292           push @DataFieldLabels, $Label;
     293           $DataFieldLabelsMap{$Label} = "PresentInSome";
     294         }
     295       }
     296     }
     297     else {
     298       # Get the initial label set and set up a map...
     299       @DataFieldLabels = GetCmpdDataHeaderLabels(\@CmpdLines);
     300       for $Label (@DataFieldLabels) {
     301         $DataFieldLabelsMap{$Label} = "PresentInAll";
     302       }
     303     }
     304   }
     305   # Identify the common data field labels...
     306   @CommonDataFieldLabels = ();
     307   for $Label (@DataFieldLabels) {
     308     if ($DataFieldLabelsMap{$Label} eq "PresentInAll") {
     309       push @CommonDataFieldLabels, $Label;
     310     }
     311   }
     312   return ($CmpdCount, \@DataFieldLabels, \@CommonDataFieldLabels);
     313 }
     314
     315 # Parse all the data header labels and return 'em as an list...
     316 #
     317 # Format:
     318 #
     319 #> Data header line
     320 #Data line(s)
     321 #Blank line
     322 #
     323 # [Data Header] (one line) precedes each item of data, starts with a greater than (>) sign, and
     324 # contains at least one of the following:
     325 #  The field name enclosed in angle brackets. For example: <melting.point>
     326 #  The field number, DTn , where n represents the number assigned to the field in a MACCS-II database
     327 #
     328 #Optional information for the data header includes:
     329 #  The compound’s external and internal registry numbers. External registry numbers must be enclosed in parentheses.
     330 #  Any combination of information
     331 #
     332 #The following are examples of valid data headers:
     333 #> <MELTING.POINT>
     334 #> 55 (MD-08974) <BOILING.POINT> DT12
     335 #> DT12 55
     336 #> (MD-0894) <BOILING.POINT> FROM ARCHIVES
     337 #
     338 #Notes: Sometimes last blank line is missing and can be just followed by $$$$
     339 #
     340 sub GetCmpdDataHeaderLabels {
     341   my($CmpdLines) = @_;
     342   my($CmpdLine, $Label, @Labels);
     343
     344   @Labels = ();
     345   CMPDLINE: for $CmpdLine (@$CmpdLines) {
     346     if ($CmpdLine !~ /^>/) {
     347       next CMPDLINE;
     348     }
     349     # Does the line contains field name enclosed in angular brackets?
     350     ($Label) = $CmpdLine =~ /<.*?>/g;
     351     if (!defined($Label)) {
     352       next CMPDLINE;
     353     }
     354     $Label =~ s/(<|>)//g;
     355     push @Labels, $Label;
     356   }
     357   return (@Labels);
     358 }
     359
     360 # Parse all the data header labels and values
     361 sub GetCmpdDataHeaderLabelsAndValues {
     362   my($CmpdLines) = @_;
     363   my($CmpdLine, $CurrentLabel, $Label, $Value, $ValueCount, $ProcessingLabelData, @Values, %DataFields);
     364
     365   %DataFields = ();
     366   $ProcessingLabelData = 0;
     367   $ValueCount = 0;
     368   CMPDLINE: for $CmpdLine (@$CmpdLines) {
     369     if ($CmpdLine =~ /^\$\$\$\$/) {
     370       last CMPDLINE;
     371     }
     372     if ($CmpdLine =~ /^>/) {
     373       # Does the line contains field name enclosed in angular brackets?
     374       ($Label) = $CmpdLine =~ /<.*?>/g;
     375       if (defined $Label) {
     376         $CurrentLabel = $Label;
     377         $CurrentLabel =~ s/(<|>)//g;
     378         $ProcessingLabelData = 0;
     379         $ValueCount = 0;
     380
     381         if ($CurrentLabel) {
     382           $ProcessingLabelData = 1;
     383           $DataFields{$CurrentLabel} = '';
     384           next CMPDLINE;
     385         }
     386       }
     387       else {
     388         if (!$ProcessingLabelData) {
     389           # Data line containing no <label> as allowed by SDF format. Just ignore it...
     390           next CMPDLINE;
     391         }
     392       }
     393     }
     394     if (!$ProcessingLabelData) {
     395       next CMPDLINE;
     396     }
     397     if (!(defined($CmpdLine) && length($CmpdLine))) {
     398       # Blank line terminates value for a label...
     399       $CurrentLabel = '';
     400       $ValueCount = 0;
     401       $ProcessingLabelData = 0;
     402       next CMPDLINE;
     403     }
     404     $ValueCount++;
     405     $Value = $CmpdLine;
     406
     407     if ($ValueCount > 1) {
     408       $DataFields{$CurrentLabel} .= "\n" . $Value;
     409     }
     410     else {
     411       $DataFields{$CurrentLabel} = $Value;
     412     }
     413   }
     414   return (%DataFields);
     415 }
     416
     417 # Return an updated compoud string after removing  data header label along with its
     418 # value from the specified compound string...
     419 #
     420 sub RemoveCmpdDataHeaderLabelAndValue {
     421   my($CmpdString, $DataHeaderLabel) = @_;
     422   my($Line, $PorcessingDataHeaderLabel, @CmpdLines);
     423
     424   @CmpdLines = ();
     425   $PorcessingDataHeaderLabel = 0;
     426
     427   CMPDLINE: for $Line (split "\n", $CmpdString) {
     428     if ($Line =~ /^>/ && $Line =~ /<$DataHeaderLabel>/i) {
     429       $PorcessingDataHeaderLabel = 1;
     430       next CMPDLINE;
     431     }
     432
     433     if ($PorcessingDataHeaderLabel) {
     434       # Blank line indicates end of fingerprints data value...
     435       if ($Line =~ /^\$\$\$\$/) {
     436         push @CmpdLines, $Line;
     437         $PorcessingDataHeaderLabel = 0;
     438       }
     439       elsif (!length($Line)) {
     440         $PorcessingDataHeaderLabel = 0;
     441       }
     442       next CMPDLINE;
     443     }
     444
     445     # Track compound lines without fingerprints data...
     446     push @CmpdLines, $Line;
     447   }
     448
     449   return join "\n", @CmpdLines;
     450 }
     451
     452 #
     453 # Using bond blocks, figure out the number of disconnected fragments  and
     454 # return their values along with the atom numbers in a string delimited by new
     455 # line character.
     456 #
     457 sub GetCmpdFragments {
     458   my($CmpdLines) = @_;
     459   my($AtomCount, $BondCount, $FirstAtomNum, $SecondAtomNum, @AtomConnections, $BondType, $FragmentString, $FragmentCount, $LineIndex, $Index, $AtomNum, $NbrAtomNum, @ProcessedAtoms, $ProcessedAtomCount, $ProcessAtomNum, @ProcessingAtoms, @ConnectedAtoms, %Fragments, $FragmentNum, $AFragmentString);
     460
     461   # Setup the connection table for each atom...
     462   @AtomConnections = ();
     463   ($AtomCount, $BondCount) = ParseCmpdCountsLine(@$CmpdLines[3]);
     464   for $AtomNum (1 .. $AtomCount) {
     465     %{$AtomConnections[$AtomNum]} = ();
     466   }
     467   for ($LineIndex = 4 + $AtomCount; $LineIndex < (4 + $AtomCount + $BondCount); $LineIndex++) {
     468     ($FirstAtomNum, $SecondAtomNum, $BondType) = ParseCmpdBondLine(@$CmpdLines[$LineIndex]);
     469     if (!$AtomConnections[$FirstAtomNum]{$SecondAtomNum}) {
     470       $AtomConnections[$FirstAtomNum]{$SecondAtomNum} = $BondType;
     471     }
     472     if (!$AtomConnections[$SecondAtomNum]{$FirstAtomNum}) {
     473       $AtomConnections[$SecondAtomNum]{$FirstAtomNum} = $BondType;
     474     }
     475   }
     476
     477   #Get set to count fragments...
     478   $ProcessedAtomCount = 0;
     479   $FragmentNum = 0;
     480   %Fragments = ();
     481   @ProcessedAtoms = ();
     482   for $AtomNum (1 .. $AtomCount) {
     483     $ProcessedAtoms[$AtomNum] = 0;
     484   }
     485   while ($ProcessedAtomCount < $AtomCount) {
     486     @ProcessingAtoms = ();
     487     @ConnectedAtoms = ();
     488     ATOMNUM: for $AtomNum (1 .. $AtomCount) {
     489       if (!$ProcessedAtoms[$AtomNum]) {
     490         $ProcessedAtomCount++;
     491         $ProcessedAtoms[$AtomNum] = 1;
     492         push @ProcessingAtoms, $AtomNum;
     493         $FragmentNum++;
     494         @{$Fragments{$FragmentNum} } = ();
     495         push @{$Fragments{$FragmentNum} }, $AtomNum;
     496         last ATOMNUM;
     497       }
     498     }
     499
     500     # Go over the neighbors and follow the connection trail while collecting the
     501     # atoms numbers present in the connected fragment...
     502     while (@ProcessingAtoms) {
     503       for ($Index = 0; $Index < @ProcessingAtoms; $Index++) {
     504         $ProcessAtomNum = $ProcessingAtoms[$Index];
     505         for $NbrAtomNum (keys %{$AtomConnections[$ProcessAtomNum]})  {
     506           if (!$ProcessedAtoms[$NbrAtomNum]) {
     507             $ProcessedAtomCount++;
     508             $ProcessedAtoms[$NbrAtomNum] = 1;
     509             push @ConnectedAtoms, $NbrAtomNum;
     510             push @{ $Fragments{$FragmentNum} }, $NbrAtomNum;
     511           }
     512         }
     513       }
     514       @ProcessingAtoms = ();
     515       @ProcessingAtoms = @ConnectedAtoms;
     516       @ConnectedAtoms = ();
     517     }
     518   }
     519   $FragmentCount = $FragmentNum;
     520   $FragmentString = "";
     521
     522   # Sort out the fragments by size...
     523   for $FragmentNum (sort { @{$Fragments{$b}} <=> @{$Fragments{$a}}  } keys %Fragments ) {
     524     # Sort the atoms in a fragment by their numbers...
     525     $AFragmentString = join " ", sort { $a <=> $b } @{ $Fragments{$FragmentNum} };
     526     if ($FragmentString) {
     527       $FragmentString .=  "\n" . $AFragmentString;
     528     }
     529     else {
     530       $FragmentString = $AFragmentString;
     531     }
     532   }
     533   return ($FragmentCount, $FragmentString);
     534 }
     535
     536 # Count number of lines present in between 4th and line containg "M END"
     537 sub GetCtabLinesCount {
     538   my($CmpdLines) = @_;
     539   my($LineIndex, $CtabLinesCount);
     540
     541   $CtabLinesCount = 0;
     542  LINE: for ($LineIndex = 4; $LineIndex < @$CmpdLines; $LineIndex++) {
     543     #
     544     # Any line after atom and bond data starting with anything other than space or
     545     # a digit indicates end of Ctab atom/bond data block...
     546     #
     547     if (@$CmpdLines[$LineIndex] !~ /^[0-9 ]/) {
     548       $CtabLinesCount = $LineIndex - 4;
     549       last LINE;
     550     }
     551   }
     552   return $CtabLinesCount;
     553 }
     554
     555 # Using atom blocks, count the number of atoms which contain special element
     556 # symbols not present in the periodic table.
     557 sub GetUnknownAtoms {
     558   my($CmpdLines) = @_;
     559   my($UnknownAtomCount, $UnknownAtoms, $UnknownAtomLines, $LineIndex, $AtomCount, $AtomSymbol);
     560
     561   $UnknownAtomCount = 0;
     562   $UnknownAtoms = "";
     563   $UnknownAtomLines = "";
     564   ($AtomCount) = ParseCmpdCountsLine(@$CmpdLines[3]);
     565   for ($LineIndex = 4; $LineIndex < (4 + $AtomCount); $LineIndex++) {
     566     ($AtomSymbol) = ParseCmpdAtomLine(@$CmpdLines[$LineIndex]);
     567     if (!IsElement($AtomSymbol)) {
     568       $UnknownAtomCount++;
     569       $UnknownAtoms .= " $AtomSymbol";
     570       if ($UnknownAtomLines) {
     571         $UnknownAtomLines .= "\n" . @$CmpdLines[$LineIndex];
     572       }
     573       else {
     574         $UnknownAtomLines = @$CmpdLines[$LineIndex];
     575       }
     576     }
     577   }
     578   return ($UnknownAtomCount, $UnknownAtoms, $UnknownAtomLines);
     579 }
     580
     581 # Check z coordinates of all atoms to see whether any of them is non-zero
     582 # which makes the compound geometry three dimensional...
     583 #
     584 sub IsCmpd3D {
     585   my($CmpdLines) = @_;
     586   my($LineIndex, $AtomCount, $AtomSymbol, $AtomX, $AtomY, $AtomZ);
     587
     588   ($AtomCount) = ParseCmpdCountsLine(@$CmpdLines[3]);
     589   for ($LineIndex = 4; $LineIndex < (4 + $AtomCount); $LineIndex++) {
     590     ($AtomSymbol, $AtomX, $AtomY, $AtomZ) = ParseCmpdAtomLine(@$CmpdLines[$LineIndex]);
     591     if ($AtomZ != 0) {
     592       return 1;
     593     }
     594   }
     595   return 0;
     596 }
     597
     598 # Check whether it's a 2D compound...
     599 #
     600 sub IsCmpd2D {
     601   my($CmpdLines) = @_;
     602
     603   return IsCmpd3D($CmpdLines) ? 0 : 1;
     604 }
     605
     606 # Using bond blocks, count the number of bond lines which contain atom numbers
     607 # greater than atom count specified in compound count line...
     608 #
     609 sub GetInvalidAtomNumbers {
     610   my($CmpdLines) = @_;
     611   my($LineIndex, $AtomCount, $BondCount, $FirstAtomNum, $SecondAtomNum, $InvalidAtomNumbersCount, $InvalidAtomNumbers, $InvalidAtomNumberLines, $Line, $InvalidAtomPropertyLine, $ValuePairIndex, $AtomNum, $Value, @ValuePairs);
     612
     613   ($AtomCount, $BondCount) = ParseCmpdCountsLine(@$CmpdLines[3]);
     614
     615   $InvalidAtomNumbersCount = 0;
     616   $InvalidAtomNumbers = "";
     617   $InvalidAtomNumberLines = "";
     618
     619   # Go over bond block lines...
     620   LINE: for ($LineIndex = 4 + $AtomCount; $LineIndex < (4 + $AtomCount + $BondCount); $LineIndex++) {
     621     ($FirstAtomNum, $SecondAtomNum) = ParseCmpdBondLine(@$CmpdLines[$LineIndex]);
     622     if ($FirstAtomNum <= $AtomCount && $SecondAtomNum <= $AtomCount) {
     623       next LINE;
     624     }
     625     if ($FirstAtomNum > $AtomCount) {
     626       $InvalidAtomNumbersCount++;
     627       $InvalidAtomNumbers .= " $FirstAtomNum";
     628     }
     629     if ($SecondAtomNum > $AtomCount) {
     630       $InvalidAtomNumbersCount++;
     631       $InvalidAtomNumbers .= " $SecondAtomNum";
     632     }
     633     if ($InvalidAtomNumberLines) {
     634       $InvalidAtomNumberLines .= "\n" . @$CmpdLines[$LineIndex];
     635     }
     636     else {
     637       $InvalidAtomNumberLines = @$CmpdLines[$LineIndex];
     638     }
     639   }
     640   # Go over property lines before M  END...
     641   #
     642   LINE: for ($LineIndex = (4 + $AtomCount + $BondCount); $LineIndex < @$CmpdLines; $LineIndex++) {
     643     $Line = @$CmpdLines[$LineIndex];
     644     @ValuePairs = ();
     645     if ($Line =~ /^M  END/i) {
     646       last LINE;
     647     }
     648     @ValuePairs = ();
     649     if ($Line =~ /^M  CHG/i) {
     650       @ValuePairs = ParseCmpdChargePropertyLine($Line);
     651     }
     652     elsif ($Line =~ /^M  RAD/i) {
     653       @ValuePairs = ParseCmpdRadicalPropertyLine($Line);
     654     }
     655     elsif ($Line =~ /^M  ISO/i) {
     656       @ValuePairs = ParseCmpdIsotopePropertyLine($Line);
     657     }
     658     elsif ($Line =~ /^A  /i) {
     659       my($NextLine);
     660       $LineIndex++;
     661       $NextLine = @$CmpdLines[$LineIndex];
     662       @ValuePairs = ParseCmpdAtomAliasPropertyLine($Line, $NextLine);
     663     }
     664     else {
     665       next LINE;
     666     }
     667
     668     $InvalidAtomPropertyLine = 0;
     669     for ($ValuePairIndex = 0; $ValuePairIndex < $#ValuePairs; $ValuePairIndex += 2) {
     670       $AtomNum = $ValuePairs[$ValuePairIndex]; $Value = $ValuePairs[$ValuePairIndex + 1];
     671       if ($AtomNum > $AtomCount) {
     672         $InvalidAtomPropertyLine = 1;
     673         $InvalidAtomNumbersCount++;
     674         $InvalidAtomNumbers .= " $AtomNum";
     675       }
     676     }
     677     if ($InvalidAtomPropertyLine) {
     678       if ($InvalidAtomNumberLines) {
     679         $InvalidAtomNumberLines .= "\n" . $Line;
     680       }
     681       else {
     682         $InvalidAtomNumberLines = $Line;
     683       }
     684     }
     685   }
     686
     687   return ($InvalidAtomNumbersCount, $InvalidAtomNumbers, $InvalidAtomNumberLines);
     688 }
     689
     690 # Ctab lines: Atom block
     691 #
     692 # Format: xxxxx.xxxxyyyyy.yyyyzzzzz.zzzz aaaddcccssshhhbbbvvvHHHrrriiimmmnnneee
     693 #         A10       A10       A10       xA3 A2A3 A3 A3 A3 A3 A3 A3 A3 A3 A3 A3
     694 # x,y,z: Atom coordinates
     695 # aaa: Atom symbol. Entry in periodic table or L for atom list, A, Q, * for unspecified
     696 #      atom, and LP for lone pair, or R# for Rgroup label
     697 # dd: Mass difference. -3, -2, -1, 0, 1, 2, 3, 4 (0 for value beyond these limits)
     698 # ccc: Charge. 0 = uncharged or value other than these, 1 = +3, 2 = +2, 3 = +1,
     699 #      4 = doublet radical, 5 = -1, 6 = -2, 7 = -3
     700 # sss: Atom stereo parity. 0 = not stereo, 1 = odd, 2 = even, 3 = either or unmarked stereo center
     701 # hhh: Hydrogen count + 1. 1 = H0, 2 = H1, 3 = H2, 4 = H3, 5 = H4
     702 # bbb: Stereo care box. 0 = ignore stereo configuration of this double bond atom, 1 = stereo
     703 #      configuration of double bond atom must match
     704 # vvv: Valence. 0 = no marking (default)(1 to 14) = (1 to 14) 15 = zero valence
     705 # HHH: H0 designator. 0 = not specified, 1 = no H atoms allowed (redundant due to hhh)
     706 # rrr: Not used
     707 # iii: Not used
     708 # mmm: Atom-atom mapping number. 1 - number of atoms
     709 # nnn: Inversion/retention flag. 0 = property not applied, 1 = configuration is inverted,
     710 #      2 = configuration is retained.
     711 # eee: Exact change flag. 0 = property not applied, 1 = change on atom must be
     712 #      exactly as shown
     713 #
     714 # Notes:
     715 #  . StereoParity: 1 - ClockwiseStereo, 2 - AntiClockwiseStereo; 3 - Either; 0 - none. These
     716 #    values determine chirailty around the chiral center; a non zero value indicates atom
     717 #    has been marked as chiral center.
     718 #
     719 sub ParseCmpdAtomLine {
     720   my($Line) = @_;
     721   my ($LineIndex, $AtomX, $AtomY, $AtomZ, $AtomSymbol, $MassDifference, $Charge, $StereoParity);
     722
     723   ($AtomX, $AtomY, $AtomZ, $AtomSymbol, $MassDifference, $Charge, $StereoParity) = ('') x 7;
     724   if (length($Line) > 31) {
     725     ($AtomX, $AtomY, $AtomZ, $AtomSymbol, $MassDifference, $Charge, $StereoParity) = unpack("A10A10A10xA3A2A3A3", $Line);
     726   }
     727   else {
     728     ($AtomX, $AtomY, $AtomZ, $AtomSymbol) = unpack("A10A10A10", $Line);
     729   }
     730   return ($AtomSymbol, $AtomX, $AtomY, $AtomZ, $MassDifference, $Charge, $StereoParity);
     731 }
     732
     733 # Map MDL charge value used in SD and MOL files to internal charge used by MayaChemTools.
     734 #
     735 sub MDLChargeToInternalCharge {
     736   my($MDLCharge) = @_;
     737   my($InternalCharge);
     738
     739   CHARGE: {
     740     if ($MDLCharge == 0) { $InternalCharge = 0; last CHARGE;}
     741     if ($MDLCharge == 1) { $InternalCharge = 3; last CHARGE;}
     742     if ($MDLCharge == 2) { $InternalCharge = 2; last CHARGE;}
     743     if ($MDLCharge == 3) { $InternalCharge = 1; last CHARGE;}
     744     if ($MDLCharge == 5) { $InternalCharge = -1; last CHARGE;}
     745     if ($MDLCharge == 6) { $InternalCharge = -2; last CHARGE;}
     746     if ($MDLCharge == 7) { $InternalCharge = -3; last CHARGE;}
     747     # All other MDL charge values, including 4 corresponding to "doublet radical",
     748     # are assigned internal value of 0.
     749     $InternalCharge = 0;
     750     if ($MDLCharge != 4) {
     751       carp "Warning: MDLChargeToInternalCharge: MDL charge value, $MDLCharge, is not supported: An internal charge value, 0, has been assigned...";
     752     }
     753   }
     754   return $InternalCharge;
     755 }
     756
     757 # Map internal charge used by MayaChemTools to MDL charge value used in SD and MOL files.
     758 #
     759 sub InternalChargeToMDLCharge {
     760   my($InternalCharge) = @_;
     761   my($MDLCharge);
     762
     763   CHARGE: {
     764     if ($InternalCharge == 3) { $MDLCharge = 1; last CHARGE;}
     765     if ($InternalCharge == 2) { $MDLCharge = 2; last CHARGE;}
     766     if ($InternalCharge == 1) { $MDLCharge = 3; last CHARGE;}
     767     if ($InternalCharge == -1) { $MDLCharge = 5; last CHARGE;}
     768     if ($InternalCharge == -2) { $MDLCharge = 6; last CHARGE;}
     769     if ($InternalCharge == -3) { $MDLCharge = 7; last CHARGE;}
     770     # All other MDL charge values, including 4 corresponding to "doublet radical",
     771     # are assigned internal value of 0.
     772     $MDLCharge = 0;
     773   }
     774   return $MDLCharge;
     775 }
     776
     777 # Ctab lines: Bond block
     778 #
     779 # Format: 111222tttsssxxxrrrccc
     780 #
     781 # 111: First atom number.
     782 # 222: Second atom number.
     783 # ttt: Bond type. 1 = Single, 2 = Double, 3 = Triple, 4 = Aromatic, 5 = Single or Double,
     784 #      6 = Single or Aromatic, 7 = Double or Aromatic, 8 = Any
     785 # sss: Bond stereo. Single bonds: 0 = not stereo, 1 = Up, 4 = Either, 6 = Down,
     786 #      Double bonds: 0 = Use x-, y-, z-coords from atom block to determine cis or trans,
     787 #      3 = Cis or trans (either) double bond
     788 # xxx: Not used
     789 # rrr: Bond topology. 0 = Either, 1 = Ring, 2 = Chain
     790 # ccc: Reacting center status. 0 = unmarked, 1 = a center, -1 = not a center,
     791 #      Additional: 2 = no change,4 = bond made/broken, 8 = bond order changes 12 = 4+8
     792 #      (both made/broken and changes); 5 = (4 + 1), 9 = (8 + 1), and 13 = (12 + 1) are also possible
     793 #
     794 sub ParseCmpdBondLine {
     795   my($Line) = @_;
     796   my($FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo);
     797
     798   ($FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo) = map {s/ //g; $_} unpack("A3A3A3A3", $Line);
     799   return ($FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo);
     800 }
     801
     802 # Map MDL bond type value used in SD and MOL files to internal bond order  and bond types
     803 # values used by MayaChemTools...
     804 #
     805 sub MDLBondTypeToInternalBondOrder {
     806   my($MDLBondType) = @_;
     807   my($InternalBondOrder, $InternalBondType);
     808
     809   $InternalBondType = '';
     810
     811   BONDTYPE: {
     812     if ($MDLBondType == 1) { $InternalBondOrder = 1; $InternalBondType = 'Single'; last BONDTYPE;}
     813     if ($MDLBondType == 2) { $InternalBondOrder = 2; $InternalBondType = 'Double'; last BONDTYPE;}
     814     if ($MDLBondType == 3) { $InternalBondOrder = 3; $InternalBondType = 'Triple'; last BONDTYPE;}
     815     if ($MDLBondType == 4) { $InternalBondOrder = 1.5; $InternalBondType = 'Aromatic'; last BONDTYPE;} # Aromatic
     816     if ($MDLBondType == 5) { $InternalBondOrder = 1; $InternalBondType = 'SingleOrDouble'; last BONDTYPE;} # Aromatic
     817     if ($MDLBondType == 6) { $InternalBondOrder = 1; $InternalBondType = 'SingleOrAromatic'; last BONDTYPE;} # Aromatic
     818     if ($MDLBondType == 7) { $InternalBondOrder = 2; $InternalBondType = 'DoubleOrAromatic'; last BONDTYPE;} # Aromatic
     819     if ($MDLBondType == 8) { $InternalBondOrder = 1; $InternalBondType = 'Any'; last BONDTYPE;} # Aromatic
     820     #
     821     # Although MDL aromatic bond values are used for query only and explicit Kekule bond order
     822     # values must be assigned, internal value of 1.5 is allowed to indicate aromatic bond orders.
     823     #
     824     # All other MDL bond type values -  5 = Single or Double, 6 = Single or Aromatic, 7 = Double or Aromatic,
     825     # 8 = Any - are also assigned appropriate internal value of 1: These are meant to be used for
     826     # structure queries by MDL products.
     827     #
     828     $InternalBondOrder = 1;
     829     $InternalBondType = 'Single';
     830
     831     carp "Warning: MDLBondTypeToInternalBondOrder: MDL bond type value, $MDLBondType, is not supported: An internal bond order value, 0, has been assigned...";
     832   }
     833   return ($InternalBondOrder, $InternalBondType);
     834 }
     835
     836 # Map internal bond order  and bond type values used by MayaChemTools to MDL bond type value used
     837 # in SD and MOL files...
     838 #
     839 sub InternalBondOrderToMDLBondType {
     840   my($InternalBondOrder, $InternalBondType) = @_;
     841   my($MDLBondType);
     842
     843   BONDTYPE: {
     844     if ($InternalBondOrder == 1) {
     845       if ($InternalBondType =~ /^SingleOrDouble$/i) {
     846         $MDLBondType = 5;
     847       }
     848       elsif ($InternalBondType =~ /^SingleOrAromatic$/i) {
     849         $MDLBondType = 6;
     850       }
     851       elsif ($InternalBondType =~ /^Any$/i) {
     852         $MDLBondType = 8;
     853       }
     854       else {
     855         $MDLBondType = 1;
     856       }
     857       $MDLBondType = 1;
     858       last BONDTYPE;
     859     }
     860     if ($InternalBondOrder == 2) {
     861       if ($InternalBondType =~ /^DoubleOrAromatic$/i) {
     862         $MDLBondType = 7;
     863       }
     864       else {
     865         $MDLBondType = 2;
     866       }
     867       last BONDTYPE;
     868     }
     869     if ($InternalBondOrder == 3) { $MDLBondType = 3; last BONDTYPE;}
     870     if ($InternalBondOrder == 1.5) { $MDLBondType = 4; last BONDTYPE;}
     871     if ($InternalBondType =~ /^Any$/i) { $MDLBondType = 8; last BONDTYPE;}
     872
     873     $MDLBondType = 1;
     874
     875     carp "Warning: InternalBondOrderToMDLBondType: Internal bond order and type values, $InternalBondOrder and $InternalBondType, don't match any valid MDL bond type: MDL bond type value, 1, has been assigned...";
     876   }
     877   return $MDLBondType;
     878 }
     879
     880 # Third line: Comments - A blank line is also allowed.
     881 sub ParseCmpdCommentsLine {
     882   my($Line) = @_;
     883   my($Comments);
     884
     885   $Comments = unpack("A80", $Line);
     886
     887   return ($Comments);
     888 }
     889
     890 # Map MDL bond stereo value used in SD and MOL files to internal bond stereochemistry values used by MayaChemTools...
     891 #
     892 sub MDLBondStereoToInternalBondStereochemistry {
     893   my($MDLBondStereo) = @_;
     894   my($InternalBondStereo);
     895
     896   $InternalBondStereo = '';
     897
     898   BONDSTEREO: {
     899     if ($MDLBondStereo == 1) { $InternalBondStereo = 'Up'; last BONDSTEREO;}
     900     if ($MDLBondStereo == 4) { $InternalBondStereo = 'UpOrDown'; last BONDSTEREO;}
     901     if ($MDLBondStereo == 6) { $InternalBondStereo = 'Down'; last BONDSTEREO;}
     902     if ($MDLBondStereo == 3) { $InternalBondStereo = 'CisOrTrans'; last BONDSTEREO;}
     903     if ($MDLBondStereo == 0) { $InternalBondStereo = 'None'; last BONDSTEREO;}
     904
     905     $InternalBondStereo = '';
     906     carp "Warning: MDLBondStereoToInternalBondType: MDL bond stereo value, $MDLBondStereo, is not supported: It has been ignored and bond order would be used to determine bond type...";
     907   }
     908   return $InternalBondStereo;
     909 }
     910
     911 # Map internal bond stereochemistry values used by MayaChemTools to MDL bond stereo value used in SD and MOL files...
     912 #
     913 sub InternalBondStereochemistryToMDLBondStereo {
     914   my($InternalBondStereo) = @_;
     915   my($MDLBondStereo);
     916
     917   $MDLBondStereo = 0;
     918
     919   BONDSTEREO: {
     920     if ($InternalBondStereo =~ /^Up$/i) { $MDLBondStereo = 1; last BONDSTEREO;}
     921     if ($InternalBondStereo =~ /^UpOrDown$/i) { $MDLBondStereo = 4; last BONDSTEREO;}
     922     if ($InternalBondStereo =~ /^Down$/) { $MDLBondStereo = 6; last BONDSTEREO;}
     923     if ($InternalBondStereo =~ /^CisOrTrans$/) { $MDLBondStereo = 3; last BONDSTEREO;}
     924
     925     $MDLBondStereo = 0;
     926   }
     927   return $MDLBondStereo;
     928 }
     929
     930 # Fourth line: Counts
     931 #
     932 # Format: aaabbblllfffcccsssxxxrrrpppiiimmmvvvvvv
     933 #
     934 # aaa: number of atoms; bbb: number of bonds; lll: number of atom lists; fff: (obsolete)
     935 # ccc: chiral flag: 0=not chiral, 1=chiral; sss: number of stext entries; xxx,rrr,ppp,iii:
     936 # (obsolete); mmm: number of lines of additional properties, including the M END line, No
     937 # longer supported, default is set to 999; vvvvvv: version
     938
     939 sub ParseCmpdCountsLine {
     940   my($Line) = @_;
     941   my($AtomCount, $BondCount, $ChiralFlag, $PropertyCount, $Version);
     942
     943   if (length($Line) >= 39) {
     944     ($AtomCount, $BondCount, $ChiralFlag, $PropertyCount, $Version) = unpack("A3A3x3x3A3x3x3x3x3x3A3A6", $Line);
     945   }
     946   elsif (length($Line) >= 15) {
     947     ($PropertyCount, $Version) = ("999", "v2000");
     948     ($AtomCount, $BondCount, $ChiralFlag) = unpack("A3A3x3x3A3", $Line);
     949   }
     950   else {
     951     ($ChiralFlag, $PropertyCount, $Version) = ("0", "999", "v2000");
     952     ($AtomCount, $BondCount) = unpack("A3A3", $Line);
     953   }
     954
     955   if ($Version =~ /V3000/i) {
     956     # Current version of MayaChemTools modules and classes for processing MDL MOL and SD don't support
     957     # V3000. So instead of relying on callers, just exit with an error to disable any processing of V3000
     958     # format.
     959     croak "Error: SDFileUtil::ParseCmpdCountsLine: The Extended Connection Table (V3000) format in MDL MOL and SD files is not supported by the current release of MayaChemTools...";
     960   }
     961
     962   return ($AtomCount, $BondCount, $ChiralFlag, $PropertyCount, $Version);
     963 }
     964
     965 # Second line: Misc info
     966 #
     967 # Format: IIPPPPPPPPMMDDYYHHmmddSSssssssssssEEEEEEEEEEEERRRRRR
     968 #         A2A8      A10       A2I2A10       A12         A6
     969 # User's first and last initials (I), program name (P), date/time (M/D/Y,H:m),
     970 # dimensional codes - 2D or 3D (d),scaling factors (S, s), energy (E) if modeling program input,
     971 # internal registry number (R) if input through MDL form. A blank line is also allowed.
     972 sub ParseCmpdMiscInfoLine {
     973   my($Line) = @_;
     974   my($UserInitial, $ProgramName, $Date, $Code, $ScalingFactor1, $ScalingFactor2, $Energy, $RegistryNum);
     975
     976   ($UserInitial, $ProgramName, $Date, $Code, $ScalingFactor1, $ScalingFactor2, $Energy, $RegistryNum) = unpack("A2A8A10A2A2A10A12A6", $Line);
     977   return ($UserInitial, $ProgramName, $Date, $Code, $ScalingFactor1, $ScalingFactor2, $Energy, $RegistryNum);
     978 }
     979
     980 # First line: Molecule name. This line is unformatted, but like all other lines in a
     981 # molfile may not extend beyond column 80. A blank line is also allowed.
     982 sub ParseCmpdMolNameLine {
     983   my($Line) = @_;
     984   my($MolName);
     985
     986   $MolName = unpack("A80", $Line);
     987
     988   return ($MolName);
     989 }
     990
     991 # Parse atom alias property line in CTAB generic properties block.
     992 #
     993 # Atom alias property line format:
     994 #
     995 # A  aaa
     996 # x...
     997 #
     998 #    aaa: Atom number
     999 #    x: Atom alias in next line
    1000 #
    1001 sub ParseCmpdAtomAliasPropertyLine {
    1002   my($Line, $NextLine) = @_;
    1003   my($Label, $AtomNumber, $AtomAlias);
    1004
    1005   ($Label, $AtomNumber) = split(' ', $Line);
    1006   $AtomAlias = $NextLine;
    1007
    1008   if (!$AtomAlias) {
    1009     carp "Warning: _ParseCmpdAtomAliasPropertyLine: No atom alias value specified on the line following atom alias property line...";
    1010   }
    1011
    1012   return ($AtomNumber, $AtomAlias);
    1013 }
    1014
    1015 # Parse charge property line in CTAB generic properties block.
    1016 #
    1017 # Charge property line format:
    1018 #
    1019 # M  CHGnn8 aaa vvv ...
    1020 #
    1021 #    nn8: Number of value pairs. Maximum of 8 pairs allowed.
    1022 #    aaa: Atom number
    1023 #    vvv: -15 to +15. Default of 0 = uncharged atom. When present, this property supersedes
    1024 #    all charge and radical values in the atom block, forcing a 0 charge on all atoms not
    1025 #    listed in an M  CHG or M  RAD line.
    1026 #
    1027 sub ParseCmpdChargePropertyLine {
    1028   my($Line) = @_;
    1029
    1030   return _ParseCmpdGenericPropertyLine('Charge', $Line);
    1031 }
    1032
    1033
    1034 # Parse isotope property line in CTAB generic properties block.
    1035 #
    1036 # Isoptope property line format:
    1037 #
    1038 # M  ISOnn8 aaa vvv ...
    1039 #
    1040 #    nn8: Number of value paris. Maximum of 8 pairs allowed.
    1041 #    aaa: Atom number
    1042 #    vvv: Absolute mass of the atom isotope as a positive integer. When present, this property
    1043 #    supersedes all isotope values in the atom block. Default (no entry) means natural
    1044 #    abundance. The difference between this absolute mass value and the natural
    1045 #    abundance value specified in the PTABLE.DAT file must be within the range of -18
    1046 #    to +12
    1047 #
    1048 # Notes:
    1049 #  . Values correspond to mass numbers...
    1050 #
    1051 sub ParseCmpdIsotopePropertyLine {
    1052   my($Line) = @_;
    1053
    1054   return _ParseCmpdGenericPropertyLine('Isotope', $Line);
    1055 }
    1056
    1057 # Parse radical property line in CTAB generic properties block.
    1058 #
    1059 # Radical property line format:
    1060 #
    1061 # M  RADnn8 aaa vvv ...
    1062 #
    1063 #    nn8: Number of value paris. Maximum of 8 pairs allowed.
    1064 #    aaa: Atom number
    1065 #    vvv: Default of 0 = no radical, 1 = singlet, 2 = doublet, 3 = triplet . When
    1066 #    present, this property supersedes all charge and radical values in the atom block,
    1067 #    forcing a 0 (zero) charge and radical on all atoms not listed in an M  CHG or
    1068 #    M  RAD line.
    1069 #
    1070 sub ParseCmpdRadicalPropertyLine {
    1071   my($Line) = @_;
    1072
    1073   return _ParseCmpdGenericPropertyLine('Radical', $Line);
    1074 }
    1075
    1076 # Map MDL radical stereo value used in SD and MOL files to internal spin multiplicity values used by MayaChemTools...
    1077 #
    1078 sub MDLRadicalToInternalSpinMultiplicity {
    1079   my($MDLRadical) = @_;
    1080   my($InternalSpinMultiplicity);
    1081
    1082   $InternalSpinMultiplicity = '';
    1083
    1084   SPINMULTIPLICITY: {
    1085     if ($MDLRadical == 0) { $InternalSpinMultiplicity = 0; last SPINMULTIPLICITY;}
    1086     if ($MDLRadical == 1) { $InternalSpinMultiplicity = 1; last SPINMULTIPLICITY;}
    1087     if ($MDLRadical == 2) { $InternalSpinMultiplicity = 2; last SPINMULTIPLICITY;}
    1088     if ($MDLRadical == 3) { $InternalSpinMultiplicity = 3; last SPINMULTIPLICITY;}
    1089     $InternalSpinMultiplicity = '';
    1090     carp "Warning: MDLRadicalToInternalSpinMultiplicity: MDL radical value, $MDLRadical, specifed on line M  RAD is not supported...";
    1091   }
    1092   return $InternalSpinMultiplicity;
    1093 }
    1094
    1095 # Map internal spin multiplicity values used by MayaChemTools to MDL radical stereo value used in SD and MOL files...
    1096 #
    1097 sub InternalSpinMultiplicityToMDLRadical {
    1098   my($InternalSpinMultiplicity) = @_;
    1099   my($MDLRadical);
    1100
    1101   $MDLRadical = 0;
    1102
    1103   SPINMULTIPLICITY: {
    1104     if ($InternalSpinMultiplicity == 1) { $MDLRadical = 1; last SPINMULTIPLICITY;}
    1105     if ($InternalSpinMultiplicity == 2) { $MDLRadical = 2; last SPINMULTIPLICITY;}
    1106     if ($InternalSpinMultiplicity == 3) { $MDLRadical = 3; last SPINMULTIPLICITY;}
    1107     $MDLRadical = 0;
    1108   }
    1109   return $MDLRadical;
    1110 }
    1111
    1112 # Process generic CTAB property line...
    1113 sub _ParseCmpdGenericPropertyLine {
    1114   my($PropertyName, $Line) = @_;
    1115
    1116   my($Label, $PropertyLabel, $ValuesCount, $ValuePairsCount, @ValuePairs);
    1117
    1118   @ValuePairs = ();
    1119   ($Label, $PropertyLabel, $ValuesCount, @ValuePairs) = split(' ', $Line);
    1120   $ValuePairsCount = (scalar @ValuePairs)/2;
    1121   if ($ValuesCount != $ValuePairsCount) {
    1122     carp "Warning: _ParseCmpdGenericPropertyLine: Number of atom number and $PropertyName value paris specified on $Label $PropertyLabel property line, $ValuePairsCount, does not match expected value of $ValuesCount...";
    1123   }
    1124
    1125   return (@ValuePairs);
    1126 }
    1127
    1128 # Generic CTAB property lines for charge, istope and radical properties...
    1129 #
    1130 sub _GenerateCmpdGenericPropertyLines {
    1131   my($PropertyName, $PropertyValuePairsRef) = @_;
    1132   my($Index, $PropertyLabel, $Line, $PropertyCount, $AtomNum, $PropertyValue, @PropertyLines);
    1133
    1134   @PropertyLines = ();
    1135   NAME: {
    1136     if ($PropertyName =~ /^Charge$/i) { $PropertyLabel = "M  CHG"; last NAME; }
    1137     if ($PropertyName =~ /^Isotope$/i) { $PropertyLabel = "M  ISO"; last NAME; }
    1138     if ($PropertyName =~ /^Radical$/i) { $PropertyLabel = "M  RAD"; last NAME; }
    1139     carp "Warning: _GenerateCmpdGenericPropertyLines: Unknown property name, $PropertyName, specified...";
    1140     return @PropertyLines;
    1141   }
    1142
    1143   # A maximum of 8 property pair values allowed per line...
    1144   $PropertyCount = 0;
    1145   $Line = '';
    1146   for ($Index = 0; $Index < $#{$PropertyValuePairsRef}; $Index += 2) {
    1147     if ($PropertyCount > 8) {
    1148       # Setup property line...
    1149       $Line = "${PropertyLabel}  8${Line}";
    1150       push @PropertyLines, $Line;
    1151
    1152       $PropertyCount = 0;
    1153       $Line = '';
    1154     }
    1155     $PropertyCount++;
    1156     $AtomNum = $PropertyValuePairsRef->[$Index];
    1157     $PropertyValue = $PropertyValuePairsRef->[$Index + 1];
    1158     $Line .= sprintf " %3i %3i", $AtomNum, $PropertyValue;
    1159   }
    1160   if ($Line) {
    1161     $Line = "${PropertyLabel}  ${PropertyCount}${Line}";
    1162     push @PropertyLines, $Line;
    1163   }
    1164   return @PropertyLines;
    1165 }
    1166
    1167 #
    1168 # Read compound data into a string and return its value
    1169 sub ReadCmpdString {
    1170   my($SDFileRef) = @_;
    1171   my($CmpdString);
    1172
    1173   $CmpdString = "";
    1174   LINE: while (defined($_ = <$SDFileRef>)) {
    1175     # Change Windows and Mac new line char to UNIX...
    1176     s/(\r\n)|(\r)/\n/g;
    1177
    1178     if (/^\$\$\$\$/) {
    1179       # Take out any new line char at the end by explicitly removing it instead of using
    1180       # chomp, which might not always work correctly on files generated on a system
    1181       # with a value of input line separator different from the current system...
    1182       s/\n$//g;
    1183
    1184       # Doesn't hurt to chomp...
    1185       chomp;
    1186
    1187       $CmpdString .=  $_;
    1188       last LINE;
    1189     }
    1190     else {
    1191       $CmpdString .=  $_;
    1192     }
    1193   }
    1194   return $CmpdString;
    1195 }
    1196
    1197 # Find out the number of fragements in the compounds. And for the compound with
    1198 # more than one fragment, remove all the others besides the largest one.
    1199 sub WashCmpd {
    1200   my($CmpdLines) = @_;
    1201   my($WashedCmpdString, $FragmentCount, $Fragments);
    1202
    1203   $WashedCmpdString = "";
    1204   ($FragmentCount, $Fragments) = GetCmpdFragments($CmpdLines);
    1205   if ($FragmentCount > 1) {
    1206     # Go over the compound data for the largest fragment including property
    1207     # data...
    1208     my (@AllFragments, @LargestFragment, %LargestFragmentAtoms, @WashedCmpdLines, $Index, $LineIndex, $AtomCount, $BondCount, $NewAtomCount, $NewBondCount, $FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo, $FirstNewAtomNum, $SecondNewAtomNum, $AtomNum, $ChiralFlag, $BondLine, $MENDLineIndex, $Line, $Value, @ValuePairs, @NewValuePairs, $ValuePairIndex, $NewAtomNum, @NewPropertyLines);
    1209
    1210     @AllFragments = (); @LargestFragment = ();
    1211     %LargestFragmentAtoms = ();
    1212     @AllFragments = split "\n", $Fragments;
    1213     @LargestFragment = split " ", $AllFragments[0];
    1214     for $Index (0 .. $#LargestFragment) {
    1215       # Map old atom numbers to new atom numbers as the fragment atom numbers are sorted
    1216       # from lowest to highest old atom numbers...
    1217       $LargestFragmentAtoms{$LargestFragment[$Index]} = $Index + 1;
    1218     }
    1219     @WashedCmpdLines = ();
    1220     push @WashedCmpdLines, @$CmpdLines[0], @$CmpdLines[1], @$CmpdLines[2], @$CmpdLines[3];
    1221     ($AtomCount, $BondCount, $ChiralFlag) = ParseCmpdCountsLine(@$CmpdLines[3]);
    1222     $NewAtomCount = @LargestFragment;
    1223     $NewBondCount = 0;
    1224     $AtomNum = 0;
    1225     # Retrieve the largest fragment atom lines...
    1226     for ($LineIndex = 4; $LineIndex < (4 + $AtomCount); $LineIndex++) {
    1227       $AtomNum++;
    1228       if ($LargestFragmentAtoms{$AtomNum}) {
    1229         push @WashedCmpdLines, @$CmpdLines[$LineIndex];
    1230       }
    1231     }
    1232     # Retrieve the largest fragment bond lines...
    1233     for ($LineIndex = 4 + $AtomCount; $LineIndex < (4 + $AtomCount + $BondCount); $LineIndex++) {
    1234       ($FirstAtomNum, $SecondAtomNum, $BondType, $BondStereo) = ParseCmpdBondLine(@$CmpdLines[$LineIndex]);
    1235       if ($LargestFragmentAtoms{$FirstAtomNum} && $LargestFragmentAtoms{$SecondAtomNum}) {
    1236         $NewBondCount++;
    1237         # Set up bond line with new atom number mapping...
    1238         $FirstNewAtomNum =  $LargestFragmentAtoms{$FirstAtomNum};
    1239         $SecondNewAtomNum =  $LargestFragmentAtoms{$SecondAtomNum};
    1240         $BondLine = GenerateCmpdBondLine($FirstNewAtomNum, $SecondNewAtomNum, $BondType, $BondStereo);
    1241         push @WashedCmpdLines, $BondLine;
    1242       }
    1243     }
    1244     # Get property lines for CHG, ISO and RAD label and map the old atom numbers to new
    1245     # atom numners; Others, property lines before M  END line are skipped as atom numbers for
    1246     # other properties might not valid anymore...
    1247     #
    1248     $MENDLineIndex = $LineIndex;
    1249     LINE: for ($LineIndex = (4 + $AtomCount + $BondCount); $LineIndex < @$CmpdLines; $LineIndex++) {
    1250       $Line = @$CmpdLines[$LineIndex];
    1251       if ($Line =~ /^M  END/i) {
    1252         push @WashedCmpdLines, "M  END";
    1253         $MENDLineIndex = $LineIndex;
    1254         last LINE;
    1255       }
    1256
    1257       @ValuePairs = ();
    1258       if ($Line =~ /^M  CHG/i) {
    1259         @ValuePairs = ParseCmpdChargePropertyLine($Line);
    1260       }
    1261       elsif ($Line =~ /^M  RAD/i) {
    1262         @ValuePairs = ParseCmpdRadicalPropertyLine($Line);
    1263       }
    1264       elsif ($Line =~ /^M  ISO/i) {
    1265         @ValuePairs = ParseCmpdIsotopePropertyLine($Line);
    1266       }
    1267       elsif ($Line =~ /^A  /i) {
    1268         my($NextLine);
    1269         $LineIndex++;
    1270         $NextLine = @$CmpdLines[$LineIndex];
    1271         @ValuePairs = ParseCmpdAtomAliasPropertyLine($Line, $NextLine);
    1272       }
    1273       else {
    1274         next LINE;
    1275       }
    1276
    1277       if (!@ValuePairs) {
    1278         next LINE;
    1279       }
    1280
    1281       # Collect values for valid atom numbers with mapping to new atom numbers...
    1282       @NewValuePairs = ();
    1283       VALUEINDEX: for ($ValuePairIndex = 0; $ValuePairIndex < $#ValuePairs; $ValuePairIndex += 2) {
    1284         $AtomNum = $ValuePairs[$ValuePairIndex]; $Value = $ValuePairs[$ValuePairIndex + 1];
    1285         if (!exists $LargestFragmentAtoms{$AtomNum}) {
    1286           next VALUEINDEX;
    1287         }
    1288         $NewAtomNum = $LargestFragmentAtoms{$AtomNum};
    1289         push @NewValuePairs, ($NewAtomNum, $Value)
    1290       }
    1291       if (!@NewValuePairs) {
    1292         next LINE;
    1293       }
    1294       @NewPropertyLines = ();
    1295       if ($Line =~ /^M  CHG/i) {
    1296         @NewPropertyLines = GenerateCmpdChargePropertyLines(\@NewValuePairs);
    1297       }
    1298       elsif ($Line =~ /^M  RAD/i) {
    1299         @NewPropertyLines = GenerateCmpdRadicalPropertyLines(\@NewValuePairs);
    1300       }
    1301       elsif ($Line =~ /^M  ISO/i) {
    1302         @NewPropertyLines = GenerateCmpdIsotopePropertyLines(\@NewValuePairs);
    1303       }
    1304       elsif ($Line =~ /^A  /i) {
    1305         @NewPropertyLines = GenerateCmpdAtomAliasPropertyLines(\@NewValuePairs);
    1306       }
    1307       push @WashedCmpdLines, @NewPropertyLines;
    1308     }
    1309
    1310     # Retrieve rest of the data label and value property data...
    1311     for ($LineIndex = (1 + $MENDLineIndex); $LineIndex < @$CmpdLines; $LineIndex++) {
    1312       push @WashedCmpdLines, @$CmpdLines[$LineIndex];
    1313     }
    1314     # Update atom and bond count line...
    1315     $WashedCmpdLines[3] = GenerateCmpdCountsLine($NewAtomCount, $NewBondCount, $ChiralFlag);
    1316
    1317     $WashedCmpdString = join "\n", @WashedCmpdLines;
    1318   }
    1319   return ($FragmentCount, $Fragments, $WashedCmpdString);
    1320 }
    1321