Version Compatibility

Reference examples tested with: RDKit 2024.09+, pharmer / pharmit (web service), PharmIT 1.1+, plip 2.4+ (interaction analysis).

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show rdkit then help(rdkit.Chem.Pharm3D) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Pharmacophore Modeling

Build 3D pharmacophore queries that capture the essential interaction features of a ligand-target binding event. A pharmacophore is the spatial arrangement of pharmacophore features (donor, acceptor, hydrophobe, aromatic, charged) sufficient for activity, abstracted from any specific chemotype. Used for scaffold-hopping (find compounds with different scaffold but matching pharmacophore), virtual screening (faster than docking), and cross-target SAR transfer. Modern best practice: derive pharmacophore from co-crystal structure if available (receptor-based; apo2ph4 workflow of Heider et al 2022/2023 J Chem Inf Model 63:147-158) or align actives if no crystal (ligand-based). Diffusion-based generation (PharmacoForge, Flynn et al 2025 Front Bioinform) lets pharmacophore drive de novo design.

For 2D scaffold-based searches, see chemoinformatics/scaffold-analysis. For 3D shape similarity, see chemoinformatics/shape-similarity. For protein-ligand interaction analysis, see chemoinformatics/virtual-screening.

Pharmacophore Feature Types

Feature	RDKit code	Definition	Geometric tolerance
H-bond donor	D	-OH, -NH	1.0-1.5 Å
H-bond acceptor	A	sp2 O / N (lone pair)	1.0-1.5 Å
Hydrophobe	H	sp3 C / aromatic ring centroid	1.5-2.0 Å
Aromatic ring	R	Aromatic ring centroid + normal	1.0-1.5 Å
Positive ionizable	P	-NH3+, -NR3+	1.0-1.5 Å
Negative ionizable	N	-COO-, -SO3-	1.0-1.5 Å
Halogen	X	Cl, Br, I (halogen bond donor)	1.0-1.5 Å
Metal coordination	M	sp/sp2 N/O near metal	0.5-1.0 Å

Tolerances are pharmacophore-feature distance windows in the search. Tighter tolerances = fewer hits but more specific.

Method Taxonomy

Method	Origin	Use case	Fails when
Ligand-based (LBP)	Catalyst, MOE, RDKit Pharm3D	Multiple actives, no crystal	<3 actives; flexible actives
Receptor-based (RBP)	apo2ph4, LigandScout	Co-crystal available	Apo structure (use AlphaFold3 or Boltz)
Common pharmacophore	Pharm3D EmbedPharmacophore	Consensus from active set	Diverse actives confound alignment
Diffusion-based (PharmacoForge)	Flynn et al 2025 Front Bioinform	De novo generation with pharmacophore prior	Pretrained model required
Active learning pharmacophore	Catalyst variant	Iterative refinement	Custom; not standard

Decision Tree by Scenario

Scenario	Method	Tools
Co-crystal structure available	Receptor-based	apo2ph4 + Pharmer/Pharmit
Multiple active compounds, no crystal	Ligand-based common pharmacophore	RDKit Pharm3D EmbedPharmacophore
Single active compound	Single-conformer pharmacophore	RDKit Pharm3D from bioactive conformer
Scaffold hopping prospective	Receptor-based + shape filter	apo2ph4 + ROCS
Cross-target SAR transfer	Common pharmacophore across targets	Manual + LigandScout
De novo design with pharmacophore	PharmacoForge	Diffusion-based generation
Library pre-filtering	Pharmacophore screen	Pharmit search

Ligand-Based Pharmacophore (RDKit Pharm3D)

Goal: Extract a common pharmacophore from a set of bioactive compounds.

Approach: Align actives by maximum-common-substructure or shape; extract conserved features; derive a pharmacophore signature.

from rdkit import Chem
from rdkit.Chem import AllChem, ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
from rdkit.RDPaths import RDDataDir
import os

fdef_file = os.path.join(RDDataDir, 'BaseFeatures.fdef')
factory = ChemicalFeatures.BuildFeatureFactory(fdef_file)

active_smiles = ['CC(C)c1ccc(C(=O)NCc2ccccn2)cc1', 'CCC(C)c1ccc(C(=O)NCc2ccccn2)cc1']
active_mols = [Chem.AddHs(Chem.MolFromSmiles(s)) for s in active_smiles]
for m in active_mols:
    AllChem.EmbedMolecule(m, AllChem.ETKDGv3())
    AllChem.MMFFOptimizeMolecule(m)

# Extract pharmacophore features (family/type/atom-ids/3D-pos) per molecule
feature_lists = [factory.GetFeaturesForMol(m) for m in active_mols]

