Tutorial 07 — Ligands and parameter files#

Load ligand chemistry, inspect a .tmol parameter bundle, score a protein–ligand complex, and refine its pocket. Assumes scoring, packing, and minimization.

Setup#

In Colab, select T4 GPU, then Run all. For local execution, follow the installation guide. Setup installs TMol and downloads the fixtures on first use.

The tutorial uses checked-in structures and parameter files. Colab downloads these fixtures during setup; the calculations do not query a live structure database. It keeps two chemistry paths separate: a pinned ADA/LG1 protein–ligand complex for scoring and refinement, and an independent AMPC MOL2 preparation example for inspecting generated parameters.

[1]:

try:
    import google.colab  # noqa: F401
except ImportError:
    IN_COLAB = False
else:
    IN_COLAB = True

if IN_COLAB:
    from urllib.request import urlopen

    exec(
        urlopen(
            "https://raw.githubusercontent.com/uw-ipd/tmol/"
            "master/docs/tutorial/colab_setup.py"
        ).read(),
        globals(),
    )
    setup_colab(
        [
            "tmol/tests/data/protein_ligand_test/ada.xtal-lig.mmff94.tmol",
            "tmol/tests/data/protein_ligand_test/ada.tmol.nomin.cif",
            "tmol/tests/data/ligand_test/ligand_ground_truth/mol2/ampc_1.mol2",
        ]
    )

[2]:

from contextlib import redirect_stderr, redirect_stdout
from io import StringIO
from pathlib import Path
import tempfile
import warnings

import attrs
import biotite.structure as struc
import biotite.structure.io
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import torch
import yaml
from IPython.display import display

import tmol
from tmol.database import ParameterDatabase
from tmol.io import build_context_from_biotite, pose_stack_from_biotite
from tmol.ligand import (
    inject_params_file,
    load_params_file,
    nonstandard_residue_info_from_mol2,
    prepare_single_ligand,
    write_params_file,
)
from tmol.ops import (
    build_coord_mask_for_mask_and_interacting_atoms,
    calculate_block_pair_ddg,
    compute_block_adjacency,
    compute_block_centroids_and_furthest_dist,
)
from tmol.score import beta2016_score_function

SEED = 20260807
np.random.seed(SEED)
torch.manual_seed(SEED)
warnings.filterwarnings(
    "ignore", message=r"Sparse invariant checks are implicitly disabled.*"
)
device = (
    torch.device("cuda", torch.cuda.current_device())
    if torch.cuda.is_available()
    else torch.device("cpu")
)
LIGAND_RES_NAME = "LG1"


def show_table(frame):
    try:
        from itables import show
    except ImportError:
        return display(frame)
    return show(frame)


def ligand_block_mask(pose_stack):
    mask = torch.zeros_like(pose_stack.block_type_ind, dtype=torch.bool)
    for pose_i in range(pose_stack.n_poses):
        for block_i in range(pose_stack.max_n_blocks):
            type_i = int(pose_stack.block_type_ind[pose_i, block_i])
            if type_i < 0:
                continue
            block_type = pose_stack.packed_block_types.active_block_types[type_i]
            mask[pose_i, block_i] = block_type.name3 == LIGAND_RES_NAME
    if not bool(mask.any()):
        raise RuntimeError(f"No {LIGAND_RES_NAME} ligand block was built")
    return mask

Environment variable CCD_MIRROR_PATH not set. Will not be able to use function requiring this variable. To set it you may:
  (1) add the line 'export VAR_NAME=path/to/variable' to your .bashrc or .zshrc file
  (2) set it in your current shell with 'export VAR_NAME=path/to/variable'
  (3) write it to a .env file in the root of the atomworks.io repository
Environment variable PDB_MIRROR_PATH not set. Will not be able to use function requiring this variable. To set it you may:
  (1) add the line 'export VAR_NAME=path/to/variable' to your .bashrc or .zshrc file
  (2) set it in your current shell with 'export VAR_NAME=path/to/variable'
  (3) write it to a .env file in the root of the atomworks.io repository

The `.tmol` parameter bundle#

