Tutorial 03 — Scoring and Analysis#

This tutorial creates a beta2016 score function, scores a pose, and inspects score terms and block-pair contributions. It also demonstrates coordinate gradients and ligand-fragment accounting.

Learning objectives#

Compare whole-pose, per-term, and block-pair scores.
Inspect one block pair and one fragmented ligand.
Differentiate a score with PyTorch autograd.

Before you begin#

Prerequisites: 01 — Working with TMol and 02 — GPU batching.
Curriculum: This completes the 01 → 02 → 03 foundation. Next, take 04 — Packing and 05 — Minimization in either order; complete both before 06.
Related: Scoring workflow · Scoring API

TMol reports score units, not kcal/mol. This tutorial uses beta2016_score_function(); TMol does not provide ref2015, centroid residue types, or centroid/full-atom switching. The examples run on CPU; CUDA is useful for larger batches.

Setup#

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}")

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

Whole-pose weighted and unweighted terms#

render_whole_pose_scoring_module() binds the score function to this pose layout. With the current scorer API, sum_terms=False retains the score-term axis and apply_weights=False returns raw term values. The default call sums weighted terms into one score per pose. These weighted values are TMol score units: beta2016-weighted term sums, not kcal/mol, thermodynamic free energies, or values guaranteed to match Rosetta score units numerically.