# Build a target pharmacophore from one active's features; the bounds matrix
# encodes per-feature-pair distance ranges across all conformer/active variability.
# Below uses 2 features (Aromatic + Donor) with 2.5-3.5 A and 5.0-6.0 A bands.
feature_types = ['Aromatic', 'Donor']
# bounds is symmetric 2x2 distance-range matrix; off-diagonal = (lo, hi)
import numpy as np
bounds_matrix = np.array([[0.0, 5.0], [3.5, 0.0]])  # upper / lower bounds
pharmacophore = Pharmacophore.Pharmacophore(feature_types)
pharmacophore.setLowerBound(0, 1, 3.5)
pharmacophore.setUpperBound(0, 1, 5.0)

# To search a target molecule for matches, use Pharm3D.EmbedLib.EmbedPharmacophore
# (see chemoinformatics/conformer-generation for 3D embedding fundamentals)

BaseFeatures.fdef (RDKit-shipped) defines feature SMARTS. For drug-like pharmacophores, this is the standard starting point.

Receptor-Based Pharmacophore (apo2ph4 workflow)

Goal: Derive a pharmacophore from a protein binding-pocket structure without requiring a bound ligand.

Approach: Identify hot-spots (donor / acceptor / hydrophobe regions) from protein geometry; assemble into a pharmacophore. The apo2ph4 workflow of Heider J, Kilian J, Garifulina A, Hering S, Langer T, Seidel T (2022/2023 J Chem Inf Model 63:147-158) describes the conceptual pipeline; the example below is illustrative — verify the exact CLI invocation against the published code/release before running.

# Conceptual apo2ph4-style workflow; flags shown are illustrative.
apo2ph4 -pdb receptor.pdb -site_residues 'A:100,A:101,A:104,A:108' \
        -output pharmacophore.ph4

apo2ph4 outputs a .ph4-style pharmacophore compatible with Phase, MOE, and Pharmer.

When a co-crystal ligand is available, derive pharmacophore directly from the ligand binding pose: each ligand feature in contact with a complementary protein residue is part of the pharmacophore.

from plip.basic import config
from plip.structure.preparation import PDBComplex

mol_complex = PDBComplex()
mol_complex.load_pdb('complex.pdb')
mol_complex.analyze()

for site in mol_complex.interaction_sets.values():
    for interaction in site.all_itypes:
        # H-bond donor / acceptor / pi-stacking / hydrophobic / salt-bridge
        feature_type = interaction.type
        feature_atom_coords = interaction.ligatom.coords

PLIP outputs interaction types per ligand atom; each interaction maps to a pharmacophore feature.

Pharmacophore Search (Pharmit / Pharmer)

For library screening, Pharmer/Pharmit are the standard tools.

# Pharmer command line (offline alternative)
pharmer search -q pharmacophore.ph4 -dbdir zinc_db -out hits.sdf