Prefer mmCIF/CIF for structure input, especially for protein–ligand systems: CIF chemical-component and bond tables can preserve connectivity and bond order that PDB coordinate records do not reliably encode. Retain MOL2 or prepared .tmol chemistry when it is the authoritative ligand source; a coordinate-only conversion cannot recover missing bond orders.

TMol’s portable .tmol file is versioned YAML with three top-level payloads:

chemical: residue atoms, bonds, internal coordinates, torsions, and properties;
elec: per-atom partial charges; and
cartbonded: residue-specific bonded parameters.

The same schemas are split across TMol’s canonical database files at tmol/database/default/chemical/chemical.yaml, tmol/database/default/scoring/elec.yaml, and tmol/database/default/scoring/cartbonded.yaml. A portable .tmol file bundles ligand additions to those three domains; chemical.yaml alone is not a complete scoring parameter set.

[3]:

repo_root = Path.cwd()
if not (
    repo_root / "tmol/tests/data/protein_ligand_test/ada.tmol.nomin.cif"
).exists():
    repo_root = Path(tmol.__file__).resolve().parents[1]
data_dir = repo_root / "tmol" / "tests" / "data" / "protein_ligand_test"
reference_tmol = data_dir / "ada.xtal-lig.mmff94.tmol"
complex_cif = data_dir / "ada.tmol.nomin.cif"
ligand_mol2 = (
    repo_root
    / "tmol"
    / "tests"
    / "data"
    / "ligand_test"
    / "ligand_ground_truth"
    / "mol2"
    / "ampc_1.mol2"
)

# This is the pinned ADA/LG1 preparation used by the complex below. The later
# AMPC MOL2 preparation is an independent pipeline demonstration.
preparation = load_params_file(reference_tmol)[0]
work_dir = Path(tempfile.mkdtemp(prefix="tmol-ligand-tutorial-"))
tmol_path = work_dir / "LG1.tmol"
write_params_file(preparation, tmol_path)

print("Wrote:", tmol_path)

Wrote: /tmp/tmol-ligand-tutorial-x0evhkzd/LG1.tmol

Verify the parameter round trip#

The reloaded .tmol bundle must reproduce atom names, types, partial charges, and all three payloads. TMol supports this parameter format only.

[4]:

tmol_document = yaml.safe_load(tmol_path.read_text())

print("TMol .tmol excerpt")
print("\n".join(tmol_path.read_text().splitlines()[:18]))
print("\nTop-level payloads:", sorted(tmol_document))
print("Format version:", tmol_document["version"])

reloaded = load_params_file(tmol_path)[0]
ada_atom_frame = pd.DataFrame(
    [
        {
            "atom": atom.name,
            "TMol_atom_type": atom.atom_type,
            "pinned_partial_charge": preparation.partial_charges[atom.name],
            "reloaded_partial_charge": reloaded.partial_charges[atom.name],
            "provenance": "pinned ada.xtal-lig.mmff94.tmol",
        }
        for atom in preparation.residue_type.atoms
    ]
)
show_table(ada_atom_frame)

assert reloaded.residue_type.atoms == preparation.residue_type.atoms
assert reloaded.residue_type.bonds == preparation.residue_type.bonds
assert reloaded.partial_charges == preparation.partial_charges
print("\nRound trip preserved atoms, bonds, and charges.")

TMol .tmol excerpt
version: '5.2'
chemical:
  residues:
  - name: LG1
    base_name: LG1
    name3: LG1
    io_equiv_class: LG1
    atoms:
    - {name: C1, atom_type: CR}
    - {name: N1, atom_type: Nad3}
    - {name: C2, atom_type: CR}
    - {name: N2, atom_type: Nim}
    - {name: C3, atom_type: CR}
    - {name: C4, atom_type: CDp}
    - {name: O1, atom_type: Oad}
    - {name: N3, atom_type: Nad}
    - {name: C5, atom_type: CS1}
    - {name: C6, atom_type: CSp}

Top-level payloads: ['cartbonded', 'chemical', 'elec', 'version']
Format version: 5.2

