Joseph - this is a really great application of molecule substructure replacement. I especially like the use of Graph to visualize the relationships.
One tip I would give is that "ReplaceSubstructureList" can also take replacement rules in the form "SMARTS"->"SMILES", which can make the code less verbose and a bit easier to read for simple substructures. So for the replacement
Molecule[{"O"}, {}] -> Molecule[{"S"}, {}]
you can write "O"->"S" and the result is the same. Likewise
MoleculePattern[{"C", "C"}, {Bond[{1, 2}, "Double"]}] -> Molecule["Ethane"]
could be written "C=C"->"CC". For more complicated substructures, the full form is needed though.
Also note that there is an operator form of MoleculeModify, so you could write an operator that converts a "CH" into a "CC" like
addOneMoreCarbon = MoleculeModify[
{"ReplaceSubstructureList", Atom["C", "HydrogenCount" -> GreaterThan[0]] -> "CC"}
]
I would consider it a bug in MoleculeModify that you have to call the "MakeHydrogensExplicit" method to get consistent results. For simulating large networks like these, keeping the hydrogen atoms implicit saves both time and memory, but the outcome should be the same whether or not hydrogens are included as explicit vertices in the molecular graph. I've filed a bug report for this and hopefully it can be fixed for the next release.
If you run into other issues where "ReplaceSubstructure" is giving wrong results please don't hesitate to report it directly.
