PhysDock

Physics-aware diffusion co-folding and conformational sampling for protein–ligand interactions, with a KRAS G12C case study.

PhysDock runs three engines over the same complex and compares what they say. DiffDock-L samples ligand poses in a rigid receptor. Boltz-2 co-folds protein and ligand together and emits a predicted affinity. OpenMM relaxes the resulting geometry under a molecular-mechanics force field and reports how far it had to move things. Everything is scored against deposited crystal structures and, where labels exist, against experimental pChEMBL values.

The pipeline hard-fails invalid chemistry before any GPU work starts, measures cross-engine agreement rather than asserting it (stage 06b), screens every reported pose through PoseBusters, and attaches bootstrap confidence intervals to every rate and correlation.

Read SCOPE_AND_LIMITATIONS.md before interpreting any number below. The KRAS G12C set is a positive control on structures both models saw in training, not evidence of generalisation, and the physics stage has been corrected since the run that produced these numbers — see CHANGELOG.md.

Runs end-to-end on a single 24 GB GPU (tested on an AWS A10G). The analytical core is CPU-only and runs anywhere.

Architecture

Tier	Engine	What it contributes
A — geometry	DiffDock-L	Diffusion sampling of ligand poses in a rigid receptor; no docking box, no scoring-function search
B — induced fit and affinity	Boltz-2	AlphaFold3-class joint co-folding from sequence + SMILES; the pocket can move; auxiliary affinity head
C — physics	OpenMM	Restrained relaxation under Amber14/GAFF2; reports pose drift and an interaction-energy proxy

Stages

Stage	Script	Module	Does
00	00_setup_check.py	—	Dependency check plus a CPU-only smoke test
01	01_prepare_target.py	receptor.py	Fetches PDBs, cleans receptor, extracts reference poses and SMILES
02	02_chem_gate.py	chem.py	Hard-fails invalid valence; flags PAINS / SA / property advisories
03	03_run_diffdock.py	docking_diffdock.py	Diffusion docking
04	04_run_boltz.py	cofold_boltz.py	Co-folding and affinity prediction
05	05_physics_rescore.py	physics_openmm.py	Relaxation, pose drift, interaction energy
06	06_ensemble_analysis.py	ensemble.py	In-frame pose scatter and clustering across all poses. Static scatter, not an MD surrogate
06b	06b_consensus_validity.py	consensus.py, validity.py	Frame-aligned DiffDock↔Boltz agreement; PoseBusters battery
07	07_evaluate_and_report.py	evaluate.py, metrics.py, provenance.py	RMSD vs crystal, top-k/oracle, Spearman vs experiment with CIs, run manifest

Results — KRAS G12C

These numbers come from the last clean end-to-end run on a single 24 GB A10G, executed before the corrections in CHANGELOG.md. The pose RMSDs are unaffected by those corrections — the docking stage did not change. The interaction energies are affected: they were computed in vacuum and are superseded by the implicit-solvent decomposition that is now the default. Those values are pending a re-run and are not reproduced here as if current.

All four ligands are PDB-deposited complexes inside the training corpora of both generative engines. This is a positive control, not a generalisation result.

The pipeline was run end-to-end on four covalent KRAS G12C inhibitors — sotorasib, adagrasib, ARS-1620, ARS-853 — against their deposited co-crystals (PDB 6OIM, 6UT0, 5V9U, 5F2E).

Pose accuracy

Symmetry-aware heavy-atom RMSD between the top-ranked generative pose and the crystal pose, computed with spyrmsd:

Ligand	RMSD (Å)	Outcome
sotorasib	1.19	docking success (sub-2 Å)
adagrasib	2.92	near-native, correct pocket
ARS-1620	4.79	partially displaced
ARS-853	7.37	misplaced

sotorasib (1.19 Å)	adagrasib (2.92 Å)

Crystal pose in grey, predicted pose in teal.