^ⓘatom	TMol_atom_type	pinned_partial_charge	reloaded_partial_charge	provenance
C1	CR	-0.302	-0.302	pinned ada.xtal-lig.mmff94.tmol
N1	Nad3	0.048	0.048	pinned ada.xtal-lig.mmff94.tmol
C2	CR	0.036	0.036	pinned ada.xtal-lig.mmff94.tmol
N2	Nim	-0.565	-0.565	pinned ada.xtal-lig.mmff94.tmol
C3	CR	0.141	0.141	pinned ada.xtal-lig.mmff94.tmol
C4	CDp	0.716	0.716	pinned ada.xtal-lig.mmff94.tmol
O1	Oad	-0.570	-0.570	pinned ada.xtal-lig.mmff94.tmol
N3	Nad	-0.800	-0.800	pinned ada.xtal-lig.mmff94.tmol
C5	CS1	0.256	0.256	pinned ada.xtal-lig.mmff94.tmol
C6	CSp	0.280	0.280	pinned ada.xtal-lig.mmff94.tmol
(26 more rows not shown)


Round trip preserved atoms, bonds, and charges.

Prepare AMPC from MOL2#

Read the checked-in AMPC MOL2, containing explicit hydrogens, Tripos bond types, formal-charge annotations, and MMFF94 partial charges. res_name="AMP" names the new residue. This separate example neither regenerates nor scores the unrelated ADA/LG1 fixture. It needs no network request or stochastic conformer generation.

AMPC’s formal charge is −1; its partial charges should sum to approximately −1 at MOL2 precision. properties.protonation.protonation_state is packing metadata, not a formal-charge calculation. A "neutral" value must not override the source chemistry.

Coordinates alone do not establish reliable bond orders or charges. Preserve authoritative chemistry throughout preparation.

[5]:

mol2_text = ligand_mol2.read_text()
assert "MMFF94_CHARGES" in mol2_text.splitlines()[:10]
assert "charge -1" in mol2_text

mol2_info = nonstandard_residue_info_from_mol2(ligand_mol2, res_name="AMP")
if mol2_info.partial_charges is None:
    raise RuntimeError("AMPC MOL2 did not expose authoritative partial charges")

# Preserve the authoritative MOL2 records before preparation. Hydrogen names may
# be regenerated, while charge transfer remains index-based.
source_atom_names = tuple(mol2_info.atom_names)
source_elements = tuple(mol2_info.elements)
source_subtypes = tuple(mol2_info.atom_array.tmol_source_subtype)
source_partial_charges = dict(mol2_info.partial_charges)

generated_preparation = prepare_single_ligand(mol2_info)
generated_tmol_path = work_dir / "AMPC.generated.tmol"
write_params_file(generated_preparation, generated_tmol_path)

def mol2_unity_atom_formal_charges(text, atom_names):
    """Read explicit per-atom formal charges from a TRIPOS MOL2 block."""
    lines = [line.strip() for line in text.splitlines()]
    try:
        line_index = lines.index("@<TRIPOS>UNITY_ATOM_ATTR") + 1
    except ValueError as error:
        raise RuntimeError("MOL2 has no UNITY_ATOM_ATTR formal charges") from error
    charges = {name: 0 for name in atom_names}
    while line_index < len(lines) and not lines[line_index].startswith("@<TRIPOS>"):
        if not lines[line_index]:
            line_index += 1
            continue
        atom_id, attribute_count = map(
            int, lines[line_index].split()[:2]
        )
        line_index += 1
        for _ in range(attribute_count):
            attribute, value = lines[line_index].split()[:2]
            line_index += 1
            if attribute == "charge":
                charges[atom_names[atom_id - 1]] = int(value)
    return charges


source_formal_charges = mol2_unity_atom_formal_charges(
    mol2_text, source_atom_names
)
source_formal_charge = sum(source_formal_charges.values())
source_partial_charge_sum = float(sum(source_partial_charges.values()))
generated_partial_charge_sum = float(
    sum(generated_preparation.partial_charges.values())
)
generic_protonation_metadata = (
    generated_preparation.residue_type.properties.protonation.protonation_state
)
assert source_formal_charge == -1
assert np.isclose(source_partial_charge_sum, -1.0, atol=5e-4)
assert np.isclose(
    generated_partial_charge_sum, source_partial_charge_sum, atol=1e-8
)