Or use Pharmit web service (https://pharmit.csb.pitt.edu) for browser-based or REST search.

import requests

url = 'https://pharmit.csb.pitt.edu/...'  # specific endpoint
response = requests.post(url, json={'pharmacophore': pharmacophore_dict})
hits = response.json()

Pharmacophore screen is orders of magnitude faster than docking; typical 100M-compound search in minutes.

Pharmacophore Quality Validation

Evaluate a pharmacophore by:

1. Retrospective enrichment: AUC-like metric on actives vs decoys (DUD-E, COCONUT)
2. Geometric tightness: feature distance variance across actives
3. Selectivity: false positives in inactive set should be low
4. Specific consistency: pharmacophore matches each active's bioactive conformer

def pharmacophore_enrichment(query_pharmacophore, actives, inactives):
    n_active_match = sum(matches_pharmacophore(a, query_pharmacophore) for a in actives)
    n_inactive_match = sum(matches_pharmacophore(d, query_pharmacophore) for d in inactives)
    enrichment = (n_active_match / len(actives)) / (n_inactive_match / len(inactives))
    return enrichment

A good pharmacophore: enrichment >= 5x relative to random.

Diffusion-Based Pharmacophore Design (PharmacoForge)

PharmacoForge (Flynn, Shah, Dunn et al 2025 Front Bioinform) generates molecules conditioned on a pharmacophore query using a diffusion model:

# Pseudo-code; depends on PharmacoForge installation
# from pharmacoforge import PharmacophoreDiffusionGenerator
# gen = PharmacophoreDiffusionGenerator()
# generated = gen.generate(pharmacophore_query, n_molecules=100)

Trade-off: PharmacoForge produces de novo molecules that satisfy pharmacophore; lower drug-likeness than REINVENT but better novelty.

Pharmacophore vs Shape vs 2D Fingerprint

Method	Captures	Best for
ECFP4 Tanimoto	Local atom environments	Lead optimization (same series)
FCFP4 Tanimoto	Pharmacophore-equivalent atoms	Loose similarity in series
Shape similarity (ROCS)	3D shape volume	Scaffold hopping by shape
Pharmacophore	Discrete features in space	Scaffold hopping with feature specificity
Combined (Tanimoto + shape)	Multi-objective	Production VS

Pharmacophore is more interpretable than shape: a hit explains why it matched (donor at position X, hydrophobe at position Y).

Per-Tool Failure Modes

Ligand-based -- diverse actives confound

Trigger: Active set spans multiple scaffolds with different bound conformations.

Mechanism: No common pharmacophore exists; algorithm forces non-consensus features.

Symptom: Pharmacophore matches no actives in retrospective.

Fix: Cluster actives by scaffold first; derive per-cluster pharmacophore.

Receptor-based -- apo structure

Trigger: Protein in apo form (no bound ligand).

Mechanism: Side-chain rotamers differ between apo and holo; "binding site" geometry is wrong.

Symptom: Pharmacophore inferred from apo doesn't match holo experimental data.

Fix: Use AlphaFold3 / Boltz-1 to predict holo conformation; derive pharmacophore from predicted holo.

Pharmacophore -- single conformer bias

Trigger: Active aligned to its first generated conformer, not bioactive conformer.

Mechanism: Crystal structure not available; generated conformer may not be the bound one.

Symptom: Pharmacophore inconsistent across runs (different starting conformer chosen).

Fix: Use conformer ensemble; align all to common scaffold; choose conformer most consistent with other actives.

Tolerance too tight

Trigger: Default geometric tolerance < 0.5 Å.

Mechanism: Real bioactive conformers have flexibility; rigid pharmacophore filters most molecules out.

Symptom: Search returns zero hits.

Fix: Use tolerance 1.0-1.5 Å for drug-like; up to 2 Å for flexible peptide-like.

Pharmacophore search misses bioisostere

Trigger: Bioisostere replacement (e.g., -COOH replaced by tetrazole).

Mechanism: Tetrazole functions as acid bioisostere but RDKit features may not classify identically.

Symptom: Known bioisosteric active not found.

Fix: Use ChemAxon-style bioisosteric feature equivalence; or pharmacophore feature class expansion (acid generic vs -COOH specific).

PLIP -- water-mediated interaction missed

Trigger: Bridging water between ligand donor and protein acceptor.

Mechanism: PLIP doesn't include explicit water.

Symptom: Pharmacophore missing critical H-bond feature.

Fix: Add water-mediated interaction manually based on crystal water positions; or use PoseView for full interaction view.

Reconciliation: Ligand-Based vs Receptor-Based

Aspect	Ligand-based	Receptor-based
Data needed	Multiple actives	Co-crystal or predicted holo
Bias	Toward known chemotype	None
Hit set	Similar to known actives	More diverse
Discovery potential	Limited	Higher
Pharmacophore confidence	Higher (validated against multiple actives)	Lower (single ligand)

For prospective scaffold-hopping, receptor-based is preferred. For ranking analogs, ligand-based is sufficient.

Common Errors

Symptom	Cause	Fix
Pharm3D.EmbedPharmacophore fails	Bounds matrix infeasible	Loosen tolerances; increase max_attempts
Pharmacophore matches everything	Too few features	Add features; tighten tolerances
Pharmacophore matches nothing	Too many features or tight bounds	Reduce feature count; loosen tolerances
BaseFeatures.fdef not found	RDKit installation issue	Check from rdkit.RDPaths import RDDataDir
Pharmacophore-conformer mismatch	Wrong conformer used	Use bioactive conformer from crystal
Pharmit search timeout	Library too large	Pre-filter by 2D fingerprint Tanimoto
apo2ph4 produces empty	No druggable hot-spots	Lower thresholds; use shape-only filter

References

• Wolber & Langer, J. Chem. Inf. Model. 45:160 (2005) -- LigandScout pharmacophore.
• Heider J, Kilian J, Garifulina A, Hering S, Langer T, Seidel T 2022/2023 J Chem Inf Model 63(1):147-158 -- apo2ph4 receptor-based pharmacophore workflow (DOI 10.1021/acs.jcim.2c00814).
• Flynn JR, Shah RH, Dunn SD et al 2025 Front Bioinform -- PharmacoForge diffusion pharmacophore.
• RDKit Chem.Pharm3D framework -- documentation at rdkit.org (the unconfirmed "Stiefl 2021" citation has been removed pending verification).
• Adasme et al., Nucleic Acids Res. 49:W530 (2021) -- PLIP interaction profiler.
• Koes et al., Nucleic Acids Res. 44:W436 (2016) -- Pharmit interactive search.

Related Skills

• chemoinformatics/molecular-io - Parse molecules
• chemoinformatics/conformer-generation - Generate 3D for pharmacophore
• chemoinformatics/shape-similarity - 3D shape adjacent to pharmacophore
• chemoinformatics/virtual-screening - Pharmacophore as docking pre-filter
• chemoinformatics/scaffold-analysis - 2D scaffold-hopping context
• chemoinformatics/generative-design - PharmacoForge for de novo
• structural-biology/structure-io - PDB handling