Generating a function which outputs possible chemical reactions

Question

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[[2]] == 0, {t}, ({t[[1]], First[#], Last[#], t[[4]]} &) /@divisors[t[[2]]]]

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

 multiply[term_] := 
 Map[Function[termvar, Select[termvar, (Not[#[[2]] == 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?