The gradient across the four ligands is a characterisation, not a failure. DiffDock-L docks non-covalently and applies no constraint tethering the electrophilic warhead to Cys12, the residue all four of these inhibitors modify. What the run maps is the applicability domain: near-native for part of the set, progressively displaced for the larger and more flexible binders. Run-to-run RMSDs move by a few tenths of an Å — DiffDock-L is stochastic — but the qualitative spread reproduces. Covalent modelling is open work; see the roadmap.

Supporting signals

• Co-folding. Boltz-2 reports high interface confidence for all four complexes (ipTM 0.98).
• Predicted affinity. Reported but underpowered: n = 4 with experimental labels unfilled, so no Spearman ρ is claimed. Boltz-2 ranks adagrasib strongest and sotorasib weakest. The latter is a miss on a potent clinical drug; it is recorded rather than dropped.
• Physics. All four poses gave interaction-energy proxies in the range −278 to −476 kJ/mol with post-minimisation drift of 1.0–1.5 Å, so the generative poses are stable under a molecular-mechanics force field. This is a relaxation and plausibility proxy, not a binding free energy. These values were computed in vacuum and are superseded.
• Conformational consistency. Where more than one pose cleared the confidence gate (ARS-1620), the survivors fell into a single binding-mode cluster (spread ≈ 1.6 Å) — internally self-consistent rather than scattering hypotheses.

Not claimed: no wet-lab validation, all signals are in silico; the OpenMM term is a proxy and not MM-GBSA or FEP free energy; the affinity correlation is underpowered until experimental labels are added; the reported confidences are model self-estimates.

Design notes

Cross-engine agreement is measured, not assumed. Two independent generative architectures produce two independent placements. Stage 06b superposes the Boltz complex onto the prepared receptor by Cα Kabsch alignment, applies the same transform to the Boltz ligand, and only then measures a symmetry-aware ligand–ligand RMSD. Without the frame alignment the two poses do not live on the same grid and the comparison is meaningless. Disagreement is flagged rather than averaged away.

Pose scatter is not dynamics. Multi-seed diffusion output is cheap and it is tempting to call it a conformational ensemble. It is a static scatter of a generative model's output distribution, which measures model uncertainty, not thermal fluctuation. Stage 06 reports it as scatter and the docs say so in both places.

Physics is inspectable. Internal strain, steric clash counts and relaxation drift are separate reported columns, not folded into one opaque score.

Success is defined against external ground truth. Geometry against X-ray crystallography (RMSD ≤ 2.0 Å), function against wet-lab assays (Spearman rank correlation). Underpowered statistics are flagged as underpowered rather than reported as findings.

Quickstart

The analytical core (RDKit, pandas, Biopython) runs on any CPU machine for data prep and reporting; the ML and physics engines need a GPU box. DiffDock-L and the PhysDock/Boltz stack require mutually incompatible PyTorch builds, so they live in two conda environments. The main environment is captured in env_physdock.yml. DiffDock-L's install is order-dependent — PyTorch from the CUDA-12.1 index, then matched PyG wheels, then ProDy --no-deps — so it is shipped as setup_diffdock.sh rather than a flat YAML that cannot rebuild.

# analytical core (CPU, any machine)
pip install -e .
python scripts/00_setup_check.py

# main environment (physics + Boltz)
conda env create -f env_physdock.yml
conda activate physdock
pip install torch==2.6.0            # separate, to match the box CUDA

# DiffDock-L, separate environment
conda create -y -n diffdock python=3.9
bash setup_diffdock.sh
git clone https://github.com/gcorso/DiffDock ~/DiffDock
export PHYSDOCK_DIFFDOCK_DIR=~/DiffDock

Two environment requirements that only surface at runtime:

• Boltz on PyTorch ≥ 2.6 needs TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD=1 in the physdock environment. The official checkpoints pickle OmegaConf objects that the safe unpickler rejects. Set it once, for example via a sitecustomize.py in the env.
• The config disables the cuEquivariance triangle kernels (boltz.no_kernels: true) and pins DiffDock to --batch_size 1. Both avoid hard crashes on current CUDA/torch builds. Leave them as shipped unless the exact legacy stack is pinned.

