Positive force-constant constraint for TS force fields: status & rationale

[ericchansen]

Background

Q2MM builds transition-state force fields by inverting the QM Hessian's
negative reaction-coordinate eigenvalue before Seminario projection
(Limé & Norrby 2015 method C, invert_ts_curvature=True), yielding
all-positive starting force constants where the TS geometry is a proper
minimum on the MM surface. But L-BFGS-B does not preserve non-negativity the
way MacroModel's original pipeline did, so bond/angle force constants could
drift negative during optimization — producing unphysical saddle-shaped
harmonic terms in the final force field.

What's already in place (as of #296)

Two mechanisms now guard against negative force constants:

1. Non-negative default bounds — ForceField.DEFAULT_BOUNDS
   (q2mm/models/forcefield.py):

   • bond_k: (0, 3600)
   • angle_k: (0, 720)
   • ub_k: (0, 500)

   torsion_k and sb_k remain signed (they are Fourier coefficients and
   legitimately span both signs).

2. Method E2 Approxn near-zero replacement (q2mm/workflows/method_e2.py):
   when Round 1 plus the non-negative bounds parks a force constant at exactly
   0.0 (a Hooke's-law spring with no restoring force), the workflow
   substitutes the Farrugia 2025 Q2MM Approxn standard before freezing for
   Round 2 — bond_k/ub_k = 5 mdyn/Å, angle_k = 0.5 mdyn·Å/rad²
   (APPROXN_DEFAULTS).

Decision: no strict epsilon floor

We considered lowering the bounds from 0.0 to a small positive epsilon to
enforce a strictly positive floor. Decided against it. An epsilon floor
risks pinning parameters at the boundary and fighting the optimizer, whereas
the (0, max) bounds + Approxn replacement already prevent negative force
constants and handle the zero-drift case with a physically motivated value.
This issue records that decision so it is not silently revisited.

Remaining / monitoring

• Confirm across all 6 systems (via the Method E2 convergence re-run) that
  no negative force constants survive into the final force fields.
• Revisit only if a system is found where the Approxn replacement is
  insufficient.

References

• feat(phase-9e): non-negative FC bounds, Approxn replacement, shared benchmark_runner #296 — non-negative bounds + Approxn near-zero replacement (merged)
• validation: published FF should beat Seminario baseline (Check 1 gate) #255 — backend parity / MM3* engine gap
• MM3 stretch-bend: OpenMM/Tinker implementation + wildcard atom type matching #261 — MM3 stretch-bend cross terms
• Limé & Norrby, J. Chem. Theory Comput. 2015, 11, 3696
• Farrugia, Helquist, Norrby & Wiest, J. Chem. Theory Comput. 2025, 22, 469 (QFUERZA / Approxn)