Skip to content

Commit 4249a11

Browse files
perf: use AMBER99SB + OBC (native GBSAOBCForce) for relaxation
Revert the relaxation force field from ff14SB + GBn2 back to AMBER99SB + amber99_obc. Every amber14/ff14SB implicit-solvent model (gbn2/obc2/obc1) falls back to OpenMM's generic CustomGBForce, which benchmarks several-fold slower than the native, hand-tuned GBSAOBCForce that amber99_obc.xml uses (~4x on a 26-residue test). GBn2's accuracy edge over OBC was marginal for this pipeline, so the speedup is the better trade for high-throughput runs. Keeps FASPR packing, OBC implicit solvent (no vacuum collapse), the two-stage minimisation and PROPKA unchanged.
1 parent 10817c2 commit 4249a11

2 files changed

Lines changed: 15 additions & 13 deletions

File tree

CLAUDE.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ The entire application lives in a single file: `esmfold.py`. It is a CLI script
3232

3333
1. **`compute_features()`** — runs ESM3 forward pass to compute per-residue log-likelihoods and mean-pooled embeddings.
3434
2. **`predict_structure()`** — runs ESM3 iterative structure generation (diffusion-style, configurable steps/temperature/schedule/strategy).
35-
3. **`relax_structure()`** (optional, `--relax`) — two-stage OpenMM energy minimisation using the ff14SB force field with GBn2 generalised-Born implicit solvent. Since ESM3 emits backbone only, sidechains are rebuilt before minimisation: when the **FASPR** binary is available (default), `pack_sidechains()` repacks them with a combinatorial backbone-dependent rotamer search; otherwise it falls back to `optimize_sidechains()` (a local GBn2 minimisation). Either way PROPKA then assigns protonation states from the packed geometry. Stage 1 uses strong Cα restraints to settle sidechains/hydrogens; stage 2 uses weak restraints to allow limited backbone movement. Implicit solvent (rather than vacuum) keeps surface charges/salt bridges from collapsing inward.
35+
3. **`relax_structure()`** (optional, `--relax`) — two-stage OpenMM energy minimisation using the AMBER99SB force field with OBC generalised-Born implicit solvent (`amber99_obc.xml`, OpenMM's native fast `GBSAOBCForce`; the newer ff14SB/GBn2 implicit models fall back to a several-fold slower `CustomGBForce`). Since ESM3 emits backbone only, sidechains are rebuilt before minimisation: when the **FASPR** binary is available (default), `pack_sidechains()` repacks them with a combinatorial backbone-dependent rotamer search; otherwise it falls back to `optimize_sidechains()` (a local GB-solvent minimisation). Either way PROPKA then assigns protonation states from the packed geometry. Stage 1 uses strong Cα restraints to settle sidechains/hydrogens; stage 2 uses weak restraints to allow limited backbone movement. Implicit solvent (rather than vacuum) keeps surface charges/salt bridges from collapsing inward.
3636
- FASPR is optional with graceful fallback. Disable with `--no-faspr`; point at a custom binary with `--faspr-bin PATH` (default: `$FASPR_BIN`, else `FASPR` on `PATH`). The Docker image builds and bundles FASPR automatically. For local runs, build it once: `git clone --depth 1 https://github.com/tommyhuangthu/FASPR && g++ -O3 --fast-math -o FASPR/FASPR FASPR/src/*.cpp`, then set `FASPR_BIN=$PWD/FASPR/FASPR` (the `dun2010bbdep.bin` rotamer library must sit next to the binary).
3737

3838
**Key design decisions:**

esmfold.py

Lines changed: 14 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ class RelaxResult(TypedDict):
6767
class OpenMMContext(NamedTuple):
6868
"""OpenMM/pdbfixer modules and objects, imported and built once for the run."""
6969

70-
forcefield: Any # openmm.app.ForceField (ff14SB + GBn2 implicit solvent)
70+
forcefield: Any # openmm.app.ForceField (AMBER99SB + OBC implicit solvent)
7171
platform: Any # openmm.Platform | None
7272
mm: Any # the openmm module
7373
app: Any # the openmm.app module
@@ -238,18 +238,20 @@ def load_openmm() -> OpenMMContext:
238238
each relax_structure() call. Everything is imported lazily (not at module top) so
239239
the cost is only paid when --relax is used.
240240
241-
The force field is ff14SB with the GBn2 generalised-Born implicit solvent. Both the
242-
sidechain pre-optimisation and the production minimisation run in GB solvent:
241+
The force field is AMBER99SB with the OBC generalised-Born implicit solvent. Both
242+
the sidechain pre-optimisation and the production minimisation run in GB solvent:
243243
minimising ESM3's rebuilt sidechains in vacuum over-packs and buries polar groups,
244244
collapsing surface charges/salt bridges and wrecking the pKa prediction, whereas GB
245-
solvent penalises burial so the geometry (and PROPKA input) stays realistic.
245+
solvent penalises burial so the geometry (and PROPKA input) stays realistic. OBC via
246+
``amber99_obc.xml`` uses OpenMM's native (fast) GBSAOBCForce; the newer ff14SB/GBn2
247+
implicit models fall back to a generic CustomGBForce that is several-fold slower.
246248
"""
247249
import openmm
248250
import openmm.app
249251
import openmm.unit
250252
import pdbfixer
251253

252-
forcefield = openmm.app.ForceField("amber14/protein.ff14SB.xml", "implicit/gbn2.xml")
254+
forcefield = openmm.app.ForceField("amber99sb.xml", "amber99_obc.xml")
253255

254256
# Prefer CUDA → OpenCL → CPU
255257
platform = None
@@ -334,7 +336,7 @@ def build_variants(topology, pka, ph: float, disulfide_residues: set) -> list:
334336
335337
Returns one entry per ``topology.residues()`` (positional — Modeller maps the
336338
list by index): a variant name, or None to use Modeller's pH default. Tyr/Arg
337-
have no neutral/deprotonated template in ff14SB and are left to default.
339+
have no neutral/deprotonated template in amber99sb and are left to default.
338340
"""
339341
variants = []
340342
for res in topology.residues():
@@ -419,7 +421,7 @@ def optimize_sidechains(
419421
sidechains first: it protonates at ``ph`` (only disulphide cysteines are pinned to
420422
CYX, since their SG–SG bonds are already in the topology and would otherwise clash
421423
with an added HG), restrains all backbone heavy atoms (N, CA, C, O) with a stiff
422-
harmonic force so only the sidechains and hydrogens move, and minimises **in GBn2
424+
harmonic force so only the sidechains and hydrogens move, and minimises **in OBC
423425
implicit solvent** — a vacuum minimisation instead over-packs sidechains and buries
424426
ionizable groups, which drives PROPKA to spurious (often non-titratable) pKa.
425427
@@ -492,7 +494,7 @@ def relax_structure(
492494
use_faspr: bool = True,
493495
faspr_bin: str | None = None,
494496
) -> RelaxResult:
495-
"""Two-stage restrained energy minimisation with the ff14SB force field in GBn2
497+
"""Two-stage restrained energy minimisation with the AMBER99SB force field in OBC
496498
implicit solvent.
497499
498500
Uses pdbfixer to reconstruct missing heavy atoms (O, sidechains) and
@@ -502,13 +504,13 @@ def relax_structure(
502504
rotamers. When ``use_faspr`` is set and the FASPR binary is available, sidechains
503505
are repacked with FASPR (a combinatorial backbone-dependent rotamer search) so
504506
the pKa prediction and final structure rest on realistic packing. If FASPR is not
505-
available (or fails), it falls back to optimize_sidechains() — a local GBn2
506-
minimisation — which de-clashes but cannot cross rotamer barriers.
507+
available (or fails), it falls back to optimize_sidechains() — a local implicit-
508+
solvent minimisation — which de-clashes but cannot cross rotamer barriers.
507509
508510
Stage 1 applies strong Cα positional restraints (stage1_k kcal/mol/Ų) so
509511
sidechains and hydrogens can relax without disturbing the backbone. Stage 2
510512
applies weak restraints (stage2_k kcal/mol/Ų) to allow limited backbone
511-
movement. max_iterations=0 runs each stage until convergence. The GBn2 implicit
513+
movement. max_iterations=0 runs each stage until convergence. The OBC implicit
512514
solvent keeps surface charges and salt bridges from collapsing inward as they
513515
would in vacuum.
514516
@@ -762,7 +764,7 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
762764
parser.add_argument(
763765
"--no-faspr",
764766
action="store_true",
765-
help="Disable FASPR sidechain repacking; use local GBn2 minimisation instead.",
767+
help="Disable FASPR sidechain repacking; use local GB-solvent minimisation instead.",
766768
)
767769
parser.add_argument(
768770
"--faspr-bin",

0 commit comments

Comments
 (0)