Tutorial 03 — Scoring and analysis#

Inspect whole-pose, per-term, and block-pair scores; differentiate coordinates; and decompose ligand-fragment interactions. Assumes structure I/O and batching.

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.

This setup fixes seeds, loads a 30-residue 1UBQ slice for the general scoring examples, and defines a table helper that uses itables when it is installed. The fragment section later uses a second checked-in fixture: a prepared ACE protein–inhibitor complex in CIF format plus its matching deterministic ligand parameter file.

[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/cif/1UBQ.cif",
            "tmol/tests/data/protein_ligand_test/ace.tmol.nomin.cif",
            "tmol/tests/data/protein_ligand_test/ace.xtal-lig.mmff94.tmol",
        ]
    )

[2]:

from collections import deque
from contextlib import redirect_stderr, redirect_stdout
from io import StringIO
from pathlib import Path

import biotite.structure as struc
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import torch
from IPython.display import display
from biotite.structure.io import load_structure

import tmol
from tmol.database import ParameterDatabase
from tmol.io import atom_records_from_pose_stack, pose_stack_from_biotite
from tmol.ligand import FRAGMENT_ID_ANNOTATION, load_params_file
from tmol.score import (
    ScoreFunction,
    ScoreType,
    beta2016_score_function,
    calculate_fragment_interactions,
)

SEED = 20260807
np.random.seed(SEED)
torch.manual_seed(SEED)
if torch.cuda.is_available():
    torch.cuda.manual_seed_all(SEED)

device = (
    torch.device("cuda", torch.cuda.current_device())
    if torch.cuda.is_available()
    else torch.device("cpu")
)
repo_root = Path.cwd()
if not (repo_root / "tmol/tests/data/cif/1UBQ.cif").exists():
    repo_root = Path(tmol.__file__).resolve().parents[1]
cif_path = repo_root / "tmol/tests/data/cif/1UBQ.cif"
atom_array = load_structure(str(cif_path), model=1, include_bonds=True)
protein_slice = atom_array[(atom_array.chain_id == "A") & (atom_array.res_id <= 30)]
pose_diagnostics = StringIO()
try:
    with redirect_stdout(pose_diagnostics), redirect_stderr(pose_diagnostics):
        # Optimize hydrogen coordinates before interpreting all-atom scores.
        pose_stack = pose_stack_from_biotite(protein_slice, device, no_optH=False)
except Exception:
    print(pose_diagnostics.getvalue())
    raise
score_function = beta2016_score_function(device)


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

print(f"device={device}; input={cif_path.name}; blocks={pose_stack.max_n_blocks}")

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

device=cpu; input=1UBQ.cif; blocks=30

Whole-pose weighted and unweighted terms#

render_whole_pose_scoring_module() binds a scorer to the pose layout. The default call returns one weighted total per pose. Set sum_terms=False to retain the term axis and apply_weights=False for raw term values.

[3]:

whole_scorer = score_function.render_whole_pose_scoring_module(pose_stack)
weighted_terms = whole_scorer(
    pose_stack.coords, sum_terms=False, apply_weights=True
)
unweighted_terms = whole_scorer(
    pose_stack.coords, sum_terms=False, apply_weights=False
)
total = whole_scorer(pose_stack.coords)
score_types = score_function.all_score_types()
weights = score_function.weights_tensor().detach().cpu().numpy()

term_frame = pd.DataFrame(
    {
        "term": [score_type.name for score_type in score_types],
        "weight": weights,
        "unweighted": unweighted_terms[:, 0].detach().cpu().numpy(),
        "weighted": weighted_terms[:, 0].detach().cpu().numpy(),
    }
)
term_frame["abs_weighted"] = term_frame["weighted"].abs()
term_frame = term_frame.sort_values("abs_weighted", ascending=False).drop(
    columns="abs_weighted"
)
show_table(term_frame)
print(f"weighted term sum: {weighted_terms[:, 0].sum().item():.4f}")
print(f"whole-pose total: {total[0].item():.4f}")

^ⓘ	term	weight	unweighted	weighted
0	fa_ljatr	1.00	-126.766144	-126.766144
2	fa_lk	1.00	99.107513	99.107513
14	dunbrack_rotdev	0.69	106.813515	73.701324
16	lk_ball_iso	-0.38	137.684601	-52.320148
17	lk_ball	0.92	56.536362	52.013454
3	fa_elec	1.00	-41.324898	-41.324898
6	cart_angles	0.50	74.373482	37.186741
15	dunbrack_semirot	0.78	27.911692	21.771118
4	hbond	1.00	-19.705750	-19.705750
1	fa_ljrep	0.55	35.223633	19.372999
(15 more rows not shown)

weighted term sum: 94.1082
whole-pose total: 94.1082

Build a focused score function#

A ScoreFunction starts empty. A nonzero weight loads the corresponding term implementation; setting it to zero removes that contribution without guaranteeing the implementation is unloaded.

This example scores repulsion and hydrogen bonding, then zeros the hydrogen-bond weight and rescores. Use such subsets for diagnostics, not as substitutes for a validated full score function.

[4]:

focused_score_function = ScoreFunction(ParameterDatabase.get_default(), device)
focused_score_function.set_weight(ScoreType.fa_ljrep, 0.55)
focused_score_function.set_weight(ScoreType.hbond, 1.0)
focused_scorer = focused_score_function.render_whole_pose_scoring_module(
    pose_stack
)
focused_score = float(focused_scorer(pose_stack.coords).detach().cpu()[0])
active_before = [term.name for term in focused_score_function.all_score_types()]

