Representing Molecules in the Wolfram Language
​
​[This notebook is a record of my first practice with graph-related and molecule-related functions in the Wolfram language. The list of chemical reaction functions at the end of this notebook are continued and improved in my second notebook “Organic Synthesis Project”. --Leonardo Cabana]​
​​
​Modeling Organic Chemistry Synthesis as Graphs in the Wolfram Language
This will be my first project using the Wolfram language (WL). First, look at the various ways to represent simple organic molecules as graphs in WL.
I can create the graph of a molecule from scratch, here for formaldehyde:
GraphPlot3D[{12,12,13,14},VertexLabels{1"C",2"O",3"H",4"H"}]
Out[]=
But many molecules are already directly available in WL, and the information comprising them is already formatted for processing using other WL functions:
MoleculePlot[Molecule["Formaldehyde"]]
Out[]=
BondList[Molecule["Formaldehyde"]]
Out[]=
{Bond[{1,2},Double],Bond[{1,3},Single],Bond[{1,4},Single]}
​
Define the molecule for formaldehyde -- does my defined instance perform as well as the “built-in” molecule?
formald=Molecule[{"C","O","H","H"},{Bond[{1,2},"Double"],Bond[{1,3}],Bond[{1,4}]}];​​MoleculePlot[formald]​​
Out[]=
​
I didn't get "Graph Minors" function to work.
Meanwhile, try to define a rule for Hydrogenation of Aldehydes:
First, add singly-bonded H atom to the C and O atoms, then delete one of the C=O bonds:
q=EdgeAdd[{12,23,24,24},{25,46}]
Out[]=
r=EdgeDelete[q,{24}]
Out[]=
GraphPlot3D[r,VertexLabels{2"C",4"O",1"H",3"H",5"H",6"H"}]
Out[]=
This is the correct connectivity, but the bond angles are incorrect.
Next try applying this rule to a more elaborate aldehyde substrate in which the aldehyde functional group is embedded.
acetald=GraphPlot3D[{12,13,13,14,25,26,27},VertexLabels{1"C",2"C",3"O",4"H",5"H",6"H",7"H"}]
Out[]=
Now apply the rule (after to taking into account that this is an aldehyde that has a C substituent.
It's better to delete old edges before adding new ones.
s=EdgeDelete[{12,13,13,14,25,26,27},{13}]
Out[]=
t=EdgeAdd[s,{18,39}]
Out[]=
GraphPlot3D[t,VertexLabels{1"C",2"C",3"O",4"H",5"H",6"H",7"H",8"H",9"H"}]
Out[]=
​
This is correct connectivity, but I used the same atom-numberings in the starting graph and rule. Will the aldehyde functional group connectivity pattern still be recognized if the numbering is different?
acetaldAltNumb=Graph[{3132,3133,3133,3134,3235,3236,3237}]
Out[]=
EdgeDelete[{3132,3133,3133,3134,3235,3236,3237},{13}]
EdgeDelete
:The argument {13} in
EdgeDelete[{31  32, 31  33, 31  33, 31  34, 32  35, 32  36, 32  37}, {1  3}]
is not a valid edge.
Out[]=
EdgeDelete[{3132,3133,3133,3134,3235,3236,3237},{13}]
The error message immediately above indicates that my above attempts are too crude, and cannot recognize a pattern independent of node labels

​
​Built-In Molecule Functions in the Wolfram Language

​
I will now try out the new experimental functionality of MoleculeGraph, starting from a slightly more complex molecule, caffeine:
​
caf=Molecule["caffeine"]
Out[]=
Molecule
Formula:
C
8
H
10
N
4
O
2
Atoms: 24 Bonds: 25

In[]:=
MoleculeGraph[caf]
​​
In[]:=
VertexDegree[MoleculeGraph[caf]]
{4,3,3,1,3,3,2,3,3,4,3,4,3,1,1,1,1,1,1,1,1,1,1,1}​​
Another built-in molecule:

​
​Finding Substructures within a Molecule

See FindMoleculeSubstructure, MoleculePattern, MoleculeContainsQ.
​
Try to find a alkene C=CH2 group within vinylbenzene:
​
Success!
Find the same pattern in a di-ene to see if internal and terminal double-bonds can be distinguished:
​
Success again -- and there was no false-positive caused by the internal C=C.
​
​
Now try finding carbonyl C=O patterns within molecules:
​
​
​
Here is another molecule (from the WL documentation), in which a substructure search is performed for the S=O bond pattern.
​
Success!

​
​Representing Chemical Reactions as Transformations of Graphs

​
Now that substructures can be located within a molecule, and the corresponding indices recovered, try a transformation (chemical reaction) on this organic molecule. For example, catalytic hydrogenation, addition of H2 across the (non-aromatic) sulfur-oxygen double-bond (S=O).
​
​
Putting all the steps together:
​
​
Here is the transformation for catalytic hydrogenation of carbonyl C=O groups:
What happens if the starting molecule (reactant) contains multiple instances of the functional group that is undergoing reaction?
​
Only one of the C=O groups was reduced, but we want all of them reduced). In the C=O pattern, C is labelled 1, O is labeled 2, but substrate C=O numbering is reversed.
The default FindMoleculeSubstructure function stops looking after finding one pattern match. Try changing that parameter to “All”:
​
The above crude approach failed, so now I’m looking into how to address each association within the list of associations:
Success. This is just a 2-dimensional array, where the array keys are listed on the right. The 1st array key specifies the association, and the 2nd array key specifies the association key.
Here is a work-around for a set containing 2 associations:
​
​
The above successfully handled 2 and 3 instances of the target functional group pattern. Now generalize this approach to an arbitrary number of associations (functional group matches).
​
Next, modify my existing reaction function definitions to automatically delete inorganic byproducts, showing only the organic products (or for the unchanged reactant if the reactant lacked the target pattern).
​
Next, Anti-Markovnikoff Hydrobromination of Alkenes (REACTION ?) (taking into account regiochemistry; stereochemistry is ?).
​
Next, Markovnikoff Hydration of Alkenes (REACTION ?) (taking into account regiochemistry; stereochemistry is racemic).
​
Next, Anti-Markovnikoff Hydration of Alkenes (REACTION ?) (taking into account regiochemistry; stereochemistry is ?).
​
Next, Zaitsev Elimination on Alkyl Bromides (REACTION ?) (taking into account regiochemistry, but ignoring E/Z stereochemistry of the alkene product).
​
Next, Hofmann Elimination on Alkyl Bromides (REACTION ?) (taking into account regiochemistry, but ignoring E/Z stereochemistry of the alkene product).
​
​