generated_atoms = tuple(generated_preparation.residue_type.atoms)
assert len(generated_atoms) == len(source_atom_names)
source_charge_vector = [source_partial_charges[name] for name in source_atom_names]
generated_charge_vector = [
    generated_preparation.partial_charges[atom.name] for atom in generated_atoms
]
assert np.allclose(generated_charge_vector, source_charge_vector, atol=1e-8)
ampc_atom_frame = pd.DataFrame(
    [
        {
            "source_atom": source_name,
            "generated_atom": generated_atom.name,
            "element": element,
            "source_Tripos_subtype": subtype,
            "source_formal_charge": source_formal_charges[source_name],
            "source_MMFF94_partial_charge": source_partial_charges[source_name],
            "generated_TMol_atom_type": generated_atom.atom_type,
            "generated_partial_charge": generated_preparation.partial_charges[
                generated_atom.name
            ],
        }
        for source_name, element, subtype, generated_atom in zip(
            source_atom_names,
            source_elements,
            source_subtypes,
            generated_atoms,
        )
    ]
)
show_table(ampc_atom_frame)

preparation_frame = pd.DataFrame(
    [
        {
            "pipeline_role": "independent AMPC MOL2 preparation demo",
            "residue_name": generated_preparation.residue_type.name3,
            "atoms": len(generated_preparation.residue_type.atoms),
            "bonds": len(generated_preparation.residue_type.bonds),
            "torsions": len(generated_preparation.residue_type.torsions),
            "molecular_formal_charge_from_source_graph": source_formal_charge,
            "formal_charge_provenance": (
                "MOL2 UNITY_ATOM_ATTR charge -1 / parsed source graph"
            ),
            "source_partial_charge_sum": source_partial_charge_sum,
            "generated_partial_charge_sum": generated_partial_charge_sum,
            "partial_charge_model": "MMFF94_CHARGES",
            "partial_charge_provenance": "checked-in ampc_1.mol2 atom records",
            "generic_protonation_state_metadata": generic_protonation_metadata,
            "metadata_interpretation": "not a molecular formal charge",
        }
    ]
)
show_table(preparation_frame)
print("Generated AMPC portable parameters (not used for ADA):", generated_tmol_path)

^ⓘsource_atom	generated_atom	element	source_Tripos_subtype	source_formal_charge	source_MMFF94_partial_charge	generated_TMol_atom_type	generated_partial_charge
O1	O1	O	2	0	-0.8167	OG2	-0.8167
S1	S1	S	3	0	1.3448	SG3	1.3448
O2	O2	O	2	0	-0.8167	OG2	-0.8167
O3	O3	O	co2	-1	-0.8167	OG2	-0.8167
C1	C1	C	3	0	0.1052	CS2	0.1052
C2	C2	C	3	0	0.0000	CS2	0.0000
C3	C3	C	3	0	0.2300	CS2	0.2300
S2	S2	S	3	0	-0.2860	Ssl	-0.2860
C4	C4	C	ar	0	0.4341	CR	0.4341
N1	N1	N	ar	0	-0.5653	Nim	-0.5653
(17 more rows not shown)

^ⓘpipeline_role	residue_name	atoms	bonds	torsions	molecular_formal_charge_from_source_graph	formal_charge_provenance	source_partial_charge_sum	generated_partial_charge_sum	partial_charge_model	partial_charge_provenance	generic_protonation_state_metadata	metadata_interpretation
independent AMPC MOL2 preparation demo	AMP	27	28	5	-1	MOL2 UNITY_ATOM_ATTR charge -1 / parsed source graph	-1.0001	-1.0001	MMFF94_CHARGES	checked-in ampc_1.mol2 atom records	neutral	not a molecular formal charge

Generated AMPC portable parameters (not used for ADA): /tmp/tmol-ligand-tutorial-x0evhkzd/AMPC.generated.tmol

Build a reusable context#

Parameter injection returns a new immutable database. Derive a BiotitePoseBuildContext from it and use that same database for scoring. Checks confirm the default database is unchanged and the context holds the extended database. An unknown LG1 usually indicates a mismatched parameter file or context.

