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.

Important

The reactmap package is part of the AMS python stack. You should use the $AMSBIN/amspython command to run the examples in this document. See this page for more info.

Useful examples/templates

There are a number of other template scripts and demonstrations of advanced reactmap functionality in the examples section of this documentation. Users who want to use custom settings, explore active atoms/bonds options, run multiple jobs, etc. will find useful examples here.

reactmap contents

Classes

The following are the main classes required to solve reaction mapping problems.

Molecule([atoms, bonds, smiles, graph, …])	The main class for inputing/using reactants and products.
Reaction(reactant, product[, name])	The main class for inputting (reactant, product) pairs to specify a reaction and investigating the results of a mapping.
Map(reaction[, settings, mapping])	The class that solves the mapping problem on a Reaction instance.
Settings(**kwargs)	This is the main class for inputting settings for solving the problem(s).
RMLogger([print_level, log_level, …])	This is the logger class for reactmap.

Submodules

The package also has a heuristics submodule. This submodule will be called during the course of a normal mapping given certain problem settings, but it is also exposed so the interested user can try heuristic solutions to (hard) problems.

heuristics

heuristics is a submodule for using heuristic techniques to generate upper and lower bounds for reaction mapping problems.

Index

Index