​
​Loading Multiple Reactant Molecule Patterns onto the Same Graph

For chemical reactions that create new C-C or C-hetero-atom bonds (hetero-atoms are N, O, S, P etc), the patterns of both reactants should be present on the same graph (but with no overlap before reacting). Since each reactant could be arbitrarily complex -- i.e., not a simple standard reagent. But how to get two (or more) molecules on the same graph?
​
​
This almost works, but the atom indices (inside of the Bond[] function) on one of the two reactant molecules graphs must be changed to avoid overlap between the the indices of the two molecules.
​
This works, after I manually renumbered the indices (inside Bond[]) of one of the molecules to avoid overlap between the two molecules.
​
​
​? But how do I automate the renumbering of indices inside of Bond[] ?
​
But each of these types or reactions involving two or more complex (“non-standard”) reactants should have a defined function with an argument for each reactant:
​
​
​
Now exploring how to deal with the orbital hybridization of carbon atoms:
​
Success, but this needs to be incorporated into a defined function.
​
Success -- the above is his is REACTION 4 (SUPERCEDED) (but passing the other reactant).​
​
​
Next: create function for Grignard Reaction with Aldehydes and Ketones. The Grignard reagent is derived directly from an organic halides (C-X, where C is any SP, SP2, or SP3 hybridized C). Later: grignard reaction CO2, esters, and epoxides.
Success!! ... This is REACTION 7 (SUPERCEDED) Grignard reagent (from organic bromide) reacts with aldehyde/ketone substrate. This reaction handles any aldehyde (including formaldehyde) or ketone. Next I will added third argument for this function: default setting passes aldehyde/ketone through the function if nor reaction occurred (because at least one required molecular pattern was missing on the reactants).

​
​Organic Chemical Reaction Functions Completed So Far​
​

Reaction 4: Alkylation of Alkynes [SUPERCEDED -- see 2nd notebook]
The function above accepts either a primary halide or methyl halide. Don’t include option for intramolecular reaction.
​
​
​
​Reaction 7: Grignard Reaction of Aldehydes and Ketones [SUPERCEDED -- see 2nd notebook]
Don’t include option for intramolecular reaction.

​
Reactions that Form C-O, C-N, and C-S Bonds
​
Reaction 8: Williamson Ether Synthesis (Alkylation of Alcohols) [SUPERCEDED -- see 2nd notebook]
​
​
Reaction 9: Gabriel Synthesis of Primary Amines from Primary Alkyl Halides [SUPERCEDED -- see 2nd notebook]
**The above halide can be either a primary or methyl halide.
​
​
Reaction 10: Fischer Ester Synthesis from Alcohols and Carboxylic Acids ****UNDER CONSTRUCTION**** [SUPERCEDED -- see 2nd notebook]
​
​Reaction 1: H2 Reduction of Aldehydes/Ketones to Alcohols [SUPERCEDED -- see 2nd notebook]
​
​
​Reaction 2: Bromination of Alcohols [SUPERCEDED -- see 2nd notebook]
​
​
​Reaction 3: Oxidation of Alcohols to Aldehyde/Ketones [SUPERCEDED -- see 2nd notebook]
​
​
​Reaction 5: H2 Reduction of Alkyne to Alkene [SUPERCEDED -- see 2nd notebook]
​
​
​Reaction 6: H2 Reduction of Alkenes and Alkynes to Alkanes [SUPERCEDED -- see 2nd notebook]
​
​
​Reaction 11: Alkene Hydration (to Alcohol) [SUPERCEDED -- see 2nd notebook]
​
​
Reaction 12: Alkene HydroHalogenation (to Alkyl Halide) [SUPERCEDED -- see 2nd notebook]
​
​
END OF THIS NOTEBOOK. This research is continued in my second notebook, “Organic Synthesis Project”.