[6]:

base_database = ParameterDatabase.get_default()
extended_database = inject_params_file(base_database, tmol_path)

complex_array = biotite.structure.io.load_structure(
    str(complex_cif), model=1, include_bonds=True
)
if isinstance(complex_array, struc.AtomArrayStack):
    complex_array = complex_array[0]

pose_diagnostics = StringIO()
try:
    with redirect_stdout(pose_diagnostics), redirect_stderr(pose_diagnostics):
        context = build_context_from_biotite(
            complex_array,
            device,
            param_db=extended_database,
            prepare_ligands=False,
        )
        pose_stack = pose_stack_from_biotite(
            complex_array,
            device,
            context=context,
            no_optH=True,
        )
except Exception:
    print(pose_diagnostics.getvalue())
    raise
score_function = beta2016_score_function(
    device, param_db=context.parameter_database
)

print("default residues:", len(base_database.chemical.residues))
print("extended residues:", len(extended_database.chemical.residues))
print("context reuses extended database:", context.parameter_database is extended_database)

default residues: 312
extended residues: 313
context reuses extended database: True

Select ligand and pocket atoms#

Select LG1 through Biotite’s NumPy annotations. An empty selection usually indicates that CIF and .tmol residue names differ.

Select complete protein residues with a heavy atom within 4.5 Å of the ligand. This defines the compact pocket; packing and minimization use the separate 5 Å rules below. Keep input atom masks separate from TMol block masks.

[7]:

ligand_query = complex_array.res_name == LIGAND_RES_NAME
ligand_atoms = complex_array[ligand_query]
if ligand_atoms.array_length() == 0:
    raise RuntimeError(
        f"No {LIGAND_RES_NAME} atoms found; check CIF and parameter residue names"
    )

ligand_heavy_query = ligand_query & (complex_array.element != "H")
ligand_heavy_coords = complex_array.coord[ligand_heavy_query]
all_to_ligand = complex_array.coord[:, None, :] - ligand_heavy_coords[None, :, :]
nearest_ligand_distance = np.linalg.norm(all_to_ligand, axis=-1).min(axis=1)
pocket_contact_query = (
    (~ligand_query)
    & (complex_array.element != "H")
    & (nearest_ligand_distance <= 4.5)
)
pocket_residues = sorted(
    {
        (str(chain), int(resid), str(name))
        for chain, resid, name in zip(
            complex_array.chain_id[pocket_contact_query],
            complex_array.res_id[pocket_contact_query],
            complex_array.res_name[pocket_contact_query],
        )
    }
)
pocket_query = np.zeros(complex_array.array_length(), dtype=bool)
for chain, resid, name in pocket_residues:
    pocket_query |= (
        (complex_array.chain_id == chain)
        & (complex_array.res_id == resid)
        & (complex_array.res_name == name)
    )
pocket_atoms = complex_array[pocket_query]
selection_table = pd.DataFrame(
    [
        {"selection": "ligand", "atoms": ligand_atoms.array_length(), "residues": 1},
        {"selection": "4.5 Å pocket", "atoms": pocket_atoms.array_length(), "residues": len(pocket_residues)},
    ]
)
show_table(selection_table)
display(
    tmol.selection_gallery(
        complex_array,
        {"ligand": ligand_query, "4.5 Å pocket": pocket_query},
    )
)
pocket_residues[:10]

^ⓘselection	atoms	residues
ligand	36	1
4.5 Å pocket	270	17

/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

pick a selection · drag to rotate · scroll to zoom · click a highlighted atom to label it

[7]:

[('', 17, 'HIS'),
 ('', 19, 'ASP'),
 ('', 58, 'LEU'),
 ('', 61, 'PHE'),
 ('', 62, 'LEU'),
 ('', 65, 'PHE'),
 ('', 66, 'ASP'),
 ('', 101, 'ARG'),
 ('', 102, 'TYR'),
 ('', 103, 'SER')]

Score ligand–protein interactions#

Sum both ligand→protein and protein→ligand entries of the weighted block-pair matrix. This measures interactions within the bound complex, without separated states, reorganization, or solvent correction.