# CPU stages
python scripts/01_prepare_target.py --config configs/kras_g12c.yaml
python scripts/02_chem_gate.py       --config configs/kras_g12c.yaml
# GPU stages
python scripts/03_run_diffdock.py    --config configs/kras_g12c.yaml
python scripts/04_run_boltz.py       --config configs/kras_g12c.yaml --max-ligands 4
# CPU / mixed
python scripts/05_physics_rescore.py --config configs/kras_g12c.yaml --top-per-ligand 3
python scripts/06_ensemble_analysis.py --config configs/kras_g12c.yaml
python scripts/06b_consensus_validity.py --config configs/kras_g12c.yaml
python scripts/07_evaluate_and_report.py --config configs/kras_g12c.yaml

# or the whole thing
bash scripts/run_all.sh configs/kras_g12c.yaml

Outputs land in results/report/report.md alongside CSV ledgers and figures.

Data integrity

data/ligands/kras_g12c_ligands.csv maps each ligand to its PDB co-crystal identifier — for example sotorasib to 6OIM. Before a full run:

1. Confirm the pdb_id / ligand_resname pairing on RCSB and set verify=1 in the manifest.
2. Populate the pchembl column with scripts/fetch_chembl_affinities.py, otherwise the functional correlation has nothing to correlate against.

Stage 01 extracts SMILES and coordinates directly from the deposited structures rather than from hand-typed strings. An incorrect PDB pairing fails loudly by design.

Layout

PhysDock/
├── configs/kras_g12c.yaml              # run parameters and thresholds
├── data/ligands/kras_g12c_ligands.csv  # co-crystal and affinity manifest
├── physdock/
│   ├── chem.py                         # cheminformatics triage (valence, SA, PAINS)
│   ├── receptor.py                     # PDB parsing, reference-pose extraction
│   ├── docking_diffdock.py             # tier A wrapper and output parser
│   ├── cofold_boltz.py                 # tier B wrapper and output parser
│   ├── physics_openmm.py               # relaxation, strain, interaction energy
│   ├── ensemble.py                     # in-frame pose scatter and clustering
│   ├── consensus.py                    # frame-aligned DiffDock↔Boltz agreement
│   ├── validity.py                     # PoseBusters battery, covalent-aware
│   ├── metrics.py                      # top-k/oracle, bootstrap CIs, calibration
│   ├── provenance.py                   # input hashing, versions, seed → run manifest
│   └── evaluate.py / report.py         # correlation and report generation
├── scripts/
│   ├── 00_setup_check.py … 07_evaluate_and_report.py
│   ├── run_all.sh
│   └── fetch_chembl_affinities.py
├── env_physdock.yml
├── setup_diffdock.sh
├── results/report/figures/
├── notebooks/01_quickstart.ipynb
├── pyproject.toml
└── requirements.txt

Open work

Ordered by what most limits the scientific claim.

1. Generalisation test. A temporal hold-out (complexes deposited after both engines' training cutoffs) or a scaffold split, plus a Gnina/Vina baseline. Until this exists, no accuracy claim here generalises — the current set is a positive control on memorised structures.
2. Covalent modelling. Explicitly model the Cys12–warhead adduct. This is the largest physical approximation in the project and a plausible driver of the RMSD gradient.
3. Affinity that means something. MM-GBSA over short restrained trajectories (scaffolded, off by default), then FEP on a shortlist. The current single-point energy is an enthalpic proxy, not ΔG.
4. Cofactor-aware physics. Parameterise the retained GDP/Mg²⁺ into the force field so the switch-II pocket is scored in its real nucleotide state.
5. Induced fit. Backbone-only restraints or restraint annealing, so pocket side chains can respond instead of being frozen at 1000 kJ/mol/nm².
6. Real dynamics. A learned conformational generator as an actual Boltzmann-ensemble surrogate, replacing pose scatter.
7. Reproducibility closure. Snapshot and version-pin the MSA so use_msa_server stops being a nondeterminism hole.
8. Closed-loop active learning. Only meaningful once 1–3 land.

License

MIT. See LICENSE.