diff --git a/scripts/discovery/overlay_to_pharmacophore/ReadMe.md b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md
new file mode 100644
index 0000000..80ee9ce
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/ReadMe.md
@@ -0,0 +1,46 @@
+# Pharmacophore Query Generator
+
+This tool provides the user with the ability to create pharmacophore queries from the results of a ligand overlay.
+
+The pharmacophore queries are produced to be used with CrossMiner.
+
+## Requirements
+
+- [CSD Python API](https://downloads.ccdc.cam.ac.uk/documentation/API/) installed.
+- Access to CSD CrossMiner and the feature definitions (`.cpf`) files.
+- Access to the CSD Ligand Overlay Tool.
+  
+## Licensing Requirements 
+
+CSD-Discovery, CSD-Enterprise and Research Partner suites would all be sufficient.
+
+## Instructions on Running
+
+### Feature Definitions
+
+The CrossMiner feature definition (`.cpf`) files are **not** shipped with this repo.
+Supply the location of the feature definitions from your CrossMiner installation with
+`-f`/`--feature_definitions`; this should be the directory containing the definition files
+(either directly, or in `any`/`protein`/`small_molecule` subdirectories).
+Usually, the location is `C:\users\<username>\CCDC\ccdc-software\csd-crossminer\feature_definitions`
+
+```
+python main.py -i <overlay_folder> -o <output_folder> -f <feature_definitions_folder>
+```
+
+The output folder (`-o`/`--output_folder`) is optional; if it is not supplied, the queries are
+written to a `queries` folder created in the current directory.
+
+### Options for Ligand Overlay Output
+
+* `cluster`: Cluster the similar pharmacophore features based on proximity
+* `projected`: Treat pharmacophore features as projected when appropriate e.g. acceptors and donors
+
+There is also the option to specify a specific Ligand Overlay from all the results. If this is not specified, all overlays are used.
+When all overlays are used, the pharmacophore query will be a union of all the features from all the overlays.
+These features are then clustered based on proximity AND prevalence.
+
+### Using the Queries Generated
+
+If you would like to use the queries generated with this tool, they can be opened in CrossMiner to run a search. 
+A file `crossminer_search.py` has also been provided which contains a Python function for the most simply kind of CrossMiner search.
diff --git a/scripts/discovery/overlay_to_pharmacophore/cluster.py b/scripts/discovery/overlay_to_pharmacophore/cluster.py
new file mode 100644
index 0000000..e49add2
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/cluster.py
@@ -0,0 +1,96 @@
+from collections import defaultdict
+
+import numpy as np
+
+from datastructures import PharmFeaturePoint
+
+
+def cluster_features(features: list[PharmFeaturePoint]) -> list[PharmFeaturePoint]:
+    """Cluster all features where relevant."""
+    groups = defaultdict(list)
+    for feature in features:
+        groups[feature.label].append(feature)
+
+    clustered_features = []
+    for label, group in groups.items():
+        if label == 'hydrophobe':
+            clustered_features.extend(_cluster_similar_features(group))
+        elif label in ('acceptor_projected', 'donor_projected'):
+            clustered_features.extend(_cluster_similar_features(group, check_vp=True))
+        else:
+            clustered_features.extend(group)
+    return clustered_features
+
+
+def _cluster_similar_features(
+        features: list[PharmFeaturePoint], cluster_radius: float = 2.0, check_vp: bool = False
+) -> list[PharmFeaturePoint]:
+    """
+    Cluster specific features that are close to each other.
+    For projected features, also check virtual point distances.
+
+    Args:
+        features: List of features to cluster
+        cluster_radius: Distance between features to be included in clustering (Ang)
+        check_vp: Whether to also require virtual points to be within the radius
+    """
+    if len(features) < 2:
+        return features
+
+    clusters = _connected_components(
+        features,
+        lambda a, b: _is_close(a, b, cluster_radius, check_vp),
+    )
+    return [_merge_cluster(cluster, check_vp) for cluster in clusters]
+
+
+def _is_close(
+        a: PharmFeaturePoint, b: PharmFeaturePoint, cluster_radius: float, check_vp: bool
+) -> bool:
+    """Whether two features are within ``cluster_radius`` (and, if ``check_vp``, their virtual points too)."""
+    if np.linalg.norm(a - b) >= cluster_radius:
+        return False
+    if check_vp and np.linalg.norm(a.virtual_point - b.virtual_point) >= cluster_radius:
+        return False
+    return True
+
+
+def _connected_components(
+        features: list[PharmFeaturePoint], is_linked
+) -> list[list[PharmFeaturePoint]]:
+    """
+    Single-linkage grouping: features are placed in the same cluster if a chain of
+    ``is_linked`` neighbours connects them.
+    """
+    unclustered = features.copy()
+    clusters = []
+
+    while unclustered:
+        cluster = [unclustered.pop(0)]
+        # Grow the cluster breadth-first: any unclustered feature linked to a member joins it.
+        i = 0
+        while i < len(cluster):
+            member = cluster[i]
+            remaining = []
+            for feature in unclustered:
+                if is_linked(feature, member):
+                    cluster.append(feature)
+                else:
+                    remaining.append(feature)
+            unclustered = remaining
+            i += 1
+        clusters.append(cluster)
+
+    return clusters
+
+
+def _merge_cluster(cluster: list[PharmFeaturePoint], check_vp: bool) -> PharmFeaturePoint:
+    """Collapse a cluster into a single feature at its centroid (singletons are returned unchanged)."""
+    if len(cluster) == 1:
+        return cluster[0]
+
+    centroid = np.mean([c.coordinates for c in cluster], axis=0).round(4)
+    vp_centroid = (
+        np.mean([c.virtual_point for c in cluster], axis=0).round(4) if check_vp else None
+    )
+    return PharmFeaturePoint(centroid, label=cluster[0].label, virtual_point=vp_centroid)
diff --git a/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py
new file mode 100644
index 0000000..962462e
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/crossminer_search.py
@@ -0,0 +1,18 @@
+from pathlib import Path
+
+from ccdc.pharmacophore import Pharmacophore
+
+def search(query_file: Path, database_file: Path):
+    settings = Pharmacophore.Search.Settings()
+    settings.max_hit_structures = 20
+    settings.max_hits_per_structure = 1
+    settings.max_hit_rmsd = 1.0
+    searcher = Pharmacophore.Search(settings)
+    feature_db = Pharmacophore.FeatureDatabase.from_file(database_file)
+    query = Pharmacophore.Query.from_file(str(query_file))
+    hits = searcher.search(
+        model=query,
+        database=feature_db,
+        verbose=True,
+    )
+    return hits
diff --git a/scripts/discovery/overlay_to_pharmacophore/datastructures.py b/scripts/discovery/overlay_to_pharmacophore/datastructures.py
new file mode 100644
index 0000000..0f451c2
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/datastructures.py
@@ -0,0 +1,114 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Iterable, Optional
+
+import numpy as np
+
+
+@dataclass
+class FeatureTolerances:
+    """
+    Allowed feature types and their tolerances.
+
+    Each field is a valid feature label; its value is the tolerance:
+    a single weight, or a (parent, virtual_point) pair for projected features.
+    """
+    acceptor: float = 1.0
+    acceptor_projected: tuple[float, float] = (0.8, 0.8)
+    donor_projected: tuple[float, float] = (1.0, 1.0)
+    hydrophobe: float = 1.0
+    ring_planar_projected: tuple[float, float] = (1.0, 1.0)
+    ring_non_planar: tuple[float, float] = (1.0, 1.0)
+    halogen: float = 1.0
+
+    def __getitem__(self, key: str) -> float | tuple[float, float]:
+        return getattr(self, key)
+
+
+class PharmFeaturePoint(np.ndarray):
+    def __new__(
+        cls,
+        *coordinates: float | Iterable[float],
+        label: Optional[str] = None,
+        virtual_point: Optional[np.ndarray] = None,
+    ):
+        if len(coordinates) == 1:
+            arr = np.asarray(coordinates[0], dtype=float)
+        elif len(coordinates) == 3:
+            arr = np.asarray(coordinates, dtype=float)
+        else:
+            raise TypeError("Coordinates must be either an iterable of length three, or three floats")
+        if arr.shape != (3,):
+            raise ValueError("Coordinates must be a 3-element array")
+        obj = arr.view(cls)
+        obj.label = label
+        obj.virtual_point = virtual_point
+        return obj
+
+    def __array_finalize__(self, obj):
+        if obj is None:
+            return
+        self.label = getattr(obj, 'label', None)
+        self.virtual_point = getattr(obj, 'virtual_point', None)
+
+    def __repr__(self):
+        return (
+            f"PharmFeaturePoint({self.x}, {self.y}, {self.z}, "
+            f"label={self.label}, virtual_points={self.virtual_point}), "
+        )
+
+    def __str__(self):
+        return self.__repr__()
+
+    @property
+    def coordinates(self) -> np.ndarray:
+        return np.asarray(self)
+
+    @property
+    def x(self):
+        return self[0]
+
+    @property
+    def y(self):
+        return self[1]
+
+    @property
+    def z(self):
+        return self[2]
+
+    @property
+    def tolerance(self) -> float | tuple[float, float]:
+        """
+        Get the tolerance for the feature based on its label.
+        Returns:
+            A single tolerance for features with a single tolerance or a tuple for those with two tolerances.
+        """
+        if self.label is None:
+            raise LookupError("Feature label must be set to determine tolerances.")
+        return FeatureTolerances()[self.label]
+
+    @property
+    def weight_parent(self) -> float:
+        if isinstance(self.tolerance, float):
+            return self.tolerance
+        elif isinstance(self.tolerance, tuple):
+            return self.tolerance[0]
+        else:
+            raise ValueError("Incorrect tolerances loaded for feature point.")
+
+    @property
+    def weight_vp(self) -> float:
+        return self.tolerance[1]
+
+
+
+@dataclass
+class OverlayData:
+    input_folder: Path
+    output_folder: Path
+    # Pharmacophore file from the pharmacophores folder
+    pharm_file: Path
+    # Overlay solution file (the chosen one from the many produced)
+    overlay_file: Path
diff --git a/scripts/discovery/overlay_to_pharmacophore/main.py b/scripts/discovery/overlay_to_pharmacophore/main.py
new file mode 100644
index 0000000..1e4954b
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/main.py
@@ -0,0 +1,76 @@
+import argparse
+from pathlib import Path
+
+from cluster import cluster_features
+from datastructures import OverlayData
+from overlay import OverlayToPharmFeatures
+from write_query import FeaturesToCrossMinerQuery
+
+
+def str_to_bool(value: str) -> bool:
+    return value.lower() in {'t', 'true', '1', 'yes', 'y'}
+
+
+def parse_args():
+    parser = argparse.ArgumentParser(
+        description="Create Pharmacophore Features from a Ligand Overlay"
+    )
+    parser.add_argument('-i', '--input_folder', type=str, help='Input file(s) path.')
+    parser.add_argument('-o', '--output_folder', type=str, default=None,
+                        help="Output folder path. Defaults to a 'queries' folder in the current directory.")
+    parser.add_argument('-f', '--feature_definitions', type=str, required=True,
+                        help='Path to the directory containing the CrossMiner feature definition (.cpf) files.')
+    parser.add_argument('-c', '--cluster', type=str_to_bool, default=False,
+                        help='Cluster features if they are close together or common across multiple inputs.')
+    parser.add_argument('-p', '--projected', type=str_to_bool, default=False,
+                        help='Use projected acceptor features or point features.')
+    parser.add_argument('-id', '--overlay_id', type=int, default=0,
+                        help='Overlay ID to process. If 0 or not specified, all overlays will be processed.')
+
+    return parser.parse_args()
+
+
+def main():
+    args = parse_args()
+
+    input_folder = Path(args.input_folder)
+    if not input_folder.exists():
+        raise FileNotFoundError(f"Input folder {input_folder} does not exist.")
+
+    feature_definitions = Path(args.feature_definitions)
+    if not feature_definitions.is_dir():
+        raise FileNotFoundError(f"Feature definitions folder {feature_definitions} does not exist.")
+
+    output_folder = Path(args.output_folder) if args.output_folder else Path('queries')
+    output_folder.mkdir(parents=True, exist_ok=True)
+
+    if (args.overlay_id == 0) or (args.overlay_id is None):
+        overlay_files = sorted(input_folder.glob('solution_*.mol2'))
+        pharm_files = sorted(input_folder.glob('pharmacophores/solution_pharm_*.mol2'))
+    else:
+        overlay_files = [input_folder / f'solution_{args.overlay_id:02}.mol2']
+        pharm_files = [input_folder / f'pharmacophores/solution_pharm_{args.overlay_id:02}.mol2']
+    feature_sets = []
+    for pharm_file, overlay_file in zip(pharm_files, overlay_files):
+        overlay_data = OverlayData(
+            input_folder=input_folder,
+            output_folder=output_folder,
+            pharm_file=pharm_file,
+            overlay_file=overlay_file
+        )
+        feature_sets.append(OverlayToPharmFeatures(overlay_data, projected=args.projected).features)
+
+    for i, feature_set in enumerate(feature_sets, 1):
+        if args.cluster:
+            feature_set = cluster_features(feature_set)
+
+        query = FeaturesToCrossMinerQuery(
+            pharm_feature_points=feature_set,
+            feature_definitions=feature_definitions,
+            output_file=output_folder / f'features_{i}.cm',
+        )
+        query.write_feature_file()
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/discovery/overlay_to_pharmacophore/overlay.py b/scripts/discovery/overlay_to_pharmacophore/overlay.py
new file mode 100644
index 0000000..1bc51bd
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/overlay.py
@@ -0,0 +1,105 @@
+from functools import cached_property
+
+import numpy as np
+
+from datastructures import OverlayData, PharmFeaturePoint
+
+
+class OverlayToPharmFeatures:
+    def __init__(self, overlay_data: OverlayData, projected: bool = True):
+        """
+        Args:
+            overlay_data: Pharmacophore and overlay data for a single ligand overlay.
+            projected: Whether to use projected acceptor features
+                (since acceptors are flexible and vary across overlays,
+                it may be better to simply use point acceptors).
+        """
+        self.overlay_data = overlay_data
+        self.projected = projected
+
+    @cached_property
+    def features(self) -> list[PharmFeaturePoint]:
+        return self._feature_centres_from_pharm()
+
+    def _read_atom_lines_from_pharm(self):
+        """
+        Generator that yields lines from the ATOM block of a pharmacophore file.
+        This block is where the pharmacophore feature points are defined.
+        """
+        reading_atoms = False
+        with open(self.overlay_data.pharm_file) as pharm_file:
+            for line in pharm_file:
+                line = line.strip()
+
+                if line == "@<TRIPOS>ATOM":
+                    reading_atoms = True
+                    continue
+                elif line.startswith("@<TRIPOS>") and reading_atoms:
+                    break  # End of ATOM block
+
+                if reading_atoms and line:
+                    yield line.split()
+
+    def _feature_centres_from_pharm(self) -> list[PharmFeaturePoint]:
+        """
+        Parse the pharmacophore mol2 file and extract feature points centres.
+
+        returns:  List of PharmFeaturePoints, each labelled with the feature name
+            (e.g. 'donor_projected', 'acceptor', 'hydrophobe'), and with associated virtual points if present.
+        """
+        feature_centres = []
+        last_feature_point = None
+
+        for parts in self._read_atom_lines_from_pharm():
+            _feature_no, feature_name, x, y, z = parts[:5]
+            # Extract first 2/3 chars which describe the feature type
+            prefix = ''.join(c for c in feature_name if c.isalpha())
+
+            if prefix in ('DON', 'ACC', 'HY'):
+                label = self._prefix_to_label(prefix)
+                last_feature_point = PharmFeaturePoint(
+                    (float(x), float(y), float(z)),
+                    label=label,
+                )
+                feature_centres.append(last_feature_point)
+            elif prefix == 'VP':
+                if last_feature_point is None:
+                    raise ValueError("Virtual point found before its parent feature point.")
+                # Acceptor feature with a virtual point is an acceptor projected feature.
+                if last_feature_point.label == "acceptor":
+                    # If preferring to use point acceptors, ignore the virtual point.
+                    if not self.projected:
+                        continue
+                    last_feature_point.label = "acceptor_projected"
+                vp_coords = np.array([x, y, z], dtype=float)
+                # If the last feature has no virtual points, assign the virtual point to it.
+                if last_feature_point.virtual_point is None:
+                    last_feature_point.virtual_point = vp_coords
+                # If the last feature does have virtual points, create a new feature.
+                else:
+                    new_feature = last_feature_point.copy()
+                    new_feature.virtual_point = vp_coords
+                    feature_centres.append(new_feature)
+            elif prefix == 'NL':
+                if last_feature_point is None:
+                    raise ValueError("Virtual point found before its parent feature point.")
+                last_feature_point.virtual_point = np.array([x, y, z], dtype=float)
+                last_feature_point.label = 'ring_planar_projected'
+            else:
+                raise ValueError(f"Unexpected feature type: {feature_name}")
+        return feature_centres
+
+    @staticmethod
+    def _prefix_to_label(prefix: str) -> str:
+        """
+        Convert a feature prefix to its corresponding label.
+        Acceptor may later be converted to acceptor_projected if it has a virtual point.
+        Args:
+            prefix: Prefix to convert to label
+        returns: Label for the feature type, e.g. 'acceptor', 'donor_projected', 'hydrophobe'
+        """
+        return {
+            'ACC': 'acceptor',
+            'DON': 'donor_projected',
+            'HY': 'hydrophobe',
+        }[prefix]
diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2 b/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2
new file mode 100644
index 0000000..71f5d16
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/tests/job1/pharmacophores/solution_pharm_01.mol2
@@ -0,0 +1,63 @@
+@<TRIPOS>MOLECULE
+pharmacophore
+    20    12     1     0     0
+SMALL
+NO_CHARGES
+****
+Generated from the CSD
+
+@<TRIPOS>ATOM
+     1 DON8.0_2    4.0507   8.2244  -1.3153   Ar        1 RES1   0.0000
+     2 VP1.0    3.7119   6.2366  -3.5091   F         1 RES1   0.0000
+     3 VP1.0    3.4830  11.1446  -0.9341   F         1 RES1   0.0000
+     4 VP1.0    5.0321   6.4941  -3.5384   F         1 RES1   0.0000
+     5 VP1.0    2.4551   6.6287  -3.1080   F         1 RES1   0.0000
+     6 ACC8.0_2    1.5469   6.6465  -0.4394   B         1 RES1   0.0000
+     7 VP1.0   -0.6255   8.6951  -0.6938   F         1 RES1   0.0000
+     8 VP1.0   -1.2681   7.3618  -1.1537   F         1 RES1   0.0000
+     9 VP1.0   -1.1566   5.8367  -1.4214   F         1 RES1   0.0000
+    10 ACC8.0_2    4.0507   8.2244  -1.3153   B         1 RES1   0.0000
+    11 VP1.0    3.7119   6.2366  -3.5091   F         1 RES1   0.0000
+    12 VP1.0    3.4830  11.1446  -0.9341   F         1 RES1   0.0000
+    13 VP1.0    5.0321   6.4941  -3.5384   F         1 RES1   0.0000
+    14 VP1.0    2.4551   6.6287  -3.1080   F         1 RES1   0.0000
+    15 HY4.0_2    1.6670   2.6310  -0.9036   C.3       1 RES1   0.0000
+    16 HY4.0_2    1.2974   1.8244  -0.9805   C.3       1 RES1   0.0000
+    17 HY4.0_2    2.0589   2.5883  -0.5021   C.3       1 RES1   0.0000
+    18 HY4.0_2    3.5562   7.3509   0.3815   C.3       1 RES1   0.0000
+    19 HY11.9_2    2.4272   9.3712   1.8815   C.2       1 RES1   0.0000
+    20 NL       3.2267   9.9618   1.7718   F         1 RES1   0.0000
+@<TRIPOS>BOND
+     1     1     2    1
+     2     1     3    1
+     3     1     4    1
+     4     1     5    1
+     5     6     7    1
+     6     6     8    1
+     7     6     9    1
+     8    10    11    1
+     9    10    12    1
+    10    10    13    1
+    11    10    14    1
+    12    19    20    1
+@<TRIPOS>SUBSTRUCTURE
+     1 RES1        1 GROUP             0 ****  ****    0  
+# ATOMS1 1 17 2 25
+# COORDS1 4.0479 4.0535 8.2010 8.2477 -1.5472 -1.0834
+# ATOMS6 1 14 2 13
+# COORDS6 1.4558 1.6380 6.6269 6.6661 -0.5341 -0.3446
+# ATOMS10 1 17 2 25
+# COORDS10 4.0479 4.0535 8.2010 8.2477 -1.5472 -1.0834
+# ATOMS15 1 5 1 6 1 7 1 8 1 9 1 10 1 11 2 3 2 4 2 5
+# ATOMS15 2 6 2 7 2 8 2 9
+# COORDS15 0.1104 3.2954 1.6196 4.0590 -2.3001 0.4906
+# ATOMS16 1 4 1 5 1 6 1 7 1 8 2 2 2 3 2 7 2 8 2 9
+# COORDS16 0.1104 2.4723 1.1327 2.2710 -2.3001 0.1771
+# ATOMS17 1 4 1 5 1 8 1 9 1 10 1 11 2 2 2 3 2 4 2 5
+# ATOMS17 2 6 2 7
+# COORDS17 0.7974 3.2954 1.1327 4.0590 -1.4998 0.4906
+# ATOMS18 1 13 1 15 1 16 1 18 2 12 2 14 2 15 2 20
+# COORDS18 2.6111 4.8913 6.3101 8.4335 -0.4078 1.3029
+# ATOMS19 1 18 1 19 1 20 1 21 1 22 1 23 2 20 2 21 2 22 2 23
+# ATOMS19 2 24
+# COORDS19 1.5791 3.2261 8.4084 10.5289 0.5379 3.2697
diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2 b/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2
new file mode 100644
index 0000000..3401d35
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/tests/job1/solution_01.mol2
@@ -0,0 +1,234 @@
+@<TRIPOS>MOLECULE
+ipratroponium
+    54    56     1     0     0
+SMALL
+USER_CHARGES
+****
+Generated from the CSD
+
+@<TRIPOS>ATOM
+     1 C1       0.4143  -0.6927   0.8235   C.3       1 RES1   0.0000
+     2 C2       1.2898  -0.2841  -0.3556   C.3       1 RES1   0.0000
+     3 C3       2.3978  -1.3166  -0.5276   C.3       1 RES1   0.0000
+     4 N4       1.8801   1.1327  -0.1432   N.4       1 RES1   1.0000
+     5 C5       0.7974   2.1776   0.0874   C.3       1 RES1   0.0000
+     6 C6       0.1175   2.2710  -1.2830   C.3       1 RES1   0.0000
+     7 C7       1.2245   1.8950  -2.3001   C.3       1 RES1   0.0000
+     8 C8       2.4723   1.6196  -1.4584   C.3       1 RES1   0.0000
+     9 C9       3.2954   2.8815  -1.1841   C.3       1 RES1   0.0000
+    10 C10      2.4141   3.9967  -0.6094   C.3       1 RES1   0.0000
+    11 C11      1.4899   3.4733   0.4906   C.3       1 RES1   0.0000
+    12 O12      3.2402   5.0698  -0.0540   O.3       1 RES1   0.0000
+    13 C13      2.7385   6.3101   0.0254   C.2       1 RES1   0.0000
+    14 O14      1.6380   6.6269  -0.3446   O.2       1 RES1   0.0000
+    15 C15      3.7487   7.2538   0.6507   C.3       1 RES1   0.0000
+    16 C16      4.7309   7.7403  -0.4078   C.3       1 RES1   0.0000
+    17 O17      4.0479   8.2477  -1.5472   O.3       1 RES1   0.0000
+    18 C18      3.0043   8.4084   1.3029   C.ar      1 RES1   0.0000
+    19 C19      2.2545   9.2919   0.5379   C.ar      1 RES1   0.0000
+    20 C20      1.5791  10.3453   1.1324   C.ar      1 RES1   0.0000
+    21 C21      1.6455  10.5289   2.4981   C.ar      1 RES1   0.0000
+    22 C22      2.3880   9.6590   3.2697   C.ar      1 RES1   0.0000
+    23 C23      3.0638   8.6053   2.6762   C.ar      1 RES1   0.0000
+    24 C24      2.9214   1.0709   0.9623   C.3       1 RES1   0.0000
+    25 H25     -0.4701  -0.0620   0.8842   H         1 RES1   0.0000
+    26 H26      0.9588  -0.6057   1.7611   H         1 RES1   0.0000
+    27 H27      0.0871  -1.7251   0.7191   H         1 RES1   0.0000
+    28 H28      0.6647  -0.3122  -1.2603   H         1 RES1   0.0000
+    29 H29      3.0775  -1.0241  -1.3253   H         1 RES1   0.0000
+    30 H30      1.9805  -2.2895  -0.7786   H         1 RES1   0.0000
+    31 H31      2.9749  -1.4245   0.3881   H         1 RES1   0.0000
+    32 H32      0.0947   1.8529   0.8533   H         1 RES1   0.0000
+    33 H33     -0.2428   3.2826  -1.4645   H         1 RES1   0.0000
+    34 H34     -0.7153   1.5722  -1.3464   H         1 RES1   0.0000
+    35 H35      1.4040   2.7184  -2.9896   H         1 RES1   0.0000
+    36 H36      0.9400   1.0065  -2.8618   H         1 RES1   0.0000
+    37 H37      3.0875   0.8468  -1.9170   H         1 RES1   0.0000
+    38 H38      3.7422   3.2249  -2.1159   H         1 RES1   0.0000
+    39 H39      4.0824   2.6452  -0.4695   H         1 RES1   0.0000
+    40 H40      1.8033   4.4069  -1.4121   H         1 RES1   0.0000
+    41 H41      2.0794   3.2920   1.3880   H         1 RES1   0.0000
+    42 H42      0.7317   4.2263   0.7003   H         1 RES1   0.0000
+    43 H43      4.3077   6.7110   1.4285   H         1 RES1   0.0000
+    44 H44      5.3796   6.9097  -0.7126   H         1 RES1   0.0000
+    45 H45      5.3695   8.5244   0.0174   H         1 RES1   0.0000
+    46 H46      4.7082   8.5543  -2.2167   H         1 RES1   0.0000
+    47 H47      2.1970   9.1558  -0.5370   H         1 RES1   0.0000
+    48 H48      0.9964  11.0276   0.5223   H         1 RES1   0.0000
+    49 H49      1.1161  11.3539   2.9633   H         1 RES1   0.0000
+    50 H50      2.4428   9.8005   4.3440   H         1 RES1   0.0000
+    51 H51      3.6457   7.9267   3.2909   H         1 RES1   0.0000
+    52 H52      2.4749   0.7263   1.8923   H         1 RES1   0.0000
+    53 H53      3.3699   2.0423   1.1516   H         1 RES1   0.0000
+    54 H54      3.7219   0.3833   0.6985   H         1 RES1   0.0000
+@<TRIPOS>BOND
+     1     1     2    1
+     2     2     3    1
+     3     2     4    1
+     4     4     5    1
+     5     5     6    1
+     6     6     7    1
+     7     7     8    1
+     8     4     8    1
+     9     8     9    1
+    10     9    10    1
+    11    10    11    1
+    12     5    11    1
+    13    10    12    1
+    14    12    13    1
+    15    13    14    2
+    16    13    15    1
+    17    15    16    1
+    18    16    17    1
+    19    15    18    1
+    20    18    19   ar
+    21    19    20   ar
+    22    20    21   ar
+    23    21    22   ar
+    24    22    23   ar
+    25    18    23   ar
+    26     4    24    1
+    27     1    25    1
+    28     1    26    1
+    29     1    27    1
+    30     2    28    1
+    31     3    29    1
+    32     3    30    1
+    33     3    31    1
+    34     5    32    1
+    35     6    33    1
+    36     6    34    1
+    37     7    35    1
+    38     7    36    1
+    39     8    37    1
+    40     9    38    1
+    41     9    39    1
+    42    10    40    1
+    43    11    41    1
+    44    11    42    1
+    45    15    43    1
+    46    16    44    1
+    47    16    45    1
+    48    17    46    1
+    49    19    47    1
+    50    20    48    1
+    51    21    49    1
+    52    22    50    1
+    53    23    51    1
+    54    24    52    1
+    55    24    53    1
+    56    24    54    1
+@<TRIPOS>SUBSTRUCTURE
+     1 RES1        1 GROUP             0 ****  ****    0  
+@<TRIPOS>MOLECULE
+tiotroponium
+    48    52     1     0     0
+SMALL
+USER_CHARGES
+****
+Generated from the CSD
+
+@<TRIPOS>ATOM
+     1 C1       1.5214  -0.1775  -0.0873   C.3       1 RES1   0.0000
+     2 N2       2.0030   1.2552  -0.0998   N.4       1 RES1   1.0000
+     3 C3       0.8928   2.2223   0.1771   C.3       1 RES1   0.0000
+     4 C4       1.5247   3.5898   0.4313   C.3       1 RES1   0.0000
+     5 C5       2.3168   4.0590  -0.8009   C.3       1 RES1   0.0000
+     6 C6       3.2057   2.9780  -1.4159   C.3       1 RES1   0.0000
+     7 C7       2.4140   1.6737  -1.4998   C.3       1 RES1   0.0000
+     8 C8       1.0622   1.8176  -2.1727   C.3       1 RES1   0.0000
+     9 C9       0.1104   2.1792  -1.1126   C.3       1 RES1   0.0000
+    10 O10      0.3365   3.0608  -2.2197   O.3       1 RES1   0.0000
+    11 O11      3.1602   5.1905  -0.4171   O.3       1 RES1   0.0000
+    12 C12      2.6111   6.4096  -0.3315   C.2       1 RES1   0.0000
+    13 O13      1.4558   6.6661  -0.5341   O.2       1 RES1   0.0000
+    14 C14      3.6619   7.4673   0.0699   C.3       1 RES1   0.0000
+    15 C15      4.8913   6.7838   0.6568   C.2       1 RES1   0.0000
+    16 C16      4.9473   5.8556   1.6681   C.2       1 RES1   0.0000
+    17 C17      6.2767   5.4539   1.9307   C.2       1 RES1   0.0000
+    18 C18      7.1812   6.0616   1.1413   C.2       1 RES1   0.0000
+    19 S19      6.4499   7.1454   0.0489   S.3       1 RES1   0.0000
+    20 C20      3.0632   8.4335   1.0857   C.2       1 RES1   0.0000
+    21 C21      3.2261   8.4489   2.4496   C.2       1 RES1   0.0000
+    22 C22      2.5043   9.5073   3.0467   C.2       1 RES1   0.0000
+    23 C23      1.8260  10.2557   2.1576   C.2       1 RES1   0.0000
+    24 S24      2.0400   9.7010   0.5610   S.3       1 RES1   0.0000
+    25 O25      4.0535   8.2010  -1.0834   O.3       1 RES1   0.0000
+    26 C26      3.1383   1.2852   0.8931   C.3       1 RES1   0.0000
+    27 H27      2.3197  -0.8530  -0.3875   H         1 RES1   0.0000
+    28 H28      0.6854  -0.3205  -0.7683   H         1 RES1   0.0000
+    29 H29      1.1955  -0.4635   0.9106   H         1 RES1   0.0000
+    30 H30      0.2859   1.9006   1.0222   H         1 RES1   0.0000
+    31 H31      2.1982   3.5196   1.2842   H         1 RES1   0.0000
+    32 H32      0.7389   4.3121   0.6478   H         1 RES1   0.0000
+    33 H33      1.6102   4.3964  -1.5576   H         1 RES1   0.0000
+    34 H34      4.0868   2.8312  -0.7929   H         1 RES1   0.0000
+    35 H35      3.5142   3.2832  -2.4147   H         1 RES1   0.0000
+    36 H36      3.0035   0.8936  -1.9793   H         1 RES1   0.0000
+    37 H37      0.7633   1.0360  -2.8696   H         1 RES1   0.0000
+    38 H38     -0.8366   1.6417  -1.0965   H         1 RES1   0.0000
+    39 H39      4.0794   5.4783   2.1988   H         1 RES1   0.0000
+    40 H40      6.5414   4.7281   2.6925   H         1 RES1   0.0000
+    41 H41      8.2501   5.8827   1.1942   H         1 RES1   0.0000
+    42 H42      3.8331   7.7374   2.9999   H         1 RES1   0.0000
+    43 H43      2.4978   9.6973   4.1150   H         1 RES1   0.0000
+    44 H44      1.2151  11.1113   2.4258   H         1 RES1   0.0000
+    45 H45      4.4434   7.5826  -1.7496   H         1 RES1   0.0000
+    46 H46      2.7864   1.0180   1.8873   H         1 RES1   0.0000
+    47 H47      3.5853   2.2755   0.9482   H         1 RES1   0.0000
+    48 H48      3.9160   0.5793   0.6091   H         1 RES1   0.0000
+@<TRIPOS>BOND
+     1     1     2    1
+     2     2     3    1
+     3     3     4    1
+     4     4     5    1
+     5     5     6    1
+     6     6     7    1
+     7     2     7    1
+     8     7     8    1
+     9     8     9    1
+    10     3     9    1
+    11     9    10    1
+    12     8    10    1
+    13     5    11    1
+    14    11    12    1
+    15    12    13    2
+    16    12    14    1
+    17    14    15    1
+    18    15    16    2
+    19    16    17    1
+    20    17    18    2
+    21    18    19    1
+    22    15    19    1
+    23    14    20    1
+    24    20    21    2
+    25    21    22    1
+    26    22    23    2
+    27    23    24    1
+    28    20    24    1
+    29    14    25    1
+    30     2    26    1
+    31     1    27    1
+    32     1    28    1
+    33     1    29    1
+    34     3    30    1
+    35     4    31    1
+    36     4    32    1
+    37     5    33    1
+    38     6    34    1
+    39     6    35    1
+    40     7    36    1
+    41     8    37    1
+    42     9    38    1
+    43    16    39    1
+    44    17    40    1
+    45    18    41    1
+    46    21    42    1
+    47    22    43    1
+    48    23    44    1
+    49    25    45    1
+    50    26    46    1
+    51    26    47    1
+    52    26    48    1
+@<TRIPOS>SUBSTRUCTURE
+     1 RES1        1 GROUP             0 ****  ****    0  
diff --git a/scripts/discovery/overlay_to_pharmacophore/tests/tests.py b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py
new file mode 100644
index 0000000..2406ea5
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/tests/tests.py
@@ -0,0 +1,105 @@
+from pathlib import Path
+
+import pytest
+
+from cluster import cluster_features
+from datastructures import OverlayData, PharmFeaturePoint
+from overlay import OverlayToPharmFeatures
+from write_query import FeaturesToCrossMinerQuery
+
+TESTS_DIR = Path(__file__).parent
+REPO_ROOT = TESTS_DIR.parent
+JOB_FOLDER = TESTS_DIR / "job1"
+
+
+@pytest.fixture(autouse=True)
+def run_from_repo_root(monkeypatch):
+    """Run every test from the repository root so ``feature_definitions/`` resolves."""
+    monkeypatch.chdir(REPO_ROOT)
+
+
+@pytest.fixture
+def overlay_data(tmp_path) -> OverlayData:
+    """Overlay/pharmacophore file locations for the bundled ``job1`` example."""
+    return OverlayData(
+        input_folder=JOB_FOLDER,
+        output_folder=tmp_path,
+        pharm_file=JOB_FOLDER / "pharmacophores/solution_pharm_01.mol2",
+        overlay_file=JOB_FOLDER / "solution_01.mol2",
+    )
+
+
+@pytest.fixture
+def features(overlay_data) -> list[PharmFeaturePoint]:
+    """Raw (unclustered) features extracted from the example pharmacophore."""
+    return OverlayToPharmFeatures(overlay_data).features
+
+
+@pytest.fixture
+def feature_definitions(features, tmp_path) -> Path:
+    """
+    Stub CrossMiner feature definition directory.
+
+    The real ``.cpf`` files ship with CrossMiner rather than this repo, so minimal
+    placeholders are generated for each label present in the example pharmacophore.
+    """
+    definitions = tmp_path / "feature_definitions" / "any"
+    definitions.mkdir(parents=True)
+    for label in {f.label for f in features}:
+        (definitions / f"features_{label}.cpf").write_text(f"SUBSTRUCTURE {label}\n")
+    return definitions.parent
+
+
+def test_features_extracted_from_pharmacophore(features):
+    """Parsing the pharmacophore mol2 yields labelled feature points."""
+    assert len(features) > 0, "No features extracted from pharmacophore."
+    assert all(feature.label is not None for feature in features)
+
+
+def test_projected_features_have_virtual_points(features):
+    """Every feature labelled as projected carries a virtual point."""
+    projected = [f for f in features if f.label.endswith("projected")]
+    assert projected, "No projected features extracted from pharmacophore."
+    assert all(f.virtual_point is not None for f in projected)
+
+
+def test_unprojected_acceptors_have_no_virtual_points(overlay_data):
+    """With ``projected=False`` acceptors are kept as plain points."""
+    features = OverlayToPharmFeatures(overlay_data, projected=False).features
+    acceptors = [f for f in features if f.label.startswith("acceptor")]
+    assert acceptors, "No acceptor features extracted from pharmacophore."
+    assert all(f.label == "acceptor" and f.virtual_point is None for f in acceptors)
+
+
+def test_clustering_reduces_features(features):
+    """Clustering merges nearby features without dropping any feature type."""
+    clustered = cluster_features(features)
+    assert 0 < len(clustered) <= len(features)
+    assert {f.label for f in clustered} == {f.label for f in features}
+
+
+def test_query_file_written(features, feature_definitions, tmp_path):
+    """A CrossMiner query file is written containing the clustered features."""
+    output_file = tmp_path / "query.cm"
+    clustered = cluster_features(features)
+
+    FeaturesToCrossMinerQuery(
+        clustered, feature_definitions=feature_definitions, output_file=output_file
+    ).write_feature_file()
+
+    assert output_file.exists()
+    contents = output_file.read_text()
+    assert "FEATURE_LIBRARY_START" in contents
+    assert "FEATURE_LIBRARY_END" in contents
+    assert contents.count("PHARMACOPHORE_FEATURE ") == len(clustered)
+
+
+def test_missing_feature_definitions_raises(features, tmp_path):
+    """An error is raised when the supplied feature definitions directory is absent."""
+    with pytest.raises(FileNotFoundError):
+        FeaturesToCrossMinerQuery(
+            features,
+            feature_definitions=tmp_path / "does_not_exist",
+            output_file=tmp_path / "query.cm",
+        )
+
diff --git a/scripts/discovery/overlay_to_pharmacophore/write_query.py b/scripts/discovery/overlay_to_pharmacophore/write_query.py
new file mode 100644
index 0000000..4d231b4
--- /dev/null
+++ b/scripts/discovery/overlay_to_pharmacophore/write_query.py
@@ -0,0 +1,110 @@
+from pathlib import Path
+
+from datastructures import PharmFeaturePoint
+
+
+class FeaturesToCrossMinerQuery:
+    """Main class for converting a set of features to a CrossMiner feature file."""
+
+    TEMPLATE_FEATURES_FROM_FILE = """
+    FEATURE_SUBSTRUCTURE_START
+
+        {}
+    FEATURE_SUBSTRUCTURE_END
+    """
+    def __init__(
+            self,
+            pharm_feature_points: list[PharmFeaturePoint],
+            feature_definitions: Path,
+            output_file: Path = Path('test.cm'),
+    ):
+        """
+        Args:
+            pharm_feature_points: List of PharmFeaturePoints to be used in the CrossMiner query.
+            feature_definitions: Directory containing the CrossMiner feature definition (``.cpf``) files.
+                These are not shipped with this repo; supply the location of your CrossMiner installation's
+                feature definitions (the directory containing the ``any``, ``protein`` and
+                ``small_molecule`` subdirectories).
+            output_file: Output file path where the features will be saved.
+        """
+        self.pharm_feature_points: list[PharmFeaturePoint] = pharm_feature_points
+        self.feature_definitions = Path(feature_definitions)
+        if not self.feature_definitions.is_dir():
+            raise FileNotFoundError(
+                f"Feature definitions directory {self.feature_definitions} does not exist."
+            )
+        self.output_file = output_file
+
+    def _feature_definition_file(self, feature: PharmFeaturePoint) -> Path:
+        """
+        Locate the ``.cpf`` definition file for a feature within the supplied feature definitions directory.
+
+        Args:
+            feature: PharmFeaturePoint whose label determines the definition file name.
+        returns: Path to the matching ``.cpf`` file.
+        """
+        if feature.label is None:
+            raise LookupError("Feature label must be set to locate its feature definition file.")
+
+        filename = f'features_{feature.label}.cpf'
+        matches = sorted(self.feature_definitions.glob(f'*/{filename}'))
+        if not matches:
+            # Also allow a flat directory of definition files.
+            flat = self.feature_definitions / filename
+            if flat.is_file():
+                return flat
+            raise FileNotFoundError(
+                f"No feature definition file '{filename}' found in {self.feature_definitions}."
+            )
+        return matches[0]
+
+
+    @staticmethod
+    def _create_feature_vector_text(feature: PharmFeaturePoint) -> str:
+        """
+        Create a formatted string for a feature vector based on the PharmFeaturePoint object.
+        Args:
+            feature: PharmFeaturePoint containing coordinates, label, and virtual point if applicable.
+        """
+
+        result = (
+            f"PHARMACOPHORE_FEATURE {feature.label}\n"
+            f"PHARMACOPHORE_SPHERE {' '.join(str(i) for i in feature.coordinates)} {feature.weight_parent}\n"
+        )
+        if feature.virtual_point is not None:
+            result += f"PHARMACOPHORE_SPHERE {' '.join(str(i) for i in feature.virtual_point)} {feature.weight_vp}\n"
+
+        result += (
+            f"PHARMACOPHORE_FEATURE_SMALL_MOLECULE\n"
+            f"PHARMACOPHORE_FEATURE_DESCRIPTION {feature.label}\n\n"
+        )
+
+        return result
+
+    def write_feature_file(self):
+        """
+        Create a feature file with proper whitespace formatting.
+        """
+        output_data = 'FEATURE_LIBRARY_START\n'
+
+        # Collect unique feature types and their corresponding files
+        feature_files = dict()
+        for feature in self.pharm_feature_points:
+            if feature.label in feature_files:
+                continue
+            feature_file = self._feature_definition_file(feature)
+            feature_files[feature.label] = feature_file
+
+            # Append the feature data to the output stream
+            with open(feature_file) as ff:
+                output_data += self.TEMPLATE_FEATURES_FROM_FILE.format(
+                    "        ".join(ff.readlines())
+                )
+
+        output_data += '\nFEATURE_LIBRARY_END\n\n'
+
+        for feature in self.pharm_feature_points:
+            output_data += self._create_feature_vector_text(feature)
+
+        with open(self.output_file, 'w') as output_file:
+            output_file.write(output_data)