diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index 6dc9e06..4fbbf94 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -26,15 +26,16 @@ jobs:
         os: [macOS-latest, ubuntu-latest, windows-latest]
         python-version: ["3.9", "3.10", "3.11", "3.12"]
         chemlib: [obabel, rdkit]
-        graphlib: [nx, gt, rx, all]
+        graphlib: [nogt, all]
         exclude:
-        # graph-tools does not work on Windows
-        - {os: "windows-latest", graphlib: "gt"}
-        - {os: "windows-latest", graphlib: "all"}
-        - {graphlib: "all", chemlib: "obabel"}
+          # graph-tool does not work on Windows
+          - {os: "windows-latest", graphlib: "all"}
+          # Don't run without graph-tool on Ubuntu and macOS
+          - {os: "ubuntu-latest", graphlib: "nogt"}
+          - {os: "macOS-latest", graphlib: "nogt"}
         include:
-        - {os: "macOS-14", graphlib: "gt", chemlib: "obabel", python-version: "3.12"}
-        - {os: "macOS-14", graphlib: "nx", chemlib: "rdkit", python-version: "3.12"}
+          - {os: "macOS-14", graphlib: "all", chemlib: "obabel", python-version: "3.12"}
+          - {os: "macOS-14", graphlib: "all", chemlib: "rdkit", python-version: "3.12"}
 
     steps:
     - uses: actions/checkout@v3
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 829dc7f..a6111ad 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,10 +8,24 @@
 Date:            XX/YY/ZZZZ
 Contributors:    @RMeli
 
