fix(qmmm): make embedding=mechanical runnable (zero embedding field)

[karmachoi]

The bug

[qmmm] embedding = mechanical aborts on the first SCF iteration:

qmmm_espf::add_potqm_contributions: Record `OQP::POTQM` not found!

scf.F90:1207 calls add_potqm_contributions on every SCF iteration whenever control%qmmm_flag is set, and that routine requires OQP::SM, OQP::ESPF_CORR and OQP::POTQM with WITH_ABORT. But compute_force only builds the embedding arrays for the electrostatic/split/ESPF paths:

potmm = potqm = None
if self.Embedding in ("electrostatic", "split") or self.espf_full:
    potmm, potqm = self.electrostatic_potential()

so mechanical leaves potmm = potqm = None, the OQP::POTQM record is never created, and the SCF dies. The unconditional oqp.grad_esp_qmmm() call needs OQP::POTMM for the same reason.

This is the same family as the known single-point QM/MM gap (runtype=energy + qmmm_flag=true aborting on a missing OQP::ESPF_CORR): scf.F90 assumes qmmm_flag implies the ESPF records exist, which only the ESPF path guarantees.

The fix

Give mechanical embedding a zero field rather than an absent one. That is precisely what mechanical embedding means for the QM subsystem: the QM Hamiltonian is unperturbed by the MM charges (a gas-phase SCF), and QM–MM electrostatics is left to OpenMM via the QM ESP charges in forces_mm().

A zero field is exact, not an approximation — it adds 0 to hcore, contributes 0 to the ESPF gradient, and leaves eqm untouched, so the existing if potmm is not None guards remain numerically no-ops.

Note potqm is already zeroed unconditionally in forces_qm_openqp (the #205 M1a fix), so add_potqm_contributions is a no-op in every code path today — the record only has to exist. A Fortran-side early return on absent records would be an equally valid fix and would avoid the (small) wasted work of building the ESPF grid; happy to switch if maintainers prefer that. espf_op_corr is needed regardless, since form_esp_charges supplies the ESP charges the mechanical path hands to OpenMM.

Verification

Test system: methanol with QM = OH, MM = CH3. The partition cuts the C–O bond, so one hydrogen link atom caps the O (g = 0.7152). HF/6-31G, NVE + Verlet, dt = 0.25 fs.

• Before: aborts on the first SCF iteration.
• After: completes 400 NVE steps.

Scope caveat — this makes mechanical run, not exact

forces_mm() writes geometry-dependent QM ESP charges into OpenMM, whose force is evaluated at fixed charges, so the charge-response term Σ_A (∂E_mm/∂q_A)(dq_A/dR) is never computed. Finite differences on the methanol system:

quantity	O_x	HO_x	C_x	H1_x
-gmm (OpenMM)	91.78	11.06	−88.48	161.50
FD, frozen charges	91.76	11.06	−88.46	161.50
FD, charges follow geometry	−243.97	282.86	−24.51	161.50

OpenMM reproduces the frozen-charge derivative to 0.02%, but the true derivative differs by up to ~336 kJ/mol/nm. The methyl hydrogens, which carry no QM-dependent charge, are exact — as expected.

This is a pre-existing property of the split/mechanical design, independent of link atoms, and is left untouched here. Mechanical embedding should not be treated as gradient-exact for dynamics on the strength of this PR; the ESPF path avoids the issue by zeroing the QM charges in OpenMM and handling charge fluctuation through grad_esp_qmmm.

🤖 Generated with Claude Code

[chatgpt-codex-connector]

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 2bc31d0bad

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

• Open a pull request for review
• Mark a draft as ready
• Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

[chatgpt-codex-connector]

Preserve DFTB mechanical embedding path

When method=dftb is used with Embedding='mechanical', this new branch now passes a non-None zero potmm into forces_qm_openqp(), which immediately dispatches to _forces_qm_dftb(). That helper explicitly treats mechanical embedding as potmm is None and raises NotImplementedError for any non-None potential when self.espf_full is false, so DFTB QM/MM mechanical jobs that the input checker still advertises as supported will now fail before the gradient call instead of running gas-phase QM plus MM coupling.

Useful? React with 👍 / 👎.

[chatgpt-codex-connector]

Reject unknown embeddings instead of treating them as mechanical

This fallback runs for every value that is not recognized by the electrostatic/split check, not just for Embedding='mechanical'. Because the QMMM config accepts user-provided embedding strings, a misspelling such as embedding=electrostatc or an unnormalized legacy value now silently runs an unembedded mechanical QM calculation rather than failing, which can produce plausible but physically different trajectories; guard this branch with an explicit mechanical check and raise for unknown values.

Useful? React with 👍 / 👎.