CITE THIS NOTEBOOK: Isomer generation using SURGE based paclet by Theodore Mollano. Wolfram Community DEC 15 2022.
If you were given the chemical formula
C
3
H
4
, how many molecules do you know of that fit this description? This kind of constraint problem is known as an isomer generation problem, and has been solved repeatedly since the 1960s. Today, one of the most efficient isomer generation algorithms is known as SURGE, and is based on the lightning-fast nauty graph isomorphism package in C. There is a new Wolfram Language paclet that uses SURGE to generate isomers. ​
​Solving the isomer generation problem using a few lines of Wolfram Language code is a great exercise to demonstrate the computational power per line of the Wolfram Language. It is also a great introduction to the way of thinking that powers the rest of the Wolfram Language Chemistry suite. We can solve the isomer generation problem in just 16 lines of code with WL.

Generating a list of Isomers

Background

Suppose we are given a chemical formula in the following form. What chemical isomers does it represent?
In[]:=
{{"C",3},{"H",4}}(*3carbons,and4hydrogens*)
Out[]=
{{C,3},{H,4}}
Here is an example output we might expect:
Out[]=
Isomers of:
C
3
H
4
This panel displays Cyclopropene, Propyne, and Propa-1,2-diene that represent isomers of
C
3
H
4
. By convention, non-important carbon atoms are represented as empty nodes. Single lines represent single bonds between atoms, and double lines between atoms represent double bonds between a pair of atoms, etc. ​​
Enumerating each column from left to right on the periodic table (without counting light green columns), we get the group numbers of each element. For example, H has group number 1, and He has group number 8.
To generate a list of isomers we need a way to connect atoms with pairwise bonds, while being able to check for each structure we build, that structure is indeed a valid molecule. Recall that a single bond between atom A and atom B requires each atom to donate a single electron to the bond, profiting each atom with +1 electrons.
According to Lewis Bonding Theory, the most important predictor of whether an atom will form bonds with other atoms is if that atom can with more bonds populate its outer-electron shell (valence shell). Each atom typically has
k
electrons in its valence shell in a neutral state, where
k
is its group number. To achieve stability, each atom needs 2, 8, or in general
2
2
n
electrons in its valence shell to be energetically stable. Atoms in row one need 2 electrons in their valence shell, while atoms in rows two and three need 8 electrons in their valence shell.
Let us take cyclopropene, for instance, and let us investigate the valence shells around each atom. We can do this by turning cyclopropene into a graph in the following way:
In[]:=
cyclopropeneGraph=Graph[{1,2,3,4,5,6,7},{12,12,23,31,14,25,36,37},VertexLabels->{1->"C (1)",2->"C (2)",3->"C (3)",4->"H (4)",5->"H (5)",6->"H (6)",7->"H (7)"}]
Out[]=
Great! Now we can see that each atom is represented by a node, with edges connecting a pair of nodes if that pair has a bond between them, with multiplicity. The rightmost upper carbon in the cyclopropene makes four bonds, two to one carbon, and one to a different carbon, and one to a hydrogen. Given carbon is in group 4, it has 4 valence electrons, and since it makes 4 single bonds, it gains +4 valence electrons. Therefore, carbon has
4+1*48
valence electrons, as necessary.
We just need to make graphs such that the vertex degree around each node gets us the number of electrons the corresponding atom needs.
Given a list of atoms how many graphs g satisfy the following: (1) a graph has vertex list in a bijection with the list of atoms (2) any given vertex in the graph has degree
8-k
if the vertex corresponds to an atom in rows 2 or 3 of the Periodic Table, where
k
is the group number of that atom. Any given vertex may have a degree of 1 or 0 if the vertex corresponds to an atom of hydrogen or helium, respectively.
We can visualise this a bit better with an adjacency matrix for each graph, where (i,j) = #edges between atom i and atom j
In[]:=
AdjacencyMatrix[cyclopropeneGraph]//MatrixForm
Out[]//MatrixForm=
AdjacencyMatrix[cyclopropeneGraph]
For example, atom #1, carbon, forms a single bond with atom #3, carbon, and atom #1 forms two bonds with atom #2, another carbon. Atom #1, also forms a single bond with atom #4, hydrogen. We can that the sum of each row in the adjacency matrix corresponding to a vertex is equal to the degree of the vertex.
Problem solved! We just construct graphs which represent atoms, and we check that for each atom, the sum of each row in the adjacency matrix corresponding to that atom, is equal to the number of bonds that atom needs to fill it’s valence.
This corresponds to finding
A
in the linear system
A

