diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 4a11101..bb625bb 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -1,118 +1,117 @@
 name: Tests
 on:
-  push:
-    branches:
-      - master
-  pull_request:
-  workflow_dispatch:
+    push:
+        branches:
+            - master
+    pull_request:
+    workflow_dispatch:
 jobs:
-  ruff:
-    runs-on: ubuntu-22.04
-    steps:
-      - uses: actions/checkout@v3
-      - uses: actions/setup-python@v4
-        with:
-          python-version: "3.11"
-          cache: "pip"
-      - run: "pip install '.[dev]'"
-      - run: ruff check .
+    ruff:
+        runs-on: ubuntu-22.04
+        steps:
+            - uses: actions/checkout@v3
+            - uses: actions/setup-python@v4
+              with:
+                  python-version: "3.11"
+                  cache: "pip"
+            - run: "pip install '.[dev]'"
+            - run: ruff check .
 
-  mypy:
-    runs-on: ubuntu-22.04
-    steps:
-      - uses: actions/checkout@v3
-      - uses: actions/setup-python@v4
-        with:
-          python-version: "3.11"
-          cache: "pip"
-      - run: "pip install -e '.[dev]'"
-      - run: mypy src examples
+    mypy:
+        runs-on: ubuntu-22.04
+        steps:
+            - uses: actions/checkout@v3
+            - uses: actions/setup-python@v4
+              with:
+                  python-version: "3.11"
+                  cache: "pip"
+            - run: "pip install -e '.[dev]'"
+            - run: mypy src examples
 
-  ruff-format:
-    runs-on: ubuntu-22.04
-    steps:
-      - uses: actions/checkout@v3
-      - uses: actions/setup-python@v4
-        with:
-          python-version: "3.11"
-          cache: "pip"
-      - run: "pip install -e '.[dev]'"
-      - run: ruff format --check .
+    ruff-format:
+        runs-on: ubuntu-22.04
+        steps:
+            - uses: actions/checkout@v3
+            - uses: actions/setup-python@v4
+              with:
+                  python-version: "3.11"
+                  cache: "pip"
+            - run: "pip install -e '.[dev]'"
+            - run: ruff format --check .
 
-  pytest:
-    # https://ericmjl.github.io/blog/2021/12/30/better-conda-environments-on-github-actions/
-    runs-on: ubuntu-22.04
-    defaults:
-      run:
-        shell: bash -l {0}
-    steps:
-      - uses: actions/checkout@v3
-      - uses: conda-incubator/setup-miniconda@v2
-        with:
-          auto-update-conda: true
-          miniforge-variant: Mambaforge
-          channels: conda-forge,defaults
-          python-version: 3.11
-          activate-environment: pytest
-          use-mamba: true
-      - name: Build environment
-        run: |
-          conda activate pytest
-          conda install -c conda-forge openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
-          conda install -c conda-forge mdanalysis openbabel
-          python -m pip install --upgrade pip
-          python -m pip install --upgrade setuptools
-          python -m pip install -e '.[dev]'
-      - run: pytest
+    pytest:
+        # https://ericmjl.github.io/blog/2021/12/30/better-conda-environments-on-github-actions/
+        runs-on: ubuntu-22.04
+        defaults:
+            run:
+                shell: bash -l {0}
 
-  example-tests:
-    # https://ericmjl.github.io/blog/2021/12/30/better-conda-environments-on-github-actions/
-    runs-on: ubuntu-22.04
-    defaults:
-      run:
-        shell: bash -l {0}
-    steps:
-      - uses: actions/checkout@v3
-      - uses: conda-incubator/setup-miniconda@v2
-        with:
-          auto-update-conda: true
-          miniforge-variant: Mambaforge
-          channels: conda-forge,defaults
-          python-version: 3.11
-          activate-environment: pytest
-          use-mamba: true
-      - name: Build environment
-        run: |
-          conda activate pytest
-          conda install -c conda-forge openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
-          conda install -c conda-forge mdanalysis openbabel
-          python -m pip install --upgrade pip
-          python -m pip install --upgrade setuptools
-          python -m pip install -e '.[dev]'
-      - name: Run script
-        run: python examples/testable_example.py
+        steps:
+            - uses: actions/checkout@v4
+            - uses: conda-incubator/setup-miniconda@v3
 