+### Added
+
+* Warnings filter for tests of multiple backends [PR #118 | @RMeli]
+* Parametrized fixture to run tests with all available backends [PR #118 | @RMeli]
+
 ### Improved
 
 * Test IDs [PR #117 | @RMeli]
 
+### Changed
+
+* Location of backend tests to standalone file [PR #118 | @RMeli]
+
+### Removed
+
+* Many CI configurations in favour of running tests for all available backends [PR #118 | @RMeli]
+* `tests/molecule.py` in favour of fixtures [PR #118 | @RMeli]
+
 ## Version 0.7.0
 
 Date:            05/04/2024
diff --git a/devtools/conda-envs/spyrmsd-test-obabel-gt.yaml b/devtools/conda-envs/spyrmsd-test-obabel-all.yaml
similarity index 87%
rename from devtools/conda-envs/spyrmsd-test-obabel-gt.yaml
rename to devtools/conda-envs/spyrmsd-test-obabel-all.yaml
index 47476cb..d067db0 100644
--- a/devtools/conda-envs/spyrmsd-test-obabel-gt.yaml
+++ b/devtools/conda-envs/spyrmsd-test-obabel-all.yaml
@@ -4,12 +4,15 @@ channels:
 dependencies:
   # Base
   - python
+  - pip
   - setuptools
 
   # Maths
   - numpy
   - scipy
   - graph-tool
+  - networkx>=2
+  - rustworkx
 
   # Chemistry
   - openbabel
diff --git a/devtools/conda-envs/spyrmsd-test-obabel-nx.yaml b/devtools/conda-envs/spyrmsd-test-obabel-nogt.yaml
similarity index 92%
rename from devtools/conda-envs/spyrmsd-test-obabel-nx.yaml
rename to devtools/conda-envs/spyrmsd-test-obabel-nogt.yaml
index f16ceaa..bca1f0c 100644
--- a/devtools/conda-envs/spyrmsd-test-obabel-nx.yaml
+++ b/devtools/conda-envs/spyrmsd-test-obabel-nogt.yaml
@@ -4,12 +4,14 @@ channels:
 dependencies:
   # Base
   - python
+  - pip
   - setuptools
 
   # Maths
   - numpy
   - scipy
   - networkx>=2
+  - rustworkx
 
   # Chemistry
   - openbabel
diff --git a/devtools/conda-envs/spyrmsd-test-obabel-rx.yaml b/devtools/conda-envs/spyrmsd-test-obabel-rx.yaml
deleted file mode 100644
index af3cf83..0000000
--- a/devtools/conda-envs/spyrmsd-test-obabel-rx.yaml
+++ /dev/null
@@ -1,26 +0,0 @@
-name: spyrmsd
-channels:
-  - conda-forge
-dependencies:
-  # Base
-  - python
-  - setuptools
-
-  # Maths
-  - numpy
-  - scipy
-  - rustworkx
-
-  # Chemistry
-  - openbabel
-
-  # Testing
-  - pytest
-  - pytest-cov
-  - pytest-benchmark
-
-  # Dev
-  - mypy
-  - flake8
-  - black
-  - codecov
diff --git a/devtools/conda-envs/spyrmsd-test-rdkit-all.yaml b/devtools/conda-envs/spyrmsd-test-rdkit-all.yaml
index d1d7ced..b2f2b9b 100644
--- a/devtools/conda-envs/spyrmsd-test-rdkit-all.yaml
+++ b/devtools/conda-envs/spyrmsd-test-rdkit-all.yaml
@@ -5,6 +5,7 @@ channels:
 dependencies:
   # Base
   - python
+  - pip
   - setuptools
 
   # Maths
diff --git a/devtools/conda-envs/spyrmsd-test-rdkit-gt.yaml b/devtools/conda-envs/spyrmsd-test-rdkit-gt.yaml
deleted file mode 100644
index 06e9dd9..0000000
--- a/devtools/conda-envs/spyrmsd-test-rdkit-gt.yaml
+++ /dev/null
@@ -1,27 +0,0 @@
-name: spyrmsd
-channels:
-  - conda-forge
-  - rdkit
-dependencies:
-  # Base
-  - python
-  - setuptools
-
-  # Maths
-  - numpy
-  - scipy
-  - graph-tool
-
-  # Chemistry
-  - rdkit
-
-  # Testing
-  - pytest
-  - pytest-cov
-  - pytest-benchmark
-
-  # Dev
-  - mypy
-  - flake8
-  - black
-  - codecov
diff --git a/devtools/conda-envs/spyrmsd-test-rdkit-nx.yaml b/devtools/conda-envs/spyrmsd-test-rdkit-nogt.yaml
similarity index 92%
rename from devtools/conda-envs/spyrmsd-test-rdkit-nx.yaml
rename to devtools/conda-envs/spyrmsd-test-rdkit-nogt.yaml
index e35ae1b..72b5429 100644
--- a/devtools/conda-envs/spyrmsd-test-rdkit-nx.yaml
+++ b/devtools/conda-envs/spyrmsd-test-rdkit-nogt.yaml
@@ -5,12 +5,14 @@ channels:
 dependencies:
   # Base
   - python
+  - pip
   - setuptools
 
   # Maths
   - numpy
   - scipy
   - networkx>=2
+  - rustworkx
 
   # Chemistry
   - rdkit
diff --git a/devtools/conda-envs/spyrmsd-test-rdkit-rx.yaml b/devtools/conda-envs/spyrmsd-test-rdkit-rx.yaml
deleted file mode 100644
index e752f97..0000000
--- a/devtools/conda-envs/spyrmsd-test-rdkit-rx.yaml
+++ /dev/null
@@ -1,27 +0,0 @@
-name: spyrmsd
-channels:
-  - conda-forge
-  - rdkit
-dependencies:
-  # Base
-  - python
-  - setuptools
-
-  # Maths
-  - numpy
-  - scipy
-  - rustworkx
-
-  # Chemistry
-  - rdkit
-
-  # Testing
-  - pytest
-  - pytest-cov
-  - pytest-benchmark
-
-  # Dev
-  - mypy
-  - flake8
-  - black
-  - codecov
diff --git a/tests/conftest.py b/tests/conftest.py
index 9c04244..ed7949b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -5,9 +5,18 @@
     https://docs.pytest.org/en/latest/example/simple.html
 """
 
+import os
+import warnings
+from collections import namedtuple
+
 import numpy as np
 import pytest
 
+import spyrmsd
+from spyrmsd import io
+
+Mol = namedtuple("Mol", ["mol", "name", "n_atoms", "n_bonds", "n_h"])
+
 
 def pytest_addoption(parser):
     parser.addoption(
@@ -69,3 +78,91 @@ def pytest_generate_tests(metafunc):
         n = metafunc.config.getoption("--n-tests")
 
         metafunc.parametrize("idx", np.random.randint(0, pytest.n_systems, size=n))
+
+
+@pytest.fixture(autouse=True, params=spyrmsd.available_backends)
+def set_backend(request):
+    # Capture warning when trying to switch to the same backend
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        spyrmsd.set_backend(request.param)
+
+
+@pytest.fixture(scope="session")
+def molpath():
+    fdir = os.path.dirname(os.path.abspath(__file__))
+    return os.path.join(fdir, f"data{os.sep}molecules")
+
+
+@pytest.fixture
+def benzene(molpath):
+    mol = io.loadmol(os.path.join(molpath, "benzene.sdf"))
+    return Mol(mol, "benzene", 12, 12, 6)
+
+
+@pytest.fixture
+def pyridine(molpath):
+    mol = io.loadmol(os.path.join(molpath, "pyridine.sdf"))
+    return Mol(mol, "pyridine", 11, 11, 5)
+
+
+@pytest.fixture(
+    params=[
+        # (name, n_atoms, n_bonds, n_h)
+        ("benzene", 12, 12, 6),
+        ("ethanol", 9, 8, 6),
+        ("pyridine", 11, 11, 5),
+        ("dialanine", 23, 22, 12),
+    ]
+)
+def mol(request, molpath):
+    """
+    Load molecule as sPyRMSD molecule.
+    """
+
+    name, n_atoms, n_bonds, n_h = request.param
+
+    mol = io.loadmol(os.path.join(molpath, f"{name}.sdf"))
+
+    return Mol(mol, name, n_atoms, n_bonds, n_h)
+
+
+@pytest.fixture
+def rawmol(mol, molpath):
+    """
+    Load molecule as a molecule of the current molecular I/O library.
+    """
+
+    RawMol = namedtuple(
+        "RawMol", ["mol", "rawmol", "name", "n_atoms", "n_bonds", "n_h"]
+    )
+
+    rawmol = io.load(os.path.join(molpath, f"{mol.name}.sdf"))
+
+    return RawMol(mol.mol, rawmol, mol.name, mol.n_atoms, mol.n_bonds, mol.n_h)
+
+
+@pytest.fixture
+def trps(molpath):
+    trp_list = []
+    for i in range(6):
+        trp_list.append(io.loadmol(os.path.join(molpath, f"trp{i}.pdb")))
+
+    return trp_list
+
+
+@pytest.fixture
+def docking_2viz(molpath):
+    mols = {}  # Dictionary (pose, molecule)
+    for i in [1, 2, 3]:
+        mols[i] = io.loadmol(os.path.join(molpath, f"2viz_{i}.sdf"))
+
+    return mols
+
+
+@pytest.fixture
+def docking_1cbr(molpath):
+    return [
+        io.loadmol(os.path.join(molpath, "1cbr_ligand.mol2")),
+        *io.loadallmols(os.path.join(molpath, "1cbr_docking.sdf")),
+    ]
diff --git a/tests/molecules.py b/tests/molecules.py
deleted file mode 100644
index 45cd20d..0000000
--- a/tests/molecules.py
+++ /dev/null
@@ -1,80 +0,0 @@
-import os
-from typing import Any, List, Tuple
-
-from spyrmsd import io, molecule
-
-fdir = os.path.dirname(os.path.abspath(__file__))
-molpath = os.path.join(fdir, "data/molecules/")
-
-
-def load(fname: str) -> Tuple[Any, molecule.Molecule]:
-    """
-    Load molecule from file.
-
-    Parameters
-    ----------
-    fname: str
-        Input file name
-
-    Returns
-    -------
-    Tuple[Any, molecule.Molecule]
-        Loaded molecule as `pybel.Molecule` or `rdkit.Chem.rdkem.Mol` and
-        `pyrmsd.molecule.Molecule`
-    """
-
-    fname = os.path.join(molpath, fname)
-
-    m = io.load(fname)
-
-    mol = io.to_molecule(m, adjacency=True)
-
-    return m, mol
-
-
-def loadall(fname: str) -> Tuple[List[Any], List[molecule.Molecule]]:
-    """
-    Load all molecule from file.
-
-    Parameters
-    ----------
-    fname: str
-        Input file name
-
-    Returns
-    -------
-    Tuple[List[Any], List[molecule.Molecule]]
-        Loaded molecule as `pybel.Molecule` or `rdkit.Chem.rdchem.Mol` and
-        `pyrmsd.molecule.Molecule`
-    """
-
-    fname = os.path.join(molpath, fname)
-
-    ms = io.loadall(fname)
-
-    mols = [io.to_molecule(m, adjacency=True) for m in ms]
-
-    return ms, mols
-
-
-obbenzene, benzene = load("benzene.sdf")
-obpyridine, pyridine = load("pyridine.sdf")
-obethanol, ethanol = load("ethanol.sdf")
-obdialanine, dialanine = load("dialanine.sdf")
-obsdf = [obbenzene, obpyridine, obethanol, obdialanine]
-sdf = [benzene, pyridine, ethanol, dialanine]
-
-allmolecules = sdf
-allobmolecules = obsdf
-
-obdocking_2viz, docking_2viz = {}, {}
-for i in [1, 2, 3]:
-    obdocking_2viz[i], docking_2viz[i] = load(f"2viz_{i}.sdf")
-
-obdocking_1cbr = [load("1cbr_ligand.mol2")[0], *loadall("1cbr_docking.sdf")[0]]
-docking_1cbr = [load("1cbr_ligand.mol2")[1], *loadall("1cbr_docking.sdf")[1]]
-
-intrp, trp = [], []
-for i in range(6):
-    intrp.append(load(f"trp{i}.pdb")[0])
-    trp.append(load(f"trp{i}.pdb")[1])
diff --git a/tests/test_backends.py b/tests/test_backends.py
new file mode 100644
index 0000000..fa77e36
--- /dev/null
+++ b/tests/test_backends.py
@@ -0,0 +1,103 @@
+import numpy as np
+import pytest
+
+import spyrmsd
+from spyrmsd import graph
+
+# TODO: Run even with two backends installed
+
+
+@pytest.mark.filterwarnings(
+    "ignore::UserWarning"
+)  # Silence "The backend is already" warning
+@pytest.mark.skipif(
+    # Run test if all supported backends are installed
+    not set(spyrmsd.graph._supported_backends) <= set(spyrmsd.available_backends),
+    reason="Not all of the required backends are installed",
+)
+def test_set_backend() -> None:
+    import graph_tool as gt
+    import networkx as nx
+    import rustworkx as rx
+
+    A = np.array([[0, 1, 1], [1, 0, 0], [1, 0, 1]])
+
+    spyrmsd.set_backend("networkx")
+    assert spyrmsd.get_backend() == "networkx"
+
+    Gnx = graph.graph_from_adjacency_matrix(A)
+    assert isinstance(Gnx, nx.Graph)
+
+    spyrmsd.set_backend("graph-tool")
+    assert spyrmsd.get_backend() == "graph_tool"
+
+    Ggt = graph.graph_from_adjacency_matrix(A)
+    assert isinstance(Ggt, gt.Graph)
+
+    spyrmsd.set_backend("rustworkx")
+    assert spyrmsd.get_backend() == "rustworkx"
+
+    Grx = graph.graph_from_adjacency_matrix(A)
+    assert isinstance(Grx, rx.PyGraph)
+
+
+def test_set_backend_unknown():
+    with pytest.raises(ValueError, match="backend is not recognized or supported"):
+        spyrmsd.set_backend("unknown")
+
+
+def test_set_backend_same():
+    current_backend = spyrmsd.get_backend()
+    with pytest.warns(UserWarning, match=f"The backend is already {current_backend}."):
+        spyrmsd.set_backend(current_backend)
+
+
+@pytest.mark.filterwarnings(
+    "ignore::UserWarning"
+)  # Silence "The backend is already" warning
+@pytest.mark.skipif(
+    # Run test if all supported backends are installed
+    not set(spyrmsd.graph._supported_backends) <= set(spyrmsd.available_backends),
+    reason="Not all of the required backends are installed",
+)
+def test_molecule_graph_cache(mol) -> None:
+    import graph_tool as gt
+    import networkx as nx
+    import rustworkx as rx
+
+    m = mol.mol
+
+    # Check molecules is in a clean state
+    assert len(m.G.items()) == 0
+    assert not m.stripped
+
+    spyrmsd.set_backend("networkx")
+    m.to_graph()
+
+    assert "networkx" in m.G.keys()
+    assert "graph_tool" not in m.G.keys()
+    assert "rustworkx" not in m.G.keys()
+
+    spyrmsd.set_backend("graph-tool")
+    m.to_graph()
+
+    assert "networkx" in m.G.keys()
+    assert "graph_tool" in m.G.keys()
+    assert "rustworkx" not in m.G.keys()
+
+    spyrmsd.set_backend("rustworkx")
+    m.to_graph()
+
+    assert "networkx" in m.G.keys()
+    assert "graph_tool" in m.G.keys()
+    assert "rustworkx" in m.G.keys()
+
+    # Make sure all backends (still) have a cache
+    assert isinstance(m.G["networkx"], nx.Graph)
+    assert isinstance(m.G["graph_tool"], gt.Graph)
+    assert isinstance(m.G["rustworkx"], rx.PyGraph)
+
+    # Strip molecule to ensure the cache is reset
+    m.strip()
+
+    assert len(m.G.items()) == 0
diff --git a/tests/test_graph.py b/tests/test_graph.py
index 1d68cef..c633a70 100644
--- a/tests/test_graph.py
+++ b/tests/test_graph.py
@@ -2,10 +2,9 @@
 import pytest
 
 import spyrmsd
-from spyrmsd import constants, graph, io, molecule
+from spyrmsd import constants, graph, io
 from spyrmsd.exceptions import NonIsomorphicGraphs
 from spyrmsd.graphs import _common as gc
-from tests import molecules
 
 
 def test_adjacency_matrix_from_atomic_coordinates_distance() -> None:
@@ -28,43 +27,42 @@ def test_adjacency_matrix_from_atomic_coordinates_distance() -> None:
     assert graph.num_edges(G) == 1
 
 
-@pytest.mark.parametrize(
-    "mol, n_bonds",
-    [(molecules.benzene, 12), (molecules.ethanol, 8), (molecules.dialanine, 22)],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_adjacency_matrix_from_atomic_coordinates(
-    mol: molecule.Molecule, n_bonds: int
-) -> None:
-    A = graph.adjacency_matrix_from_atomic_coordinates(mol.atomicnums, mol.coordinates)
+def test_adjacency_matrix_from_atomic_coordinates(mol) -> None:
+    A = graph.adjacency_matrix_from_atomic_coordinates(
+        mol.mol.atomicnums, mol.mol.coordinates
+    )
 
     G = graph.graph_from_adjacency_matrix(A)
 
-    assert graph.num_vertices(G) == len(mol)
-    assert graph.num_edges(G) == n_bonds
+    assert graph.num_vertices(G) == mol.n_atoms
+    assert graph.num_edges(G) == mol.n_bonds
 
 
-@pytest.mark.parametrize("mol", molecules.allobmolecules)
-def test_adjacency_matrix_from_mol(mol) -> None:
-    natoms = io.numatoms(mol)
-    nbonds = io.numbonds(mol)
+def test_adjacency_matrix_from_mol(rawmol) -> None:
+    natoms = io.numatoms(rawmol.rawmol)
+    nbonds = io.numbonds(rawmol.rawmol)
 
-    A = io.adjacency_matrix(mol)
+    assert natoms == rawmol.n_atoms
+    assert nbonds == rawmol.n_bonds
+
+    A = io.adjacency_matrix(rawmol.rawmol)
 
     assert A.shape == (natoms, natoms)
     assert np.all(A == A.T)
     assert np.sum(A) == nbonds * 2
 
-    for i, j in io.bonds(mol):
+    for i, j in io.bonds(rawmol.rawmol):
         assert A[i, j] == 1
 
 
-@pytest.mark.parametrize("mol", molecules.allobmolecules)
-def test_graph_from_adjacency_matrix(mol) -> None:
-    natoms = io.numatoms(mol)
-    nbonds = io.numbonds(mol)
+def test_graph_from_adjacency_matrix(rawmol) -> None:
+    natoms = io.numatoms(rawmol.rawmol)
+    nbonds = io.numbonds(rawmol.rawmol)
+
+    assert natoms == rawmol.n_atoms
+    assert nbonds == rawmol.n_bonds
 
-    A = io.adjacency_matrix(mol)
+    A = io.adjacency_matrix(rawmol.rawmol)
 
     assert A.shape == (natoms, natoms)
     assert np.all(A == A.T)
@@ -76,14 +74,13 @@ def test_graph_from_adjacency_matrix(mol) -> None:
     assert graph.num_edges(G) == nbonds
 
 
-@pytest.mark.parametrize(
-    "rawmol, mol", zip(molecules.allobmolecules, molecules.allmolecules)
-)
-def test_graph_from_adjacency_matrix_atomicnums(rawmol, mol) -> None:
-    natoms = io.numatoms(rawmol)
-    nbonds = io.numbonds(rawmol)
+def test_graph_from_adjacency_matrix_atomicnums(rawmol) -> None:
+    mol = rawmol.mol
+
+    natoms = io.numatoms(rawmol.rawmol)
+    nbonds = io.numbonds(rawmol.rawmol)
 
-    A = io.adjacency_matrix(rawmol)
+    A = io.adjacency_matrix(rawmol.rawmol)
 
     assert len(mol) == natoms
     assert mol.adjacency_matrix.shape == (natoms, natoms)
@@ -99,28 +96,44 @@ def test_graph_from_adjacency_matrix_atomicnums(rawmol, mol) -> None:
         assert graph.vertex_property(G, "aprops", idx) == atomicnum
 
 
-@pytest.mark.parametrize(
-    "G1, G2",
-    [
-        *[(graph.lattice(n, n), graph.lattice(n, n)) for n in range(2, 5)],
-        *[(graph.cycle(n), graph.cycle(n)) for n in range(2, 5)],
-    ],
-)
-def test_match_graphs_isomorphic(G1, G2) -> None:
+@pytest.mark.parametrize("n", list(range(2, 5)))
+def test_match_graphs_isomorphic_lattice(n) -> None:
+    G1 = graph.lattice(n, n)
+    G2 = graph.lattice(n, n)
+
     with pytest.warns(UserWarning, match=gc.warn_no_atomic_properties):
         isomorphisms = graph.match_graphs(G1, G2)
 
     assert len(isomorphisms) != 0
 
 
-@pytest.mark.parametrize(
-    "G1, G2",
-    [
-        *[(graph.lattice(n, n), graph.lattice(n + 1, n)) for n in range(2, 5)],
-        *[(graph.cycle(n), graph.cycle(n + 1)) for n in range(1, 5)],
-    ],
-)
-def test_match_graphs_not_isomorphic(G1, G2) -> None:
+@pytest.mark.parametrize("n", list(range(2, 5)))
+def test_match_graphs_isomorphic_cycle(n) -> None:
+    G1 = graph.cycle(n)
+    G2 = graph.cycle(n)
+
+    with pytest.warns(UserWarning, match=gc.warn_no_atomic_properties):
+        isomorphisms = graph.match_graphs(G1, G2)
+
+    assert len(isomorphisms) != 0
+
+
+@pytest.mark.parametrize("n", list(range(2, 5)))
+def test_match_graphs_not_isomorphic_lattice(n) -> None:
+    G1 = graph.lattice(n, n)
+    G2 = graph.lattice(n + 1, n)
+
+    with pytest.raises(
+        NonIsomorphicGraphs, match=gc.error_non_isomorphic_graphs
+    ), pytest.warns(UserWarning, match=gc.warn_no_atomic_properties):
+        graph.match_graphs(G1, G2)
+
+
+@pytest.mark.parametrize("n", range(2, 5))
+def test_match_graphs_not_isomorphic_cycle(n) -> None:
+    G1 = graph.cycle(n)
+    G2 = graph.cycle(n + 1)
+
     with pytest.raises(
         NonIsomorphicGraphs, match=gc.error_non_isomorphic_graphs
     ), pytest.warns(UserWarning, match=gc.warn_no_atomic_properties):
@@ -144,48 +157,15 @@ def test_build_graph_node_features(property) -> None:
     assert graph.num_edges(G) == 3
 
 
-@pytest.mark.skipif(
-    spyrmsd.get_backend() != "graph_tool",
-    reason="NetworkX supports all Python objects as node properties.",
-)
 def test_build_graph_node_features_unsupported() -> None:
+    if spyrmsd.get_backend() != "graph-tool":
+        pytest.skip(
+            "NetworkX and RustworkX support all Python objects as node properties."
+        )
+
     A = np.array([[0, 1, 1], [1, 0, 0], [1, 0, 1]])
 
     property = [True, False, True]
 
     with pytest.raises(ValueError, match="Unsupported property type:"):
         _ = graph.graph_from_adjacency_matrix(A, property)
-
-
-@pytest.mark.skipif(
-    # Run test if all supported backends are installed
-    not set(spyrmsd.graph._supported_backends) <= set(spyrmsd.available_backends),
-    reason="Not all of the required backends are installed",
-)
-def test_set_backend() -> None:
-    import graph_tool as gt
-    import networkx as nx
-    import rustworkx as rx
-
-    A = np.array([[0, 1, 1], [1, 0, 0], [1, 0, 1]])
-
-    spyrmsd.set_backend("networkx")
-    assert spyrmsd.get_backend() == "networkx"
-
-    Gnx = graph.graph_from_adjacency_matrix(A)
-    assert isinstance(Gnx, nx.Graph)
-
-    spyrmsd.set_backend("graph-tool")
-    assert spyrmsd.get_backend() == "graph_tool"
-
-    Ggt = graph.graph_from_adjacency_matrix(A)
-    assert isinstance(Ggt, gt.Graph)
-
-    spyrmsd.set_backend("rustworkx")
-    assert spyrmsd.get_backend() == "rustworkx"
-
-    Grx = graph.graph_from_adjacency_matrix(A)
-    assert isinstance(Grx, rx.PyGraph)
-
-    with pytest.raises(ValueError, match="backend is not recognized or supported"):
-        spyrmsd.set_backend("unknown")
diff --git a/tests/test_hungarian.py b/tests/test_hungarian.py
index 0532f1f..cc83c99 100644
--- a/tests/test_hungarian.py
+++ b/tests/test_hungarian.py
@@ -3,14 +3,12 @@
 import numpy as np
 import pytest
 
-from spyrmsd import hungarian, molecule
-from tests import molecules
+from spyrmsd import hungarian
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_cost_mtx(mol: molecule.Molecule):
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_cost_mtx(mol):
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     M = hungarian.cost_mtx(mol1.coordinates, mol2.coordinates)
 
@@ -27,10 +25,9 @@ def test_cost_mtx(mol: molecule.Molecule):
             assert M[i, j] == pytest.approx(np.dot(ab, ab))
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_optimal_assignement_same_molecule(mol: molecule.Molecule):
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_optimal_assignement_same_molecule(mol):
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     assert len(mol1) == len(mol2)
 
diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index b0f6064..c216ee6 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -1,107 +1,90 @@
 import copy
 import os
 from collections import defaultdict
-from typing import DefaultDict, List, Tuple
+from typing import DefaultDict
 
 import numpy as np
 import pytest
 
-import spyrmsd
 from spyrmsd import constants, graph, io, molecule, utils
-from tests import molecules
 
 
-# atoms is a list of atomic numbers and atom counts
-@pytest.mark.parametrize(
-    "mol, atoms",
-    [
-        (molecules.benzene, [(1, 6), (6, 6)]),
-        (molecules.ethanol, [(1, 6), (6, 2), (8, 1)]),
-        (molecules.dialanine, [(1, 12), (6, 6), (7, 2), (8, 3)]),
-    ],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_load(mol: molecule.Molecule, atoms: List[Tuple[int, int]]) -> None:
-    n = sum([n_atoms for _, n_atoms in atoms])
+def test_load(mol) -> None:
+    # Atoms for each type
+    atoms = {
+        "benzene": [(1, 6), (6, 6)],
+        "ethanol": [(1, 6), (6, 2), (8, 1)],
+        "pyridine": [(1, 5), (6, 5), (7, 1)],
+        "dialanine": [(1, 12), (6, 6), (7, 2), (8, 3)],
+    }
 
-    assert len(mol) == n
-    assert mol.atomicnums.shape == (n,)
-    assert mol.coordinates.shape == (n, 3)
+    n = sum([n_atoms for _, n_atoms in atoms[mol.name]])
+
+    assert mol.n_atoms == n
+    assert mol.mol.atomicnums.shape == (n,)
+    assert mol.mol.coordinates.shape == (n, 3)
 
     # Count number of atoms of different elements
     atomcount: DefaultDict[int, int] = defaultdict(int)
-    for atomicnum in mol.atomicnums:
+    for atomicnum in mol.mol.atomicnums:
         atomcount[atomicnum] += 1
 
-    assert len(atomcount) == len(atoms)
+    assert len(atomcount) == len(atoms[mol.name])
 
-    for Z, n_atoms in atoms:
+    for Z, n_atoms in atoms[mol.name]:
         assert atomcount[Z] == n_atoms
 
 
-def test_loadall() -> None:
-    path = os.path.join(molecules.molpath, "1cbr_docking.sdf")
+def test_loadall(molpath) -> None:
+    path = os.path.join(molpath, "1cbr_docking.sdf")
 
     mols = io.loadall(path)
 
     assert len(mols) == 10
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_molecule_translate(mol: molecule.Molecule) -> None:
-    mt = copy.deepcopy(mol)
+def test_molecule_translate(mol) -> None:
+    mt = copy.deepcopy(mol.mol)
 
     t = np.array([0.5, 1.1, -0.1])
     mt.translate(t)
 
-    for tcoord, coord in zip(mt.coordinates, mol.coordinates):
+    for tcoord, coord in zip(mt.coordinates, mol.mol.coordinates):
         assert np.allclose(tcoord - t, coord)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_molecule_rotate_z(mol: molecule.Molecule) -> None:
+def test_molecule_rotate_z(mol) -> None:
     z_axis = np.array([0, 0, 1])
 
     for angle in [0, 45, 90]:
-        rotated = np.zeros((len(mol), 3))
-        for i, coord in enumerate(mol.coordinates):
+        rotated = np.zeros((mol.n_atoms, 3))
+        for i, coord in enumerate(mol.mol.coordinates):
             rotated[i] = utils.rotate(coord, angle, z_axis, units="deg")
 
-        mol.rotate(angle, z_axis, units="deg")
-
-        assert np.allclose(mol.coordinates, rotated)
+        mol.mol.rotate(angle, z_axis, units="deg")
 
-        # Reset
-        mol.rotate(-angle, z_axis, units="deg")
+        assert np.allclose(mol.mol.coordinates, rotated)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_molecule_rotate(mol: molecule.Molecule) -> None:
+def test_molecule_rotate(mol) -> None:
     axis = np.random.rand(3)
 
     for angle in np.random.rand(10) * 180:
-        rotated = np.zeros((len(mol), 3))
-        for i, coord in enumerate(mol.coordinates):
+        rotated = np.zeros((mol.n_atoms, 3))
+        for i, coord in enumerate(mol.mol.coordinates):
             rotated[i] = utils.rotate(coord, angle, axis, units="deg")
 
-        mol.rotate(angle, axis, units="deg")
+        mol.mol.rotate(angle, axis, units="deg")
 
-        assert np.allclose(mol.coordinates, rotated)
+        assert np.allclose(mol.mol.coordinates, rotated)
 
-        # Reset
-        mol.rotate(-angle, axis, units="deg")
 
+def test_molecule_center_of_geometry_benzene(benzene) -> None:
+    assert np.allclose(benzene.mol.center_of_geometry(), np.zeros(3))
 
-def test_molecule_center_of_geometry_benzene() -> None:
-    mol = molecules.benzene
 
-    assert np.allclose(mol.center_of_geometry(), np.zeros(3))
-
-
-def test_molecule_center_of_mass_benzene() -> None:
-    mol = molecules.benzene
-
-    assert np.allclose(mol.center_of_mass(), np.zeros(3))
+def test_molecule_center_of_mass_benzene(benzene) -> None:
+    assert np.allclose(benzene.mol.center_of_mass(), np.zeros(3))
 
 
 def test_molecule_center_of_mass_H2() -> None:
@@ -127,49 +110,28 @@ def test_molecule_center_of_mass_HF() -> None:
     assert np.allclose(mol.center_of_mass(), np.array([0, 0, z_com]))
 
 
-@pytest.mark.parametrize(
-    "mol, n_atoms, stripped",
-    [
-        (molecules.benzene, 12, 6),
-        (molecules.ethanol, 9, 6),
-        (molecules.dialanine, 23, 12),
-    ],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_molecule_strip(mol: molecule.Molecule, n_atoms: int, stripped: int) -> None:
-    m = copy.deepcopy(mol)
+def test_molecule_strip(mol) -> None:
+    m = copy.deepcopy(mol.mol)
 
-    assert len(m) == n_atoms
+    assert len(m) == mol.n_atoms
 
     m.strip()
 
-    assert len(m) == n_atoms - stripped
+    assert len(m) == mol.n_atoms - mol.n_h
 
 
-@pytest.mark.parametrize(
-    "mol, n_bonds",
-    [(molecules.benzene, 12), (molecules.ethanol, 8), (molecules.dialanine, 22)],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_graph_from_adjacency_matrix(mol: molecule.Molecule, n_bonds: int) -> None:
-    G = mol.to_graph()
+def test_graph_from_adjacency_matrix(mol) -> None:
+    G = mol.mol.to_graph()
 
-    assert graph.num_vertices(G) == len(mol)
-    assert graph.num_edges(G) == n_bonds
+    assert graph.num_vertices(G) == mol.n_atoms
+    assert graph.num_edges(G) == mol.n_bonds
 
-    for idx, atomicnum in enumerate(mol.atomicnums):
+    for idx, atomicnum in enumerate(mol.mol.atomicnums):
         assert graph.vertex_property(G, "aprops", idx) == atomicnum
 
 
-@pytest.mark.parametrize(
-    "mol, n_bonds",
-    [(molecules.benzene, 12), (molecules.ethanol, 8), (molecules.dialanine, 22)],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_graph_from_atomic_coordinates_perception(
-    mol: molecule.Molecule, n_bonds: int
-) -> None:
-    m = copy.deepcopy(mol)
+def test_graph_from_atomic_coordinates_perception(mol) -> None:
+    m = copy.deepcopy(mol.mol)
 
     delattr(m, "adjacency_matrix")
     m.G = {}
@@ -178,10 +140,10 @@ def test_graph_from_atomic_coordinates_perception(
         # Uses automatic bond perception
         G = m.to_graph()
 
-        assert graph.num_vertices(G) == len(m)
-        assert graph.num_edges(G) == n_bonds
+        assert graph.num_vertices(G) == mol.n_atoms
+        assert graph.num_edges(G) == mol.n_bonds
 
-        for idx, atomicnum in enumerate(mol.atomicnums):
+        for idx, atomicnum in enumerate(mol.mol.atomicnums):
             assert graph.vertex_property(G, "aprops", idx) == atomicnum
 
 
@@ -190,13 +152,13 @@ def test_graph_from_atomic_coordinates_perception(
     [True, False],
     ids=["adjacency", "no_adjacency"],
 )
-def test_from_obmol(adjacency):
+def test_from_obmol(molpath, adjacency):
     pytest.importorskip("openbabel")
 
     from spyrmsd.optional import obabel as ob
 
     # Load molecules with OpenBabel
-    path = os.path.join(molecules.molpath, "1cbr_docking.sdf")
+    path = os.path.join(molpath, "1cbr_docking.sdf")
     mols = ob.loadall(path)
 
     # Convert OpenBabel molecules to spyrmsd molecules
@@ -220,13 +182,13 @@ def test_from_obmol(adjacency):
     [True, False],
     ids=["adjacency", "no_adjacency"],
 )
-def test_from_rdmol(adjacency):
+def test_from_rdmol(molpath, adjacency):
     pytest.importorskip("rdkit")
 
     from spyrmsd.optional import rdkit as rd
 
     # Load molecules with RDKit
-    path = os.path.join(molecules.molpath, "1cbr_docking.sdf")
+    path = os.path.join(molpath, "1cbr_docking.sdf")
     mols = rd.loadall(path)
 
     # Convert OpenBabel molecules to spyrmsd molecules
@@ -243,43 +205,3 @@ def test_from_rdmol(adjacency):
             with pytest.raises(AttributeError):
                 # No adjacency_matrix attribute
                 mol.adjacency_matrix
-
-
-@pytest.mark.skipif(
-    # Run test if all supported backends are installed
-    not set(spyrmsd.graph._supported_backends) <= set(spyrmsd.available_backends),
-    reason="Not all of the required backends are installed",
-)
-@pytest.mark.parametrize(
-    "mol",
-    [(molecules.benzene), (molecules.ethanol), (molecules.dialanine)],
-    ids=["benzene", "ethanol", "dialanine"],
-)
-def test_molecule_graph_cache(mol) -> None:
-    import graph_tool as gt
-    import networkx as nx
-    import rustworkx as rx
-
-    ## Graph cache persists from previous tests, manually reset them
-    mol.G = {}
-    spyrmsd.set_backend("networkx")
-    mol.to_graph()
-
-    assert isinstance(mol.G["networkx"], nx.Graph)
-    assert "graph_tool" not in mol.G.keys()
-
-    spyrmsd.set_backend("graph-tool")
-    mol.to_graph()
-
-    spyrmsd.set_backend("rustworkx")
-    mol.to_graph()
-
-    ## Make sure all backends (still) have a cache
-    assert isinstance(mol.G["networkx"], nx.Graph)
-    assert isinstance(mol.G["graph_tool"], gt.Graph)
-    assert isinstance(mol.G["rustworkx"], rx.PyGraph)
-
-    ## Strip the molecule to ensure the cache is reset
-    mol.strip()
-
-    assert len(mol.G.items()) == 0
diff --git a/tests/test_qcp.py b/tests/test_qcp.py
index 3556eb0..2b72d93 100644
--- a/tests/test_qcp.py
+++ b/tests/test_qcp.py
@@ -5,14 +5,12 @@
 import numpy as np
 import pytest
 
-from spyrmsd import molecule, qcp
-from tests import molecules
+from spyrmsd import qcp
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_M_mtx(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_M_mtx(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     # Build rotated coordinate set
     mol2.rotate(10, np.random.rand(3))
@@ -29,10 +27,9 @@ def S(i, j):
         assert M[i, j] == pytest.approx(S(i, j))
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_K_mtx(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_K_mtx(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     # Build rotated coordinate set
     mol2.rotate(10, np.random.rand(3))
@@ -63,10 +60,9 @@ def test_lambda_max(
     assert qcp.lambda_max(*input) == pytest.approx(result)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_lambda_max_eig(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_lambda_max_eig(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     # Build rotated coordinate set
     mol2.rotate(10, np.random.rand(3))
diff --git a/tests/test_rmsd.py b/tests/test_rmsd.py
index f442d73..57e45c7 100644
--- a/tests/test_rmsd.py
+++ b/tests/test_rmsd.py
@@ -4,8 +4,7 @@
 import numpy as np
 import pytest
 
-from spyrmsd import molecule, qcp, rmsd
-from tests import molecules
+from spyrmsd import qcp, rmsd
 
 
 @pytest.fixture(
@@ -42,9 +41,9 @@ def lambda_max_failure(Ga, Gb, c2, c1, c0):
 @pytest.mark.parametrize(
     "t, RMSD", [(0.0, 0.0), (1.0, 1.0), (2.0, 2.0)], ids=["t0", "t1", "t2"]
 )
-def test_rmsd_benzene(t: float, RMSD: float) -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+def test_rmsd_benzene(benzene, t: float, RMSD: float) -> None:
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     mol2.translate(np.array([0, 0, t]))
 
@@ -60,9 +59,9 @@ def test_rmsd_benzene(t: float, RMSD: float) -> None:
     [(1, 2, 2.60065218), (1, 3, 9.94411523), (2, 3, 9.4091711)],
     ids=["1-2", "1-3", "2-3"],
 )
-def test_rmsd_2viz(i: int, j: int, result: float) -> None:
-    moli = copy.deepcopy(molecules.docking_2viz[i])
-    molj = copy.deepcopy(molecules.docking_2viz[j])
+def test_rmsd_2viz(docking_2viz, i: int, j: int, result: float) -> None:
+    moli = copy.deepcopy(docking_2viz[i])
+    molj = copy.deepcopy(docking_2viz[j])
 
     assert rmsd.rmsd(
         moli.coordinates, molj.coordinates, moli.atomicnums, molj.atomicnums
@@ -76,9 +75,9 @@ def test_rmsd_2viz(i: int, j: int, result: float) -> None:
     [(1, 2, 2.65327362), (1, 3, 10.11099065), (2, 3, 9.57099612)],
     ids=["1-2", "1-3", "2-3"],
 )
-def test_rmsd_2viz_stripped(i: int, j: int, result: float) -> None:
-    moli = copy.deepcopy(molecules.docking_2viz[i])
-    molj = copy.deepcopy(molecules.docking_2viz[j])
+def test_rmsd_2viz_stripped(docking_2viz, i: int, j: int, result: float) -> None:
+    moli = copy.deepcopy(docking_2viz[i])
+    molj = copy.deepcopy(docking_2viz[j])
 
     moli.strip()
     molj.strip()
@@ -88,9 +87,9 @@ def test_rmsd_2viz_stripped(i: int, j: int, result: float) -> None:
     ) == pytest.approx(result)
 
 
-def test_rmsd_centred_benzene() -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+def test_rmsd_centred_benzene(benzene) -> None:
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     mol2.translate(np.array([0, 0, 1]))
 
@@ -107,10 +106,9 @@ def test_rmsd_centred_benzene() -> None:
     ) == pytest.approx(0)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_rmsd_minimize(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_rmsd_minimize(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     assert rmsd.rmsd(
         mol1.coordinates, mol2.coordinates, mol1.atomicnums, mol2.atomicnums
@@ -151,9 +149,9 @@ def test_rmsd_minimize(mol: molecule.Molecule) -> None:
     [(1, 2, 1.95277757), (1, 3, 3.11801105), (2, 3, 2.98609758)],
     ids=["1-2", "1-3", "2-3"],
 )
-def test_rmsd_qcp_2viz(i: int, j: int, result: float) -> None:
-    moli = copy.deepcopy(molecules.docking_2viz[i])
-    molj = copy.deepcopy(molecules.docking_2viz[j])
+def test_rmsd_qcp_2viz(docking_2viz, i: int, j: int, result: float) -> None:
+    moli = copy.deepcopy(docking_2viz[i])
+    molj = copy.deepcopy(docking_2viz[j])
 
     assert rmsd.rmsd(
         moli.coordinates,
@@ -171,9 +169,9 @@ def test_rmsd_qcp_2viz(i: int, j: int, result: float) -> None:
     [(1, 2, 1.98171656), (1, 3, 3.01799306), (2, 3, 2.82917355)],
     ids=["1-2", "1-3", "2-3"],
 )
-def test_rmsd_qcp_2viz_stripped(i: int, j: int, result: float) -> None:
-    moli = copy.deepcopy(molecules.docking_2viz[i])
-    molj = copy.deepcopy(molecules.docking_2viz[j])
+def test_rmsd_qcp_2viz_stripped(docking_2viz, i: int, j: int, result: float) -> None:
+    moli = copy.deepcopy(docking_2viz[i])
+    molj = copy.deepcopy(docking_2viz[j])
 
     # Strip hydrogen atoms
     moli.strip()
@@ -208,9 +206,9 @@ def test_rmsd_qcp_2viz_stripped(i: int, j: int, result: float) -> None:
     ],
     ids=["1", "2", "3", "4", "5"],
 )
-def test_rmsd_qcp_protein(i: int, rmsd_dummy: float, rmsd_min: float):
-    mol0 = copy.deepcopy(molecules.trp[0])
-    mol = copy.deepcopy(molecules.trp[i])
+def test_rmsd_qcp_protein(trps, i: int, rmsd_dummy: float, rmsd_min: float):
+    mol0 = copy.deepcopy(trps[0])
+    mol = copy.deepcopy(trps[i])
 
     assert rmsd.rmsd(
         mol0.coordinates, mol.coordinates, mol0.atomicnums, mol.atomicnums
@@ -230,9 +228,9 @@ def test_rmsd_qcp_protein(i: int, rmsd_dummy: float, rmsd_min: float):
     [(60, 1e-4), (120, 1e-4), (180, 1e-4), (240, 1e-4), (300, 1e-4)],
     ids=["60", "120", "180", "240", "300"],
 )
-def test_rmsd_hungarian_benzene_rotated(angle: float, tol: float) -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+def test_rmsd_hungarian_benzene_rotated(benzene, angle: float, tol: float) -> None:
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     assert rmsd.rmsd(
         mol1.coordinates, mol2.coordinates, mol1.atomicnums, mol2.atomicnums
@@ -263,10 +261,10 @@ def test_rmsd_hungarian_benzene_rotated(angle: float, tol: float) -> None:
     ids=["60", "120", "180", "240", "300"],
 )
 def test_rmsd_hungarian_benzene_shifted_rotated(
-    d: float, angle: float, tol: float
+    benzene, d: float, angle: float, tol: float
 ) -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     mol2.translate([0, 0, d])
 
@@ -289,10 +287,9 @@ def test_rmsd_hungarian_benzene_shifted_rotated(
     ) == pytest.approx(abs(d), abs=tol)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_rmsd_hungarian_centred(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_rmsd_hungarian_centred(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     mol2.translate(np.random.rand(3))
 
@@ -310,10 +307,9 @@ def test_rmsd_hungarian_centred(mol: molecule.Molecule) -> None:
     ) == pytest.approx(0)
 
 
-@pytest.mark.parametrize("mol", molecules.allmolecules)
-def test_symmrmsd_centred(mol: molecule.Molecule) -> None:
-    mol1 = copy.deepcopy(mol)
-    mol2 = copy.deepcopy(mol)
+def test_symmrmsd_centred(mol) -> None:
+    mol1 = copy.deepcopy(mol.mol)
+    mol2 = copy.deepcopy(mol.mol)
 
     mol2.translate(np.random.rand(3))
 
@@ -341,9 +337,9 @@ def test_symmrmsd_centred(mol: molecule.Molecule) -> None:
 
 
 @pytest.mark.parametrize("angle", [60, 120, 180, 240, 300, 360])
-def test_symmrmsd_rotated_benzene(angle: float) -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+def test_symmrmsd_rotated_benzene(benzene, angle: float) -> None:
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     mol2.rotate(angle, np.array([0, 0, 1]), units="deg")
 
@@ -375,9 +371,9 @@ def test_symmrmsd_rotated_benzene(angle: float) -> None:
 
 
 @pytest.mark.parametrize("angle", [60, 120, 180, 240, 300, 360])
-def test_symmrmsd_rotated_benzene_stripped(angle: float) -> None:
-    mol1 = copy.deepcopy(molecules.benzene)
-    mol2 = copy.deepcopy(molecules.benzene)
+def test_symmrmsd_rotated_benzene_stripped(benzene, angle: float) -> None:
+    mol1 = copy.deepcopy(benzene.mol)
+    mol2 = copy.deepcopy(benzene.mol)
 
     mol2.rotate(angle, np.array([0, 0, 1]), units="deg")
 
@@ -411,9 +407,9 @@ def test_symmrmsd_rotated_benzene_stripped(angle: float) -> None:
     ) == pytest.approx(0, abs=1e-4)
 
 
-def test_symmrmsd_atomicnums_matching_pyridine_stripped() -> None:
-    mol1 = copy.deepcopy(molecules.pyridine)
-    mol2 = copy.deepcopy(molecules.pyridine)
+def test_symmrmsd_atomicnums_matching_pyridine_stripped(pyridine) -> None:
+    mol1 = copy.deepcopy(pyridine.mol)
+    mol2 = copy.deepcopy(pyridine.mol)
 
     mol2.rotate(60, np.array([0, 0, 1]), units="deg")
 
@@ -485,9 +481,9 @@ def test_symmrmsd_atomicnums_matching_pyridine_stripped() -> None:
         "10-minimize",
     ],
 )
-def test_rmsd_symmrmsd(index: int, RMSD: float, minimize: bool) -> None:
-    molc = copy.deepcopy(molecules.docking_1cbr[0])
-    mol = copy.deepcopy(molecules.docking_1cbr[index])
+def test_rmsd_symmrmsd(docking_1cbr, index: int, RMSD: float, minimize: bool) -> None:
+    molc = copy.deepcopy(docking_1cbr[0])
+    mol = copy.deepcopy(docking_1cbr[index])
 
     molc.strip()
     mol.strip()
@@ -551,9 +547,11 @@ def test_rmsd_symmrmsd_disconnected_node() -> None:
     ],
     ids=["no_minimize", "minimize"],
 )
-def test_multi_spyrmsd(minimize: bool, referenceRMSDs: List[float]) -> None:
-    molc = copy.deepcopy(molecules.docking_1cbr[0])
-    mols = [copy.deepcopy(mol) for mol in molecules.docking_1cbr[1:]]
+def test_multi_spyrmsd(
+    docking_1cbr, minimize: bool, referenceRMSDs: List[float]
+) -> None:
+    molc = copy.deepcopy(docking_1cbr[0])
+    mols = [copy.deepcopy(mol) for mol in docking_1cbr[1:]]
 
     molc.strip()
 
@@ -611,9 +609,11 @@ def test_multi_spyrmsd(minimize: bool, referenceRMSDs: List[float]) -> None:
     ],
     ids=["no_minimize", "minimize"],
 )
-def test_symmrmsd_cache(minimize: bool, referenceRMSDs: List[float]) -> None:
-    molc = copy.deepcopy(molecules.docking_1cbr[0])
-    mols = [copy.deepcopy(mol) for mol in molecules.docking_1cbr[1:]]
+def test_symmrmsd_cache(
+    docking_1cbr, minimize: bool, referenceRMSDs: List[float]
+) -> None:
+    molc = copy.deepcopy(docking_1cbr[0])
+    mols = [copy.deepcopy(mol) for mol in docking_1cbr[1:]]
 
     molc.strip()
 
@@ -756,7 +756,7 @@ def test_issue_35_2():
     [(1, 2, 1.95277757), (1, 3, 3.11801105), (2, 3, 2.98609758)],
     ids=["1-2", "1-3", "2-3"],
 )
-def test_rmsd_atol(i: int, j: int, result: float):
+def test_rmsd_atol(docking_2viz, i: int, j: int, result: float):
     """
     Test usage of the :code:`atol` parameter for the QCP method.
 
@@ -764,8 +764,8 @@ def test_rmsd_atol(i: int, j: int, result: float):
     (https://github.com/RMeli/spyrmsd/issues/35)
     """
 
-    moli = copy.deepcopy(molecules.docking_2viz[i])
-    molj = copy.deepcopy(molecules.docking_2viz[j])
+    moli = copy.deepcopy(docking_2viz[i])
+    molj = copy.deepcopy(docking_2viz[j])
 
     # Check results are different from 0.0
     assert not result == pytest.approx(0.0)
@@ -792,9 +792,9 @@ def test_rmsd_atol(i: int, j: int, result: float):
 @pytest.mark.parametrize(
     "i, reference", [(1, 0.476858), (2, 1.68089), (3, 1.50267)], ids=["1", "2", "3"]
 )
-def test_symmrmsd_atol(i: bool, reference: float) -> None:
-    moli = copy.deepcopy(molecules.docking_1cbr[0])
-    molj = copy.deepcopy(molecules.docking_1cbr[i])
+def test_symmrmsd_atol(docking_1cbr, i: bool, reference: float) -> None:
+    moli = copy.deepcopy(docking_1cbr[0])
+    molj = copy.deepcopy(docking_1cbr[i])
 
     moli.strip()
     molj.strip()
@@ -824,11 +824,11 @@ def test_symmrmsd_atol(i: bool, reference: float) -> None:
     ) == pytest.approx(0.0)
 
 
-def test_symmrmsd_atol_multi() -> None:
+def test_symmrmsd_atol_multi(docking_1cbr) -> None:
     references = [0.476858, 1.68089, 1.50267]
 
-    molc = copy.deepcopy(molecules.docking_1cbr[0])
-    mols = [copy.deepcopy(mol) for mol in molecules.docking_1cbr[1:4]]
+    molc = copy.deepcopy(docking_1cbr[0])
+    mols = [copy.deepcopy(mol) for mol in docking_1cbr[1:4]]
 
     molc.strip()
 
@@ -901,9 +901,9 @@ def test_symmrmsd_atol_multi() -> None:
     ],
     ids=["no_minimize", "minimize"],
 )
-def test_rmsdwrapper_nosymm_protein(minimize: bool, referenceRMSDs: List[float]):
-    mol0 = copy.deepcopy(molecules.trp[0])
-    mols = [copy.deepcopy(mol) for mol in molecules.trp[1:]]
+def test_rmsdwrapper_nosymm_protein(trps, minimize: bool, referenceRMSDs: List[float]):
+    mol0 = copy.deepcopy(trps[0])
+    mols = [copy.deepcopy(mol) for mol in trps[1:]]
 
     RMSDs = rmsd.rmsdwrapper(mol0, mols, symmetry=False, minimize=minimize, strip=False)
 
@@ -948,9 +948,11 @@ def test_rmsdwrapper_nosymm_protein(minimize: bool, referenceRMSDs: List[float])
     ],
     ids=["minimize", "no_minimize"],
 )
-def test_rmsdwrapper_isomorphic(minimize: bool, referenceRMSDs: List[float]) -> None:
-    molref = copy.deepcopy(molecules.docking_1cbr[0])
-    mols = [copy.deepcopy(mol) for mol in molecules.docking_1cbr[1:]]
+def test_rmsdwrapper_isomorphic(
+    docking_1cbr, minimize: bool, referenceRMSDs: List[float]
+) -> None:
+    molref = copy.deepcopy(docking_1cbr[0])
+    mols = [copy.deepcopy(mol) for mol in docking_1cbr[1:]]
 
     RMSDs = rmsd.rmsdwrapper(molref, mols, minimize=minimize, strip=True)
 
@@ -964,9 +966,11 @@ def test_rmsdwrapper_isomorphic(minimize: bool, referenceRMSDs: List[float]) ->
     [(True, 0.476858), (False, 0.592256)],
     ids=["minimize", "no_minimize"],
 )
-def test_rmsdwrapper_single_molecule(minimize: bool, referenceRMSD: float) -> None:
-    molref = copy.deepcopy(molecules.docking_1cbr[0])
-    mols = copy.deepcopy(molecules.docking_1cbr[1])
+def test_rmsdwrapper_single_molecule(
+    docking_1cbr, minimize: bool, referenceRMSD: float
+) -> None:
+    molref = copy.deepcopy(docking_1cbr[0])
+    mols = copy.deepcopy(docking_1cbr[1])
 
     RMSD = rmsd.rmsdwrapper(molref, mols, minimize=minimize, strip=True)