reactmap¶

The reactmap python package determines an optimal atom mapping between reactants and products in a chemical reaction. An atom mapping is an assignment of every atom in the reactants to a unique atom in the products. More simply, an atom mapping is an answer to the question where does each atom in the reactant(s) end up in the product(s)?. The optimal atom mapping as defined in reactmap is the mapping that minimizes the number of required (induced) bond changes (bond breakages in the reactant and bond formations in the product).

Fig. 1 An optimal atom mapping for a hydrolysis reaction. Note the methyl groups are mapped as a unit to reduce the figure complexity. The red bonds indicate required bond breakages in the reactants, and the green bonds indicate required bond formations in the product. The 2 bond formations and 2 bond breakages result in a cost of 4 for this mapping.

Example:

Let’s take a look at the hydrolysis example above using reactmap. To solve this reaction mapping problem, we’ll first input the reactant and product molecules as SMILES strings and then create a reaction instance for the problem. This reaction instance will be given to the Map class, which will determine the optimal mapping. All results from the mapping are stored in the reaction instance. Below, we’ll take a look at 2 important values: the cost of the mapping (the number of required bond changes); and the mapping itself (the assignment of reactant atoms to product atoms). Note that the indices of the mapping list correspond to the atom indices of the reactants. The entries are the atom indices in the products. Additionally, we’ll print some more descriptive output and generate some useful images illustrating the mapping.

>>> import scm.reactmap as rm
>>> reactant = rm.Molecule(smiles="CC(=O)OC.O")
>>> product  = rm.Molecule(smiles="CC(=O)O.CO")
>>> reaction = rm.Reaction(reactant = reactant, product = product, name = "hydrolysis")
>>> rm.Map(reaction)
>>> print(reaction.cost)
4
>>> print(reaction.mapping)
[0, 1, 3, 2, 4, 5, 6, 7, 8, 10, 11, 12, 13, 9]
>>> reaction.show_results()
Results for hydrolysis
-> Problem details:
        Number of atoms           : 14
        Atom counts               : C:3, O:3, H:8
        Number of reactant bonds  : 12 (12 active)
        Number of product bonds   : 12 (12 active)
-> Atom mapping details:
        Reactant      Product   Atom  Atom Status
        index         index     type  (A/I/D)
        ─────────────────────────────────────────
               0  ->       0    C     A
               1  ->       1    C     A
               2  ->       3    O     A
               3  ->       2    O     A
               4  ->       4    C     A
               5  ->       5    O     A
               6  ->       6    H     A
               7  ->       7    H     A
               8  ->       8    H     A
               9  ->      10    H     A
              10  ->      11    H     A
              11  ->      12    H     A
              12  ->      13    H     A
              13  ->       9    H     A
        ─────────────────────────────────────────
-> Bond change details:
        Broken: O-C bond between atoms 3 and 4 in the Reactants
        Broken: O-H bond between atoms 5 and 13 in the Reactants
        Formed: C-O bond between atoms 4 and 5 in the Products
        Formed: O-H bond between atoms 3 and 9 in the Products
        (Total broken: 2, Total formed: 2)
-> Solution details:
        Found solution          : True
        Time to find solution   : 0.022364 s
        Cost                    : 4

We can additionally call the function output_image to produce 2 files corresponding to the reactant(s) and product(s). These images show the optimal mapping between reactant and product atoms (atoms with the same index in the images are mapped to one another; these indices correspond to the original reactant indices but need not correspond to the original product indices). These images also show required bond breakages in red and required bond formations in green.

>>> reaction.output_image("example")

This function produces two output files (example_reactants.svg and example_products.svg). These files are shown below.

example_reactants.svg

example_products.svg

This solution does not represent the true mechanism of this reaction, but it is still optimal in terms of the required bond changes (it is degenerate with a more chemically correct mapping). The mapping that is more consistent with the hydrolysis mechanism can be achieved using active_bonds, as shown in this example.

Index¶

Index