# Generating a function which outputs possible chemical reactions

I want to make a list of chemical reactions and I write them down in a $\require{mhchem}\LaTeX$ format. They are of the following form

$$NA_n^i+MB_m^j \rightarrow \hat NA_{\hat n}^{\hat i}+\hat MA_{\hat m}^{\hat j},$$

the relevant part is

$$NA_n^i+MB_m^j$$

which depends on up to 8 parameters $$A,B,N,M,n,m,i,j,$$ i.e. atomar species $A$, quantity $N$, number of elements in the molecule $n$ and an index $i$ for modi like excitation $A^*$ or ionization $A^+,A^-$.

Like e.g.

$$\cf{O}^*+\cf{2Cl2^-}.$$

I want to create an "object" of function

reactions[A,B,N,M,n,m,i,j,k]


with values $k$, where

$k=1$ evaluates to the TeX-form $NA_n^i+MB_m^j$,

$k=2$ is some string with information,

and the tricky part

$k=3$

gives the list of permissible products of the reaction.

This amounts to implement conservation of atomar species, but as $\ce{2O->O2}$ is legal, I think $n+N=const.$ is the bound. Then there is charge conservation.

I know that there is a ToTeXForm function so, while I'm open for optimization, this problem should be easy to handle. But how would I go about solving this, I guess, combinatorial problem? Does Mathematica want me to write some subroutine, which gets called by my function. Or do you do this all at once? Is there a method to add more values to $k$ afterwards, i.e. if I want to add information later. And does updating of the $k=2$ string just amount to redefining the string all reactions[...,3]:="new information"? You see, I'm thinking about this from an object oriented point of view, but this might be wrong.

The last step would be to also take into account that there are object like $\cf{CH4}$, which can decompose into e.g. $\cf{C}+\cf{2H2}$, but enough for now.

My progress

 tex[{A_, N_, n_, q_}] :=
If[Head[A] === List,
tex[N] <> " + " <> tex[A] <> " + " <> tex[n] <> " + " <> tex[q],
ToString[N] <> " " <> ToString[A] <> "_" <> ToString[n] <> "^" <>
ToString[q]]
tex[{a_, b_, c_}] := tex[a] <> " + " <> tex[b] <> " + " <> tex[c]
tex[{Q_, e_}] :=
If[Head[Q] === List, tex[Q] <> " + " <> tex[e],
ToString[Q] <> " " <> ToString[e]]
tex[{a_}] := tex[a]


Example

 tex@{{CH, 2, 4, 0}}
tex@{{H, 2, 2, 0}, {O, 1, 2, 1}} norm[{{A_, N_, n_, q_}}] := {{A, N*n, 1, q}}
norm[{A_, N_, n_, q_}] := norm[{{A, N, n, q}}]
norm[{{A_, N_, n_, q_}, {B_, M_, m_, p_}}] :=
If[A === B, {{A, N*n + M*m, 1, q + p}}, {{A, N*n, 1, q + p}, {B, M*m,
1, 0}}]
norm[{{A_, N_, n_, q_}, {Q_, e_}}] := {{A, N*n, 1, q - Q}}
norm[{{A_, N_, n_, q_}, {B_, M_, m_, p_}, {Q_, e_}}] :=
If[A === B, {{A, N*n + M*m, 1, q + p - Q}}, {{A, N*n, 1,
q + p - Q}, {B, M*m, 1, 0}}]

summands[int_] := Table[{int - index, index}, {index, 0, Floor[int/2]}]
divisors[int_] := {int/#, #} & /@ Divisors[int]

split[t_] := ({{t[[1, 1]], First[#], 1, t[[1, 4]]}, {t[[1, 1]],
Last[#], 1, t[[1, 4]]}} &) /@ summands[t[[1, 2]]]

shift[t_] :=
If[t[] == 0, {t}, ({t[], First[#], Last[#], t[]} &) /@divisors[t[]]]

multiplyTemp[term_] :=
Flatten[
Map[Table[{#[[1, ind1]], #[[2, ind2]]}, {ind1, 1,
Length[#[] ]}, {ind2, 1, Length[#[] ] }] &,
Map[shift, split[norm[term]], {2}]],2]

multiply[term_] :=
Map[Function[termvar, Select[termvar, (Not[#[] == 0] &)]],
DeleteDuplicates[multiplyTemp[term], (#1 === Reverse[#2] &)]]

products[t_] :=
Which[
Length[norm[t]] == 1, multiply[t],
Length[norm[t]] == 2,
Flatten[Table[
Flatten[{#[[1, ind1]], #[[2, ind2]]}, 1], {ind1, 1,
Length[First[#]]}, {ind2, 1, Length[Last[#]]}] &@
(multiply /@ t), 1]]


Example

 myTerm := {{CH, 2, 3, 0}}
(tex@myTerm " \[LongRightArrow] " <> # &) /@
tex /@ (products@myTerm)    // TableForm myTerm := {{H, 2, 2, 0}, {O, 1, 1, 1}}
(tex@myTerm " \[LongRightArrow] " <> # &) /@
tex /@ (products@myTerm)    // TableForm elaboration

So I realized that if I'm going to generate all possible outcomes of a reaction, based on atomic species number and charge conservation, the general expression

$$NA_n^i+Qe$$

can be normalized to

$$(N\cdot n)A_1^{i-Q}$$

before the generation. This generalizes to

$$NA_n^i+MA_m^i+Qe$$

and

$$NA_n^i+MB_m^i+Qe,$$

with $A\ne B$. There are some If-clauses in the code which result from special treatment of $A\ne B$.

So firstly I normalize (norm), then generate all the possible decompositions (split)

$$(N\cdot n)A_1^i \Rightarrow (N\cdot n-k)A_1^i+kA_1^i$$

and then I generate all the combinations $KA_k^i, \ \ \ \ \ \ K\cdot k==N\cdot n$. These are the functions shift and multiply, where multiplyTemp is just the result pre filter (deletion of expressions $0A_k^i$ and dublicates).

I'm open for suggestions regarding the above code up to now.

I stick to reactions with integer coefficients here, but in case one wants to generalize that, one would just have to change the shift function to do something else.

Next step are the charge conservation shift. The total charge of an expression is easy to read off. I guess it makes sense to define a global variable ionizationlevel (which we will usually set to 1 or 2, since higher ionizations aren't that common) and then take every list from products and generate a couple of new lists consisting of the same species, but with charges shifted around between the species as well as possible freed electrons $Q e$.

Them manybe comes the question of putting species together. This will give a lot of new options of course. Any ideas how to implement it?

• I did what I can, but just so you know: MathJax supports the mhchem extension; you might want to use it the next time you edit... – J. M.'s torpor Nov 19 '12 at 11:10
• @Nick Kidman I don't understand the reasoning for wanting to use an integer to switch the behavior of a function. Why not just create different functions? like creating one reactionsTexForm[] and one reactionsInformation[] instead of having the last integer argument completely switch what the function is doing (and switching both type and meaning of output arguments even). – jVincent Nov 19 '12 at 11:23
• @jVincent: Yeah, okay, makes sense. The function which generates reaction is really the big part. I wanted to start doing this as a Java object, when I decided that all the list generation would maybe be better suited for Mathematica, this is why it looks like that (something which returns output for k=1,2,3) in my brain atm. – Nikolaj-K Nov 19 '12 at 11:26
• A single function might work well if, for all three cases of k, the other arguments always have the same pattern. Should that hold, I would implement k as a option. Something like: output -> "TeX", output -> "Information", output -> "Products". – m_goldberg Nov 19 '12 at 17:10