Check for finite scores. The viewer uses PDB coordinates; CIF/.tmol remain authoritative for ligand chemistry.

[8]:

ligand_blocks = ligand_block_mask(pose_stack)
real_blocks = pose_stack.block_type_ind >= 0
protein_blocks = real_blocks & ~ligand_blocks

block_pair_scorer = score_function.render_block_pair_scoring_module(pose_stack)
block_pair_scores = block_pair_scorer(pose_stack.coords)
interaction_score = (
    block_pair_scores[0]
    * ligand_blocks[0, :, None]
    * protein_blocks[0, None, :]
).sum() + (
    block_pair_scores[0]
    * protein_blocks[0, :, None]
    * ligand_blocks[0, None, :]
).sum()

show_table(
    pd.DataFrame(
        [
            {
                "metric": "weighted ligand–protein block-pair interaction",
                "value": float(interaction_score.detach().cpu()),
            }
        ]
    )
)

try:
    viewer = tmol.view(pose_stack, zoom_to={"resn": LIGAND_RES_NAME})
    viewer.setStyle(
        {"resn": LIGAND_RES_NAME},
        {"stick": {"colorscheme": "cyanCarbon", "radius": 0.22}},
    )
    viewer.show()
except ImportError:
    print("Install py3Dmol for the interactive ligand-pocket view.")

^ⓘmetric	value
weighted ligand–protein block-pair interaction	12.211267

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Refine the pocket#

Packing uses Dunbrack, fixed-amino-acid-chi, and input-conformation samplers. No ligand-chi sampler is registered, so LG1 stays fixed during packing and moves only during Cartesian minimization.

The three selection rules differ:

Operation	Scope
Build the compact pocket	Complete protein residues with a heavy atom within 4.5 Å of the ligand
Pack	Blocks separated by less than the sum of their enclosing radii plus 5 Å
Minimize	All ligand atoms and protein side-chain atoms within 5 Å of the ligand

The purple preview uses input coordinates. The helper recomputes minimization membership after packing, so atoms can cross the cutoff.

calculate_block_pair_ddg() sums cross-mask interactions. The tables show the total and per-term changes after refinement.

[9]:

refinement_residues = set(pocket_residues)
protein_refinement_query = np.zeros(complex_array.array_length(), dtype=bool)
for chain, resid, name in refinement_residues:
    protein_refinement_query |= (
        (complex_array.chain_id == chain)
        & (complex_array.res_id == resid)
        & (complex_array.res_name == name)
    )
refinement_array = complex_array[ligand_query | protein_refinement_query]
refinement_diagnostics = StringIO()
try:
    with redirect_stdout(refinement_diagnostics), redirect_stderr(refinement_diagnostics):
        refinement_context = build_context_from_biotite(
            refinement_array,
            device,
            param_db=extended_database,
            prepare_ligands=False,
        )
        refinement_pose = pose_stack_from_biotite(
            refinement_array,
            device,
            context=refinement_context,
            no_optH=True,
        )
except Exception:
    print(refinement_diagnostics.getvalue())
    raise
refinement_sfxn = beta2016_score_function(
    device, param_db=refinement_context.parameter_database
)
refinement_ligand_mask = ligand_block_mask(refinement_pose)
movable_coord_mask = build_coord_mask_for_mask_and_interacting_atoms(
    refinement_pose, refinement_ligand_mask
)
refinement_centroids, refinement_radii = compute_block_centroids_and_furthest_dist(
    refinement_pose
)
packing_adjacency = compute_block_adjacency(
    refinement_centroids, refinement_radii, constant=5.0
)
packing_blocks = refinement_ligand_mask | (
    refinement_ligand_mask.unsqueeze(2) & packing_adjacency
).any(dim=1)
movable_real_atom_count = int(
    (movable_coord_mask & refinement_pose.real_atoms).sum().item()
)
ligand_atom_count = int(
    refinement_pose.n_ats_per_block[refinement_ligand_mask].sum().item()
)
show_table(
    pd.DataFrame(
        [
            {
                "complete_protein_residues_selected_at_4.5A": len(refinement_residues),
                "blocks_in_5A_expanded_sphere_pack_scope": int(packing_blocks.sum()),
                "movable_ligand_atoms": ligand_atom_count,
                "protein_sidechain_atoms_in_input_pose_5A_min_shell": (
                    movable_real_atom_count - ligand_atom_count
                ),
                "input_pose_min_shell_atoms_total": movable_real_atom_count,
            }
        ]
    )
)
try:
    movable_viewer = tmol.view(
        refinement_pose,
        highlighted=movable_coord_mask[refinement_pose.real_atoms],
        highlight_color="#7b2cbf",
    )
    movable_viewer.show()