focused_score_function.set_weight(ScoreType.hbond, 0.0)
active_after = [term.name for term in focused_score_function.all_score_types()]
repulsion_only_scorer = focused_score_function.render_whole_pose_scoring_module(
    pose_stack
)
repulsion_only_score = float(
    repulsion_only_scorer(pose_stack.coords).detach().cpu()[0]
)
assert float(focused_score_function.get_weight(ScoreType.hbond)) == 0.0

show_table(
    pd.DataFrame(
        [
            {
                "stage": "fa_ljrep + hbond",
                "active_score_types": ", ".join(active_before),
                "score_units": focused_score,
            },
            {
                "stage": "hbond deactivated",
                "active_score_types": ", ".join(active_after),
                "score_units": repulsion_only_score,
            },
        ]
    )
)

^ⓘstage	active_score_types	score_units
fa_ljrep + hbond	fa_ljatr, fa_ljrep, fa_lk, hbond	-0.332750
hbond deactivated	fa_ljatr, fa_ljrep, fa_lk	19.372999

Weighted terms must sum to the whole-pose result within floating-point tolerance. Their scales differ and contributions can cancel. Hydrogens were optimized during construction.

Block-pair accounting#

render_block_pair_scoring_module() returns [n_poses, n_blocks, n_blocks]; sum_terms=False adds a leading term axis.

The matrix need not be symmetric. One-body and intra-block terms occupy the diagonal. Each inter-block contribution is stored once, usually in the upper triangle. Sum M[i, j] + M[j, i] for an unordered pair.

The ranking uses weighted fa_ljatr, fa_ljrep, fa_lk, fa_elec, and hbond, excluding pairs in inter_residue_connections. The all-term matrix remains available for total-score accounting. The per-block profile splits each off-diagonal contribution equally between partners.

[5]:

block_scorer = score_function.render_block_pair_scoring_module(pose_stack)
block_pair_by_term = block_scorer(
    pose_stack.coords, sum_terms=False, apply_weights=True
)
block_pair_total = block_scorer(
    pose_stack.coords, sum_terms=True, apply_weights=True
)

matrix = block_pair_total[0].detach().cpu().numpy()
chain_labels = np.asarray(pose_stack.pdb_info.chain_labels[0]).astype(str)
residue_labels = np.asarray(pose_stack.pdb_info.residue_labels[0]).astype(str)
insertion_codes = np.asarray(
    pose_stack.pdb_info.residue_insertion_codes[0]
).astype(str)
block_labels = np.asarray(
    [
        f"{chain}:{residue}{insertion}"
        for chain, residue, insertion in zip(
            chain_labels, residue_labels, insertion_codes
        )
    ]
)