x
b
, where

x

1
⋮

,
b

v
i
⋮

, where
v
i
is the number of bonds necessary to achieve the full valence for each atom.
A
takes the form
0
y
z
…
y
0
c
…
z
c
0
…
⋮
⋮
⋮
0
, where all of its entries are in the range
{0,1,2,3,4}
. Because
A
is N
×
N, we know that
A
admits a finite number of solutions. One way of understanding
A
is that we are performing stars and bars at each row. For example, for the first atom’s row, we are partitioning the number of necessary bonds
v
among variables
x,y,z…
. From the theory of generating functions, the number of ways we can do this is the
v
x
coefficient of
N-1
(1+x+
2
x
+
3
x
+
4
x
)
, for the
N
atom scenario. A weak upper bound on the number of solutions to the system would then be
∏
v
[
v
x
]
N-1
(1+x+
2
x
+
3
x
+
4
x
)
for each
v
.
The stars and bars approach shows that there is a lot of redundancy in this problem. We can effectively permute row
i↔j
and column
i↔j
so long as
i
and
j
represent atoms that need the same number of valence electrons, and the solution to the system will still be the same. Therefore, once we find a solution
A
, we will have to check for all pairs
i,j
that we will only take a single solution out of all the redundant solutions.

Coding it up

Now, we can code up the problem. Let us solve the problem with
C
3
H
4
, which has plenty of redundancy and bonding flexibility from carbon C
​
So, we want to have each carbon atom be connected to at most 4 atoms, and each hydrogen atom connected to one other atom. From above, we can parametrize this problem as an adjacency matrix problem, where atoms
i
and
j
are connected by
n
bonds if and only if there is a
n
at location
(i,j)
in the adjacency matrix.
Now, we would like to make an adjacency matrix out of the following chemical formula:
In[]:=
formula={{"C",3},{"H",4}}
Out[]=
{{C,3},{H,4}}
We can create a list of the atoms from formula in the following way:
In[]:=
atomSymbolList=ConstantArray[#1,#2]&@@@formula//Catenate
Out[]=
{C,C,C,H,H,H,H}
This line of code applies the anonymous function to each element of
formula
(@@@), and takes #1 to be the atom, and #2 to be the frequency, which makes a constant array of an expression, with it’s corresponding frequency. Catenate joins the resulting list of lists into a single list.
Each row in the adjacency matrix needs to be equal to the number of bonds the corresponding atom needs. We can find this list - designated
b
above - by calling
ElementData[]
on each element. We denote this list sumVal
In[]:=
sumVal=(If[MatchQ[#1,"H"|"He"],2,8]-ElementData[#,"ValenceElectronCount"])&/@atomSymbolList​​(*8||2-ElementData[,"ValenceElectronCount"]getsnumberofnecessarybondsoneachelement*)
Out[]=
{4,4,4,1,1,1,1}
Given we have 7 atoms, we can now make the 7

7 adjacency matrix with the additional two constraints that:​​(1) a atom cannot be bonded to itself and ​(2) that bonding is symmetric, so if atom i is bonded to atom j, then atom j is bonded to atom i. ​(1) means that
u[i_,i_]->0
in the matrix. (2) corresponds to a symmetry
u
ij

u
ji
in the adjacency matrix. We can put these ideas into code by creating a symbolic matrix and modifying it with
/.
the
ReplaceAll
and
Rule
. This rule takes elements
u
ji
from the array, and if
j>i
replaces the entry with
u
ij
, as well as replaces any elements with
u
ii
with 0.
In[]:=
arr=Array[u[#1,#2]&,{Length[sumVal],Length[sumVal]}]/.{u[j_,i_]/;j>i:>u[i,j],u[i_,i_]->0};​​​​TableForm[Table[Append[arr[[i]],sumVal[[i]]],{i,Length[arr]}],TableHeadings->{atomSymbolList,Join[atomSymbolList,{"Sum #"}]}]//Framed
Out[]=
C
C
C
H
H
H
H
Sum #
C
0
u[1,2]
u[1,3]
u[1,4]
u[1,5]
u[1,6]
u[1,7]
4
C
u[1,2]
0
u[2,3]
u[2,4]
u[2,5]
u[2,6]
u[2,7]
4
C
u[1,3]
u[2,3]
0
u[3,4]
u[3,5]
u[3,6]
u[3,7]
4
H
u[1,4]
u[2,4]
u[3,4]
0
u[4,5]
u[4,6]
u[4,7]
1
H
u[1,5]
u[2,5]
u[3,5]
u[4,5]
0
u[5,6]
u[5,7]
1
H
u[1,6]
u[2,6]
u[3,6]
u[4,6]
u[5,6]
0
u[6,7]
1
H
u[1,7]
u[2,7]
u[3,7]
u[4,7]
u[5,7]
u[6,7]
0
1
We can solve for
A
, here being the variables
u[1,2],...,u[1,7]
in this case by using the
Solve[eqns,vars,domain]
function
◼
  • Eqns:
    Total[arr[[#]]]&/@Range[Length[atomSymbolList]]==sumVal
    , which gives the sum of every row of the matrix, and it’s corresponding sum rule value.
    • ◼
  • Vars: is a list of the variables
    u[1,2],...,u[1,7]
    , extracted from the matrix
    • ◼
  • Domain:
    NonNegativeIntegers
    specfies that the solutions must be nonnegative integers, so that no half-bonds, or zero-bonds exist.
    • There are a total number of solutions of:
    The second solution to the problem is:
    This solution corresponds to a graph system and therefore, a molecule.

    Deleting Equivalent Permutations

    Next, we want to create a list of complete permutations of all elements, deleting the cases where there are fixed points.
    Note there are 7 solutions, without deleting molecules which are not connected. It is not clear yet that some molecules can still be disconnected.

    Rendering Adjacency Matrix Solutions

    We can visualise this solution by calling MoleculePlot on its output.

    generateIsomersFromFormula

    Calculating Hydrogen Count with DoU

    Another important thing to do is to generate the number of hydrogens automatically based on the other elements in the molecule. We can do this if we are given the number of double-bond equivalents, this is a measure of the number of double bonds or equivalent ring structures as defined here. This is also known as the number of degrees of unsaturation (DoU).

    getHydrogenCount

    Examples and Timing

    The number of hydrogens for three carbon with two double bonds, or a double bond and a ring:
    Below are the execution times and the outputs of this algorithm, and SURGE, respectively.
    We can sort the results by isomeric SMILES, the linear notation for chemical structures, and check that they produce identical outputs
    Now, let us do this with a few more complicated molecules.
    The isomer generation above is interesting in that it shows the weaknesses of modeling programs such as SURGE, and from this algorithm. Many of these molecules that are generated are not realizable, due to instability in O-O bonds. Additionally, both programs can not account for dative bonds between N and O, which can change the stable valence number of N, P and S. For example, the output does not include nitromethane, which would have a nitrogen with 4 bonds, due to an electron donation to O. See here

    This./SURGE comparison

    We may first create a list of potential molecules with C, N, O, Cl. For each valid molecule, we want to compare the of this algorithm vs. SURGE in determining its isomers. Note: we constrain i,j,k,w to keep runtimes generally low.
    We see there is only one isomer for this molecule.
    Now, we filter out non-molecules from complexityComp, based on whether they throw an error. For each list of atoms that possesses contains a list of isomers, we return back a pair. The first element in the pair is the length of that list and the second element is the corresponding runtime of this algorithm/SURGE.

    Conclusion

    This little project illustrates the ability of Wolfram Language to handle problems in chemistry by casting them as graph theory problems. However, this small algorithm does not handle all bonding situations. The program could handle dative bonds by considering “extra-valencies,” but cannot account for specific interactions between specific atoms, for example, the interactions between phosphorus and oxygen in phosphate. In contrast, SURGE builds non-isomorphic graphs without hydrogen, and then renames vertices to atoms and builds bonds based on these specific interactions.
    ​
    Despite this, this program described here relies on only a few key lines of code at runtime, and is a good project for understanding the basics of WL, chemistry, and graph theory. This little project was an extension of my WSS 22 project on graph template problems. Thank you sincerely to Dr. Robert Nachbar for all of the help and support with this project.