except ImportError:
    print("Install py3Dmol to inspect the movable refinement atoms.")

torch.manual_seed(SEED)
interaction_before = calculate_block_pair_ddg(
    refinement_pose,
    refinement_ligand_mask,
    sfxn=refinement_sfxn,
    minimize=False,
    pack=False,
    database=extended_database,
)
interaction_after, refined_pose = calculate_block_pair_ddg(
    refinement_pose,
    refinement_ligand_mask,
    sfxn=refinement_sfxn,
    minimize=True,
    pack=True,
    database=extended_database,
    return_pose_stack=True,
)


def ligand_protein_scores_by_term(pose):
    scorer = refinement_sfxn.render_block_pair_scoring_module(pose)
    matrices = scorer(pose.coords, sum_terms=False, apply_weights=True)[:, 0]
    ligand = refinement_ligand_mask[0]
    protein = (pose.block_type_ind64[0] >= 0) & ~ligand
    return matrices[:, ligand][:, :, protein].sum((1, 2)) + matrices[:, protein][
        :, :, ligand
    ].sum((1, 2))


before_by_term = ligand_protein_scores_by_term(refinement_pose)
after_by_term = ligand_protein_scores_by_term(refined_pose)
torch.testing.assert_close(before_by_term.sum().reshape(1), interaction_before)
torch.testing.assert_close(after_by_term.sum().reshape(1), interaction_after)
term_frame = pd.DataFrame(
    {
        "term": [score_type.name for score_type in refinement_sfxn.all_score_types()],
        "before": before_by_term.detach().cpu().numpy(),
        "after": after_by_term.detach().cpu().numpy(),
        "delta": (after_by_term - before_by_term).detach().cpu().numpy(),
    }
)
term_frame = term_frame.loc[term_frame["delta"].abs().sort_values(ascending=False).index]
show_table(term_frame.head(10))
plot_terms = term_frame.head(8).sort_values("delta")
fig, ax = plt.subplots(figsize=(8, 4))
ax.barh(plot_terms["term"], plot_terms["delta"], color="#4c78a8")
ax.axvline(0, color="black", linewidth=0.8)
ax.set(
    xlabel="weighted interaction-score change",
    title="Largest ligand–protein term changes after refinement",
)
plt.tight_layout()
plt.show()

refinement_total_scorer = refinement_sfxn.render_whole_pose_scoring_module(
    refinement_pose
)
refined_total_scorer = refinement_sfxn.render_whole_pose_scoring_module(refined_pose)
total_before = float(refinement_total_scorer(refinement_pose.coords).detach().cpu()[0])
total_after = float(refined_total_scorer(refined_pose.coords).detach().cpu()[0])
real_atoms = refinement_pose.real_atoms & refined_pose.real_atoms
coordinate_delta = refined_pose.coords[real_atoms] - refinement_pose.coords[real_atoms]
all_real_atom_rms = float(
    torch.sqrt(torch.mean(torch.sum(coordinate_delta.square(), dim=-1))).detach().cpu()
)
movable_atoms = movable_coord_mask & real_atoms
movable_coordinate_delta = (
    refined_pose.coords[movable_atoms] - refinement_pose.coords[movable_atoms]
)
movable_atom_rms = float(
    torch.sqrt(torch.mean(torch.sum(movable_coordinate_delta.square(), dim=-1)))
    .detach()
    .cpu()
)