fig, ax = plt.subplots(figsize=(7, 6))
image = ax.imshow(matrix, cmap="coolwarm", aspect="equal")
tick_step = max(1, len(block_labels) // 10)
tick_positions = np.arange(0, len(block_labels), tick_step)
ax.set_xticks(tick_positions, block_labels[tick_positions], rotation=90)
ax.set_yticks(tick_positions, block_labels[tick_positions])
ax.set(
    title="Directed beta2016 block-pair storage (all terms)",
    xlabel="author chain:residue label",
    ylabel="author chain:residue label",
)
fig.colorbar(image, ax=ax, label="weighted beta2016 score units")
plt.tight_layout()
plt.show()

print("per-term block-pair shape:", tuple(block_pair_by_term.shape))
print("summed block-pair shape:", tuple(block_pair_total.shape))
print("matrix sum:", matrix.sum(), "whole-pose score:", total[0].item())
np.testing.assert_allclose(matrix.sum(), total[0].item(), rtol=1e-4, atol=1e-3)

# Rank only direct nonbonded terms, combining both directed storage locations.
nonbonded_score_types = [
    ScoreType.fa_ljatr,
    ScoreType.fa_ljrep,
    ScoreType.fa_lk,
    ScoreType.fa_elec,
    ScoreType.hbond,
]
nonbonded_term_indices = [score_types.index(term) for term in nonbonded_score_types]
nonbonded_directed = (
    block_pair_by_term[nonbonded_term_indices, 0].detach().cpu().numpy()
)
nonbonded_two_orientation = nonbonded_directed + nonbonded_directed.transpose(0, 2, 1)
nonbonded_pair_matrix = nonbonded_two_orientation.sum(axis=0)

# Exclude directly covalently connected blocks using the pose's chemical graph.
n_blocks = matrix.shape[0]
covalently_connected = np.zeros((n_blocks, n_blocks), dtype=bool)
connection_partners = (
    pose_stack.inter_residue_connections64[0, :, :, 0].detach().cpu().numpy()
)
for block_index, partners in enumerate(connection_partners):
    for partner in partners:
        if 0 <= partner < n_blocks:
            covalently_connected[block_index, partner] = True
            covalently_connected[partner, block_index] = True
candidate_mask = np.triu(np.ones_like(covalently_connected), k=1)
candidate_mask &= ~covalently_connected
ranked_nonbonded = np.where(candidate_mask, nonbonded_pair_matrix, np.inf)
pair_i, pair_j = np.unravel_index(
    np.argmin(ranked_nonbonded), ranked_nonbonded.shape
)

pair_term_breakdown = pd.DataFrame(
    {
        "term": [term.name for term in nonbonded_score_types],
        "stored_i_to_j": nonbonded_directed[:, pair_i, pair_j],
        "stored_j_to_i": nonbonded_directed[:, pair_j, pair_i],
        "two_orientation_sum": nonbonded_two_orientation[:, pair_i, pair_j],
    }
)
show_table(pair_term_breakdown)

# Optional analytical profile: diagonal once, every all-term off-diagonal pair half each.
all_term_two_orientation = matrix + matrix.T
np.fill_diagonal(all_term_two_orientation, 0.0)
equal_split_profile = np.diag(matrix) + 0.5 * all_term_two_orientation.sum(axis=1)
np.testing.assert_allclose(equal_split_profile.sum(), matrix.sum(), atol=1e-3)
equal_split_frame = pd.DataFrame(
    {
        "block": np.arange(n_blocks),
        "pdb_label": block_labels,
        "equal_split_all_term_score_units": equal_split_profile,
    }
).sort_values("equal_split_all_term_score_units")
show_table(equal_split_frame)


def deposited_mask_for_block(block_index):
    mask = (
        (protein_slice.chain_id.astype(str) == chain_labels[block_index])
        & (protein_slice.res_id.astype(str) == residue_labels[block_index])
    )
    insertion = insertion_codes[block_index]
    if insertion:
        mask &= protein_slice.ins_code.astype(str) == insertion
    return mask


pair_masks = {
    f"block {pair_i} / {block_labels[pair_i]}": deposited_mask_for_block(pair_i),
    f"block {pair_j} / {block_labels[pair_j]}": deposited_mask_for_block(pair_j),
    "selected nonbonded pair": deposited_mask_for_block(pair_i)
    | deposited_mask_for_block(pair_j),
}
print(
    f"Most favorable noncovalently connected pair by the five selected terms: "
    f"blocks {pair_i}, {pair_j} ({block_labels[pair_i]}, {block_labels[pair_j]}); "
    f"two-orientation weighted score "
    f"{nonbonded_pair_matrix[pair_i, pair_j]:.3f} score units"
)
try:
    display(tmol.selection_gallery(protein_slice, pair_masks))
except ImportError as exc:
    print("Interactive residue-pair gallery unavailable:", exc)

../_images/tutorial_03_scoring_and_analysis_9_0.png

per-term block-pair shape: (25, 1, 30, 30)
summed block-pair shape: (1, 30, 30)
matrix sum: 94.10825 whole-pose score: 94.10823822021484

^ⓘterm	stored_i_to_j	two_orientation_sum
fa_ljatr	-3.301455	-3.301455
fa_ljrep	0.089763	0.089763
fa_lk	3.824763	3.824763
fa_elec	-4.535543	-4.535543
hbond	-2.508232	-2.508232

^ⓘ	block	pdb_label	equal_split_all_term_score_units
4	4	A:5	-1.889499
16	16	A:17	-1.736070
2	2	A:3	-1.693172
20	20	A:21	-1.591118
18	18	A:19	-1.341876
25	25	A:26	-1.169288
21	21	A:22	-0.683446
5	5	A:6	-0.446275
11	11	A:12	-0.022984
9	9	A:10	-0.009457
(20 more rows not shown)

Most favorable noncovalently connected pair by the five selected terms: blocks 4, 12 (A:5, A:13); two-orientation weighted score -6.431 score units

/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

Both the full matrix sum and the equal-split profile reproduce the whole-pose score. A row sum alone misses contributions stored in the opposite orientation. The profile is an attribution convention, not a cached residue energy.

The selected pair table combines both orientations of the five nonbonded terms and excludes directly covalently connected blocks.

Ligand-fragment interactions in one connected complex#

The next workflow uses the checked-in ace.tmol.nomin.cif protein–small-molecule complex and its matching ace.xtal-lig.mmff94.tmol parameters. This is the same deterministic CIF-plus-parameters path exercised by TMol’s ligand regression tests: no stochastic conformer generation is needed. Tutorial 07 covers ligand preparation; this section assumes the checked-in CIF and .tmol chemistry are already authoritative.

Concrete API and input requirements:

Load a Biotite AtomArray with bond information. The ligand atom/residue names must match a prepared ligand definition (here, residue LG1 in the checked-in .tmol file).
Before pose construction, add one integer tmol_fragment_id per atom. Use 0 outside fragmented residues and positive IDs for every atom of the ligand. At least two IDs must occur.
Each fragment must be connected, contain at least three heavy atoms, and have at most four inter-fragment connections. No atom may participate in two cuts, no four-atom bonded path may cross two cuts, and cuts through hbond/lk-ball acceptor frame geometry are rejected.
Build with prepare_ligands=True, the matching parameter file, and return_context=True. The annotation causes the connected ligand to become fragment block types with explicit inter-block connections.
Construct beta2016_score_function from context.parameter_database, the same ligand-extended ParameterDatabase used for the pose. The default database does not contain the generated fragment scoring parameters.
Supply calculate_fragment_interactions a boolean partner mask with shape [n_poses, max_n_blocks], on the pose device, that excludes every ligand fragment block. Here it selects amino-acid polymer blocks only.

The result decomposes fragment–partner interactions in the connected complex. Fragment–fragment energies remain separate in the block-pair tensor.

[6]:

ligand_data_dir = repo_root / "tmol/tests/data/protein_ligand_test"
complex_cif_path = ligand_data_dir / "ace.tmol.nomin.cif"
ligand_params_path = ligand_data_dir / "ace.xtal-lig.mmff94.tmol"

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

ligand_name = "LG1"
is_ligand = complex_array.res_name == ligand_name
assert is_ligand.any(), f"{ligand_name} was not found in {complex_cif_path.name}"
preparation = load_params_file(ligand_params_path)[0]
restype = preparation.residue_type

# Functional-group-sized cuts validated by the fragmented-ligand regression test.
cut_bonds = (
    ("C1", "C2"),   # amide/pyrrolidine arm from the central scaffold
    ("C3", "C4"),   # central carboxylate from its substituted carbon
    ("C5", "C9"),   # pyrrolidine carboxylate from the ring
    ("C11", "C12"), # terminal aminoethyl group from the alkyl linker
    ("C15", "C16"), # phenyl ring from the ethyl linker
)


def components_after_cuts(residue_type, cuts):
    """Connected components of prepared ligand chemistry after conceptual cuts."""
    removed = {frozenset(cut) for cut in cuts}
    adjacency = {atom.name: set() for atom in residue_type.atoms}
    for atom_a, atom_b, *_ in residue_type.bonds:
        if frozenset((atom_a, atom_b)) not in removed:
            adjacency[atom_a].add(atom_b)
            adjacency[atom_b].add(atom_a)

    components = []
    unseen = set(adjacency)
    while unseen:
        queue = deque([next(iter(unseen))])
        component = set()
        while queue:
            atom_name = queue.popleft()
            if atom_name in component:
                continue
            component.add(atom_name)
            queue.extend(adjacency[atom_name] - component)
        unseen -= component
        components.append(component)

    atom_order = {atom.name: index for index, atom in enumerate(residue_type.atoms)}
    components.sort(key=lambda component: min(atom_order[name] for name in component))
    return components


components = components_after_cuts(restype, cut_bonds)
fragment_for_atom = {
    atom_name: fragment_id
    for fragment_id, component in enumerate(components, start=1)
    for atom_name in component
}
fragment_region = {
    1: "amide/pyrrolidine core",
    2: "central carboxylate",
    3: "pyrrolidine carboxylate",
    4: "central linker scaffold",
    5: "terminal aminoethyl group",
    6: "phenyl ring",
}
assert set(fragment_region) == set(range(1, len(components) + 1))

fragment_ids = np.zeros(complex_array.array_length(), dtype=np.int32)
for atom_index in np.flatnonzero(is_ligand):
    atom_name = str(complex_array.atom_name[atom_index])
    fragment_ids[atom_index] = fragment_for_atom[atom_name]
assert np.all(fragment_ids[is_ligand] > 0)

annotated_complex = complex_array.copy()
annotated_complex.set_annotation(FRAGMENT_ID_ANNOTATION, fragment_ids)

atom_type_by_name = {atom.name: atom.atom_type for atom in restype.atoms}
fragment_definition_frame = pd.DataFrame(
    [
        {
            "fragment": fragment_id,
            "region": fragment_region[fragment_id],
            "heavy atoms": ", ".join(
                name
                for name in sorted(component)
                if not atom_type_by_name[name].upper().startswith("H")
            ),
            "heavy-atom count": sum(
                not atom_type_by_name[name].upper().startswith("H")
                for name in component
            ),
        }
        for fragment_id, component in enumerate(components, start=1)
    ]
)
cut_frame = pd.DataFrame(
    [
        {
            "cut bond": f"{atom_a}–{atom_b}",
            "fragment A": fragment_for_atom[atom_a],
            "fragment B": fragment_for_atom[atom_b],
        }
        for atom_a, atom_b in cut_bonds
    ]
)
show_table(fragment_definition_frame)
show_table(cut_frame)

^ⓘfragment	region	heavy atoms	heavy-atom count
1	amide/pyrrolidine core	C1, C5, C6, C7, C8, N2, O1	7
2	central carboxylate	C3, O2, O3	3
3	pyrrolidine carboxylate	C9, O4, O5	3
4	central linker scaffold	C10, C11, C14, C15, C2, C4, N1	7
5	terminal aminoethyl group	C12, C13, N3	3
6	phenyl ring	C16, C17, C18, C19, C20, C21	6

^ⓘcut bond	fragment A	fragment B
C1–C2	1	4
C3–C4	2	4
C5–C9	1	3
C11–C12	4	5
C15–C16	4	6

What is being cut?#

The five selected single bonds separate recognizable functional-group-sized regions while leaving every fragment connected and above the three-heavy-atom minimum. Conceptually, the fragment graph is:

F3 pyrrolidine carboxylate -- C9–C5 -- F1 amide/pyrrolidine -- C1–C2 -- F4 central scaffold -- C11–C12 -- F5 aminoethyl
                                                                           |                 |
                                                                        C3–C4           C15–C16
                                                                           |                 |
                                                               F2 central carboxylate   F6 phenyl ring

Thus the five fragment-graph edges are F3–F1, F1–F4, F2–F4, F4–F5, and F4–F6; in particular, the central carboxylate F2 attaches to the central scaffold F4, not to F1.

The diagram depicts conceptual partition boundaries, not broken chemistry in the scored complex. TMol prepares LG1 as one molecule, creates one fragment block per connected component, and installs paired connections across every listed bond. Bonded separation and bonded terms can therefore traverse those explicit inter-block links.

[7]:

fragment_graph = nx.Graph()
fragment_graph.add_nodes_from(fragment_region)
for atom_a, atom_b in cut_bonds:
    fragment_graph.add_edge(
        fragment_for_atom[atom_a],
        fragment_for_atom[atom_b],
        bond=f"{atom_a}–{atom_b}",
    )
fragment_positions = {
    3: (-2.0, 0.0),
    1: (-1.0, 0.0),
    4: (0.0, 0.0),
    5: (1.2, 0.0),
    2: (0.0, 1.0),
    6: (0.0, -1.0),
}
fragment_colors = plt.cm.tab10(np.linspace(0, 1, len(fragment_region)))
fig, ax = plt.subplots(figsize=(10, 5))
nx.draw_networkx_nodes(
    fragment_graph,
    fragment_positions,
    node_color=fragment_colors,
    node_size=2400,
    ax=ax,
)
nx.draw_networkx_edges(fragment_graph, fragment_positions, width=2, ax=ax)
nx.draw_networkx_labels(
    fragment_graph,
    fragment_positions,
    labels={fragment_id: f"F{fragment_id}" for fragment_id in fragment_region},
    font_weight="bold",
    ax=ax,
)
nx.draw_networkx_edge_labels(
    fragment_graph,
    fragment_positions,
    edge_labels=nx.get_edge_attributes(fragment_graph, "bond"),
    font_size=9,
    ax=ax,
)
ax.set_title("LG1 fragment connectivity retained across conceptual cut bonds")
ax.axis("off")
plt.tight_layout()
plt.show()

../_images/tutorial_03_scoring_and_analysis_14_0.png

Build the connected fragment-block pose#

tmol_fragment_id is already present on annotated_complex before this call. pose_stack_from_biotite first prepares the complete ligand from the checked-in definition, then expands it into LG1.1 through LG1.6 blocks and attaches an explicit pair of connections for each conceptual cut.

The fixture contains deposited, already prepared hydrogens matching the frozen parameter file, so this deterministic path uses no_optH=True and disables proton-chi sampling, as in the regression test. For generated or incompletely hydrogenated inputs, use an appropriately validated hydrogen-preparation/optimization protocol instead.

[8]:

fragment_pose, fragment_context = pose_stack_from_biotite(
    annotated_complex,
    device,
    param_db=ParameterDatabase.get_default(),
    prepare_ligands=True,
    ligand_params_files=[str(ligand_params_path)],
    no_optH=True,
    return_context=True,
)
fragment_mapping = fragment_pose.split_block_mapping
fragment_entries = sorted(
    (entry for entry in fragment_mapping.entries if entry.pose_ind == 0),
    key=lambda entry: entry.block_ind,
)
assert len(fragment_entries) == len(fragment_region)
fragment_id_by_block = {
    entry.block_ind: fragment_id
    for fragment_id, entry in enumerate(fragment_entries, start=1)
}

fragment_block_mask = torch.zeros_like(
    fragment_pose.block_type_ind, dtype=torch.bool
)
for entry in fragment_mapping.entries:
    fragment_block_mask[entry.pose_ind, entry.block_ind] = True

# Select protein polymer blocks explicitly; do not use the complement blindly.
protein_partner_mask = torch.zeros_like(fragment_block_mask)
for pose_index in range(fragment_pose.n_poses):
    for block_index in range(fragment_pose.max_n_blocks):
        block_type_index = int(
            fragment_pose.block_type_ind64[pose_index, block_index].item()
        )
        if block_type_index < 0:
            continue
        block_type = fragment_pose.packed_block_types.active_block_types[
            block_type_index
        ]
        polymer = block_type.properties.polymer
        protein_partner_mask[pose_index, block_index] = (
            polymer.is_polymer and polymer.polymer_type == "amino_acid"
        )

assert protein_partner_mask.dtype == torch.bool
assert protein_partner_mask.shape == fragment_pose.block_type_ind.shape
assert protein_partner_mask.device == fragment_pose.device
assert not torch.any(protein_partner_mask & fragment_block_mask)
assert torch.any(protein_partner_mask)

# Scoring must use the exact ligand-extended database returned by this build.
fragment_score_function = beta2016_score_function(
    device,
    param_db=fragment_context.parameter_database,
)

print(
    f"fragment blocks={int(fragment_block_mask.sum())}; "
    f"protein partner blocks={int(protein_partner_mask.sum())}; "
    f"annotated cut bonds={len(cut_bonds)}"
)

fragment blocks=6; protein partner blocks=574; annotated cut bonds=5

Calculate fragment–protein interactions#

calculate_fragment_interactions() sums both orientations between each fragment and the selected protein blocks. SplitBlockEntry records align columns with pose and block indices. The fragment columns sum to the connected ligand–protein cross-mask score.

[9]:

fragment_interactions = calculate_fragment_interactions(
    fragment_pose,
    protein_partner_mask,
    sfxn=fragment_score_function,
    sum_terms=False,
)
fragment_score_types = fragment_score_function.all_score_types()
fragment_term_names = [score_type.name for score_type in fragment_score_types]
fragment_score_matrix = (
    fragment_interactions.scores[:, 0, :].detach().cpu().numpy()
)
assert fragment_score_matrix.shape == (
    len(fragment_term_names),
    len(fragment_interactions.mapping),
)

fragment_rows = []
for column, record in enumerate(fragment_interactions.mapping):
    fragment_id = fragment_id_by_block[record.block_ind]
    row = {
        "fragment": f"F{fragment_id}",
        "region": fragment_region[fragment_id],
        "pose block": record.block_ind,
        "weighted total": fragment_score_matrix[:, column].sum(),
    }
    row.update(
        {
            term_name: fragment_score_matrix[term_index, column]
            for term_index, term_name in enumerate(fragment_term_names)
        }
    )
    fragment_rows.append(row)
fragment_interaction_frame = pd.DataFrame(fragment_rows)
show_table(fragment_interaction_frame.round(4))

# Independent block-pair check: all fragment↔protein entries equal the API sum.
fragment_block_scorer = (
    fragment_score_function.render_block_pair_scoring_module(fragment_pose)
)
fragment_block_pair_terms = fragment_block_scorer(
    fragment_pose.coords, sum_terms=False, apply_weights=True
)
fragment_protein_cross_mask = (
    fragment_block_mask.unsqueeze(2) & protein_partner_mask.unsqueeze(1)
) | (
    protein_partner_mask.unsqueeze(2) & fragment_block_mask.unsqueeze(1)
)
direct_ligand_protein_terms = (
    fragment_block_pair_terms
    * fragment_protein_cross_mask.unsqueeze(0)
).sum(dim=(2, 3))
torch.testing.assert_close(
    fragment_interactions.scores.sum(dim=2),
    direct_ligand_protein_terms,
    rtol=1e-5,
    atol=1e-5,
)

# Keep fragment–fragment interactions in a separate diagnostic table.
fragment_fragment_rows = []
for column_a, record_a in enumerate(fragment_interactions.mapping):
    for column_b, record_b in enumerate(fragment_interactions.mapping):
        if column_b <= column_a:
            continue
        block_a, block_b = record_a.block_ind, record_b.block_ind
        pair_terms = (
            fragment_block_pair_terms[:, 0, block_a, block_b]
            + fragment_block_pair_terms[:, 0, block_b, block_a]
        )
        fragment_a = fragment_id_by_block[block_a]
        fragment_b = fragment_id_by_block[block_b]
        fragment_fragment_rows.append(
            {
                "fragment pair": f"F{fragment_a}–F{fragment_b}",
                "separate weighted interaction": float(pair_terms.sum().detach().cpu()),
            }
        )
fragment_fragment_frame = pd.DataFrame(fragment_fragment_rows)
print("Fragment–fragment interactions (separate; not attributed above):")
show_table(fragment_fragment_frame.round(4))

^ⓘfragment	region	pose block	weighted total	fa_ljatr	fa_ljrep	fa_lk	fa_elec	hbond	lk_ball_iso	lk_ball	lk_bridge	lk_bridge_uncpl
F1	amide/pyrrolidine core	574	-8.7480	-8.6317	0.7737	2.2328	-2.2899	-1.3816	-1.5967	2.1455	0.0000	0.0000
F2	central carboxylate	575	4.7783	-3.7861	0.0801	8.0684	-0.2244	-1.6343	-3.2270	5.5018	0.0000	0.0000
F3	pyrrolidine carboxylate	576	-3.6570	-3.1002	0.0552	6.4262	-5.7848	-1.8866	-2.6292	3.3871	-0.0165	-0.1082
F4	central linker scaffold	577	-3.3711	-10.8923	1.2292	3.8025	1.7800	0.0000	-2.5877	3.2974	0.0000	0.0000
F5	terminal aminoethyl group	578	-4.1776	-2.5273	0.0643	0.7697	-2.3863	0.0000	-0.6968	0.8410	-0.0439	-0.1983
F6	phenyl ring	579	-5.3081	-3.7398	0.0858	-1.2742	-0.2641	0.0000	-0.1342	0.0184	0.0000	0.0000

Fragment–fragment interactions (separate; not attributed above):

^ⓘfragment pair	separate weighted interaction
F1–F2	-0.3553
F1–F3	7.1373
F1–F4	5.4638
F1–F5	-0.2495
F1–F6	0.0000
F2–F3	0.0000
F2–F4	12.8716
F2–F5	0.0000
F2–F6	-0.0973
F3–F4	0.0901
(5 more rows not shown)

[10]:

fragment_labels = [
    f"F{fragment_id_by_block[record.block_ind]}"
    for record in fragment_interactions.mapping
]
fragment_totals = fragment_score_matrix.sum(axis=0)
active_term_mask = np.max(np.abs(fragment_score_matrix), axis=1) > 1e-8
active_fragment_terms = fragment_score_matrix[active_term_mask]
active_fragment_term_names = np.asarray(fragment_term_names)[active_term_mask]
heatmap_limit = max(float(np.max(np.abs(active_fragment_terms))), 1e-8)

fig, (total_ax, term_ax) = plt.subplots(
    2,
    1,
    figsize=(10, 4 + 0.38 * len(active_fragment_term_names)),
    gridspec_kw={"height_ratios": [1, 2]},
    constrained_layout=True,
)
total_colors = np.where(fragment_totals >= 0, "tab:red", "tab:blue")
total_ax.bar(fragment_labels, fragment_totals, color=total_colors)
total_ax.axhline(0, color="black", linewidth=0.8)
total_ax.set(
    ylabel="weighted score units",
    title="Connected-ligand fragment interactions with the ACE protein partner",
)

term_image = term_ax.imshow(
    active_fragment_terms,
    cmap="coolwarm",
    vmin=-heatmap_limit,
    vmax=heatmap_limit,
    aspect="auto",
)
term_ax.set_xticks(np.arange(len(fragment_labels)), fragment_labels)
term_ax.set_yticks(
    np.arange(len(active_fragment_term_names)), active_fragment_term_names
)
term_ax.set(xlabel="ligand fragment", ylabel="weighted beta2016 term")
fig.colorbar(term_image, ax=term_ax, label="weighted score units", shrink=0.85)
plt.show()

# Interactive pocket viewer: one button per fragment in the bound complex, each
# captioned with that fragment's weighted interaction total with the protein
# partner -- the same one-complex interaction score plotted above.
fragment_total_by_id = {
    fragment_id_by_block[record.block_ind]: float(fragment_totals[column])
    for column, record in enumerate(fragment_interactions.mapping)
}
whole_ligand_total = float(fragment_totals.sum())
fragment_pocket_selections = {"all LG1": is_ligand}
fragment_pocket_notes = {
    "all LG1": (
        f"whole connected ligand · summed fragment interaction "
        f"{whole_ligand_total:+.2f} score units"
    )
}
for fragment_id in fragment_region:
    key = f"F{fragment_id}: {fragment_region[fragment_id]}"
    fragment_pocket_selections[key] = is_ligand & (fragment_ids == fragment_id)
    fragment_pocket_notes[key] = (
        f"weighted fragment–protein interaction "
        f"{fragment_total_by_id[fragment_id]:+.2f} score units "
        f"(one bound complex; no state subtraction)"
    )
try:
    try:
        fragment_viewer = tmol.selection_gallery(
            complex_array,
            fragment_pocket_selections,
            notes=fragment_pocket_notes,
        )
    except TypeError as exc:
        if "notes" not in str(exc):
            raise
        # The current release wheel predates selection captions but still
        # supports the interactive fragment buttons and shared 3D viewer.
        fragment_viewer = tmol.selection_gallery(
            complex_array,
            fragment_pocket_selections,
        )
    display(fragment_viewer)
except ImportError as exc:
    print("Interactive fragment pocket viewer unavailable:", exc)

../_images/tutorial_03_scoring_and_analysis_19_0.png

/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

The viewer and bar plot report the same fragment–protein scores. All LG1 shows their sum. Fragment–fragment entries are displayed separately. The cross-mask assertion checks the decomposition against the block-pair tensor.

Differentiate and reweight an interface#

Scoring modules accept coordinate tensors and support PyTorch autograd. Clone and detach coordinates to avoid modifying the original pose.

First differentiate the 1UBQ total. Then multiply both ligand–protein matrix orientations by 1.5 before summing. This changes the analytical objective, not the beta2016 parameters or the interpretation as a one-complex score.

[11]:

differentiable_coords = pose_stack.coords.detach().clone().requires_grad_(True)
differentiable_total = whole_scorer(differentiable_coords).sum()
differentiable_total.backward()
atom_gradient_norm = differentiable_coords.grad[0].norm(dim=-1)
real_gradient_norm = atom_gradient_norm[pose_stack.real_atoms[0]]

records = atom_records_from_pose_stack(pose_stack)
gradient_frame = pd.DataFrame(
    {
        "chain": records["chain"],
        "residue": records["resi"],
        "res_name": records["resn"],
        "atom": records["atomn"],
        "gradient_norm": real_gradient_norm.detach().cpu().numpy(),
    }
).sort_values("gradient_norm", ascending=False)
show_table(gradient_frame.head(20))
print("finite gradients:", torch.isfinite(real_gradient_norm).all().item())

# Reweight the chemically defined ligand↔protein block-pair interface before
# reduction. Both directed storage orientations receive the same multiplier.
interface_coords = fragment_pose.coords.detach().clone().requires_grad_(True)
interface_matrix = fragment_block_scorer(interface_coords)
interface_weights = torch.ones_like(interface_matrix)
interface_weights = torch.where(
    fragment_protein_cross_mask,
    torch.full_like(interface_weights, 1.5),
    interface_weights,
)
reweighted_interface_objective = (interface_matrix * interface_weights).sum()
reweighted_interface_objective.backward()
interface_gradient = interface_coords.grad[fragment_pose.real_atoms]
assert torch.isfinite(interface_gradient).all()
print(
    "1.5× ligand–protein block-pair objective:",
    float(reweighted_interface_objective.detach().cpu()),
)
print(
    "reweighted interface coordinate-gradient norm:",
    float(torch.linalg.vector_norm(interface_gradient).detach().cpu()),
)

^ⓘ	chain	residue	res_name	atom	gradient_norm
203	A	13	ILE	N	72.699387
207	A	13	ILE	CB	61.871479
204	A	13	ILE	CA	51.504276
38	A	3	ILE	C	37.663975
237	A	15	LEU	CA	37.211075
236	A	15	LEU	N	35.575195
77	A	5	VAL	C	33.410530
238	A	15	LEU	C	33.066227
208	A	13	ILE	CG1	32.681229
55	A	4	PHE	N	28.599651
(10 more rows not shown)

finite gradients: True
1.5× ligand–protein block-pair objective: 108.47476196289062
reweighted interface coordinate-gradient norm: 654.6614990234375

Perturb a nonlocal contact#

Find the closest heavy-atom pair separated by at least three sequence positions and move one atom by 0.10 Å along the interatomic direction. This tests local score sensitivity. The move is not chemically constrained; both bonded and nonbonded terms may respond.

[12]:

pbt = pose_stack.packed_block_types
heavy_indices_by_block = []
for block_index in range(pose_stack.max_n_blocks):
    block_type_index = int(pose_stack.block_type_ind64[0, block_index].item())
    block_type = pbt.active_block_types[block_type_index]
    offset = int(pose_stack.block_coord_offset64[0, block_index].item())
    local_heavy = torch.nonzero(
        pbt.atom_is_hydrogen[block_type_index, : len(block_type.atoms)] == 0,
        as_tuple=False,
    ).flatten()
    heavy_indices_by_block.append(local_heavy + offset)

best_distance = float("inf")
best_contact = None
for block_i in range(pose_stack.max_n_blocks):
    for block_j in range(block_i + 3, pose_stack.max_n_blocks):
        atoms_i = heavy_indices_by_block[block_i]
        atoms_j = heavy_indices_by_block[block_j]
        distances = torch.cdist(
            pose_stack.coords[0, atoms_i], pose_stack.coords[0, atoms_j]
        )
        flat_index = int(torch.argmin(distances).item())
        local_i = flat_index // distances.shape[1]
        local_j = flat_index % distances.shape[1]
        distance = float(distances[local_i, local_j].item())
        if distance < best_distance:
            best_distance = distance
            best_contact = (
                block_i,
                block_j,
                int(atoms_i[local_i].item()),
                int(atoms_j[local_j].item()),
            )

block_i, block_j, atom_i, atom_j = best_contact
perturbed_coords = pose_stack.coords.detach().clone()
direction = perturbed_coords[0, atom_j] - perturbed_coords[0, atom_i]
perturbed_coords[0, atom_j] += 0.10 * direction / torch.linalg.vector_norm(direction)

before_terms = whole_scorer(
    pose_stack.coords, sum_terms=False, apply_weights=True
)[:, 0]
after_terms = whole_scorer(
    perturbed_coords, sum_terms=False, apply_weights=True
)[:, 0]
perturbation_frame = pd.DataFrame(
    {
        "term": [score_type.name for score_type in score_types],
        "weighted_before": before_terms.detach().cpu().numpy(),
        "weighted_after": after_terms.detach().cpu().numpy(),
        "delta": (after_terms - before_terms).detach().cpu().numpy(),
    }
).sort_values("delta", key=np.abs, ascending=False)
show_table(perturbation_frame.head(15))
print(
    f"contact blocks {block_i} and {block_j}; initial heavy-atom distance "
    f"{best_distance:.3f} Å"
)

perturbed_pose = pose_stack.clone()
perturbed_pose.coords.copy_(perturbed_coords)
try:
    display(
        tmol.switchable_view(
            {"original": pose_stack, "0.10 Å perturbation": perturbed_pose},
            notes={
                "original": f"closest selected contact: {best_distance:.3f} Å",
                "0.10 Å perturbation": "One contact atom moved along the interatomic direction",
            },
        )
    )
except ImportError as exc:
    print("Interactive coordinate comparison unavailable:", exc)

^ⓘ	term	weighted_before	weighted_after	delta
5	cart_lengths	6.905261	7.648484	0.743223
3	fa_elec	-41.324898	-41.003883	0.321014
4	hbond	-19.705750	-19.877102	-0.171352
1	fa_ljrep	19.372999	19.221270	-0.151730
2	fa_lk	99.107513	99.018799	-0.088715
7	cart_torsions	8.646235	8.581419	-0.064816
17	lk_ball	52.013454	51.970871	-0.042583
16	lk_ball_iso	-52.320148	-52.284924	0.035225
6	cart_angles	37.186741	37.159916	-0.026825
0	fa_ljatr	-126.766144	-126.760635	0.005508
(5 more rows not shown)

contact blocks 0 and 16; initial heavy-atom distance 2.564 Å

Structure:

closest selected contact: 2.564 Å

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

One contact atom moved along the interatomic direction

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

The perturbation should produce finite, nonzero score changes. Their sizes depend on local geometry and do not establish chemical plausibility or a physical response.

Plot dominant weighted terms#

[13]:

dominant = term_frame.reindex(term_frame["weighted"].abs().sort_values().index).tail(12)
fig, ax = plt.subplots(figsize=(8, 5))
colors = np.where(dominant["weighted"] >= 0, "tab:red", "tab:blue")
ax.barh(dominant["term"], dominant["weighted"], color=colors)
ax.axvline(0, color="black", linewidth=0.8)
ax.set(
    xlabel="weighted beta2016 score units",
    title="Largest beta2016 term contributions",
)
plt.tight_layout()
plt.show()

../_images/tutorial_03_scoring_and_analysis_25_0.png

Next#

Continue with packing and minimization before FastRelax.

Exercises#

Confirm numerically that block_pair_by_term.sum((2, 3)) matches weighted_terms.
Build an unweighted block-pair tensor and inspect how one selected nonbonded term changes after the perturbation.
Select a submatrix by pdb_info chain, residue, and insertion-code labels rather than by assuming file order.
Compare the equal-split profile with an alternative declared attribution convention without calling either a cached residue energy.
Repeat the gradient analysis after cloning and perturbing a different contact; compare cosine similarity between gradients.

References#

Download this notebook