-  doctest:
-    runs-on: ubuntu-22.04
-    defaults:
-      run:
-        shell: bash -l {0}
-    steps:
-      - uses: actions/checkout@v3
-      - uses: conda-incubator/setup-miniconda@v2
-        with:
-          auto-update-conda: true
-          miniforge-variant: Mambaforge
-          channels: conda-forge,defaults
-          python-version: 3.11
-          activate-environment: pytest
-          use-mamba: true
-      - name: Build environment
-        run: |
-          conda activate pytest
-          conda install -c conda-forge openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
-          conda install -c conda-forge mdanalysis openbabel
-          python -m pip install --upgrade pip
-          python -m pip install --upgrade setuptools
-          python -m pip install -e '.[dev]'
-      - run: make -C docs doctest
+              with:
+                  miniforge-version: latest
+                  python-version: 3.11
+                  activate-environment: pytest
+
+            - name: Build environment
+              run: |
+                  conda activate pytest
+                  conda install -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
+                  conda install -c conda-forge mdanalysis openbabel
+                  python -m pip install --upgrade pip
+                  python -m pip install --upgrade setuptools
+                  python -m pip install -e '.[dev]'
+            - run: pytest
+
+    example-tests:
+        # https://ericmjl.github.io/blog/2021/12/30/better-conda-environments-on-github-actions/
+        runs-on: ubuntu-22.04
+        defaults:
+            run:
+                shell: bash -l {0}
+        steps:
+            - uses: actions/checkout@v4
+            - uses: conda-incubator/setup-miniconda@v3
+
+              with:
+                  miniforge-version: latest
+                  python-version: 3.11
+                  activate-environment: pytest
+
+            - name: Build environment
+              run: |
+                  conda activate pytest
+                  conda install -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
+                  conda install -c conda-forge mdanalysis openbabel
+                  python -m pip install --upgrade pip
+                  python -m pip install --upgrade setuptools
+                  python -m pip install -e '.[dev]'
+            - name: Run script
+              run: python examples/testable_example.py
+
+    doctest:
+        runs-on: ubuntu-22.04
+        defaults:
+            run:
+                shell: bash -l {0}
+
+        steps:
+            - uses: actions/checkout@v4
+            - uses: conda-incubator/setup-miniconda@v3
+
+              with:
+                  miniforge-version: latest
+                  python-version: 3.11
+                  activate-environment: pytest
+
+            - name: Build environment
+              run: |
+                  conda activate pytest
+                  conda install -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
+                  conda install -c conda-forge mdanalysis openbabel
+                  python -m pip install --upgrade pip
+                  python -m pip install --upgrade setuptools
+                  python -m pip install -e '.[dev]'
+            - run: make -C docs doctest
diff --git a/.gitignore b/.gitignore
index 5c34184..b18dc7e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,3 +22,4 @@ tests/failures.txt
 /examples/.model.pt
 
 /examples/*/md_optimisation
+.model.pt
diff --git a/README.rst b/README.rst
index 8de1d7e..2676908 100644
--- a/README.rst
+++ b/README.rst
@@ -35,7 +35,7 @@ Some optional dependencies are only available through conda:
   # for OpenMM and espaloma charge
   # note the temporary issue with rdkit versions and conda will overwrite pip
   # installed software
-  mamba install -c conda-forge openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
+  mamba install -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools dgl rdkit==2024.3.4 espaloma_charge
   # for xtb
   mamba install xtb
   # for openbabel
@@ -43,6 +43,24 @@ Some optional dependencies are only available through conda:
   # for mdanalysis
   mamba install mdanalysis
 
+Developer Setup
+---------------
+
+1. Install `just`_.
+2. In a new virtual environment run:
+
+.. code-block:: bash
+
+  just dev
+
+3. Run code checks:
+
+.. code-block:: bash
+
+  just check
+
+.. _`just`: https://github.com/casey/just
+
 Examples
 ========
 
diff --git a/docs/source/cage_analysis.rst b/docs/source/cage_analysis.rst
index a094576..7d3dd64 100644
--- a/docs/source/cage_analysis.rst
+++ b/docs/source/cage_analysis.rst
@@ -257,7 +257,7 @@ We can get measures of pore size and cage geometry.
 .. testcode:: analysing-cage
     :hide:
 
-    assert analyser.get_min_centroid_distance(apdcage) == 6.612215150137052
+    assert np.isclose(analyser.get_min_centroid_distance(apdcage), 6.612215150137052)
 
 Giving:
 
@@ -274,4 +274,4 @@ Giving:
 
     import shutil
 
-    shutil.rmtree('cage_output')
\ No newline at end of file
+    shutil.rmtree('cage_output')
diff --git a/docs/source/index.rst b/docs/source/index.rst
index d5fcbb0..d4539a0 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -23,8 +23,12 @@
 Welcome to stko's documentation!
 ================================
 
-GitHub: https://github.com/JelfsMaterialsGroup/stko
+| GitHub: https://github.com/JelfsMaterialsGroup/stko
+| Discord: https://discord.gg/zbCUzuxe2B
 
+.. tip::
+
+  ⭐ Star us on `GitHub <https://github.com/JelfsMaterialsGroup/stko>`_! ⭐
 
 Install
 =======
@@ -42,7 +46,7 @@ Some optional dependencies are only available through conda:
   # for OpenMM and espaloma charge
   # note the temporary issue with rdkit versions and conda will overwrite pip
   # installed software
-  mamba install -c conda-forge openff-toolkit openmm openmmtools dgl rdkit==2024.3.4  espaloma_charge
+  mamba install -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools dgl rdkit==2024.3.4  espaloma_charge
   # for xtb
   mamba install xtb
   # for openbabel
@@ -50,6 +54,19 @@ Some optional dependencies are only available through conda:
   # for mdanalysis
   mamba install mdanalysis
 
+Developer Setup
+---------------
+
+#. Install `just`_.
+#. In a new virtual environment run::
+
+    $ just dev
+
+#. Run code checks::
+
+    $ just check
+
+.. _`just`: https://github.com/casey/just
 
 Dependencies
 ------------
diff --git a/docs/source/molecular.rst b/docs/source/molecular.rst
index dfb1d50..9ca3a4e 100644
--- a/docs/source/molecular.rst
+++ b/docs/source/molecular.rst
@@ -16,6 +16,16 @@ Atoms, Functional Groups and Molecules
   ThreeSiteFactory <_autosummary/stko.functional_groups.ThreeSiteFactory>
 
 
+Molecular utilities
+-------------------
+
+.. toctree::
+  :maxdepth: 1
+
+  separate_molecule <_autosummary/stko.molecular_utilities.separate_molecule>
+  merge_stk_molecules <_autosummary/stko.molecular_utilities.merge_stk_molecules>
+  update_stk_from_rdkit_conformer <_autosummary/stko.molecular_utilities.update_stk_from_rdkit_conformer>
+
 Molecular Analysis
 ------------------
 
diff --git a/examples/cage_analysis_example.py b/examples/cage_analysis_example.py
index 29a6678..8ce9d67 100644
--- a/examples/cage_analysis_example.py
+++ b/examples/cage_analysis_example.py
@@ -82,20 +82,14 @@ def main() -> None:
     # Distance between binders in the organic linkers?
     print(f"avg. binder distance: {np.mean(ligand_dict['binder_dist'])}")
     # How twisted is the molecule, or aligned are the binding groups?
-    print(
-        "avg. binder binder angle: " f"{np.mean(ligand_dict['binder_binder'])}"
-    )
+    print(f"avg. binder binder angle: {np.mean(ligand_dict['binder_binder'])}")
     # How twisted is the molecule, or what is the torsion between
     # the binding groups?
-    print(
-        "avg. binder adjacent torsion: " f"{np.mean(ligand_dict['torsion'])}"
-    )
+    print(f"avg. binder adjacent torsion: {np.mean(ligand_dict['torsion'])}")
     # What is the angle made by the binders?
     print(f"avg. binder angles: {np.mean(ligand_dict['binder_angle'])}")
     # And the resultant bite-angle [caution!]?
-    print(
-        "avg. bite angle [caution]: " f"{np.mean(ligand_dict['bite_angle'])}"
-    )
+    print(f"avg. bite angle [caution]: {np.mean(ligand_dict['bite_angle'])}")
 
     # Get the centroid and atom ids of distinct building blocks.
     # This is similar to the above, but does not perform any disconnections
@@ -133,11 +127,11 @@ def main() -> None:
     # And some geometrical measures.
     print(
         f"avg. N-Pd bond length:"
-        f' {np.mean(analyser.calculate_bonds(apdcage)[("N", "Pd")])}'
+        f" {np.mean(analyser.calculate_bonds(apdcage)[('N', 'Pd')])}"
     )
     print(
         f"N-Pd-N angles: "
-        f'{analyser.calculate_angles(apdcage)[("N", "Pd", "N")]}'
+        f"{analyser.calculate_angles(apdcage)[('N', 'Pd', 'N')]}"
     )
 
 
diff --git a/examples/obabel_example.py b/examples/obabel_example.py
index 18b773c..1d295ee 100644
--- a/examples/obabel_example.py
+++ b/examples/obabel_example.py
@@ -74,7 +74,7 @@ def main() -> None:
 
     # Define spacer building block.
     bb3 = stk.BuildingBlock(
-        smiles=("C1=CC(C2=CC=C(Br)C=C2)=C" "C=C1Br"),
+        smiles=("C1=CC(C2=CC=C(Br)C=C2)=CC=C1Br"),
         functional_groups=[stk.BromoFactory()],
     )
 
diff --git a/justfile b/justfile
index b8f251e..6579277 100644
--- a/justfile
+++ b/justfile
@@ -11,7 +11,7 @@ docs:
 # Install development environment.
 dev:
   pip install -e '.[dev]'
-  mamba install -y -c conda-forge openff-toolkit openmm openmmtools rdkit==2024.3.4 dgl espaloma_charge
+  mamba install -y -c conda-forge pytorch==2.3.1 torchdata==0.7.1 openff-toolkit openmm openmmtools rdkit!=2024.3.5 dgl espaloma_charge
   mamba install -y xtb
   mamba install -y openbabel
   mamba install -y mdanalysis
diff --git a/mypy.ini b/mypy.ini
new file mode 100644
index 0000000..5233381
--- /dev/null
+++ b/mypy.ini
@@ -0,0 +1,42 @@
+[mypy]
+show_error_codes = true
+implicit_optional = false
+warn_no_return = true
+strict_optional = true
+disallow_untyped_defs = true
+disallow_incomplete_defs = true
+check_untyped_defs = true
+disallow_untyped_decorators = true
+warn_unreachable = true
+disallow_any_generics = false
+
+[mypy-rdkit.*]
+follow_imports = skip
+follow_imports_for_stubs = True
+
+[mypy-scipy.*]
+ignore_missing_imports = True
+
+[mypy-openbabel.*]
+ignore_missing_imports = True
+
+[mypy-openmm.*]
+ignore_missing_imports = True
+
+[mypy-rmsd.*]
+ignore_missing_imports = True
+
+[mypy-espaloma_charge.*]
+ignore_missing_imports = True
+
+[mypy-openff.*]
+ignore_missing_imports = True
+
+[mypy-MDAnalysis.*]
+ignore_missing_imports = True
+
+[mypy-networkx.*]
+ignore_missing_imports = True
+
+[mypy-spindry.*]
+ignore_missing_imports = True
diff --git a/pyproject.toml b/pyproject.toml
index e1cd0a8..3f25464 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,17 +5,16 @@ build-backend = "setuptools.build_meta"
 [project]
 name = "stko"
 maintainers = [
-  { name = "Andrew Tarzia", email = "andrew.tarzia@gmail.com" },
-  { name = "Lukas Turcani", email = "lukasturcani93@gmail.com" },
+    { name = "Andrew Tarzia", email = "andrew.tarzia@gmail.com" },
+    { name = "Lukas Turcani", email = "lukasturcani93@gmail.com" },
 ]
 dependencies = [
-  # Remove pin once openmm infrastructure is safely moved to >= 2.
-  "numpy < 2",
-  # Pinned while rdkit changes stk results.
-  "rdkit == 2024.3.4",
-  "stk",
-  "networkx",
-  "rmsd",
+    "numpy",
+    # Pinned while rdkit changes stk results.
+    "rdkit == 2024.3.4",
+    "stk",
+    "networkx",
+    "rmsd",
 ]
 requires-python = ">=3.11"
 dynamic = ["version"]
@@ -23,21 +22,20 @@ readme = "README.rst"
 
 [project.optional-dependencies]
 dev = [
-  "ruff",
-  "mypy",
-  "moldoc==3.0.0",
-  # TODO: Remove pin when https://github.com/TvoroG/pytest-lazy-fixture/issues/65 is resolved.
-  # pytest-lazy-fixture 0.6.0 is incompatible with pytest 8.0.0
-  "pytest<8",
-  "pytest-benchmark",
-  "pytest-datadir",
-  "pytest-lazy-fixture",
-  "pytest-cov",
-  "sphinx",
-  "sphinx-copybutton",
-  "sphinx-rtd-theme",
-  "twine",
-  "build",
+    "ruff",
+    "mypy",
+    "moldoc==3.0.0",
+    # TODO: Remove pin when https://github.com/TvoroG/pytest-lazy-fixture/issues/65 is resolved.
+    # pytest-lazy-fixture 0.6.0 is incompatible with pytest 8.0.0
+    "pytest<8",
+    "pytest-datadir",
+    "pytest-lazy-fixture",
+    "pytest-cov",
+    "sphinx",
+    "sphinx-copybutton",
+    "sphinx-rtd-theme",
+    "twine",
+    "build",
 ]
 
 [project.urls]
@@ -49,92 +47,37 @@ write_to = "src/stko/_version.py"
 
 [tool.setuptools.packages.find]
 where = [
-  # list of folders that contain the packages (["."] by default)
-  "src",
+    # list of folders that contain the packages (["."] by default)
+    "src",
 ]
 
 [tool.ruff]
 line-length = 79
 
 [tool.ruff.lint]
+ignore = ["ANN401", "COM812", "ISC001", "FBT001", "FBT002"]
 select = ["ALL"]
-ignore = [
-  "ANN401",
-  "COM812",
-  "ISC001",
-  "FBT001",
-  "FBT002",
-]
 
 [tool.ruff.lint.pydocstyle]
 convention = "google"
 
 [tool.ruff.lint.per-file-ignores]
-"examples/*" = [
-  "D100",
-  "INP001",
-]
+"examples/*" = ["D100", "INP001"]
 "tests/*" = [
-  "D100",
-  "D101",
-  "D102",
-  "D103",
-  "D104",
-  "D105",
-  "D106",
-  "D107",
-  "S101",
-  "INP001",
+    "D100",
+    "D101",
+    "D102",
+    "D103",
+    "D104",
+    "D105",
+    "D106",
+    "D107",
+    "S101",
+    "INP001",
 ]
 "docs/source/conf.py" = ["D100", "INP001"]
 
 [tool.pytest.ini_options]
-testpaths = [
-  "tests",
-]
-python_files = [
-    "test_*.py",
-    "*_test.py",
-]
-python_functions = [
-  "test_*",
-]
-
-[tool.mypy]
-show_error_codes = true
-implicit_optional = false
-warn_no_return = true
-strict_optional = true
-disallow_untyped_defs = true
-disallow_incomplete_defs = true
-check_untyped_defs = true
-disallow_untyped_decorators = true
-warn_unreachable = true
-disallow_any_generics = false
-
-[[tool.mypy.overrides]]
-module = [
-  "rdkit.*",
-  "scipy.*",
-  "pytest_lazyfixture.*",
-  "pathos.*",
-  "matplotlib.*",
-  "pandas.*",
-  "seaborn.*",
-  "mchammer.*",
-  "spindry.*",
-  "pymongo.*",
-  "vabene.*",
-  "setuptools.*",
-  "stk.*",
-  "networkx.*",
-  "openbabel.*",
-  "MDAnalysis.*",
-  "openff.*",
-  "openmm.*",
-  "openmmforcefields.*",
-  "rmsd.*",
-  "espaloma_charge.*",
-  "traitlets.*",
-]
-ignore_missing_imports = true
+testpaths = ["tests"]
+python_files = ["test_*.py", "*_test.py"]
+python_functions = ["test_*"]
diff --git a/src/stko/__init__.py b/src/stko/__init__.py
index 71c8bf2..c851b2c 100644
--- a/src/stko/__init__.py
+++ b/src/stko/__init__.py
@@ -2,7 +2,12 @@
 
 import contextlib
 
-from stko import functional_groups, molecule_analysis, topology_functions
+from stko import (
+    functional_groups,
+    molecular_utilities,
+    molecule_analysis,
+    topology_functions,
+)
 from stko._internal.calculators.extractors.orca_extractor import OrcaExtractor
 from stko._internal.calculators.extractors.xtb_extractor import XTBExtractor
 from stko._internal.calculators.open_babel_calculators import OpenBabelEnergy
@@ -37,6 +42,7 @@
     TorsionCalculator,
 )
 from stko._internal.calculators.xtb_calculators import XTBEnergy
+from stko._internal.internal_types import ConstructedMoleculeT, MoleculeT
 from stko._internal.molecular.atoms.dummy_atom import Du
 from stko._internal.molecular.atoms.positioned_atom import PositionedAtom
 from stko._internal.molecular.conversion.md_analysis import MDAnalysis
@@ -85,7 +91,6 @@
     move_generated_macromodel_files,
 )
 from stko._internal.optimizers.xtb import XTB, XTBCREST, XTBFF, XTBFFCREST
-from stko._internal.types import ConstructedMoleculeT, MoleculeT
 from stko._internal.utilities.exceptions import (
     CalculatorError,
     ConvergenceError,
@@ -119,6 +124,7 @@
     from stko._internal.calculators.openmm_calculators import OpenMMEnergy
     from stko._internal.optimizers.openmm import OpenMMForceField, OpenMMMD
 
+
 MoleculeT = MoleculeT  # noqa: PLW0127
 """Type parameter matching any :class:`stk.Molecule` or subclasses."""
 
@@ -217,6 +223,7 @@
     "get_torsion_info_angles",
     "is_valid_xtb_solvent",
     "mol_from_mae_file",
+    "molecular_utilities",
     "molecule_analysis",
     "move_generated_macromodel_files",
     "topology_functions",
diff --git a/src/stko/_internal/calculators/openmm_calculators.py b/src/stko/_internal/calculators/openmm_calculators.py
index f2889c8..79ec828 100644
--- a/src/stko/_internal/calculators/openmm_calculators.py
+++ b/src/stko/_internal/calculators/openmm_calculators.py
@@ -95,6 +95,7 @@ def calculate(self, mol: stk.Molecule) -> abc.Generator:
                     self._partial_charges_method,
                     toolkit_registry=EspalomaChargeToolkitWrapper(),
                 )
+
             openff_molecules.append(molecule)
 
         topology = Topology.from_molecules(openff_molecules)
@@ -105,7 +106,9 @@ def calculate(self, mol: stk.Molecule) -> abc.Generator:
             force_field=self._force_field,
             topology=topology,
             positions=mol.get_position_matrix() * openmm.unit.angstrom,
-            charge_from_molecules=openff_molecules,
+            # Test this to check if molecules are _eq_, which is defined in
+            # openff.
+            charge_from_molecules=list(set(openff_molecules)),
         )
         system = interchange.to_openmm_system()
 
diff --git a/src/stko/_internal/calculators/orca_calculators.py b/src/stko/_internal/calculators/orca_calculators.py
index 88b880f..17be49b 100644
--- a/src/stko/_internal/calculators/orca_calculators.py
+++ b/src/stko/_internal/calculators/orca_calculators.py
@@ -182,9 +182,7 @@ def _write_input_file(self, path: Path, xyz_file: Path) -> None:
             f"%scf\n   MaxIter 2000\nend\n\n"
         )
         # Add geometry section.
-        string += (
-            f"* xyzfile {self._charge} {self._multiplicity} " f"{xyz_file}\n"
-        )
+        string += f"* xyzfile {self._charge} {self._multiplicity} {xyz_file}\n"
 
         path.write_text(string)
 
diff --git a/src/stko/_internal/calculators/rmsd_calculators.py b/src/stko/_internal/calculators/rmsd_calculators.py
index 344ad7f..f497a99 100644
--- a/src/stko/_internal/calculators/rmsd_calculators.py
+++ b/src/stko/_internal/calculators/rmsd_calculators.py
@@ -47,7 +47,7 @@ class RmsdCalculator:
         .. testcode:: rmsd-calc
             :hide:
 
-            assert np.isclose(rmsd, 0.21, rtol=0, atol=1E-1)
+            assert np.isclose(rmsd, 0.21, rtol=0, atol=1E-2)
 
     """
 
@@ -176,7 +176,7 @@ class RmsdMappedCalculator(RmsdCalculator):
         .. testcode:: rmsd-mapped-calc
             :hide:
 
-            assert np.isclose(rmsd, 0.24, rtol=0, atol=1E-1)
+            assert np.isclose(rmsd, 0.24, rtol=0, atol=1E-2)
 
     """
 
@@ -266,7 +266,7 @@ class KabschRmsdCalculator:
         .. testcode:: rmsd-kabsch-calc
             :hide:
 
-            assert np.isclose(rmsd, 0.20, rtol=0, atol=1E-1)
+            assert np.isclose(rmsd, 0.20, rtol=0, atol=1E-2)
 
     """
 
diff --git a/src/stko/_internal/types.py b/src/stko/_internal/internal_types.py
similarity index 100%
rename from src/stko/_internal/types.py
rename to src/stko/_internal/internal_types.py
diff --git a/src/stko/_internal/molecular/atoms/positioned_atom.py b/src/stko/_internal/molecular/atoms/positioned_atom.py
index db28728..4d1935d 100644
--- a/src/stko/_internal/molecular/atoms/positioned_atom.py
+++ b/src/stko/_internal/molecular/atoms/positioned_atom.py
@@ -76,8 +76,7 @@ def __repr__(self) -> str:
             else ""
         )
         return (
-            f"{self._atom.__class__.__name__}({self._atom.get_id()}"
-            f"{charge})"
+            f"{self._atom.__class__.__name__}({self._atom.get_id()}{charge})"
         )
 
     def __str__(self) -> str:
diff --git a/src/stko/_internal/molecular/conversion/md_analysis.py b/src/stko/_internal/molecular/conversion/md_analysis.py
index 5d5a4a4..c4bb26a 100644
--- a/src/stko/_internal/molecular/conversion/md_analysis.py
+++ b/src/stko/_internal/molecular/conversion/md_analysis.py
@@ -52,12 +52,10 @@ class MDAnalysis:
 
     def __init__(self) -> None:
         if mda is None:
-            msg = (
-                "MDAnalysis is not installed; see README for " "installation."
-            )
+            msg = "MDAnalysis is not installed; see README for installation."
             raise WrapperNotInstalledError(msg)
 
-    def get_universe(self, mol: stk.Molecule) -> mda.Universe:
+    def get_universe(self, mol: stk.Molecule):  # type: ignore[no-untyped-def] # noqa: ANN202
         """Get an MDAnalysis object.
 
         Parameters:
diff --git a/src/stko/_internal/molecular/conversion/z_matrix.py b/src/stko/_internal/molecular/conversion/z_matrix.py
index babccd4..ada6dfa 100644
--- a/src/stko/_internal/molecular/conversion/z_matrix.py
+++ b/src/stko/_internal/molecular/conversion/z_matrix.py
@@ -85,8 +85,7 @@ def get_zmatrix(self, molecule: stk.Molecule) -> str:
                 distance = round(distance, 2)
                 angle = round(angle, 2)
                 zmatrix.append(
-                    f"{atom.__class__.__name__} {i} {distance}"
-                    f"{i-1} {angle}"
+                    f"{atom.__class__.__name__} {i} {distance}{i - 1} {angle}"
                 )
                 coords += 2
 
@@ -137,13 +136,13 @@ def get_zmatrix(self, molecule: stk.Molecule) -> str:
                 torsion = round(torsion, 2)
                 zmatrix.append(
                     f"{atom.__class__.__name__} {i} {distance}"
-                    f"{i-1} {angle} {i-2} {torsion}"
+                    f"{i - 1} {angle} {i - 2} {torsion}"
                 )
                 coords += 3
 
         if 3 * len(zmatrix) - 6 != coords:
             msg = (
-                f"zmatrix: There are {coords}, not {3*len(zmatrix)-6} as "
+                f"zmatrix: There are {coords}, not {3 * len(zmatrix) - 6} as "
                 "expected. Therefore, the conversion has failed."
             )
             raise ConversionError(msg)
diff --git a/src/stko/_internal/molecular/intermediates/unreacted.py b/src/stko/_internal/molecular/intermediates/unreacted.py
index ab343bb..66e2d40 100644
--- a/src/stko/_internal/molecular/intermediates/unreacted.py
+++ b/src/stko/_internal/molecular/intermediates/unreacted.py
@@ -4,11 +4,10 @@
 from dataclasses import dataclass
 from functools import partial
 
-import numpy as np
 import stk
 from rdkit import Chem
 
-from stko._internal.molecular.networkx.network import Network
+from stko._internal.molecular.molecular_utilities import separate_molecule
 
 logger = logging.getLogger(__name__)
 
@@ -19,7 +18,7 @@ class NamedIntermediate:
 
     intermediate_name: str
     molecule: stk.BuildingBlock
-    present_bbs: dict[stk.BuildingBlock, int]
+    present_bbs: dict[stk.BuildingBlock, list[int]]
     num_atoms: int
     count_bbs: int
 
@@ -36,7 +35,7 @@ def __repr__(self) -> str:
 class IntermediatePool:
     """Container of a set of intermediates."""
 
-    intermediates: abc.Sequence[NamedIntermediate]
+    intermediates: list[NamedIntermediate]
 
     def __str__(self) -> str:
         """String representation."""
@@ -155,7 +154,7 @@ def get_reaction_results(self) -> abc.Sequence[stk.ReactionResult]:
     def yield_constructed_molecules(
         self,
         n: int | None = None,
-    ) -> abc.Sequence[stk.ConstructedMolecule]:
+    ) -> abc.Iterator[stk.ConstructedMolecule]:
         """Yield constructed molecules for possible reaction combos.
 
         If `n` is None, this produces all reactions, which could be
@@ -178,56 +177,11 @@ def yield_constructed_molecules(
                     )
                 )
 
-    def separate_molecule(
-        self, molecule: stk.Molecule
-    ) -> abc.Sequence[tuple[stk.Molecule, list[int]]]:
-        """Given a molecule, it returns distinct disconnected molecules."""
-        network = Network.init_from_molecule(molecule)
-        connected = network.get_connected_components()
-        molecules = []
-        for cg in connected:
-            # Get atoms from nodes.
-            atoms = list(cg)
-            atom_ids = tuple(i.get_id() for i in atoms)
-
-            # Sort both by atom id.
-            atom_ids, atoms = zip(  # type: ignore[assignment]
-                *sorted(zip(atom_ids, atoms, strict=True)), strict=True
-            )
-
-            atom_ids_map = {atom_ids[i]: i for i in range(len(atom_ids))}
-            new_mol = stk.BuildingBlock.init(
-                atoms=(
-                    stk.Atom(
-                        id=atom_ids_map[i.get_id()],
-                        atomic_number=i.get_atomic_number(),
-                        charge=i.get_charge(),
-                    )
-                    for i in atoms
-                ),
-                bonds=(
-                    i.with_ids(id_map=atom_ids_map)
-                    for i in molecule.get_bonds()
-                    if i.get_atom1().get_id() in atom_ids
-                    and i.get_atom2().get_id() in atom_ids
-                ),
-                position_matrix=np.array(
-                    tuple(
-                        i
-                        for i in molecule.get_atomic_positions(
-                            atom_ids=atom_ids
-                        )
-                    )
-                ),
-            )
-            molecules.append((new_mol, atom_ids))
-        return molecules
-
     def get_reacted_smiles(self, n: int | None = None) -> set[str]:
         """Yield constructed molecules with n reactions performed."""
         yielded_smiles = set()
         for const_mol in self.yield_constructed_molecules(n=n):
-            distinct_molecules = self.separate_molecule(const_mol)
+            distinct_molecules = separate_molecule(const_mol)
             for dmol, _ in distinct_molecules:
                 smiles = Chem.CanonSmiles(
                     Chem.MolToSmiles(dmol.to_rdkit_mol())
@@ -247,16 +201,16 @@ def get_present_building_blocks(
         self,
         const_mol: stk.ConstructedMolecule,
         subset_ids: list[int],
-    ) -> dict[stk.BuildingBlock, int]:
+    ) -> dict[stk.BuildingBlock, list[int]]:
         """Get the building blocks present in a constructed molecule."""
-        bbs = defaultdict(list)
+        bbs: dict[stk.BuildingBlock, list[int]] = defaultdict(list)
         for atom_info in const_mol.get_atom_infos():
             if atom_info.get_atom().get_id() in subset_ids and (
                 atom_info.get_building_block_id()
-                not in bbs[atom_info.get_building_block()]
+                not in bbs[atom_info.get_building_block()]  # type: ignore[index]
             ):
-                bbs[atom_info.get_building_block()].append(
-                    atom_info.get_building_block_id()
+                bbs[atom_info.get_building_block()].append(  # type: ignore[index]
+                    atom_info.get_building_block_id()  # type: ignore[arg-type]
                 )
 
         return bbs
@@ -268,7 +222,7 @@ def get_named_intermediates(
         yielded_smiles = set()
         pool = IntermediatePool(intermediates=[])
         for const_mol in self.yield_constructed_molecules(n=n):
-            distinct_molecules = self.separate_molecule(const_mol)
+            distinct_molecules = separate_molecule(const_mol)
             for dmol, datom_ids in distinct_molecules:
                 smiles = Chem.CanonSmiles(
                     Chem.MolToSmiles(dmol.to_rdkit_mol())
diff --git a/src/stko/_internal/molecular/molecular_utilities.py b/src/stko/_internal/molecular/molecular_utilities.py
new file mode 100644
index 0000000..7c192dd
--- /dev/null
+++ b/src/stko/_internal/molecular/molecular_utilities.py
@@ -0,0 +1,111 @@
+"""Module of molecular utilities."""
+
+from collections import abc
+
+import numpy as np
+import stk
+from rdkit.Chem import AllChem as rdkit  # noqa: N813
+
+from stko._internal.molecular.networkx.network import Network
+
+
+def merge_stk_molecules(
+    molecules: abc.Sequence[stk.Molecule],
+) -> stk.BuildingBlock:
+    """Merge any list of separate molecules to be in one class."""
+    atoms: list[stk.Atom] = []
+    bonds: list[stk.Bond] = []
+    pos_mat = []
+    for molecule in molecules:
+        atom_ids_map = {}
+        for atom in molecule.get_atoms():
+            new_id = len(atoms)
+            atom_ids_map[atom.get_id()] = new_id
+            atoms.append(
+                stk.Atom(
+                    id=atom_ids_map[atom.get_id()],
+                    atomic_number=atom.get_atomic_number(),
+                    charge=atom.get_charge(),
+                )
+            )
+
+        bonds.extend(
+            i.with_ids(id_map=atom_ids_map) for i in molecule.get_bonds()
+        )
+        pos_mat.extend(list(molecule.get_position_matrix()))
+
+    return stk.BuildingBlock.init(
+        atoms=atoms,
+        bonds=bonds,
+        position_matrix=np.array(pos_mat),
+    )
+
+
+def separate_molecule(
+    molecule: stk.Molecule,
+) -> abc.Sequence[tuple[stk.BuildingBlock, list[int]]]:
+    """Given a molecule, it returns distinct disconnected molecules."""
+    network = Network.init_from_molecule(molecule)
+    connected = network.get_connected_components()
+    molecules = []
+    for cg in connected:
+        # Get atoms from nodes.
+        atoms = list(cg)
+        atom_ids = tuple(i.get_id() for i in atoms)
+
+        # Sort both by atom id.
+        atom_ids, atoms = zip(  # type: ignore[assignment]
+            *sorted(zip(atom_ids, atoms, strict=True)), strict=True
+        )
+
+        atom_ids_map = {atom_ids[i]: i for i in range(len(atom_ids))}
+        new_mol = stk.BuildingBlock.init(
+            atoms=(
+                stk.Atom(
+                    id=atom_ids_map[i.get_id()],
+                    atomic_number=i.get_atomic_number(),
+                    charge=i.get_charge(),
+                )
+                for i in atoms
+            ),
+            bonds=(
+                i.with_ids(id_map=atom_ids_map)
+                for i in molecule.get_bonds()
+                if i.get_atom1().get_id() in atom_ids
+                and i.get_atom2().get_id() in atom_ids
+            ),
+            position_matrix=np.array(
+                tuple(
+                    i for i in molecule.get_atomic_positions(atom_ids=atom_ids)
+                )
+            ),
+        )
+        molecules.append((new_mol, list(atom_ids)))
+    return molecules
+
+
+def update_stk_from_rdkit_conformer(
+    stk_mol: stk.Molecule,
+    rdk_mol: rdkit.Mol,
+    conf_id: int,
+) -> stk.Molecule:
+    """Update the structure to match `conf_id` of `mol`.
+
+    Parameters:
+        struct:
+            The molecule whoce coordinates are to be updated.
+
+        mol:
+            The :mod:`rdkit` molecule to use for the structure update.
+
+        conf_id:
+            The conformer ID of the `mol` to update from.
+
+    Returns:
+        The molecule.
+
+    TODO: Add to stko.
+
+    """
+    pos_mat = rdk_mol.GetConformer(id=conf_id).GetPositions()
+    return stk_mol.with_position_matrix(pos_mat)
diff --git a/src/stko/_internal/molecular/topology_extractor/topology_info.py b/src/stko/_internal/molecular/topology_extractor/topology_info.py
index d235ec2..ae4e2ae 100644
--- a/src/stko/_internal/molecular/topology_extractor/topology_info.py
+++ b/src/stko/_internal/molecular/topology_extractor/topology_info.py
@@ -98,7 +98,7 @@ def write(self, path: Path) -> None:
             a2 = edge[1]
             if a1 in atoms and a2 in atoms:
                 content.append(
-                    f"{conect:<6}{a1+1:>5}{a2+1:>5}               \n"
+                    f"{conect:<6}{a1 + 1:>5}{a2 + 1:>5}               \n"
                 )
 
         content.append("END\n")
diff --git a/src/stko/_internal/optimizers/aligner.py b/src/stko/_internal/optimizers/aligner.py
index 51f9e05..9cc92b4 100644
--- a/src/stko/_internal/optimizers/aligner.py
+++ b/src/stko/_internal/optimizers/aligner.py
@@ -9,8 +9,8 @@
 from scipy.spatial.distance import cdist
 
 from stko._internal.calculators.rmsd_calculators import RmsdMappedCalculator
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import Optimizer
-from stko._internal.types import MoleculeT
 
 logger = logging.getLogger(__name__)
 
diff --git a/src/stko/_internal/optimizers/collapser.py b/src/stko/_internal/optimizers/collapser.py
index f89a62e..10c94b5 100644
--- a/src/stko/_internal/optimizers/collapser.py
+++ b/src/stko/_internal/optimizers/collapser.py
@@ -11,8 +11,8 @@
 from scipy.spatial.distance import pdist
 from stk import PdbWriter
 
+from stko._internal.internal_types import ConstructedMoleculeT
 from stko._internal.molecular.periodic.unitcell import UnitCell
-from stko._internal.types import ConstructedMoleculeT
 from stko._internal.utilities.exceptions import InputError
 from stko._internal.utilities.utilities import get_atom_distance
 
@@ -209,7 +209,7 @@ def _get_bb_vectors(
             norms = {
                 i: np.linalg.norm(bb_cent_vectors[i]) for i in bb_cent_vectors
             }
-            max_distance = max(list(norms.values()))
+            max_distance = max(list(norms.values()))  # type: ignore[type-var]
             bb_cent_scales = {i: norms[i] / max_distance for i in norms}
         else:
             bb_cent_scales = {i: np.floating(1) for i in bb_cent_vectors}
@@ -305,10 +305,7 @@ def p_optimize(
 
         """
         if not isinstance(mol, stk.ConstructedMolecule):
-            msg = (
-                f"{mol} needs to be a ConstructedMolecule for "
-                f"this optimizer"
-            )
+            msg = f"{mol} needs to be a ConstructedMolecule for this optimizer"
             raise InputError(msg)
 
         output_dir = self._output_dir.resolve()
@@ -742,9 +739,7 @@ def optimize(self, mol: ConstructedMoleculeT) -> ConstructedMoleculeT:
                 "Step system_potential nonbond_potential max_dist "
                 "opt_bbs updated?\n"
             )
-            f.write(
-                f"{steps[-1]} {spots[-1]} {npots[-1]} {maxds[-1]} " "-- --\n"
-            )
+            f.write(f"{steps[-1]} {spots[-1]} {npots[-1]} {maxds[-1]} -- --\n")
             for step in range(1, self._num_steps):
                 # Randomly select a long bond.
                 lb_ids = random.choice(list(long_bond_infos.keys()))  # noqa: S311
@@ -843,7 +838,7 @@ def optimize(self, mol: ConstructedMoleculeT) -> ConstructedMoleculeT:
             f.write(
                 "Optimisation done:\n"
                 f"{len(passed)} steps passed: "
-                f"{len(passed)/self._num_steps}"
+                f"{len(passed) / self._num_steps}"
             )
 
         self._plot_progess(steps, maxds, spots, npots, output_dir)
diff --git a/src/stko/_internal/optimizers/gulp.py b/src/stko/_internal/optimizers/gulp.py
index 0419045..ef4d107 100644
--- a/src/stko/_internal/optimizers/gulp.py
+++ b/src/stko/_internal/optimizers/gulp.py
@@ -11,6 +11,7 @@
 import stk
 from rdkit.Chem import AllChem as rdkit  # noqa: N813
 
+from stko._internal.internal_types import MoleculeT
 from stko._internal.molecular.periodic.unitcell import UnitCell
 from stko._internal.optimizers.optimizers import Optimizer
 from stko._internal.optimizers.utilities import (
@@ -20,7 +21,6 @@
     has_metal_atom,
     to_rdkit_mol_without_metals,
 )
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.exceptions import (
     ExpectedMetalError,
     ForceFieldSetupError,
@@ -468,8 +468,8 @@ def _bond_section(
                 bond_type = "triple"
 
             string = (
-                f"connect {bond.get_atom1().get_id()+1} "
-                f"{bond.get_atom2().get_id()+1} {bond_type}"
+                f"connect {bond.get_atom1().get_id() + 1} "
+                f"{bond.get_atom2().get_id() + 1} {bond_type}"
             )
             bond_section += string + "\n"
 
diff --git a/src/stko/_internal/optimizers/macromodel.py b/src/stko/_internal/optimizers/macromodel.py
index 6785226..1c283e0 100644
--- a/src/stko/_internal/optimizers/macromodel.py
+++ b/src/stko/_internal/optimizers/macromodel.py
@@ -9,13 +9,13 @@
 import rdkit.Chem.AllChem as rdkit  # noqa: N813
 import stk
 
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import Optimizer
 from stko._internal.optimizers.utilities import (
     MAEExtractor,
     mol_from_mae_file,
     move_generated_macromodel_files,
 )
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.exceptions import (
     ConversionError,
     ForceFieldError,
diff --git a/src/stko/_internal/optimizers/open_babel.py b/src/stko/_internal/optimizers/open_babel.py
index 0d06827..40a93a2 100644
--- a/src/stko/_internal/optimizers/open_babel.py
+++ b/src/stko/_internal/optimizers/open_babel.py
@@ -8,8 +8,8 @@
 except ImportError:
     openbabel = None
 
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import Optimizer
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.exceptions import (
     ForceFieldSetupError,
     WrapperNotInstalledError,
@@ -75,7 +75,7 @@ def __init__(
         cg_steps: int = 50,
     ) -> None:
         if openbabel is None:
-            msg = "openbabel is not installed; see README for " "installation."
+            msg = "openbabel is not installed; see README for installation."
             raise WrapperNotInstalledError(msg)
 
         self._forcefield = forcefield
diff --git a/src/stko/_internal/optimizers/openmm.py b/src/stko/_internal/optimizers/openmm.py
index ac97227..1d3359c 100644
--- a/src/stko/_internal/optimizers/openmm.py
+++ b/src/stko/_internal/optimizers/openmm.py
@@ -11,8 +11,8 @@
 from openmm import app, openmm
 
 from stko._internal.calculators.openmm_calculators import OpenMMEnergy
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import NullOptimizer, Optimizer
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.exceptions import InputError
 from stko._internal.utilities.utilities import get_atom_distance
 
@@ -177,7 +177,9 @@ def optimize(self, mol: MoleculeT) -> MoleculeT:
             force_field=self._force_field,
             topology=topology,
             positions=mol.get_position_matrix() * openmm.unit.angstrom,
-            charge_from_molecules=openff_molecules,
+            # Test this to check if molecules are _eq_, which is defined in
+            # openff.
+            charge_from_molecules=list(set(openff_molecules)),
         )
         system = interchange.to_openmm_system()
         # Add constraints.
@@ -418,7 +420,9 @@ def optimize(self, mol: MoleculeT) -> MoleculeT:
             force_field=self._force_field,
             topology=topology,
             positions=mol.get_position_matrix() * openmm.unit.angstrom,
-            charge_from_molecules=openff_molecules,
+            # Test this to check if molecules are _eq_, which is defined in
+            # openff.
+            charge_from_molecules=list(set(openff_molecules)),
         )
         simulation = interchange.to_openmm_simulation(
             integrator=self._integrator,
diff --git a/src/stko/_internal/optimizers/optimizers.py b/src/stko/_internal/optimizers/optimizers.py
index 9341137..1cd1148 100644
--- a/src/stko/_internal/optimizers/optimizers.py
+++ b/src/stko/_internal/optimizers/optimizers.py
@@ -4,7 +4,7 @@
 
 import stk
 
-from stko._internal.types import MoleculeT
+from stko._internal.internal_types import MoleculeT
 
 logger = logging.getLogger(__name__)
 
@@ -182,7 +182,7 @@ def _suffix_from_writer(self) -> str:
             return "pdb"
         if isinstance(self._writer, stk.MolWriter):
             return "mol"
-        return None
+        raise NotImplementedError
 
     def optimize(self, mol: MoleculeT) -> MoleculeT:
         filesuffix = self._suffix_from_writer()
diff --git a/src/stko/_internal/optimizers/rdkit.py b/src/stko/_internal/optimizers/rdkit.py
index 20952f9..6ee1890 100644
--- a/src/stko/_internal/optimizers/rdkit.py
+++ b/src/stko/_internal/optimizers/rdkit.py
@@ -5,13 +5,13 @@
 import rdkit.Chem.AllChem as rdkit  # noqa: N813
 import stk
 
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import Optimizer
 from stko._internal.optimizers.utilities import (
     get_metal_atoms,
     get_metal_bonds,
     to_rdkit_mol_without_metals,
 )
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.utilities import vector_angle
 
 logger = logging.getLogger(__name__)
diff --git a/src/stko/_internal/optimizers/xtb.py b/src/stko/_internal/optimizers/xtb.py
index 3e4c79d..b881aa0 100644
--- a/src/stko/_internal/optimizers/xtb.py
+++ b/src/stko/_internal/optimizers/xtb.py
@@ -9,8 +9,8 @@
 import stk
 
 from stko._internal.calculators.extractors.xtb_extractor import XTBExtractor
+from stko._internal.internal_types import MoleculeT
 from stko._internal.optimizers.optimizers import Optimizer
-from stko._internal.types import MoleculeT
 from stko._internal.utilities.exceptions import (
     ConvergenceError,
     InvalidSolventError,
@@ -415,8 +415,8 @@ def _run_optimizations(
 
         """
         for run in range(self._max_runs):
-            xyz = f"input_structure_{run+1}.xyz"
-            out_file = f"optimization_{run+1}.output"
+            xyz = f"input_structure_{run + 1}.xyz"
+            out_file = f"optimization_{run + 1}.output"
             mol.write(xyz)
             self._write_detailed_control()
             self._run_xtb(xyz=xyz, out_file=out_file)
diff --git a/src/stko/molecular_utilities.py b/src/stko/molecular_utilities.py
new file mode 100644
index 0000000..1cb58ba
--- /dev/null
+++ b/src/stko/molecular_utilities.py
@@ -0,0 +1,13 @@
+"""Tools for molecules."""
+
+from stko._internal.molecular.molecular_utilities import (
+    merge_stk_molecules,
+    separate_molecule,
+    update_stk_from_rdkit_conformer,
+)
+
+__all__ = [
+    "merge_stk_molecules",
+    "separate_molecule",
+    "update_stk_from_rdkit_conformer",
+]
diff --git a/tests/calculators/geometry/conftest.py b/tests/calculators/geometry/conftest.py
index 2d86164..ed099ed 100644
--- a/tests/calculators/geometry/conftest.py
+++ b/tests/calculators/geometry/conftest.py
@@ -2327,9 +2327,7 @@
                             position_matrix=np.array([[0, 0, 0]]),
                         ),
                         stk.BuildingBlock(
-                            smiles=(
-                                "C1=NC=CC(C2=CC=CC(C3=C" "C=NC=C3)=C2)=C1"
-                            ),
+                            smiles=("C1=NC=CC(C2=CC=CC(C3=CC=NC=C3)=C2)=C1"),
                             functional_groups=[
                                 stk.SmartsFunctionalGroupFactory(
                                     smarts="[#6]~[#7X2]~[#6]",
diff --git a/tests/calculators/openmm/__init__.py b/tests/calculators/openmm/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/calculators/openmm/conftest.py b/tests/calculators/openmm/conftest.py
new file mode 100644
index 0000000..49e8801
--- /dev/null
+++ b/tests/calculators/openmm/conftest.py
@@ -0,0 +1,49 @@
+from dataclasses import dataclass
+
+import pytest
+import stk
+
+
+@dataclass(frozen=True, slots=True)
+class CaseData:
+    molecule: stk.Molecule
+    energy: float
+
+
+@pytest.fixture(
+    scope="session",
+    params=[
+        CaseData(
+            molecule=stk.BuildingBlock("NCCN"),
+            energy=142.99500933275297,
+        ),
+        CaseData(
+            molecule=stk.ConstructedMolecule(
+                topology_graph=stk.host_guest.Complex(
+                    host=stk.BuildingBlock("NCCN"),
+                    guests=stk.host_guest.Guest(
+                        building_block=stk.BuildingBlock("CC"),
+                        displacement=(5, 5, 5),
+                    ),
+                    optimizer=stk.Spinner(),
+                ),
+            ),
+            energy=165.25686166058674,
+        ),
+        CaseData(
+            molecule=stk.ConstructedMolecule(
+                topology_graph=stk.host_guest.Complex(
+                    host=stk.BuildingBlock("NCCN"),
+                    guests=stk.host_guest.Guest(
+                        building_block=stk.BuildingBlock("NCCN"),
+                        displacement=(4, 4, 4),
+                    ),
+                    optimizer=stk.Spinner(),
+                ),
+            ),
+            energy=285.665342151835,
+        ),
+    ],
+)
+def case_molecule(request: pytest.FixtureRequest) -> CaseData:
+    return request.param
diff --git a/tests/calculators/openmm/test_openmm_energy.py b/tests/calculators/openmm/test_openmm_energy.py
new file mode 100644
index 0000000..a5dbe36
--- /dev/null
+++ b/tests/calculators/openmm/test_openmm_energy.py
@@ -0,0 +1,21 @@
+try:
+    from openff.toolkit import ForceField
+except ImportError:
+    ForceField = None
+import numpy as np
+
+import stko
+
+from .conftest import CaseData
+
+
+def test_openmm_energy(case_molecule: CaseData) -> None:
+    if ForceField is not None:
+        energy = (
+            stko.OpenMMEnergy(
+                force_field=ForceField("openff_unconstrained-2.1.0.offxml"),
+                partial_charges_method="espaloma-am1bcc",
+            ).get_energy(case_molecule.molecule),
+        )
+
+        assert np.isclose(energy, case_molecule.energy)
diff --git a/tests/calculators/planarity/conftest.py b/tests/calculators/planarity/conftest.py
index 67fcf39..ed0f3ed 100644
--- a/tests/calculators/planarity/conftest.py
+++ b/tests/calculators/planarity/conftest.py
@@ -82,8 +82,7 @@ class CaseData:
         ),
         CaseData(
             molecule=stk.BuildingBlock(
-                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)"
-                "cc3c2c1"
+                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)cc3c2c1"
             ),
             plane_ids=None,
             deviation_ids=None,
diff --git a/tests/calculators/pore/conftest.py b/tests/calculators/pore/conftest.py
index eb2b9ab..2b35abe 100644
--- a/tests/calculators/pore/conftest.py
+++ b/tests/calculators/pore/conftest.py
@@ -85,9 +85,7 @@
                             position_matrix=np.array([[0, 0, 0]]),
                         ),
                         stk.BuildingBlock(
-                            smiles=(
-                                "C1=NC=CC(C2=CC=CC(C3=C" "C=NC=C3)=C2)=C1"
-                            ),
+                            smiles=("C1=NC=CC(C2=CC=CC(C3=CC=NC=C3)=C2)=C1"),
                             functional_groups=[
                                 stk.SmartsFunctionalGroupFactory(
                                     smarts="[#6]~[#7X2]~[#6]",
diff --git a/tests/calculators/rmsd/conftest.py b/tests/calculators/rmsd/conftest.py
index e616d8b..fc85717 100644
--- a/tests/calculators/rmsd/conftest.py
+++ b/tests/calculators/rmsd/conftest.py
@@ -248,7 +248,7 @@ def ordering_case_data(request: pytest.FixtureRequest) -> CaseData:
                 np.array((0, 0, 0)),
             )
             .with_displacement(np.array((0, 0, 1))),
-            rmsd=0.5943193981905652,
+            rmsd=0.5922202905739481,
             kabsch_rmsd=0.5943193981905652,
         ),
         CaseData(
diff --git a/tests/molecular/utilities/__init__.py b/tests/molecular/utilities/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/molecular/utilities/conftest.py b/tests/molecular/utilities/conftest.py
new file mode 100644
index 0000000..013d7c7
--- /dev/null
+++ b/tests/molecular/utilities/conftest.py
@@ -0,0 +1,139 @@
+from dataclasses import dataclass
+
+import pytest
+import stk
+
+import stko
+
+
+@dataclass(frozen=True, slots=True)
+class CaseData:
+    molecules: list[stk.Molecule]
+    name: str
+
+
+@pytest.fixture(
+    scope="session",
+    params=(
+        lambda name: CaseData(
+            molecules=[stk.BuildingBlock("C1=CC=CC=C1")],
+            name=name,
+        ),
+        lambda name: CaseData(
+            molecules=[
+                stk.BuildingBlock("C1=CC=CC=C1"),
+                stk.BuildingBlock("C1N=CC(CCC2CCOC2)N=1"),
+            ],
+            name=name,
+        ),
+        lambda name: CaseData(
+            molecules=[
+                stk.BuildingBlock("C1=CC=CC=C1"),
+                stk.BuildingBlock("C1=CC=C(C=C1)C#CC2=CN=CC=C2"),
+            ],
+            name=name,
+        ),
+        lambda name: CaseData(
+            molecules=[
+                stk.BuildingBlock("C1=CC=CC=C1"),
+                stk.BuildingBlock("C1N=CC(CCC2CCOC2)N=1"),
+                stk.BuildingBlock("C1=CC=C(C=C1)C#CC2=CN=CC=C2"),
+            ],
+            name=name,
+        ),
+        lambda name: CaseData(
+            molecules=[
+                i[0]
+                for i in stko.molecular_utilities.separate_molecule(
+                    stk.ConstructedMolecule(
+                        topology_graph=stk.host_guest.Complex(
+                            host=stk.BuildingBlock("NCCN"),
+                            guests=stk.host_guest.Guest(
+                                building_block=stk.BuildingBlock("NCCN"),
+                                displacement=(5, 5, 5),
+                            ),
+                        ),
+                    )
+                )
+            ],
+            name=name,
+        ),
+    ),
+)
+def case_data(request: pytest.FixtureRequest) -> CaseData:
+    return request.param(
+        f"{request.fixturename}{request.param_index}",
+    )
+
+
+@dataclass(frozen=True, slots=True)
+class TogetherCaseData:
+    molecule: stk.Molecule
+    name: str
+
+
+@pytest.fixture(
+    scope="session",
+    params=(
+        lambda name: TogetherCaseData(
+            molecule=stk.BuildingBlock("C1=CC=CC=C1"),
+            name=name,
+        ),
+        lambda name: TogetherCaseData(
+            molecule=stko.molecular_utilities.merge_stk_molecules(
+                [
+                    stk.BuildingBlock("C1=CC=CC=C1"),
+                    stk.BuildingBlock("C1N=CC(CCC2CCOC2)N=1"),
+                    stk.BuildingBlock("C1=CC=C(C=C1)C#CC2=CN=CC=C2"),
+                ]
+            ),
+            name=name,
+        ),
+        lambda name: TogetherCaseData(
+            molecule=stk.ConstructedMolecule(
+                topology_graph=stk.host_guest.Complex(
+                    host=stk.BuildingBlock("NCCN"),
+                    guests=stk.host_guest.Guest(
+                        building_block=stk.BuildingBlock("NCCN"),
+                        displacement=(5, 5, 5),
+                    ),
+                ),
+            ),
+            name=name,
+        ),
+        lambda name: TogetherCaseData(
+            molecule=stk.ConstructedMolecule(
+                topology_graph=stk.host_guest.Complex(
+                    host=stk.BuildingBlock("C1N=CC(CCC2CCOC2)N=1"),
+                    guests=stk.host_guest.Guest(
+                        building_block=stk.BuildingBlock("NCCN"),
+                        displacement=(5, 5, 5),
+                    ),
+                ),
+            ),
+            name=name,
+        ),
+        lambda name: TogetherCaseData(
+            molecule=stk.ConstructedMolecule(
+                topology_graph=stk.host_guest.Complex(
+                    host=stk.BuildingBlock("C1N=CC(CCC2CCOC2)N=1"),
+                    guests=(
+                        stk.host_guest.Guest(
+                            building_block=stk.BuildingBlock("NCCN"),
+                            displacement=(5, 5, 5),
+                        ),
+                        stk.host_guest.Guest(
+                            building_block=stk.BuildingBlock("C1=CC=CC=C1"),
+                            displacement=(5, 5, 5),
+                        ),
+                    ),
+                ),
+            ),
+            name=name,
+        ),
+    ),
+)
+def together_case_data(request: pytest.FixtureRequest) -> CaseData:
+    return request.param(
+        f"{request.fixturename}{request.param_index}",
+    )
diff --git a/tests/molecular/utilities/test_merge_molecules.py b/tests/molecular/utilities/test_merge_molecules.py
new file mode 100644
index 0000000..c174c05
--- /dev/null
+++ b/tests/molecular/utilities/test_merge_molecules.py
@@ -0,0 +1,13 @@
+import stko
+
+from .conftest import CaseData
+
+
+def test_merge_molecules(case_data: CaseData) -> None:
+    merged = stko.molecular_utilities.merge_stk_molecules(case_data.molecules)
+    assert merged.get_num_atoms() == sum(
+        [i.get_num_atoms() for i in case_data.molecules]
+    )
+    assert merged.get_num_bonds() == sum(
+        [i.get_num_bonds() for i in case_data.molecules]
+    )
diff --git a/tests/molecular/utilities/test_seperate_molecules.py b/tests/molecular/utilities/test_seperate_molecules.py
new file mode 100644
index 0000000..ab6b630
--- /dev/null
+++ b/tests/molecular/utilities/test_seperate_molecules.py
@@ -0,0 +1,15 @@
+import stko
+
+from .conftest import CaseData
+
+
+def test_seperate_molecules(together_case_data: CaseData) -> None:
+    distinct = stko.molecular_utilities.separate_molecule(
+        together_case_data.molecule
+    )
+    assert together_case_data.molecule.get_num_atoms() == sum(
+        [i[0].get_num_atoms() for i in distinct]
+    )
+    assert together_case_data.molecule.get_num_bonds() == sum(
+        [i[0].get_num_bonds() for i in distinct]
+    )
diff --git a/tests/molecular/z_matrix/conftest.py b/tests/molecular/z_matrix/conftest.py
index 21af5a4..5a15039 100644
--- a/tests/molecular/z_matrix/conftest.py
+++ b/tests/molecular/z_matrix/conftest.py
@@ -26,10 +26,7 @@
         CaseData(
             molecule=stk.BuildingBlock("BrC#CBr"),
             zmatrix=(
-                "Br\n"
-                "C 1 1.9\n"
-                "C 2 1.211 0.17\n"
-                "Br 3 1.892 179.14 1 -171.59"
+                "Br\nC 1 1.9\nC 2 1.211 0.17\nBr 3 1.892 179.14 1 -171.59"
             ),
         ),
     ),
diff --git a/tests/optimizers/aligner/conftest.py b/tests/optimizers/aligner/conftest.py
index 923ff92..c4f26a4 100644
--- a/tests/optimizers/aligner/conftest.py
+++ b/tests/optimizers/aligner/conftest.py
@@ -24,7 +24,7 @@ class CaseData:
                     np.array((0, 0, 0)),
                 )
             ),
-            rmsd=0.22366328852274148,
+            rmsd=0.28595531556943543,
         ),
         CaseData(
             molecule=stk.BuildingBlock("CCCCCC"),
diff --git a/tests/optimizers/aligner/test_aligner.py b/tests/optimizers/aligner/test_aligner.py
index 9207ff6..6525a8d 100644
--- a/tests/optimizers/aligner/test_aligner.py
+++ b/tests/optimizers/aligner/test_aligner.py
@@ -21,7 +21,7 @@ def test_aligner(case_molecule: CaseData) -> None:
 
     test_rmsd = calculator.get_results(opt_res).get_rmsd()
 
-    assert np.isclose(test_rmsd, case_molecule.rmsd, atol=1e-6)
+    assert np.isclose(test_rmsd, case_molecule.rmsd, atol=1e-2)
     assert test_rmsd < test_rmsd_unopt
 
 
@@ -40,5 +40,5 @@ def test_alignment_potential(case_potential: CasePotential) -> None:
         potential.compute_potential(
             supramolecule,
         ),
-        atol=1e-6,
+        atol=1e-2,
     )
diff --git a/tests/optimizers/openmm/__init__.py b/tests/optimizers/openmm/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/optimizers/openmm/conftest.py b/tests/optimizers/openmm/conftest.py
new file mode 100644
index 0000000..ae608b9
--- /dev/null
+++ b/tests/optimizers/openmm/conftest.py
@@ -0,0 +1,23 @@
+from dataclasses import dataclass
+
+import pytest
+import stk
+
+
+@dataclass(frozen=True, slots=True)
+class CaseData:
+    molecule: stk.Molecule
+    opt_energy: float
+
+
+@pytest.fixture(
+    scope="session",
+    params=[
+        CaseData(
+            molecule=stk.BuildingBlock("C"),
+            opt_energy=0.22629718073104876,
+        ),
+    ],
+)
+def case_molecule(request: pytest.FixtureRequest) -> CaseData:
+    return request.param
diff --git a/tests/optimizers/openmm/test_openmm.py b/tests/optimizers/openmm/test_openmm.py
new file mode 100644
index 0000000..8211429
--- /dev/null
+++ b/tests/optimizers/openmm/test_openmm.py
@@ -0,0 +1,30 @@
+try:
+    from openff.toolkit import ForceField
+except ImportError:
+    ForceField = None
+import numpy as np
+
+import stko
+
+from .conftest import CaseData
+
+
+def test_openmm(case_molecule: CaseData) -> None:
+    if ForceField is not None:
+        optimiser = stko.OpenMMForceField(
+            # Load the openff-2.1.0 force field appropriate for
+            # vacuum calculations (without constraints)
+            force_field=ForceField("openff_unconstrained-2.1.0.offxml"),
+            restricted=False,
+            partial_charges_method="espaloma-am1bcc",
+        )
+        opt_molecule = optimiser.optimize(case_molecule.molecule)
+
+        energy = (
+            stko.OpenMMEnergy(
+                force_field=ForceField("openff_unconstrained-2.1.0.offxml"),
+                partial_charges_method="espaloma-am1bcc",
+            ).get_energy(opt_molecule),
+        )
+
+        assert np.isclose(energy, case_molecule.opt_energy)
diff --git a/tests/optimizers/rdkit/conftest.py b/tests/optimizers/rdkit/conftest.py
index 2f69a59..c245799 100644
--- a/tests/optimizers/rdkit/conftest.py
+++ b/tests/optimizers/rdkit/conftest.py
@@ -19,8 +19,7 @@ class CaseData:
         ),
         CaseData(
             molecule=stk.BuildingBlock(
-                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)"
-                "cc3c2c1"
+                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)cc3c2c1"
             ),
             unoptimised_energy=276.0206611549808,
         ),
@@ -67,8 +66,7 @@ def case_uff_molecule(request: pytest.FixtureRequest) -> CaseData:
         ),
         CaseData(
             molecule=stk.BuildingBlock(
-                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)"
-                "cc3c2c1"
+                "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)cc3c2c1"
             ),
             unoptimised_energy=226.18914087716263,
         ),
@@ -111,7 +109,7 @@ def case_mmff_molecule(request: pytest.FixtureRequest) -> CaseData:
     params=[
         stk.BuildingBlock("NCCN"),
         stk.BuildingBlock(
-            "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)" "cc3c2c1"
+            "C(#Cc1cccc2ccncc21)c1ccc2[nH]c3ccc(C#Cc4cccc5cnccc54)cc3c2c1"
         ),
         stk.BuildingBlock("CCCCCC"),
         stk.BuildingBlock("c1ccccc1"),