# The whole-pocket total contains protein intramolecular and one-body terms as
# well as ligand–protein interactions; it is not an interaction or binding score.
refinement_frame = pd.DataFrame(
    [
        {
            "stage": "input pocket",
            "whole_pocket_total_score": total_before,
            "ligand_protein_interaction": float(interaction_before.detach().cpu()[0]),
            "RMS_displacement_all_real_atoms_A": 0.0,
            "RMS_displacement_input_pose_5A_mask_atoms_A": 0.0,
        },
        {
            "stage": "repacked + minimized",
            "whole_pocket_total_score": total_after,
            "ligand_protein_interaction": float(interaction_after.detach().cpu()[0]),
            "RMS_displacement_all_real_atoms_A": all_real_atom_rms,
            "RMS_displacement_input_pose_5A_mask_atoms_A": movable_atom_rms,
        },
    ]
)
show_table(refinement_frame)
display(
    tmol.switchable_view(
        {"input pocket": refinement_pose, "refined pocket": refined_pose},
        notes={
            "input pocket": f"interaction {float(interaction_before.detach().cpu()[0]):.3f}",
            "refined pocket": (
                f"interaction {float(interaction_after.detach().cpu()[0]):.3f}; "
                f"all-real-atom RMS motion {all_real_atom_rms:.3f} Å; "
                f"input-pose 5 Å-mask atom RMS motion {movable_atom_rms:.3f} Å"
            ),
        },
    )
)

^ⓘcomplete_protein_residues_selected_at_4.5A	blocks_in_5A_expanded_sphere_pack_scope	movable_ligand_atoms	protein_sidechain_atoms_in_input_pose_5A_min_shell	input_pose_min_shell_atoms_total
17	18	36	133	169

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

^ⓘ	term	before	after	delta
1	fa_ljrep	35.176788	1.671224	-33.505566
0	fa_ljatr	-28.486050	-21.700432	6.785618
17	lk_ball	6.559413	2.267619	-4.291794
2	fa_lk	8.610508	4.721464	-3.889044
16	lk_ball_iso	-6.101487	-4.003672	2.097815
3	fa_elec	-0.514676	-1.781459	-1.266783
4	hbond	-1.081031	0.000000	1.081031
19	lk_bridge_uncpl	-0.051930	-0.288215	-0.236285
18	lk_bridge	-0.011995	-0.040149	-0.028154
8	cart_impropers	0.000000	0.000000	0.000000

../_images/tutorial_07_ligand_and_params_17_3.png

^ⓘstage	whole_pocket_total_score	ligand_protein_interaction	RMS_displacement_all_real_atoms_A	RMS_displacement_input_pose_5A_mask_atoms_A
input pocket	58.221756	14.099542	0.000000	0.000000
repacked + minimized	-12.260399	-19.153622	1.394061	1.607037

Structure:

interaction 14.100

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

interaction -19.154; all-real-atom RMS motion 1.394 Å; input-pose 5 Å-mask atom RMS motion 1.607 Å

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Next#

Try ligand pose sensitivity or DNA/RNA examples.

Exercises#

Reload AMPC.generated.tmol and verify that AMPC atom names, partial charges, and cartbonded payload survive that .tmol round trip. Confirm that its source formal charge and partial-charge sum remain separate concepts. Do not compare it numerically with the chemically unrelated pinned ADA/LG1 artifact as a parity test.
Inspect which chi_samples and proton-chi declarations the .tmol bundle carries for LG1, and confirm they survive a write/reload cycle unchanged.
Change the 4.5 Å AtomArray pocket cutoff while leaving the two 5 Å protocol scopes unchanged, then quantify the effect on runtime and refined interaction score.
Separate ligand desolvation and protein desolvation contributions in the per-term interaction table.
Repeat pocket refinement with several packer seeds and report the interaction-score and structural spread.

References#

TMol ligand workflow

Download this notebook

Tutorial 07 — Ligands and parameter files#

Setup#

The .tmol parameter bundle#

Verify the parameter round trip#

Prepare AMPC from MOL2#

Build a reusable context#

Select ligand and pocket atoms#

Score ligand–protein interactions#

Refine the pocket#

Next#

Exercises#

References#

The `.tmol` parameter bundle#