How to compute protein–ligand interactions#

Goal#

Detect hydrogen bonds, π–π stacking, cation–π, and σ-hole interactions between a protein receptor and a small-molecule ligand.

Minimal example#

import numpy as np
from rdkit.Chem import ChemicalFeatures
from rdkit import RDConfig
import os
from moleculekit.molecule import Molecule
from moleculekit.interactions.interactions import hbonds_calculate

mol = Molecule("3PTB")
mol.templateResidueFromSmiles(
    mol.resname == "BEN",
    smiles="NC(=N)c1ccccc1",
    addHs=True,
)

# Convert the ligand to RDKit, then ask RDKit which atoms are H-bond
# donors and acceptors.
lig = mol.copy(sel="resname BEN")
rdlig = lig.toRDKitMol(sanitize=True)

fdef = os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")
factory = ChemicalFeatures.BuildFeatureFactory(fdef)
features = factory.GetFeaturesForMol(rdlig)

start_idx = int(mol.atomselect("resname BEN", indexes=True)[0])
acceptors = [idx + start_idx for f in features if f.GetFamily() == "Acceptor"
             for idx in f.GetAtomIds()]
donor_pairs = [(idx + start_idx, idx + start_idx)
               for f in features if f.GetFamily() == "Donor"
               for idx in f.GetAtomIds()]

# For non-periodic structures hbonds_calculate needs a zero box.
mol.box = np.zeros((3, mol.numFrames), dtype=np.float32)
hbonds = hbonds_calculate(mol, np.array(donor_pairs), np.array(acceptors), sel2="protein")
print(f"H-bonds in frame 0: {len(hbonds[0])}")

Parameters that matter#

Function/parameter

What it does

toRDKitMol(sanitize=True)

Convert the ligand selection to an RDKit Mol for feature detection.

RDConfig.RDDataDir / "BaseFeatures.fdef"

RDKit’s standard donor/acceptor feature definition file.

start_idx

Offset of the first ligand atom in the parent Molecule; required when the ligand does not start at atom 0.

hbonds_calculate(mol, donors, acceptors, sel2=...)

Returns a list (one entry per frame) of H-bond arrays.

Common variations#

# Visualise H-bonds in VMD (requires VMD on PATH)
from moleculekit.interactions.interactions import view_hbonds

view_hbonds(mol, hbonds)

# Full interaction set (H-bond + π–π + cation–π + σ-hole) via the high-level API
from moleculekit.interactions.interactions import (
    pipi_calculate,
    cationpi_calculate,
    sigmahole_calculate,
    get_protein_rings,
)

rings = get_protein_rings(mol)

Gotchas#

  • Reasonable hydrogens must be present — run systemPrepare() before this step if the input lacks explicit H atoms.

  • Feature detection uses RDKit through toRDKitMol(); make sure RDKit is installed.

  • start_idx is critical when the ligand atoms do not start at index 0 in the parent molecule.

  • For non-periodic structures set mol.box to a zero array before calling hbonds_calculate.

See also#