[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.778610	-126.778610
2	fa_lk	1.00	99.107491	99.107491
14	dunbrack_rotdev	0.69	106.813499	73.701317
16	lk_ball_iso	-0.38	137.684586	-52.320141
17	lk_ball	0.92	56.544212	52.020676
3	fa_elec	1.00	-41.767223	-41.767223
6	cart_angles	0.50	74.373535	37.186768
15	dunbrack_semirot	0.78	27.911692	21.771118
4	hbond	1.00	-19.797235	-19.797235
1	fa_ljrep	0.55	35.315784	19.423681
(14 more rows not shown)

weighted term sum: 87.9629
whole-pose total: 87.9629

Build a focused score function and set a weight to zero#

beta2016_score_function() is a complete preset, but a ScoreFunction can also start empty. Setting a nonzero weight loads the implementation that covers that score type. Setting the weight back to zero removes that term’s contribution from the weighted score; callers should not depend on the implementation being unloaded.

The example scores repulsion and hydrogen bonding, then sets the hydrogen-bond weight to zero and rescoring the same coordinates. This is useful for debugging and controlled experiments, not a replacement 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.373554
hbond deactivated	fa_ljatr, fa_ljrep, fa_lk	19.423681

Expected observations. Raw and weighted columns differ wherever a term’s weight is not one. The weighted column must sum to the default whole-pose result within floating-point tolerance. A large magnitude is not automatically an error: bonded, solvation, electrostatic, and reference contributions have different scales and cancellation is common. Hydrogen coordinates were optimized during pose construction because these scores are interpreted scientifically.

Directed accounting and a selected nonbonded pair#

render_block_pair_scoring_module() returns weighted block-pair values with shape [n_poses, n_blocks, n_blocks] when sum_terms=True. Keeping sum_terms=False adds a leading score-term axis.

The all-term tensor is an accounting representation, not a symmetric contact map. One-body and intra-block contributions lie on the diagonal. Inter-block two-body contributions are evaluated once per unordered block pair and stored in one directed entry (normally the upper triangle), so M[i, j] need not equal M[j, i]. Use M[i, j] + M[j, i] for an off-diagonal residue-pair interaction.

For a chemically clearer example, the ranking below uses only weighted fa_ljatr, fa_ljrep, fa_lk, fa_elec, and hbond values and excludes every block pair listed in inter_residue_connections. The full all-term matrix is still retained for total-score accounting. The optional per-block profile assigns each off-diagonal all-term pair equally to its two blocks; it is an explicit analytical convention, not a cached residue energy.

[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: (24, 1, 30, 30)
summed block-pair shape: (1, 30, 30)
matrix sum: 87.96295 whole-pose score: 87.962890625

^ⓘterm	stored_i_to_j	two_orientation_sum
fa_ljatr	-3.301456	-3.301456
fa_ljrep	0.089763	0.089763
fa_lk	3.824765	3.824765
fa_elec	-4.535557	-4.535557
hbond	-2.508230	-2.508230

^ⓘ	block	pdb_label	equal_split_all_term_score_units
18	18	A:19	-2.777947
16	16	A:17	-1.963329
21	21	A:22	-1.699825
25	25	A:26	-1.579815
2	2	A:3	-1.560023
4	4	A:5	-1.553674
20	20	A:21	-1.336094
11	11	A:12	-1.168918
6	6	A:7	-0.746334
5	5	A:6	-0.256000
(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

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

Expected observations. Summing every stored all-term matrix entry reproduces the weighted whole-pose score within numerical tolerance. Diagonal entries include one-body/self contributions; an inter-block two-body contribution appears once in one of the two off-diagonal orientations. The selected pair is ranked only after adding both orientations of the five stated nonbonded terms and removing directly covalently connected candidates. Its table makes that focused term breakdown explicit.

The equal-split profile also sums to the whole-pose score, but it is only one declared attribution rule: each off-diagonal all-term interaction contributes half to each partner while a diagonal contribution stays with its block. A naive row sum misses interactions stored in the opposite orientation, and neither profile is a cached per-residue energy table.

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 is a fragment–partner interaction score decomposition from one connected complex. It is not a binding free energy or thermodynamic ddG: no separated state, solvent correction, reorganization, packing, or minimization is evaluated. Fragment–fragment energies remain separate in the block-pair tensor and are not attributed to either fragment.

[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,
    sample_proton_chi=False,
    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 weighted fragment–protein interactions#

calculate_fragment_interactions() scores the connected pose once and sums both stored orientations between each split block and the selected protein blocks. Its SplitBlockEntry records keep each score column aligned with a pose and block index.

The table reports weighted TMol score units. Summing fragment columns recovers the connected ligand-versus-protein cross-mask score; it is still not a binding free energy.

[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.1155
F1–F4	5.4605
F1–F5	-0.2495
F1–F6	0.0000
F2–F3	0.0000
F2–F4	12.8556
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

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

The interactive pocket viewer above colors each ligand fragment inside the bound complex, and its caption reports that fragment’s weighted interaction total with the protein partner — the same one-complex score shown in the bar plot. Selecting all LG1 reports the summed whole-ligand interaction. These captions make the per-fragment attribution legible in three dimensions, but remain weighted score units from one complex, not binding free energies.

Expected observations. Different connected regions of the same bound ligand can have favorable and unfavorable weighted interactions with the selected protein blocks, and different score terms can dominate each region. The independent cross-mask assertion verifies the attribution against the underlying directed block-pair tensor.

These values answer “which fragment–protein block-pair terms contribute to this one complex?” They do not answer “how much binding free energy does this fragment provide?” Fragment–fragment entries are displayed separately and never folded into a fragment’s protein-partner column. The remainder of the notebook returns to the original 1UBQ pose_stack for differentiation and coordinate sensitivity.

Differentiate and reweight an interface#

Because whole-pose and block-pair scoring modules consume coordinate tensors, ordinary PyTorch autograd provides derivatives. Clone and detach first so the tutorial never mutates a pose’s coordinate storage.

The first calculation differentiates the ordinary 1UBQ total. The second returns to the connected ligand pose and multiplies both stored ligand↔protein orientations by 1.5 before reducing the block-pair matrix. That factor defines a deliberate analytical/training objective; it does not alter beta2016 parameters and does not turn the one-complex interaction into a binding free energy.

[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.344444
207	A	13	ILE	CB	61.870312
204	A	13	ILE	CA	51.022690
38	A	3	ILE	C	37.618317
237	A	15	LEU	CA	35.370716
236	A	15	LEU	N	33.712692
77	A	5	VAL	C	33.463181
238	A	15	LEU	C	33.092995
208	A	13	ILE	CG1	32.718456
191	A	12	THR	C	28.827141
(10 more rows not shown)

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

Perturb one nonlocal contact as a sensitivity diagnostic#

The next cell finds the closest heavy-atom pair between blocks separated by at least three sequence positions. It clones coordinates and moves one atom by only 0.10 Å along the interatomic direction. This one-sided perturbation is solely a score-sensitivity diagnostic: it is not a finite-difference estimate of a physical response and not a chemically valid conformational move. Bonded and nonbonded terms may both 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.905287	7.648530	0.743243
3	fa_elec	-41.767223	-41.446209	0.321014
4	hbond	-19.797235	-19.968584	-0.171349
1	fa_ljrep	19.423681	19.271950	-0.151731
2	fa_lk	99.107491	99.018791	-0.088699
7	cart_torsions	8.672227	8.607414	-0.064813
17	lk_ball	52.020676	51.978092	-0.042583
16	lk_ball_iso	-52.320141	-52.284916	0.035225
6	cart_angles	37.186768	37.159943	-0.026825
0	fa_ljatr	-126.778610	-126.773117	0.005493
(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.

Expected observations. A 0.10 Å change should produce finite, nonzero score deltas. The largest responding terms depend on the chosen contact and local geometry. Interpret only the local sensitivity of the score function—not chemical plausibility, relaxation behavior, or an experimentally measurable energy change.

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

Rosetta comparison#

Rosetta can switch a pose between centroid and full-atom representations and stores evaluated energies in a pose-associated Energies object. TMol’s current implementation is all-atom only. It renders a scorer for a PoseStack layout and returns tensors directly; it does not maintain Rosetta’s per-residue Energies cache.

Consequently, recompute a whole-pose or block-pair tensor when coordinates change, keep the score-function weights/options with the analysis, and do not call a block-pair matrix a cached per-residue decomposition. The Rosetta-to-TMol crosswalk collects these distinctions.

See the full official Rosetta scoring tutorial, analysis tutorial, full-atom versus centroid tutorial, PyRosetta score-function basics, and analyzing energy between residues.

Next: choose a sampling branch#

Tutorials 04 — Packing and a Small Mutation Scan and 05 — Minimization, constraints, and kinematics are parallel after this notebook: 04 changes discrete identities/conformers, while 05 changes continuous Cartesian or internal coordinates. Complete both before FastRelax composes them in tutorial 06.

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