Repository: OpenFreeEnergy/gufe (language: python)

Issue:
Move `openfe` `LigandAtomMapping` into `gufe`?
@ijpulidos and I are working on [nonequilibrium cycling](https://github.com/choderalab/perses/pull/1066) in `perses`, and were wondering if [`LigandAtomMapping`](https://github.com/OpenFreeEnergy/openfe/blob/82a6a6db4c22cc5b30ddf26253db10130a24662f/openfe/setup/atom_mapping/ligandatommapping.py#L13) from OpenFreeEnergy/openfe made more sense here in `gufe`? 

The `AtomMapping` base class is an abstract base class, so we can't directly use that, and think it best if we avoid implementing our own child class in `perses` for the mappings used there.


Discussion / hints:
I wonder if it may make more sense for the protocol to depend on openfe? My understanding is that the long term hope is that openfe acts as a gateway to packaged protocols, so it would make sense for shared classes that implements gufe abstract classes that multiple protocols can use to sit in openfe?
`gufe` is the place for data models (the "nouns" in FE). Ligand atom mappings (or more accurately, small molecule atom mappings) are an example of that. I think it is reasonable to move that into here. We have `SmallMoleculeComponent`, why not a `SmallMoleculeAtomMapping`?

OTOH, the `AtomMapper` base class (which is expected by the `openfe` network planners) is an OpenFE-specific API for generating mappings. I'm less convinced that it would belong in `gufe`.
I think the only problem I can see is that `LigandAtomMapping` has a lot of visualisation code.  So if we're ok with also including that in gufe then sure.

One problem I can see with the status quo is that if we want to package the non-eq Protocol in the openfe app, there'll be a gross circular dependency
> I think the only problem I can see is that LigandAtomMapping has a lot of visualisation code. So if we're ok with also including that in gufe then sure.

Maybe we can put the core small molecule mapping data model in `gufe` and leave the visualization as an `openfe` subclass/extension?

> One problem I can see with the status quo is that if we want to package the non-eq Protocol in the openfe app, there'll be a gross circular dependency

I had the same thought. It's a solvable issue, but yeah, kinda gross. I think `perses` would have to make `openfe` an optional dependency. For `openfe`, we're already treating perses as optional for atom mappers/scorers.
The viz stuff is something that I know other folks using atom mapping want to use (I know it was a request on alchemiscale and iirc someone mentioned in a perses issue at some point that it would be great to bring over access to the 3D viz stuff).

Fix patch:
diff --git a/gufe/__init__.py b/gufe/__init__.py
index 929ecdb0..cc023de6 100644
--- a/gufe/__init__.py
+++ b/gufe/__init__.py
@@ -6,6 +6,8 @@
 
 from . import tokenization
 
+from . import visualization
+
 from .components import (
     Component, 
     SmallMoleculeComponent, 
@@ -18,6 +20,7 @@
 from .mapping import (
     ComponentMapping,  # how individual Components relate
     AtomMapping, AtomMapper,  # more specific to atom based components
+    LigandAtomMapping,
 )
 from .protocols import (
     Protocol,  # description of a method
diff --git a/gufe/mapping/__init__.py b/gufe/mapping/__init__.py
index 8ea4f5b2..74988ac4 100644
--- a/gufe/mapping/__init__.py
+++ b/gufe/mapping/__init__.py
@@ -4,3 +4,4 @@
 from .componentmapping import ComponentMapping
 from .atom_mapping import AtomMapping
 from .atom_mapper import AtomMapper
+from .ligandatommapping import LigandAtomMapping
diff --git a/gufe/mapping/ligandatommapping.py b/gufe/mapping/ligandatommapping.py
new file mode 100644
index 00000000..c6f63ddb
--- /dev/null
+++ b/gufe/mapping/ligandatommapping.py
@@ -0,0 +1,181 @@
+# This code is part of gufe and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/gufe
+import json
+from typing import Any, Optional
+import numpy as np
+from numpy.typing import NDArray
+
+from gufe.components import SmallMoleculeComponent
+from gufe.visualization.mapping_visualization import draw_mapping
+from . import AtomMapping
+
+
+class LigandAtomMapping(AtomMapping):
+    """Simple container with the mapping between two Molecules
+
+    """
+    componentA: SmallMoleculeComponent
+    componentB: SmallMoleculeComponent
+
+    def __init__(
+        self,
+        componentA: SmallMoleculeComponent,
+        componentB: SmallMoleculeComponent,
+        componentA_to_componentB: dict[int, int],
+        annotations: Optional[dict[str, Any]] = None,
+    ):
+        """
+        Parameters
+        ----------
+        componentA, componentB : SmallMoleculeComponent
+          the ligand molecules on either end of the mapping
+        componentA_to_componentB : dict[int, int]
+          correspondence of indices of atoms between the two ligands
+        annotations : dict[str, Any]
+          Mapping of annotation identifier to annotation data. Annotations may
+          contain arbitrary JSON-serializable data. Annotation identifiers
+          starting with ``ofe-`` may have special meaning in other parts of
+          OpenFE. ``score`` is a reserved annotation identifier.
+        """
+        super().__init__(componentA, componentB)
+        self._compA_to_compB = componentA_to_componentB
+
+        if annotations is None:
+            # TODO: this should be a frozen dict
+            annotations = {}
+
+        self._annotations = annotations
+
+    @classmethod
+    def _defaults(self):
+        return {}
+
+    @property
+    def componentA_to_componentB(self) -> dict[int, int]:
+        return dict(self._compA_to_compB)
+
+    @property
+    def componentB_to_componentA(self) -> dict[int, int]:
+        return {v: k for k, v in self._compA_to_compB.items()}
+
+    @property
+    def componentA_unique(self):
+        return (i for i in range(self.componentA.to_rdkit().GetNumAtoms())
+                if i not in self._compA_to_compB)
+
+    @property
+    def componentB_unique(self):
+        return (i for i in range(self.componentB.to_rdkit().GetNumAtoms())
+                if i not in self._compA_to_compB.values())
+
+    @property
+    def annotations(self):
+        # return a copy (including copy of nested)
+        return json.loads(json.dumps(self._annotations))
+
+    def _to_dict(self):
+        """Serialize to dict"""
+        return {
+            # openff serialization doesn't go deep, so stringify at this level
+            'componentA': json.dumps(self.componentA.to_dict(), sort_keys=True),
+            'componentB': json.dumps(self.componentB.to_dict(), sort_keys=True),
+            'componentA_to_componentB': self._compA_to_compB,
+            'annotations': json.dumps(self.annotations),
+        }
+
+    @classmethod
+    def _from_dict(cls, d: dict):
+        """Deserialize from dict"""
+        # the mapping dict gets mangled sometimes
+        mapping = d['componentA_to_componentB']
+        fixed = {int(k): int(v) for k, v in mapping.items()}
+
+        return cls(
+            componentA=SmallMoleculeComponent.from_dict(
+                json.loads(d['componentA'])),
+            componentB=SmallMoleculeComponent.from_dict(
+                json.loads(d['componentB'])),
+            componentA_to_componentB=fixed,
+            annotations=json.loads(d['annotations'])
+        )
+
+    def __repr__(self):
+        return (f"{self.__class__.__name__}(componentA={self.componentA!r}, "
+                f"componentB={self.componentB!r}, "
+                f"componentA_to_componentB={self._compA_to_compB!r}, "
+                f"annotations={self.annotations!r})")
+
+    def _ipython_display_(self, d2d=None):  # pragma: no-cover
+        """
+        Visualize atom mapping in a Jupyter Notebook.
+
+        Parameters
+        ---------
+        d2d : :class:`rdkit.Chem.Draw.rdMolDraw2D.MolDraw2D`
+            If desired specify an instance of a MolDraw2D object.
+            Default ``None`` will use the MolDraw2DCairo backend.
+
+        Returns
+        -------
+        Image: IPython.core.display.Image
+            Image of the atom map
+        """
+        from IPython.display import Image, display
+
+        return display(Image(draw_mapping(self._compA_to_compB,
+                                          self.componentA.to_rdkit(),
+                                          self.componentB.to_rdkit(), d2d)))
+
+    def draw_to_file(self, fname: str, d2d=None):
+        """
+        Save atom map visualization to disk
+
+        Parameters
+        ---------
+        d2d : :class:`rdkit.Chem.Draw.rdMolDraw2D.MolDraw2D`
+            If desired specify an instance of a MolDraw2D object.
+            Default ``None`` will write a .png file using the MolDraw2DCairo
+            backend.
+
+        fname : str
+            Name of file to save atom map
+        """
+        data = draw_mapping(self._compA_to_compB, self.componentA.to_rdkit(),
+                            self.componentB.to_rdkit(), d2d)
+        if type(data) == bytes:
+            mode = "wb"
+        else:
+            mode = "w"
+
+        with open(fname, mode) as f:
+            f.write(draw_mapping(self._compA_to_compB, self.componentA.to_rdkit(),
+                                 self.componentB.to_rdkit(), d2d))
+
+    def with_annotations(self, annotations: dict[str, Any]):
+        """Create an new mapping based on this one with extra annotations.
+
+        Parameters
+        ----------
+        annotations : dict[str, Any]
+            Annotation update for this mapping. New annotation keys will be
+            added to the annotations dict; existing keys will be replaced by
+            the data provided here.
+        """
+        return self.__class__(
+            componentA=self.componentA,
+            componentB=self.componentB,
+            componentA_to_componentB=self._compA_to_compB,
+            annotations=dict(**self.annotations, **annotations)
+        )
+
+    def get_distances(self) -> NDArray[np.float64]:
+        """Return the distances between pairs of atoms in the mapping"""
+        dists = []
+        molA = self.componentA.to_rdkit().GetConformer()
+        molB = self.componentB.to_rdkit().GetConformer()
+        for i, j in self._compA_to_compB.items():
+            dA = molA.GetAtomPosition(i)
+            dB = molB.GetAtomPosition(j)
+            dists.append(dA.Distance(dB))
+
+        return np.array(dists)
diff --git a/gufe/visualization/__init__.py b/gufe/visualization/__init__.py
new file mode 100644
index 00000000..3cd3e5db
--- /dev/null
+++ b/gufe/visualization/__init__.py
@@ -0,0 +1,4 @@
+# This code is part of gufe and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/gufe
+
+from . import mapping_visualization
diff --git a/gufe/visualization/mapping_visualization.py b/gufe/visualization/mapping_visualization.py
new file mode 100644
index 00000000..18a2677d
--- /dev/null
+++ b/gufe/visualization/mapping_visualization.py
@@ -0,0 +1,284 @@
+# This code is part of gufe and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/openfe
+from typing import Any, Collection, Optional
+from itertools import chain
+
+from rdkit import Chem
+from rdkit.Chem import Draw
+from rdkit.Chem import AllChem
+
+
+def _match_elements(
+        mol1: Chem.Mol, idx1: int, mol2: Chem.Mol, idx2: int
+) -> bool:
+    """
+    Convenience method to check if elements between two molecules (molA
+    and molB) are the same.
+
+    Parameters
+    ----------
+    mol1 : RDKit.Mol
+        RDKit representation of molecule 1.
+    idx1 : int
+        Index of atom to check in molecule 1.
+    mol2 : RDKit.Mol
+        RDKit representation of molecule 2.
+    idx2 : int
+        Index of atom to check in molecule 2.
+
+    Returns
+    -------
+    bool
+        True if elements are the same, False otherwise.
+    """
+    elem_mol1 = mol1.GetAtomWithIdx(idx1).GetAtomicNum()
+    elem_mol2 = mol2.GetAtomWithIdx(idx2).GetAtomicNum()
+    return elem_mol1 == elem_mol2
+
+
+def _get_unique_bonds_and_atoms(
+        mapping: dict[int, int], mol1: Chem.Mol, mol2: Chem.Mol
+) -> dict:
+    """
+    Given an input mapping, returns new atoms, element changes, and
+    involved bonds.
+
+    Parameters
+    ----------
+    mapping : dict of int:int
+        Dictionary describing the atom mapping between molecules 1 and 2.
+    mol1 : RDKit.Mol
+        RDKit representation of molecule 1.
+    mol2 : RDKit.Mol
+        RDKit representation of molecule 2.
+
+    Returns
+    -------
+    uniques : dict
+        Dictionary containing; unique atoms ("atoms"), new elements
+        ("elements"), and bonds involved with either of the previous two
+        ("bonds") for molecule 1.
+    """
+
+    uniques: dict[str, set] = {
+        "atoms": set(),  # atoms which fully don't exist in molB
+        "elements": set(),  # atoms which exist but change elements in molB
+        "bonds": set(),  # bonds involving either unique atoms or elements
+    }
+
+    for at in mol1.GetAtoms():
+        idx = at.GetIdx()
+        if idx not in mapping:
+            uniques["atoms"].add(idx)
+        elif not _match_elements(mol1, idx, mol2, mapping[idx]):
+            uniques["elements"].add(idx)
+
+    for bond in mol1.GetBonds():
+        bond_at_idxs = [bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()]
+        for at in chain(uniques["atoms"], uniques["elements"]):
+            if at in bond_at_idxs:
+                bond_idx = bond.GetIdx()
+                uniques["bonds"].add(bond_idx)
+
+    return uniques
+
+
+def _draw_molecules(
+    d2d,
+    mols: Collection[Chem.Mol],
+    atoms_list: Collection[set[int]],
+    bonds_list: Collection[set[int]],
+    atom_colors: Collection[dict[Any, tuple[float, float, float, float]]],
+    highlight_color: tuple[float, float, float, float],
+    atom_mapping: Optional[dict[tuple[int, int], dict[int, int]]] = None,
+) -> str:
+    """
+    Internal method to visualize a molecule, possibly with mapping info
+
+    Parameters
+    ----------
+    d2d :
+        renderer to draw the molecule; currently we only support
+        rdkit.rdMolDraw2D
+    mols : Collection[Chem.Mol]
+        molecules to visualize
+    atoms_list: Collection[Set[int]]
+        iterable containing one set per molecule in ``mols``, with each set
+        containing the indices of the atoms to highlight
+    bonds_list: Collection[Set[int]]
+        iterable containing one set per molecule in ``mols``, with each set
+        containing the indices of the atoms involved in bonds to highlight
+    atom_colors: Collection[Dict[Any, Tuple[float, float, float, float]]]
+        iterable containing one dict per molecule in ``mols``, with each
+        dict containing a mapping of RDKit atom to color, expressed as an
+        RGBA tuple, for atoms that need special coloring (e.g., element
+        changes)
+    highlight_color: Tuple[float, float, float, float]
+        RGBA tuple for the default highlight color used in the mapping
+        visualization
+    atom_mapping: Dict[Tuple[int,int], Dict[int, int]], optional
+        used to align the molecules to each othter for clearer visualization.
+        The structure contains the indices of the molecules in mols as key
+        Tuple[molA, molB] and maps the atom indices via the value Dict[
+        molA_atom_idx, molB_atom_idx]
+        default None
+    """
+    # input standardization:
+    if atom_mapping is None:
+        atom_mapping = {}
+
+    if d2d is None:
+        # select default layout based on number of molecules
+        grid_x, grid_y = {1: (1, 1), 2: (2, 1), }[len(mols)]
+        d2d = Draw.rdMolDraw2D.MolDraw2DCairo(
+            grid_x * 300, grid_y * 300, 300, 300)
+
+    # squash to 2D
+    copies = [Chem.Mol(mol) for mol in mols]
+    for mol in copies:
+        AllChem.Compute2DCoords(mol)
+
+    # mol alignments if atom_mapping present
+    for (i, j), atomMap in atom_mapping.items():
+        AllChem.AlignMol(
+            copies[j], copies[i],
+            atomMap=[(k, v) for v, k in atomMap.items()]
+        )
+
+    # standard settings for our visualization
+    d2d.drawOptions().useBWAtomPalette()
+    d2d.drawOptions().continousHighlight = False
+    d2d.drawOptions().setHighlightColour(highlight_color)
+    d2d.drawOptions().addAtomIndices = True
+    d2d.DrawMolecules(
+        copies,
+        highlightAtoms=atoms_list,
+        highlightBonds=bonds_list,
+        highlightAtomColors=atom_colors,
+    )
+    d2d.FinishDrawing()
+    return d2d.GetDrawingText()
+
+
+def draw_mapping(
+        mol1_to_mol2: dict[int, int], mol1: Chem.Mol, mol2: Chem.Mol, d2d=None
+):
+    """
+    Method to visualise the atom map correspondence between two rdkit
+    molecules given an input mapping.
+
+    Legend:
+        * Red highlighted atoms: unique atoms, i.e. atoms which are not
+          mapped.
+        * Blue highlighted atoms: element changes, i.e. atoms which are
+          mapped but change elements.
+        * Red highlighted bonds: any bond which involves at least one
+          unique atom or one element change.
+
+    Parameters
+    ----------
+    mol1_to_mol2 : dict of int:int
+        Atom mapping between input molecules.
+    mol1 : RDKit.Mol
+        RDKit representation of molecule 1
+    mol2 : RDKit.Mol
+        RDKit representation of molecule 2
+    d2d : :class:`rdkit.Chem.Draw.rdMolDraw2D.MolDraw2D`
+        Optional MolDraw2D backend to use for visualisation.
+    """
+    mol1_uniques = _get_unique_bonds_and_atoms(mol1_to_mol2, mol1, mol2)
+
+    # invert map
+    mol2_to_mol1_map = {v: k for k, v in mol1_to_mol2.items()}
+    mol2_uniques = _get_unique_bonds_and_atoms(mol2_to_mol1_map, mol2, mol1)
+
+    atoms_list = [
+        mol1_uniques["atoms"] | mol1_uniques["elements"],
+        mol2_uniques["atoms"] | mol2_uniques["elements"],
+    ]
+
+    # highlight core element changes differently from unique atoms
+    # RGBA color value needs to be between 0 and 1, so divide by 255
+    red = (220 / 255, 50 / 255, 32 / 255, 1.0)
+    blue = (0.0, 90 / 255, 181 / 255, 1.0)
+
+    at1_colors = {at: blue for at in mol1_uniques["elements"]}
+    at2_colors = {at: blue for at in mol2_uniques["elements"]}
+
+    atom_colors = [at1_colors, at2_colors]
+
+    bonds_list = [mol1_uniques["bonds"], mol2_uniques["bonds"]]
+
+    return _draw_molecules(
+        d2d,
+        [mol1, mol2],
+        atoms_list=atoms_list,
+        bonds_list=bonds_list,
+        atom_colors=atom_colors,
+        highlight_color=red,
+        atom_mapping={(0, 1): mol1_to_mol2},
+    )
+
+
+def draw_one_molecule_mapping(mol1_to_mol2, mol1, mol2, d2d=None):
+    """Draw the mapping visualization for a single molecular of a mapping
+
+    This will always draw ``molA``. To draw ``molB``, switch order/invert
+    ``molA_to_molB`` mapping.
+
+    See :func:`.draw_mapping` for details on the meaning of different
+    colors.
+
+    Parameters
+    ----------
+    mol1_to_mol2 : dict of int:int
+        Atom mapping between input molecules.
+    mol1 : RDKit.Mol
+        RDKit representation of molecule 1
+    mol2 : RDKit.Mol
+        RDKit representation of molecule 2
+    d2d : :class:`rdkit.Chem.Draw.rdMolDraw2D.MolDraw2D`
+        Optional MolDraw2D backend to use for visualisation.
+
+    """
+    uniques = _get_unique_bonds_and_atoms(mol1_to_mol2, mol1, mol2)
+    atoms_list = [uniques["atoms"] | uniques["elements"]]
+    bonds_list = [uniques["bonds"]]
+    red = (220 / 255, 50 / 255, 32 / 255, 1.0)
+    blue = (0, 90 / 255, 181 / 255, 1.0)
+
+    atom_colors = [{at: blue for at in uniques["elements"]}]
+
+    return _draw_molecules(d2d,
+                           [mol1],
+                           atoms_list,
+                           bonds_list,
+                           atom_colors,
+                           red,
+                           )
+
+
+def draw_unhighlighted_molecule(mol, d2d=None):
+    """
+    Draw a molecule without any mapping information.
+
+    This uses the same mechanisms as the mapping visualizations, and so can
+    be useful in cases where you want the same visual style both for a plain
+    molecule and a molecule with the mapping highlighted.
+
+    Parameters
+    ----------
+    mol : RDKit.Mol
+        RDKit representation of the molecule
+    d2d : :class:`rdkit.Chem.Draw.rdMolDraw2D.MolDraw2D`
+        Optional MolDraw2D backend to use for visualisation.
+    """
+    red = (220 / 255, 50 / 255, 32 / 255, 1.0)
+    return _draw_molecules(
+        d2d,
+        [mol],
+        atoms_list=[[]],
+        bonds_list=[[]],
+        atom_colors=[{}],
+        highlight_color=red,
+    )


Test patch:
diff --git a/gufe/tests/conftest.py b/gufe/tests/conftest.py
index 82220d0a..2d7c42db 100644
--- a/gufe/tests/conftest.py
+++ b/gufe/tests/conftest.py
@@ -150,7 +150,7 @@ def PDBx_181L_openMMClean_path():
 
 ## RDKit molecules
 
-@pytest.fixture
+@pytest.fixture(scope='session')
 def benzene_modifications():
     with importlib.resources.path('gufe.tests.data',
                                   'benzene_modifications.sdf') as f:
@@ -161,6 +161,12 @@ def benzene_modifications():
     return {m.GetProp('_Name'): m for m in mols}
 
 
+@pytest.fixture(scope='session')
+def benzene_transforms(benzene_modifications):
+    return {k: gufe.SmallMoleculeComponent(v)
+            for k, v in benzene_modifications.items()}
+
+
 ## Components
 
 
diff --git a/gufe/tests/data/lomap_basic/1,3,7-trimethylnaphthalene.mol2 b/gufe/tests/data/lomap_basic/1,3,7-trimethylnaphthalene.mol2
new file mode 100644
index 00000000..4075dbcd
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/1,3,7-trimethylnaphthalene.mol2
@@ -0,0 +1,63 @@
+@<TRIPOS>MOLECULE
+*****
+ 27 28 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.9251    0.1494   -0.0025 C.3     1  LIG1       -0.0397
+      2 C           2.4242    0.1552    0.0465 C.ar    1  LIG1       -0.0498
+      3 C           3.1297   -1.0574    0.0810 C.ar    1  LIG1       -0.0509
+      4 C           4.5422   -1.1042    0.1427 C.ar    1  LIG1       -0.0147
+      5 C           5.2861   -2.3145    0.1727 C.ar    1  LIG1       -0.0425
+      6 C           6.6903   -2.2735    0.2331 C.ar    1  LIG1       -0.0553
+      7 C           7.3760   -1.0578    0.2620 C.ar    1  LIG1       -0.0496
+      8 C           6.6436    0.1328    0.2358 C.ar    1  LIG1       -0.0511
+      9 C           5.2413    0.1261    0.1765 C.ar    1  LIG1       -0.0174
+     10 C           4.5298    1.3362    0.1522 C.ar    1  LIG1       -0.0538
+     11 C           3.1361    1.3519    0.0914 C.ar    1  LIG1       -0.0582
+     12 C           8.8737   -1.0147    0.3281 C.3     1  LIG1       -0.0397
+     13 C           4.6094   -3.6603    0.1443 C.3     1  LIG1       -0.0391
+     14 H           0.5282    1.1476   -0.2149 H       1  LIG1        0.0278
+     15 H           0.5194   -0.1845    0.9575 H       1  LIG1        0.0278
+     16 H           0.5728   -0.5198   -0.7944 H       1  LIG1        0.0278
+     17 H           2.5577   -1.9826    0.0575 H       1  LIG1        0.0626
+     18 H           7.2540   -3.2043    0.2569 H       1  LIG1        0.0623
+     19 H           7.1781    1.0810    0.2619 H       1  LIG1        0.0626
+     20 H           5.0603    2.2854    0.1820 H       1  LIG1        0.0624
+     21 H           2.6157    2.3065    0.0777 H       1  LIG1        0.0620
+     22 H           9.3069   -2.0200    0.3015 H       1  LIG1        0.0278
+     23 H           9.1960   -0.5332    1.2567 H       1  LIG1        0.0278
+     24 H           9.2731   -0.4567   -0.5249 H       1  LIG1        0.0278
+     25 H           5.3371   -4.4789    0.1591 H       1  LIG1        0.0278
+     26 H           4.0157   -3.7691   -0.7690 H       1  LIG1        0.0278
+     27 H           3.9645   -3.7798    1.0209 H       1  LIG1        0.0278
+@<TRIPOS>BOND
+     1     6     7   ar
+     2     5     6   ar
+     3     7     8   ar
+     4     8     9   ar
+     5     9    10   ar
+     6     4     9   ar
+     7    10    11   ar
+     8     2    11   ar
+     9     2     3   ar
+    10     3     4   ar
+    11     4     5   ar
+    12     1     2    1
+    13     5    13    1
+    14     7    12    1
+    15     1    14    1
+    16     1    15    1
+    17     1    16    1
+    18     3    17    1
+    19     6    18    1
+    20     8    19    1
+    21    10    20    1
+    22    11    21    1
+    23    12    22    1
+    24    12    23    1
+    25    12    24    1
+    26    13    25    1
+    27    13    26    1
+    28    13    27    1
diff --git a/gufe/tests/data/lomap_basic/1-butyl-4-methylbenzene.mol2 b/gufe/tests/data/lomap_basic/1-butyl-4-methylbenzene.mol2
new file mode 100644
index 00000000..3c4d3de9
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/1-butyl-4-methylbenzene.mol2
@@ -0,0 +1,62 @@
+@<TRIPOS>MOLECULE
+*****
+ 27 27 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           1.9899   -1.5496    1.4095 C.3     1  LIG1       -0.0653
+      2 C           3.0612   -0.6401    0.8303 C.3     1  LIG1       -0.0556
+      3 C           3.5100    0.4139    1.8419 C.3     1  LIG1       -0.0493
+      4 C           4.5859    1.3256    1.2470 C.3     1  LIG1       -0.0277
+      5 C           5.0777    2.3433    2.2464 C.ar    1  LIG1       -0.0473
+      6 C           4.4793    3.6056    2.3316 C.ar    1  LIG1       -0.0583
+      7 C           4.9350    4.5439    3.2608 C.ar    1  LIG1       -0.0586
+      8 C           5.9858    4.2317    4.1299 C.ar    1  LIG1       -0.0504
+      9 C           6.5860    2.9711    4.0431 C.ar    1  LIG1       -0.0586
+     10 C           6.1309    2.0320    3.1146 C.ar    1  LIG1       -0.0583
+     11 C           6.4975    5.2565    5.0979 C.3     1  LIG1       -0.0397
+     12 H           1.6931   -2.3034    0.6726 H       1  LIG1        0.0230
+     13 H           1.0980   -0.9785    1.6885 H       1  LIG1        0.0230
+     14 H           2.3560   -2.0720    2.3000 H       1  LIG1        0.0230
+     15 H           2.6700   -0.1517   -0.0688 H       1  LIG1        0.0263
+     16 H           3.9190   -1.2484    0.5229 H       1  LIG1        0.0263
+     17 H           3.8957   -0.0824    2.7418 H       1  LIG1        0.0268
+     18 H           2.6473    1.0144    2.1583 H       1  LIG1        0.0268
+     19 H           4.1939    1.8426    0.3620 H       1  LIG1        0.0313
+     20 H           5.4358    0.7279    0.8940 H       1  LIG1        0.0313
+     21 H           3.6566    3.8691    1.6710 H       1  LIG1        0.0620
+     22 H           4.4634    5.5233    3.3010 H       1  LIG1        0.0620
+     23 H           7.4138    2.7097    4.6988 H       1  LIG1        0.0620
+     24 H           6.6080    1.0557    3.0700 H       1  LIG1        0.0620
+     25 H           7.0175    4.7827    5.9374 H       1  LIG1        0.0278
+     26 H           5.6710    5.8412    5.5155 H       1  LIG1        0.0278
+     27 H           7.1937    5.9336    4.5926 H       1  LIG1        0.0278
+@<TRIPOS>BOND
+     1     5     6   ar
+     2     6     7   ar
+     3     5    10   ar
+     4     9    10   ar
+     5     8     9   ar
+     6     7     8   ar
+     7     1     2    1
+     8     2     3    1
+     9     3     4    1
+    10     4     5    1
+    11     8    11    1
+    12     1    12    1
+    13     1    13    1
+    14     1    14    1
+    15     2    15    1
+    16     2    16    1
+    17     3    17    1
+    18     3    18    1
+    19     4    19    1
+    20     4    20    1
+    21     6    21    1
+    22     7    22    1
+    23     9    23    1
+    24    10    24    1
+    25    11    25    1
+    26    11    26    1
+    27    11    27    1
diff --git a/gufe/tests/data/lomap_basic/2,6-dimethylnaphthalene.mol2 b/gufe/tests/data/lomap_basic/2,6-dimethylnaphthalene.mol2
new file mode 100644
index 00000000..1ead0044
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/2,6-dimethylnaphthalene.mol2
@@ -0,0 +1,57 @@
+@<TRIPOS>MOLECULE
+*****
+ 24 25 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.8817   -0.0759    0.0737 C.3     1  LIG1       -0.0397
+      2 C           2.3811   -0.0540    0.0354 C.ar    1  LIG1       -0.0499
+      3 C           3.0834    1.1532    0.0347 C.ar    1  LIG1       -0.0582
+      4 C           4.4799    1.1582   -0.0011 C.ar    1  LIG1       -0.0538
+      5 C           5.2032   -0.0442   -0.0356 C.ar    1  LIG1       -0.0176
+      6 C           6.6087   -0.0550   -0.0772 C.ar    1  LIG1       -0.0511
+      7 C           7.3242   -1.2594   -0.1176 C.ar    1  LIG1       -0.0499
+      8 C           8.8235   -1.2393   -0.1486 C.3     1  LIG1       -0.0397
+      9 C           6.6224   -2.4667   -0.0984 C.ar    1  LIG1       -0.0582
+     10 C           5.2258   -2.4717   -0.0573 C.ar    1  LIG1       -0.0538
+     11 C           4.5024   -1.2692   -0.0283 C.ar    1  LIG1       -0.0176
+     12 C           3.0971   -1.2586    0.0081 C.ar    1  LIG1       -0.0511
+     13 H           0.4638    0.9355    0.0362 H       1  LIG1        0.0278
+     14 H           0.5354   -0.5508    0.9973 H       1  LIG1        0.0278
+     15 H           0.4903   -0.6304   -0.7856 H       1  LIG1        0.0278
+     16 H           2.5515    2.1015    0.0601 H       1  LIG1        0.0620
+     17 H           5.0016    2.1125   -0.0033 H       1  LIG1        0.0624
+     18 H           7.1539    0.8869   -0.0807 H       1  LIG1        0.0626
+     19 H           9.2330   -2.2332   -0.3578 H       1  LIG1        0.0278
+     20 H           9.2150   -0.9034    0.8167 H       1  LIG1        0.0278
+     21 H           9.1785   -0.5651   -0.9354 H       1  LIG1        0.0278
+     22 H           7.1546   -3.4149   -0.1177 H       1  LIG1        0.0620
+     23 H           4.7045   -3.4261   -0.0479 H       1  LIG1        0.0624
+     24 H           2.5515   -2.2003    0.0134 H       1  LIG1        0.0626
+@<TRIPOS>BOND
+     1     2     3   ar
+     2     3     4   ar
+     3     2    12   ar
+     4    11    12   ar
+     5    10    11   ar
+     6     5    11   ar
+     7     9    10   ar
+     8     7     9   ar
+     9     6     7   ar
+    10     5     6   ar
+    11     4     5   ar
+    12     1     2    1
+    13     7     8    1
+    14     1    13    1
+    15     1    14    1
+    16     1    15    1
+    17     3    16    1
+    18     4    17    1
+    19     6    18    1
+    20     8    19    1
+    21     8    20    1
+    22     8    21    1
+    23     9    22    1
+    24    10    23    1
+    25    12    24    1
diff --git a/gufe/tests/data/lomap_basic/2-methyl-6-propylnaphthalene.mol2 b/gufe/tests/data/lomap_basic/2-methyl-6-propylnaphthalene.mol2
new file mode 100644
index 00000000..e5595c49
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/2-methyl-6-propylnaphthalene.mol2
@@ -0,0 +1,69 @@
+@<TRIPOS>MOLECULE
+*****
+ 30 31 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.7969    1.3147   -1.5722 C.3     1  LIG1       -0.0650
+      2 C           2.3162    1.3551   -1.5843 C.3     1  LIG1       -0.0519
+      3 C           2.8874    1.3784   -0.1657 C.3     1  LIG1       -0.0279
+      4 C           4.3949    1.4239   -0.1763 C.ar    1  LIG1       -0.0468
+      5 C           5.1419    0.2376   -0.2084 C.ar    1  LIG1       -0.0508
+      6 C           6.5466    0.2630   -0.2482 C.ar    1  LIG1       -0.0176
+      7 C           7.2164    1.5056   -0.2598 C.ar    1  LIG1       -0.0176
+      8 C           6.4634    2.6893   -0.2285 C.ar    1  LIG1       -0.0538
+      9 C           5.0678    2.6494   -0.1901 C.ar    1  LIG1       -0.0580
+     10 C           8.6208    1.5306   -0.3057 C.ar    1  LIG1       -0.0511
+     11 C           9.3673    0.3447   -0.3470 C.ar    1  LIG1       -0.0499
+     12 C           8.6960   -0.8804   -0.3232 C.ar    1  LIG1       -0.0582
+     13 C           7.3000   -0.9207   -0.2780 C.ar    1  LIG1       -0.0538
+     14 C          10.8656    0.4019   -0.3796 C.3     1  LIG1       -0.0397
+     15 H           0.4106    1.2994   -2.5961 H       1  LIG1        0.0230
+     16 H           0.3854    2.1942   -1.0670 H       1  LIG1        0.0230
+     17 H           0.4325    0.4199   -1.0579 H       1  LIG1        0.0230
+     18 H           2.6477    2.2415   -2.1391 H       1  LIG1        0.0266
+     19 H           2.6945    0.4813   -2.1292 H       1  LIG1        0.0266
+     20 H           2.5534    0.4943    0.3922 H       1  LIG1        0.0313
+     21 H           2.5001    2.2454    0.3843 H       1  LIG1        0.0313
+     22 H           4.6215   -0.7184   -0.2065 H       1  LIG1        0.0626
+     23 H           6.9606    3.6568   -0.2384 H       1  LIG1        0.0624
+     24 H           4.5102    3.5829   -0.1747 H       1  LIG1        0.0620
+     25 H           9.1424    2.4858   -0.3125 H       1  LIG1        0.0626
+     26 H           9.2523   -1.8147   -0.3422 H       1  LIG1        0.0620
+     27 H           6.8033   -1.8881   -0.2657 H       1  LIG1        0.0624
+     28 H          11.2977   -0.5759   -0.6162 H       1  LIG1        0.0278
+     29 H          11.2498    0.7216    0.5938 H       1  LIG1        0.0278
+     30 H          11.2028    1.1048   -1.1483 H       1  LIG1        0.0278
+@<TRIPOS>BOND
+     1     4     9   ar
+     2     8     9   ar
+     3     4     5   ar
+     4     5     6   ar
+     5     6    13   ar
+     6     6     7   ar
+     7    12    13   ar
+     8    11    12   ar
+     9    10    11   ar
+    10     7    10   ar
+    11     7     8   ar
+    12     1     2    1
+    13     2     3    1
+    14     3     4    1
+    15    11    14    1
+    16     1    15    1
+    17     1    16    1
+    18     1    17    1
+    19     2    18    1
+    20     2    19    1
+    21     3    20    1
+    22     3    21    1
+    23     5    22    1
+    24     8    23    1
+    25     9    24    1
+    26    10    25    1
+    27    12    26    1
+    28    13    27    1
+    29    14    28    1
+    30    14    29    1
+    31    14    30    1
diff --git a/gufe/tests/data/lomap_basic/2-methylnaphthalene.mol2 b/gufe/tests/data/lomap_basic/2-methylnaphthalene.mol2
new file mode 100644
index 00000000..4f4a6ea5
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/2-methylnaphthalene.mol2
@@ -0,0 +1,51 @@
+@<TRIPOS>MOLECULE
+*****
+ 21 22 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.9750   -0.0348    0.0833 C.3     1  LIG1       -0.0397
+      2 C           2.4749   -0.0532    0.0754 C.ar    1  LIG1       -0.0499
+      3 C           3.2098    1.1400    0.0597 C.ar    1  LIG1       -0.0511
+      4 C           4.6141    1.1264   -0.0127 C.ar    1  LIG1       -0.0176
+      5 C           5.3554    2.3190   -0.0614 C.ar    1  LIG1       -0.0540
+      6 C           6.7465    2.2894   -0.1768 C.ar    1  LIG1       -0.0612
+      7 C           7.4142    1.0702   -0.2321 C.ar    1  LIG1       -0.0612
+      8 C           6.6939   -0.1236   -0.1628 C.ar    1  LIG1       -0.0540
+      9 C           5.2934   -0.1107   -0.0560 C.ar    1  LIG1       -0.0179
+     10 C           4.5528   -1.3018   -0.0057 C.ar    1  LIG1       -0.0538
+     11 C           3.1574   -1.2723    0.0551 C.ar    1  LIG1       -0.0582
+     12 H           0.5850    0.9882    0.1017 H       1  LIG1        0.0278
+     13 H           0.5915   -0.5586    0.9643 H       1  LIG1        0.0278
+     14 H           0.5924   -0.5209   -0.8204 H       1  LIG1        0.0278
+     15 H           2.6835    2.0917    0.0844 H       1  LIG1        0.0626
+     16 H           4.8512    3.2813   -0.0218 H       1  LIG1        0.0624
+     17 H           7.3059    3.2198   -0.2301 H       1  LIG1        0.0618
+     18 H           8.4962    1.0459   -0.3322 H       1  LIG1        0.0618
+     19 H           7.2331   -1.0669   -0.2063 H       1  LIG1        0.0624
+     20 H           5.0575   -2.2644   -0.0289 H       1  LIG1        0.0624
+     21 H           2.6062   -2.2099    0.0730 H       1  LIG1        0.0620
+@<TRIPOS>BOND
+     1     1     2    1
+     2     2     3   ar
+     3     3     4   ar
+     4     4     5   ar
+     5     5     6   ar
+     6     6     7   ar
+     7     7     8   ar
+     8     8     9   ar
+     9     4     9   ar
+    10     9    10   ar
+    11    10    11   ar
+    12     2    11   ar
+    13     1    12    1
+    14     1    13    1
+    15     1    14    1
+    16     3    15    1
+    17     5    16    1
+    18     6    17    1
+    19     7    18    1
+    20     8    19    1
+    21    10    20    1
+    22    11    21    1
diff --git a/gufe/tests/data/lomap_basic/2-naftanol.mol2 b/gufe/tests/data/lomap_basic/2-naftanol.mol2
new file mode 100644
index 00000000..79bfc04e
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/2-naftanol.mol2
@@ -0,0 +1,47 @@
+@<TRIPOS>MOLECULE
+*****
+ 19 20 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 O           9.9477   -4.1875    1.4569 O.3     1  LIG1       -0.5067
+      2 C           9.0676   -3.1850    1.1716 C.ar    1  LIG1        0.1175
+      3 C           7.7820   -3.1808    1.6974 C.ar    1  LIG1       -0.0196
+      4 C           6.9153   -2.1329    1.3798 C.ar    1  LIG1       -0.0506
+      5 C           7.3294   -1.0867    0.5375 C.ar    1  LIG1       -0.0177
+      6 C           6.4703   -0.0252    0.2106 C.ar    1  LIG1       -0.0540
+      7 C           6.8977    1.0066   -0.6271 C.ar    1  LIG1       -0.0612
+      8 C           8.1878    0.9906   -1.1489 C.ar    1  LIG1       -0.0612
+      9 C           9.0548   -0.0574   -0.8344 C.ar    1  LIG1       -0.0539
+     10 C           8.6411   -1.1023    0.0071 C.ar    1  LIG1       -0.0144
+     11 C           9.5008   -2.1625    0.3346 C.ar    1  LIG1       -0.0125
+     12 H           9.5144   -4.8197    2.0534 H       1  LIG1        0.2921
+     13 H           7.4363   -3.9757    2.3505 H       1  LIG1        0.0654
+     14 H           5.9102   -2.1360    1.7951 H       1  LIG1        0.0625
+     15 H           5.4583    0.0059    0.6080 H       1  LIG1        0.0624
+     16 H           6.2223    1.8225   -0.8713 H       1  LIG1        0.0618
+     17 H           8.5216    1.7931   -1.8011 H       1  LIG1        0.0618
+     18 H          10.0593   -0.0532   -1.2512 H       1  LIG1        0.0624
+     19 H          10.5134   -2.1967   -0.0603 H       1  LIG1        0.0660
+@<TRIPOS>BOND
+     1     2     3   ar
+     2     2    11   ar
+     3     3     4   ar
+     4     4     5   ar
+     5     5     6   ar
+     6     5    10   ar
+     7     6     7   ar
+     8     7     8   ar
+     9     8     9   ar
+    10     9    10   ar
+    11    10    11   ar
+    12     1     2    1
+    13     1    12    1
+    14     3    13    1
+    15     4    14    1
+    16     6    15    1
+    17     7    16    1
+    18     8    17    1
+    19     9    18    1
+    20    11    19    1
diff --git a/gufe/tests/data/lomap_basic/README.md b/gufe/tests/data/lomap_basic/README.md
new file mode 100644
index 00000000..5e4485d7
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/README.md
@@ -0,0 +1,1 @@
+Testfiles from Lomap's basic tests
diff --git a/gufe/tests/data/lomap_basic/__init__.py b/gufe/tests/data/lomap_basic/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/gufe/tests/data/lomap_basic/methylcyclohexane.mol2 b/gufe/tests/data/lomap_basic/methylcyclohexane.mol2
new file mode 100644
index 00000000..926963c2
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/methylcyclohexane.mol2
@@ -0,0 +1,50 @@
+@<TRIPOS>MOLECULE
+*****
+ 21 21 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.9730    0.1554    0.1412 C.3     1  LIG1       -0.0624
+      2 C           2.4934    0.1043   -0.0303 C.3     1  LIG1       -0.0439
+      3 C           2.9737   -1.3048   -0.4072 C.3     1  LIG1       -0.0505
+      4 C           2.6388   -1.6853   -1.8487 C.3     1  LIG1       -0.0528
+      5 C           3.1458   -0.6424   -2.8397 C.3     1  LIG1       -0.0530
+      6 C           2.6503    0.7581   -2.4959 C.3     1  LIG1       -0.0528
+      7 C           2.9839    1.1382   -1.0536 C.3     1  LIG1       -0.0505
+      8 H           0.6557    1.1542    0.4607 H       1  LIG1        0.0232
+      9 H           0.6478   -0.5562    0.9081 H       1  LIG1        0.0232
+     10 H           0.4392   -0.0840   -0.7835 H       1  LIG1        0.0232
+     11 H           2.9316    0.3596    0.9431 H       1  LIG1        0.0298
+     12 H           2.5484   -2.0480    0.2778 H       1  LIG1        0.0268
+     13 H           4.0632   -1.3509   -0.2778 H       1  LIG1        0.0268
+     14 H           3.0861   -2.6588   -2.0822 H       1  LIG1        0.0265
+     15 H           1.5550   -1.8040   -1.9616 H       1  LIG1        0.0265
+     16 H           4.2430   -0.6482   -2.8423 H       1  LIG1        0.0265
+     17 H           2.8244   -0.9091   -3.8522 H       1  LIG1        0.0265
+     18 H           3.1052    1.4851   -3.1790 H       1  LIG1        0.0265
+     19 H           1.5667    0.8148   -2.6555 H       1  LIG1        0.0265
+     20 H           4.0736    1.2348   -0.9613 H       1  LIG1        0.0268
+     21 H           2.5648    2.1260   -0.8271 H       1  LIG1        0.0268
+@<TRIPOS>BOND
+     1     2     7    1
+     2     2     3    1
+     3     6     7    1
+     4     5     6    1
+     5     4     5    1
+     6     3     4    1
+     7     1     2    1
+     8     1     8    1
+     9     1     9    1
+    10     1    10    1
+    11     2    11    1
+    12     3    12    1
+    13     3    13    1
+    14     4    14    1
+    15     4    15    1
+    16     5    16    1
+    17     5    17    1
+    18     6    18    1
+    19     6    19    1
+    20     7    20    1
+    21     7    21    1
diff --git a/gufe/tests/data/lomap_basic/toluene.mol2 b/gufe/tests/data/lomap_basic/toluene.mol2
new file mode 100644
index 00000000..472586c6
--- /dev/null
+++ b/gufe/tests/data/lomap_basic/toluene.mol2
@@ -0,0 +1,38 @@
+@<TRIPOS>MOLECULE
+*****
+ 15 15 0 0 0
+SMALL
+GASTEIGER
+
+@<TRIPOS>ATOM
+      1 C           0.9562    0.0874    0.1225 C.3     1  LIG1       -0.0397
+      2 C           2.4532    0.1169    0.0327 C.ar    1  LIG1       -0.0504
+      3 C           3.1480    1.3313   -0.0192 C.ar    1  LIG1       -0.0588
+      4 C           4.5410    1.3447   -0.1106 C.ar    1  LIG1       -0.0615
+      5 C           5.2502    0.1457   -0.1538 C.ar    1  LIG1       -0.0617
+      6 C           4.5672   -1.0681   -0.1065 C.ar    1  LIG1       -0.0615
+      7 C           3.1745   -1.0830   -0.0148 C.ar    1  LIG1       -0.0588
+      8 H           0.5405    1.0944    0.2295 H       1  LIG1        0.0278
+      9 H           0.6412   -0.4972    0.9928 H       1  LIG1        0.0278
+     10 H           0.5343   -0.3590   -0.7834 H       1  LIG1        0.0278
+     11 H           2.6097    2.2757    0.0116 H       1  LIG1        0.0620
+     12 H           5.0726    2.2915   -0.1488 H       1  LIG1        0.0618
+     13 H           6.3343    0.1573   -0.2252 H       1  LIG1        0.0618
+     14 H           5.1186   -2.0039   -0.1418 H       1  LIG1        0.0618
+     15 H           2.6530   -2.0368    0.0195 H       1  LIG1        0.0620
+@<TRIPOS>BOND
+     1     2     3   ar
+     2     3     4   ar
+     3     2     7   ar
+     4     6     7   ar
+     5     5     6   ar
+     6     4     5   ar
+     7     1     2    1
+     8     1     8    1
+     9     1     9    1
+    10     1    10    1
+    11     3    11    1
+    12     4    12    1
+    13     5    13    1
+    14     6    14    1
+    15     7    15    1
diff --git a/gufe/tests/test_ligandatommapping.py b/gufe/tests/test_ligandatommapping.py
new file mode 100644
index 00000000..158be6ae
--- /dev/null
+++ b/gufe/tests/test_ligandatommapping.py
@@ -0,0 +1,222 @@
+# This code is part of gufe and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/gufe
+import importlib
+import pytest
+import pathlib
+import json
+from rdkit import Chem
+
+import gufe
+from gufe import LigandAtomMapping, SmallMoleculeComponent
+
+
+def mol_from_smiles(smiles: str) -> Chem.Mol:
+    m = Chem.MolFromSmiles(smiles)
+    m.Compute2DCoords()
+
+    return m
+
+
+@pytest.fixture(scope='session')
+def simple_mapping():
+    """Disappearing oxygen on end
+
+    C C O
+
+    C C
+    """
+    molA = gufe.SmallMoleculeComponent(mol_from_smiles('CCO'))
+    molB = gufe.SmallMoleculeComponent(mol_from_smiles('CC'))
+
+    m = LigandAtomMapping(molA, molB, componentA_to_componentB={0: 0, 1: 1})
+
+    return m
+
+
+@pytest.fixture(scope='session')
+def other_mapping():
+    """Disappearing middle carbon
+
+    C C O
+
+    C   C
+    """
+    molA = SmallMoleculeComponent(mol_from_smiles('CCO'))
+    molB = SmallMoleculeComponent(mol_from_smiles('CC'))
+
+    m = LigandAtomMapping(molA, molB, componentA_to_componentB={0: 0, 2: 1})
+
+    return m
+
+
+@pytest.fixture
+def annotated_simple_mapping(simple_mapping):
+    mapping = LigandAtomMapping(simple_mapping.componentA,
+                                simple_mapping.componentB,
+                                simple_mapping.componentA_to_componentB,
+                                annotations={'foo': 'bar'})
+    return mapping
+
+
+@pytest.fixture(scope='session')
+def benzene_maps():
+    MAPS = {
+        'phenol': {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6,
+                   7: 7, 8: 8, 9: 9, 10: 12, 11: 11},
+        'anisole': {0: 5, 1: 6, 2: 7, 3: 8, 4: 9, 5: 10,
+                    6: 11, 7: 12, 8: 13, 9: 14, 10: 2, 11: 15}}
+    return MAPS
+
+
+@pytest.fixture(scope='session')
+def benzene_phenol_mapping(benzene_transforms, benzene_maps):
+    molA = SmallMoleculeComponent(benzene_transforms['benzene'].to_rdkit())
+    molB = SmallMoleculeComponent(benzene_transforms['phenol'].to_rdkit())
+    m = LigandAtomMapping(molA, molB, benzene_maps['phenol'])
+    return m
+
+
+@pytest.fixture(scope='session')
+def benzene_anisole_mapping(benzene_transforms, benzene_maps):
+    molA = SmallMoleculeComponent(benzene_transforms['benzene'].to_rdkit())
+    molB = SmallMoleculeComponent(benzene_transforms['anisole'].to_rdkit())
+    m = LigandAtomMapping(molA, molB, benzene_maps['anisole'])
+    return m
+
+
+@pytest.fixture(scope='session')
+def atom_mapping_basic_test_files():
+    # a dict of {filenames.strip(mol2): SmallMoleculeComponent} for a simple
+    # set of ligands
+    files = {}
+    for f in [
+        '1,3,7-trimethylnaphthalene',
+        '1-butyl-4-methylbenzene',
+        '2,6-dimethylnaphthalene',
+        '2-methyl-6-propylnaphthalene',
+        '2-methylnaphthalene',
+        '2-naftanol',
+        'methylcyclohexane',
+        'toluene']:
+        with importlib.resources.path('gufe.tests.data.lomap_basic',
+                                      f + '.mol2') as fn:
+            mol = Chem.MolFromMol2File(str(fn), removeHs=False)
+            files[f] = SmallMoleculeComponent(mol, name=f)
+
+    return files
+
+
+def test_atommapping_usage(simple_mapping):
+    assert simple_mapping.componentA_to_componentB[1] == 1
+    assert simple_mapping.componentA_to_componentB.get(2, None) is None
+    assert simple_mapping.annotations == {}
+
+    with pytest.raises(KeyError):
+        simple_mapping.componentA_to_componentB[3]
+
+
+def test_mapping_inversion(benzene_phenol_mapping):
+    assert benzene_phenol_mapping.componentB_to_componentA == {
+        0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7, 8: 8, 9: 9, 11: 11,
+        12: 10
+    }
+
+
+def test_uniques(atom_mapping_basic_test_files):
+    mapping = LigandAtomMapping(
+        componentA=atom_mapping_basic_test_files['methylcyclohexane'],
+        componentB=atom_mapping_basic_test_files['toluene'],
+        componentA_to_componentB={
+            0: 6, 1: 7, 2: 8, 3: 9, 4: 10, 5: 11, 6: 12
+        }
+    )
+
+    assert list(mapping.componentA_unique) == [7, 8, 9, 10, 11, 12, 13, 14, 15,
+                                               16, 17, 18, 19, 20]
+    assert list(mapping.componentB_unique) == [0, 1, 2, 3, 4, 5, 13, 14]
+
+
+def test_modification(benzene_phenol_mapping):
+    # check that we get a copy of the mapping and we can't modify
+    AtoB = benzene_phenol_mapping.componentA_to_componentB
+    before = len(AtoB)
+
+    AtoB.pop(10)
+
+    assert len(benzene_phenol_mapping.componentA_to_componentB) == before
+
+
+def test_atommapping_hash(simple_mapping, other_mapping):
+    # these two mappings map the same molecules, but with a different mapping
+    assert simple_mapping is not other_mapping
+
+
+def test_draw_mapping_cairo(tmpdir, simple_mapping):
+    with tmpdir.as_cwd():
+        simple_mapping.draw_to_file('test.png')
+        filed = pathlib.Path('test.png')
+        assert filed.exists()
+
+
+def test_draw_mapping_svg(tmpdir, other_mapping):
+    with tmpdir.as_cwd():
+        d2d = Chem.Draw.rdMolDraw2D.MolDraw2DSVG(600, 300, 300, 300)
+        other_mapping.draw_to_file('test.svg', d2d=d2d)
+        filed = pathlib.Path('test.svg')
+        assert filed.exists()
+
+
+class TestLigandAtomMappingSerialization:
+    def test_to_dict(self, benzene_phenol_mapping):
+        d = benzene_phenol_mapping.to_dict()
+
+        assert isinstance(d, dict)
+        assert 'componentA' in d
+        assert 'componentB' in d
+        assert 'annotations' in d
+        assert isinstance(d['componentA'], str)
+
+    def test_deserialize_roundtrip(self, benzene_phenol_mapping,
+                                   benzene_anisole_mapping):
+
+        roundtrip = LigandAtomMapping.from_dict(
+                        benzene_phenol_mapping.to_dict())
+
+        assert roundtrip == benzene_phenol_mapping
+
+        # We don't check coordinates since that's already done in guefe for
+        # SmallMoleculeComponent
+
+        assert roundtrip != benzene_anisole_mapping
+
+    def test_file_roundtrip(self, benzene_phenol_mapping, tmpdir):
+        with tmpdir.as_cwd():
+            with open('tmpfile.json', 'w') as f:
+                f.write(json.dumps(benzene_phenol_mapping.to_dict()))
+
+            with open('tmpfile.json', 'r') as f:
+                d = json.load(f)
+
+            assert isinstance(d, dict)
+            roundtrip = LigandAtomMapping.from_dict(d)
+
+            assert roundtrip == benzene_phenol_mapping
+
+
+def test_annotated_atommapping_hash_eq(simple_mapping,
+                                       annotated_simple_mapping):
+    assert annotated_simple_mapping != simple_mapping
+    assert hash(annotated_simple_mapping) != hash(simple_mapping)
+
+
+def test_annotation_immutability(annotated_simple_mapping):
+    annot1 = annotated_simple_mapping.annotations
+    annot1['foo'] = 'baz'
+    annot2 = annotated_simple_mapping.annotations
+    assert annot1 != annot2
+    assert annot2 == {'foo': 'bar'}
+
+
+def test_with_annotations(simple_mapping, annotated_simple_mapping):
+    new_annot = simple_mapping.with_annotations({'foo': 'bar'})
+    assert new_annot == annotated_simple_mapping
diff --git a/gufe/tests/test_mapping_visualization.py b/gufe/tests/test_mapping_visualization.py
new file mode 100644
index 00000000..8da2d013
--- /dev/null
+++ b/gufe/tests/test_mapping_visualization.py
@@ -0,0 +1,188 @@
+import pytest
+from unittest import mock
+import inspect
+
+from rdkit import Chem
+
+import gufe
+from gufe.visualization.mapping_visualization import (
+    _match_elements, _get_unique_bonds_and_atoms, draw_mapping,
+    draw_one_molecule_mapping, draw_unhighlighted_molecule
+)
+
+# default colors currently used
+_HIGHLIGHT_COLOR = (220/255, 50/255, 32/255, 1)
+_CHANGED_ELEMENTS_COLOR = (0, 90/255, 181/255, 1)
+
+
+def bound_args(func, args, kwargs):
+    """Return a dictionary mapping parameter name to value.
+
+    Parameters
+    ----------
+    func : Callable
+        this must be inspectable; mocks will require a spec
+    args : List
+        args list
+    kwargs : Dict
+        kwargs Dict
+
+    Returns
+    -------
+    Dict[str, Any] :
+        mapping of string name of function parameter to the value it would
+        be bound to
+    """
+    sig = inspect.Signature.from_callable(func)
+    bound = sig.bind(*args, **kwargs)
+    return bound.arguments
+
+
+@pytest.mark.parametrize("at1, idx1, at2, idx2, response", [
+    ["N", 0, "C", 0, False],
+    ["C", 0, "C", 0, True],
+    ["COC", 1, "NOC", 1, True],
+    ["COON", 2, "COC", 2, False]]
+)
+def test_match_elements(at1, idx1, at2, idx2, response):
+    mol1 = Chem.MolFromSmiles(at1)
+    mol2 = Chem.MolFromSmiles(at2)
+    retval = _match_elements(mol1, idx2, mol2, idx2)
+
+    assert retval == response
+
+
+@pytest.fixture(scope='module')
+def maps():
+    MAPS = {
+        'phenol': {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6,
+                   7: 7, 8: 8, 9: 9, 10: 12, 11: 11},
+        'anisole': {0: 5, 1: 6, 2: 7, 3: 8, 4: 9, 5: 10,
+                    6: 11, 7: 12, 8: 13, 9: 14, 10: 2, 11: 15}}
+    return MAPS
+
+
+@pytest.fixture(scope='module')
+def benzene_phenol_mapping(benzene_transforms, maps):
+    mol1 = benzene_transforms['benzene'].to_rdkit()
+    mol2 = benzene_transforms['phenol'].to_rdkit()
+    mapping = maps['phenol']
+    return mapping, mol1, mol2
+
+
+@pytest.mark.parametrize('molname, atoms, elems, bonds', [
+    ['phenol', {10, }, {12, }, {10, 12}],
+    ['anisole', {0, 1, 3, 4}, {2, }, {0, 1, 2, 3, 13}]
+])
+def test_benzene_to_phenol_uniques(molname, atoms, elems, bonds,
+                                   benzene_transforms, maps):
+    mol1 = benzene_transforms['benzene']
+    mol2 = benzene_transforms[molname]
+
+    mapping = maps[molname]
+
+    uniques = _get_unique_bonds_and_atoms(mapping,
+                                          mol1.to_rdkit(), mol2.to_rdkit())
+
+    # The benzene perturbations don't change
+    # no unique atoms in benzene
+    assert uniques['atoms'] == set()
+    # H->O
+    assert uniques['elements'] == {10, }
+    # One bond involved
+    assert uniques['bonds'] == {10, }
+
+    # invert and check the molB uniques
+    inv_map = {v: k for k, v in mapping.items()}
+
+    uniques = _get_unique_bonds_and_atoms(inv_map,
+                                          mol2.to_rdkit(), mol1.to_rdkit())
+
+    assert uniques['atoms'] == atoms
+    assert uniques['elements'] == elems
+    assert uniques['bonds'] == bonds
+
+
+@mock.patch("gufe.visualization.mapping_visualization._draw_molecules", autospec=True)
+def test_draw_mapping(mock_func, benzene_phenol_mapping):
+    # ensure that draw_mapping passes the desired parameters to our internal
+    # _draw_molecules method
+    mapping, mol1, mol2 = benzene_phenol_mapping
+    draw_mapping(mapping, mol1, mol2)
+
+    mock_func.assert_called_once()
+    args = bound_args(mock_func, mock_func.call_args.args,
+                      mock_func.call_args.kwargs)
+    assert args['mols'] == [mol1, mol2]
+    assert args['atoms_list'] == [{10}, {10, 12}]
+    assert args['bonds_list'] == [{10}, {10, 12}]
+    assert args['atom_colors'] == [{10: _CHANGED_ELEMENTS_COLOR},
+                                   {12: _CHANGED_ELEMENTS_COLOR}]
+    assert args['highlight_color'] == _HIGHLIGHT_COLOR
+
+
+@pytest.mark.parametrize('inverted', [True, False])
+@mock.patch("gufe.visualization.mapping_visualization._draw_molecules", autospec=True)
+def test_draw_one_molecule_mapping(mock_func, benzene_phenol_mapping,
+                                   inverted):
+    # ensure that draw_one_molecule_mapping passes the desired parameters to
+    # our internal _draw_molecules method
+    mapping, mol1, mol2 = benzene_phenol_mapping
+    if inverted:
+        mapping = {v: k for k, v in mapping.items()}
+        mol1, mol2 = mol2, mol1
+        atoms_list = [{10, 12}]
+        bonds_list = [{10, 12}]
+        atom_colors = [{12: _CHANGED_ELEMENTS_COLOR}]
+    else:
+        atoms_list = [{10}]
+        bonds_list = [{10}]
+        atom_colors = [{10: _CHANGED_ELEMENTS_COLOR}]
+
+    draw_one_molecule_mapping(mapping, mol1, mol2)
+
+    mock_func.assert_called_once()
+    args = bound_args(mock_func, mock_func.call_args.args,
+                      mock_func.call_args.kwargs)
+
+    assert args['mols'] == [mol1]
+    assert args['atoms_list'] == atoms_list
+    assert args['bonds_list'] == bonds_list
+    assert args['atom_colors'] == atom_colors
+    assert args['highlight_color'] == _HIGHLIGHT_COLOR
+
+
+@mock.patch("gufe.visualization.mapping_visualization._draw_molecules", autospec=True)
+def test_draw_unhighlighted_molecule(mock_func, benzene_transforms):
+    # ensure that draw_unhighlighted_molecule passes the desired parameters
+    # to our internal _draw_molecules method
+    mol = benzene_transforms['benzene'].to_rdkit()
+    draw_unhighlighted_molecule(mol)
+
+    mock_func.assert_called_once()
+    args = bound_args(mock_func, mock_func.call_args.args,
+                      mock_func.call_args.kwargs)
+    assert args['mols'] == [mol]
+    assert args['atoms_list'] == [[]]
+    assert args['bonds_list'] == [[]]
+    assert args['atom_colors'] == [{}]
+    # technically, we don't care what the highlight color is, so no
+    # assertion on that
+
+
+def test_draw_mapping_integration_smoke(benzene_phenol_mapping):
+    # integration test/smoke test to catch errors if the upstream drawing
+    # code changes
+    draw_mapping(*benzene_phenol_mapping)
+
+
+def test_draw_one_molecule_integration_smoke(benzene_phenol_mapping):
+    # integration test/smoke test to catch errors if the upstream drawing
+    # code changes
+    draw_one_molecule_mapping(*benzene_phenol_mapping)
+
+
+def test_draw_unhighlighted_molecule_integration_smoke(benzene_transforms):
+    # integration test/smoke test to catch errors if the upstream drawing
+    # code changes
+    draw_unhighlighted_molecule(benzene_transforms['benzene'].to_rdkit())