diff --git a/docs/qmmm-tb-reconciliation.md b/docs/qmmm-tb-reconciliation.md new file mode 100644 index 000000000..1ed4120db --- /dev/null +++ b/docs/qmmm-tb-reconciliation.md @@ -0,0 +1,193 @@ +# QM/MM ↔ tight-binding reconciliation (full-ESPF covalent boundary for `method=dftb` and `method=xtb`) + +This note records the merge of `origin/codex/soc-namd-options` (the correct +full-ESPF covalent-boundary QM/MM) into the TB-adapter branch +(`feat/openqp-xtb-adapter`: DFTB adapter PR #266 + the `method=xtb` adapter + the +`oqp.utils.tb_backends` dispatch generalization). The goal is a single tree where +the full-ESPF QM/MM driver internals coexist with the TB dispatch so that +`method=dftb` **and** `method=xtb` QM/MM both conserve energy at covalent +(hydrogen link-atom) boundaries. + +Merge-base: `55befcb8f` (PR #253). `soc-namd-options` was 83 commits ahead in +QM/MM; the TB branch was 13 commits ahead. The merge was performed +`soc-namd-options → feat/xtb-qmmm-fullespf`. + +## What full-ESPF brings + +The released ESPF path embedded covalent boundaries with **frontier-charge +redistribution** (`frontier_scheme = none|rcd|rc|z1`, PR #258): the M1 host +charge could be deleted/redistributed onto virtual midpoint charges. `soc-namd- +options` replaced that with a cleaner, gradient-exact **full-field ESPF** scheme +and removed the redistribution machinery entirely: + +- `OpenQpQMMM.__init__` builds a per-QM-centre MM-exclusion structure + (`_build_qm_mm_exclusions`) instead of a frontier-host redistribution map. The + embedding potential (`_full_field_potmm`) and the coupling force + (`_coupling_forces`) use the *same* per-centre view of the MM charges, so the + analytic gradient stays exact across a covalent cut. +- The full-ESPF native path adds the **nuclear-MM interaction** `Σ_A Z_A φ_A` + to complete the QM–MM electrostatic energy (commit `fac2be89c`, "nuclear-MM + energy → FD-exact covalent-boundary forces"), and drops the frontier + redistribution / virtual-site scatter. +- `espf_full` is now selected by `Embedding ∈ {espf, espf_full, electrostatic}` + and is the default. +- PBC / min-image support (`_box_lengths_bohr`, `_min_image`) is retained. + +The `redistribute_frontier_charges` / `assemble_embedding_sites` helpers, +`_frontier_hosts` / `_embedding_sites`, the `frontier_scheme` schema key and API +parameter, and the frontier tests/examples were all removed on `soc-namd-options`; +this merge takes that removal. + +## TB dispatch contract (unchanged by full-ESPF) + +For `method=dftb`/`xtb` the openqp-dftb/openqp-xtb library folds the per-QM-centre +MM potential (Hartree/e) **directly into the SCC Hamiltonian** via +`mol.dftb_external_potential = POTMM`: + +- the returned state energy is the **complete** embedded QM energy (it already + contains `E_ext = Σ_A q_A φ_A` with `q_A` the *net* atomic charge), so the + driver must **not** add the native nuclear `Σ_A Z_A φ_A` term for TB; +- the returned analytic gradient is `d(E_embedded)/dR_QM` at fixed potential (the + Pulay/charge-response coupling is already inside it), so the native + `espf_op_corr` hcore mutation and `grad_esp_qmmm(_excited)` / `OQP::ESPF_GRAD` + additions are **skipped**; +- the adapter publishes relaxed net atomic charges of the active state into + `OQP::partial_charges`; the driver's backend-agnostic `_coupling_forces` reads + those for the classical `dφ/dR` MM-field forces (the DFTB analog of + `form_esp_charges`). + +Only full-ESPF (`espf_full`) or mechanical (`potmm is None`) embedding is +supported for TB; the legacy `split` scheme is rejected (it would double-count +the coupling already inside the embedded TB energy). + +## TB-dispatch call sites re-applied on the full-ESPF driver + +Dispatch helpers live in `pyoqp/oqp/utils/tb_backends.py`: +`is_tb_method`, `tb_section_name`, `tb_config`, `make_tb_adapter`. + +`pyoqp/oqp/library/qmmm_driver.py` (`OpenQpQMMM`): + +- import `is_tb_method`; +- `_forces_qm`, **mol mode** — after `_update_mol_positions()`, before the native + `SinglePoint`/`espf_op_corr` path: `if is_tb_method(...): return + self._forces_qm_dftb(self.mol, potmm)`; +- `_forces_qm`, **config mode** — after `self.op = OPENQP(...)`, same early + dispatch on `self.op.mol`; +- new method `_forces_qm_dftb(mol, potmm)`: sets `mol.dftb_external_potential`, + runs `Gradient(mol).gradient()` (dispatches to the TB adapter), reads the active + state energy/gradient and `OQP::partial_charges`; guards `espf_full` (else + `NotImplementedError`). + +`pyoqp/oqp/library/namd.py`: + +- import `is_tb_method, make_tb_adapter, tb_section_name`; +- `NAMD_QMMM._electronic_structure` — TB branch sets `mol.dftb_external_potential` + and runs the embedded SCF/excitation without the ESPF hcore mutation + (covalent-boundary electronic step); +- `NAMD_QMMM._forces_qm` — TB early-return of the raw gradient (skips + `grad_esp_qmmm`/`OQP::ESPF_GRAD`); +- `NAMD_QMMM._total_force` — TB energy uses `mol.energies[active]` with **no** + nuclear `Σ Z_A φ_A` term; +- `NAMD_SOC`, `NAMD_SOC_MCH`, `NAMD_SOC_QMMM`, `NAMD_SOC_MCH_QMMM` — TB SOC arms: + `tb_section_name(...)['target_multiplicity']` bookkeeping, `_dftb_soc_tags` + (one-centre SOC + numpy `eigh`) in place of `oqp.soc_mrsf`, and + `_dftb_spatial_overlap` in place of `oqp.get_states_overlap`; +- module helpers `_dftb_soc_tags`, `_dftb_spatial_overlap` (TB, dftb+xtb via + `make_tb_adapter`). + +`pyoqp/oqp/library/runfunc.py` — `compute_soc` TB branch (`xtb_soc`/`dftb_soc`); +`compute_namd` (soc's QM/MM NAMD dispatch) retained. + +`pyoqp/oqp/library/single_point.py` — TB dispatch at `energy`, `reference`, +`excitation`, `gradient`, and the `dftb_overlap`/`dftb_states_overlap` state- +overlap paths (all via `is_tb_method`/`make_tb_adapter`); soc's +`import oqp.utils.qmmm` and `scf_grad` tweaks retained. + +`pyoqp/oqp/utils/input_checker.py` — NAMD gate allows +`method ∈ {tdhf, dftb, xtb}`; the `method=dftb/xtb` QM/MM validator (rejects +`split`, requires full-ESPF electrostatic, gates SOC-QMMM off as "not yet wired") +retained from the TB branch and now describes the full-ESPF scheme. + +`pyoqp/oqp/molecule/oqpdata.py` — `[dftb]`/`[xtb]` schema sections retained; +`frontier_scheme` removed. + +`pyoqp/oqp/openqp.py` — `[dftb]`/`[xtb]` builder API retained; `frontier_scheme` +argument removed from `qmmm(...)`. + +`pyoqp/oqp/utils/oqp_tester.py` — openqp-dftb/openqp-xtb "backend not found → +SKIPPED" handling and the `method=dftb`/`method=xtb` example-skip predicate +retained. + +## Conflict resolution per file + +| File | Type | Resolution | +| --- | --- | --- | +| `pyoqp/oqp/library/namd.py` | add/add | Take soc (full-ESPF), re-apply the TB delta (`674aab601..HEAD`), which applied cleanly (soc namd.py ≡ released namd.py apart from the 2-line `frontier_scheme` removal). | +| `pyoqp/oqp/library/qmmm_driver.py` | add/add | Take soc (full-ESPF); re-apply the TB delta — 3/4 hunks clean, the `tb_backends` import added manually (soc had edited the adjacent `qmmm_connectivity` import). | +| `pyoqp/oqp/library/qmmm_connectivity.py` | add/add | Take soc (frontier helpers removed). TB delta never touched this file. | +| `pyoqp/oqp/library/qmmm_md.py` | add/add | Take soc (removes `frontier_scheme` + the internal `runtype→energy` forcing). TB dispatch is in the driver, not here. | +| `pyoqp/oqp/utils/qmmm.py` | (auto) | Identical on both sides; auto-merged. | +| `pyoqp/oqp/utils/input_checker.py` | content | Take HEAD for the NAMD `method ∈ {tdhf,dftb,xtb}` gate; rest auto-merged (soc QM/MM checks + TB validator both present). | +| `pyoqp/oqp/molecule/oqpdata.py` | content | Take HEAD `timestep` (`float`, safer than soc's `int` — the MD/NAMD code casts to float and `int` truncates sub-fs steps); take soc (remove `frontier_scheme`). `[dftb]`/`[xtb]` sections auto-merged. | +| `pyoqp/oqp/openqp.py` | content | Take soc (remove `frontier_scheme` param/docstring/update). `[dftb]`/`[xtb]` API auto-merged. | +| `pyoqp/oqp/utils/oqp_tester.py` | content | Take HEAD (openqp-dftb/xtb SKIPPED handling + skip predicate). | +| `pyoqp/oqp/library/runfunc.py` | (auto) | Auto-merged: soc `compute_namd`/`QMMMOpt` + TB `compute_soc` branch both present. | +| `pyoqp/oqp/library/single_point.py` | (auto) | Auto-merged: soc `import qmmm`/`scf_grad` + TB dispatch both present. | +| `include/oqp.h`, `source/tagarray_driver.F90` | (auto) | Auto-merged; soc additions verified present. | +| `tests/test_openqp_api.py` | content | Take soc (drop `frontier_scheme` schema + `test_qmmm_frontier_scheme_sets_section_key`); TB tests auto-merged. | +| `tests/test_rohf_status_and_interface.py` | content | Take HEAD (the `tb_backends` stub bootstrap the shared `try` block depends on). | +| `README.md`, `examples/QMMM/README.md` | content/add | Take HEAD tutorial links in README; take soc for the examples README (drop the `frontier_scheme` section). | +| `tests/test_qmmm_frontier.py`, `tests/test_qmmm_frontier_openmm.py`, `examples/QMMM/ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp` | stale | **Removed** — HEAD-only artifacts of the removed `frontier_scheme` feature; they import `redistribute_frontier_charges`/`assemble_embedding_sites` (deleted on soc) and would fail to load. Absent on soc. | + +Adapter files `pyoqp/oqp/library/openqp_dftb.py`, `openqp_xtb.py`, +`pyoqp/oqp/utils/tb_backends.py` are TB-branch-only (no conflict) and kept as-is. + +## Static verification (this agent) + +- Every changed `.py` file parses (`ast.parse`). +- No conflict markers remain in the tree. +- Dependency-free unit tests (no compiled `liboqp`): the DFTB/XTB schema-hook + tests, the runfunc/interface dispatch test, and the QM/MM builder API test pass + — **54 passed, 20 skipped** (skips are OpenMM/liboqp-gated). +- Full `tests/` run vs the pre-merge HEAD baseline: **identical set of 51 + pre-existing failures** (all require the compiled `oqp` extension, which is not + built in this environment) — the merge introduces **no new failures**. The + drop in passing count equals exactly the removed frontier tests. + +## Remaining runtime-validation checklist (needs the integrated runtime) + +The compiled `_oqp`/`liboqp` extension plus the external openqp-dftb / openqp-xtb +libraries are **not** built in this agent's environment, so the numeric +covalent-boundary behaviour was not exercised. Once the integrated runtime is up: + +1. **Covalent-boundary energy-vs-force FD, `method=dftb`** — a QM/MM partition + that cuts a covalent bond (H link-atom cap), full-ESPF electrostatic + embedding, `runtype=grad`: verify `F ≈ −dE/dR` by finite differences on both + QM and frontier/host MM atoms (this is the original DFTB QM/MM bug this merge + targets). +2. **Same FD check for `method=xtb`.** +3. **Whole-molecule QM/MM regression** (no link atom) for dftb and xtb — confirm + full-ESPF still matches the pre-merge energies/forces. +4. **Ground-state QM/MM MD** (`runtype=md`, `QMMM_MD`) energy conservation for + dftb and xtb across a covalent boundary. +5. **NAMD-QMMM** (`runtype=namd`, non-SOC, `embedding=electrostatic`) for dftb — + the `NAMD_QMMM._electronic_structure`/`_forces_qm`/`_total_force` TB arms. + +### Honest flags for the runtime tests + +- **`_forces_qm_dftb` is the released method name/body applied onto soc's driver.** + Its `espf_full` guard and its reliance on `_coupling_forces` reading + `OQP::partial_charges` are consistent with soc's `_coupling_forces` (verified + by reading both), but the *numerical* interaction of the TB net-charge coupling + with soc's **new full-field exclusion structure** (`_build_qm_mm_exclusions`, + which in practice applies *no* exclusions — see its NOTE) has not been run. + Item 1/2 above is the authoritative check. +- **SOC-QMMM for TB.** `namd.py` carries TB SOC-QMMM arms + (`NAMD_SOC_QMMM`/`NAMD_SOC_MCH_QMMM`), but `input_checker` still gates + `dftb/xtb` SOC-QMMM off as "not yet wired". The code path exists but is not + reachable through the checker; do not assume it is validated. +- **`qmmm_md.py`** dropped soc's removal of the internal `runtype→energy` + forcing block. The TB QM/MM MD path routes through + `OpenQpQMMM._forces_qm` (which has the TB dispatch), so this is expected to be + fine, but a `runtype=md` dftb/xtb smoke run (item 4) should confirm the QM + sub-run accepts the config. diff --git a/examples/QMMM/README.md b/examples/QMMM/README.md index 3be13d6aa..2ead67cbc 100644 --- a/examples/QMMM/README.md +++ b/examples/QMMM/README.md @@ -39,42 +39,3 @@ contract and the compact `job.qmmm(...)` / `job.workflow.namd(...)` Python API. `ala.inp` and `2E4E_RHF-DFT-QMMM_energy.inp` are QM/MM single-point energy decks (QM selection via `[input] system = file.pdb `); `run.inp` is a ground-state OpenMM-integrator QM/MM MD deck. - -## Covalent QM/MM boundary — `[qmmm] frontier_scheme` - -When the QM/MM partition cuts a covalent bond, the dangling QM bond is capped -with a hydrogen link atom and the MM host atom (`M1`) sits ~1.5 Å from the QM -density. `[qmmm] frontier_scheme` selects how that frontier charge is treated in -the ESPF electrostatics. Covalent QM/MM boundaries are handled by the -ground-state QM/MM MD path (`QMMM_MD`); the nonadiabatic `runtype=namd` path does -not yet append link atoms to its QM molecule and raises on a covalent cut. - -| value | meaning | -| --- | --- | -| `none` (default) | Full-field: the QM density sees the complete MM charge set. This is the **validated ESPF baseline** — ESPF couples the MM potential to QM *atomic-charge operators* (`h += Σ_A φ_A Q̂_A`, Huix-Rotllant & Ferré, *JCTC* 2021, 17, 538, eq 6), which already suppresses the electron spill-out that motivates redistribution in density-based embedding, so the ESPF papers use full MM charges even at a covalent protein boundary. | -| `rcd` | Delete `M1`'s charge and redistribute it to virtual point charges at the `M1–M2` bond midpoints, conserving the **total charge and the dipole about `M1`**. Gradient-consistent (the midpoints are linear in the real atom positions). | -| `rc` | As `rcd` but conserving only the total charge. | -| `z1` | Delete `M1`'s charge (conserves neither; for comparison). | - -`rcd`/`rc`/`z1` are **optional refinements**, not the ESPF default. Enable via the -input (`[qmmm] frontier_scheme = rcd`) or the Python API -(`job.qmmm(..., frontier_scheme="rcd")`). It is a no-op for whole-molecule QM -regions (no cut bond). - -A runnable covalent-boundary deck is -`ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp` — the alanine dipeptide (ACE-ALA-NH2) with -AMBER-14, QM = the C-terminal amide so the QM/MM partition cuts the `ALA C–CA` -backbone bond, run as ground-state QM/MM MD (`runtype=md`) with -`frontier_scheme=rcd`: - -```bash -cd examples/QMMM && openqp ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp -``` - -Like the other ground-state QM/MM decks it is skipped by `openqp --run_tests all` -(covalent-boundary QM/MM is the ground-state MD path, not `runtype=namd`). The -same alanine boundary is exercised automatically — link-atom detection + -frontier-charge conservation on the real AMBER-14 charges — in -`tests/test_qmmm_frontier_openmm.py` (OpenMM-gated), and the pure redistribution -math (including a finite-difference gradient check) in -`tests/test_qmmm_frontier.py`. diff --git a/examples/QMMM/ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp b/examples/QMMM/ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp deleted file mode 100644 index 714c72a84..000000000 --- a/examples/QMMM/ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp +++ /dev/null @@ -1,37 +0,0 @@ -# Ground-state QM/MM molecular dynamics (runtype=md) across a COVALENT QM/MM -# boundary, demonstrating [qmmm] frontier_scheme. -# -# System: the alanine dipeptide (ACE-ALA-NH2, ala.pdb) with AMBER-14. QM = the -# C-terminal amide (ALA C,O + NH2), so the QM/MM partition cuts the ALA C-CA -# backbone bond: it is capped with a hydrogen link atom and the MM host-atom -# charge is redistributed by frontier_scheme=rcd (delete the M1 charge; virtual -# point charges at the M1-M2 bond midpoints conserving total charge + dipole). -# -# Covalent-boundary QM/MM is handled by this ground-state QM/MM MD path; the -# nonadiabatic surface-hopping (NAMD) path does not append link atoms and raises -# on a cut. Requires the optional OpenMM backend. Like the other ground-state -# QM/MM decks (run.inp, ala.inp) it is skipped by the full example suite; the -# frontier-charge math is regression-tested in tests/test_qmmm_frontier*.py. -[input] -runtype = md -qmmm_flag = true -method = hf -functional = bhhlyp -basis = 6-31g -charge = 0 - -[scf] -type = rhf -multiplicity = 1 - -[qmmm] -pdb_file = ala.pdb -forcefield_files = amber14-all.xml -qm_atoms = 8,9,16,17,18 -cutoff = NoCutoff -embedding = electrostatic -frontier_scheme = rcd -n_steps = 2 -timestep = 0.5 -ensemble = nve -temperature = 300.0 diff --git a/pyoqp/oqp/library/namd.py b/pyoqp/oqp/library/namd.py index cc5eb87be..a76dc61c0 100644 --- a/pyoqp/oqp/library/namd.py +++ b/pyoqp/oqp/library/namd.py @@ -392,7 +392,6 @@ def _resolve_aux(name): self.cutoff = _resolve_cutoff(str(q['cutoff']).strip()) # NoCutoff | PME | Ewald | ... self.periodic = self.cutoff is not app.NoCutoff embedding = str(q['embedding']).strip() - frontier_scheme = str(q.get('frontier_scheme', 'none')).strip() self.pdb = app.PDBFile(pdb_file) self.forcefield = app.ForceField(*ff_files) @@ -404,7 +403,6 @@ def _resolve_aux(name): mol=mol, Cutoff=self.cutoff, Embedding=embedding, - frontier_scheme=frontier_scheme, ) self.mm = self.driver.mm_systems diff --git a/pyoqp/oqp/library/qmmm_connectivity.py b/pyoqp/oqp/library/qmmm_connectivity.py index 080cae087..6743d1e81 100644 --- a/pyoqp/oqp/library/qmmm_connectivity.py +++ b/pyoqp/oqp/library/qmmm_connectivity.py @@ -173,182 +173,3 @@ def project_link_gradient(grad_link, g): """ grad_link = np.asarray(grad_link, dtype=float) return (1.0 - g) * grad_link, g * grad_link - - -# --------------------------------------------------------------------------- -# Frontier (M1) charge redistribution for electrostatic embedding -# --------------------------------------------------------------------------- -# -# A covalent QM/MM cut leaves an MM *host* atom M1 (the MM end of the severed -# bond) roughly a bond length (~1.5 A) from the QM density and from the capping -# link-atom hydrogen. Embedding M1's full force-field point charge there badly -# over-polarizes the QM density. The remedy (Lin & Truhlar, J. Phys. Chem. A -# 2005, 109, 3991) is to delete q1 from M1 and redistribute it so the embedding -# charge cloud reproduces: -# -# * the total charge of the deleted q1 (RC and RCD), and -# * its dipole about M1 (which is identically zero) (RCD only), -# -# using *virtual* point charges placed at the midpoints of the M1-M2 bonds, -# where M2 are the MM atoms bonded to M1. Let M1 have N such neighbours and let -# d_k = (R_M2k - R_M1)/2 be the displacement of midpoint k from M1: -# -# Z1 : delete q1. Conserves neither charge nor dipole. -# RC : + q1/N at each midpoint. Conserves charge; not dipole. -# RCD : + 2 q1/N at each midpoint AND Conserves charge AND dipole: -# - q1/N added to each M2. sum q = N(2 q1/N) - N(q1/N) = q1; -# dipole = sum_k [2q1/N d_k - q1/N (2 d_k)] = 0. -# -# Each virtual charge sits at a *fixed linear combination* of two real atom -# positions (a bond midpoint), so - exactly like the link-atom capping position - -# any electrostatic force on it redistributes onto M1 and M2 by the chain rule. -# This keeps the analytic QM/MM gradient exact, which naive charge *deletion* -# does not (deletion drops a chunk of the field with no matching force term). -# -# NOTE (ESPF): OpenQP's electrostatic embedding is ESPF (Huix-Rotllant & Ferre, -# JCTC 2021, 17, 538). There the MM potential phi_A couples to the QM *atomic -# charge operators* Q_A (H += sum_A phi_A Q_A, eq 6), NOT to the raw electron -# density via 1/|r-R_M| point-charge integrals. That coarse-grained coupling -# already suppresses the electron spill-out / over-polarization that plagues -# density-based embedding, which is why the ESPF papers use the FULL MM charges -# with "no scaling of the QM-MM electrostatic interactions" even at a covalent -# protein boundary (cryptochrome). So for ESPF the DEFAULT is full-field -# ('none'); RCD/RC/Z1 are OPTIONAL refinements (physically motivated - they keep -# the MM charge/dipole faithful near the cut - but a departure from the validated -# ESPF baseline, so opt-in). - - -@dataclass(frozen=True) -class VirtualCharge: - """A redistributed (virtual) embedding point charge at an M1-M2 bond - midpoint. - - Attributes - ---------- - charge : float - Point charge (elementary charge units). - hosts : tuple[int, int] - Absolute atom indices ``(m1_index, m2_index)`` whose positions define - the virtual site. - weights : tuple[float, float] - Barycentric weights so that ``r_virtual = w0*R[m1] + w1*R[m2]`` (``0.5, - 0.5`` for a bond midpoint). A force ``f`` on the virtual charge maps to - ``w0*f`` on ``m1`` and ``w1*f`` on ``m2``. - """ - - charge: float - hosts: tuple - weights: tuple - - -def redistribute_frontier_charges(frontier, charge_of, scheme="rcd"): - """Redistribute MM frontier-host (M1) charges for electrostatic embedding. - - Parameters - ---------- - frontier : iterable of (int, sequence[int]) - One ``(m1_index, [m2_index, ...])`` per *unique* MM host atom cut by the - QM/MM boundary, with its MM neighbours M2. - charge_of : callable(int) -> float - MM force-field charge (e) of an absolute atom index. - scheme : {'rcd', 'rc', 'z1', 'none'} - Redistribution scheme (see module notes above). ``'none'`` disables it - (full-field embedding, the pre-RCD behaviour). - - Returns - ------- - deleted : set[int] - MM atoms whose embedding charge is removed (the M1 hosts). - delta_q : dict[int, float] - Charge added to a *real* MM atom's embedding charge (the M2 corrections). - virtuals : list[VirtualCharge] - Virtual midpoint charges to add to the embedding set. - """ - scheme = (scheme or "none").lower() - if scheme not in ("rcd", "rc", "z1", "none"): - raise ValueError(f"Unknown frontier charge scheme: {scheme!r}") - - deleted, delta_q, virtuals = set(), {}, [] - if scheme == "none": - return deleted, delta_q, virtuals - - for m1, m2s in frontier: - m1 = int(m1) - m2s = [int(m) for m in m2s] - q1 = float(charge_of(m1)) - deleted.add(m1) - # Nothing to redistribute onto (isolated host) -> fall back to deletion. - if scheme == "z1" or not m2s: - continue - n = len(m2s) - for m2 in m2s: - if scheme == "rc": - virtuals.append(VirtualCharge(q1 / n, (m1, m2), (0.5, 0.5))) - else: # rcd - virtuals.append(VirtualCharge(2.0 * q1 / n, (m1, m2), (0.5, 0.5))) - delta_q[m2] = delta_q.get(m2, 0.0) - q1 / n - return deleted, delta_q, virtuals - - -def assemble_embedding_sites(mm_idx, mm_charges, mm_positions, - deleted, delta_q, virtuals): - """Assemble the embedding point-charge set the QM density is embedded in. - - Combines the raw MM atoms with the frontier redistribution from - :func:`redistribute_frontier_charges`. - - Parameters - ---------- - mm_idx : sequence[int] - Absolute atom indices of the real MM atoms, in array order. - mm_charges : sequence[float] - Force-field charges (e), aligned with ``mm_idx``. - mm_positions : (M, 3) array - Positions (any length unit), aligned with ``mm_idx``. - deleted, delta_q, virtuals : - Output of :func:`redistribute_frontier_charges`. - - Returns - ------- - charges : (S,) ndarray - Embedding charges. - positions : (S, 3) ndarray - Embedding positions (same unit as ``mm_positions``). - scatter : list[tuple[tuple[int, float], ...]] - For each site, the ``(absolute_atom_index, weight)`` contributions that - build its position, so a force ``f`` on site ``s`` adds ``weight*f`` to - each contributing real atom (identity for a real atom, 0.5/0.5 for a - midpoint). - """ - mm_positions = np.asarray(mm_positions, dtype=float) - pos_of = {int(mm_idx[k]): mm_positions[k] for k in range(len(mm_idx))} - - charges, positions, scatter = [], [], [] - applied = set() - for k in range(len(mm_idx)): - a = int(mm_idx[k]) - if a in deleted: - continue - charges.append(float(mm_charges[k]) + delta_q.get(a, 0.0)) - positions.append(mm_positions[k]) - scatter.append(((a, 1.0),)) - applied.add(a) - # A delta_q correction whose target atom is itself a deleted host (two - # adjacent cuts / mutually-bonded frontier hosts) would be dropped by the - # loop above. Keep it as a charge-only residual site at that atom so total - # charge and the dipole about each host stay conserved regardless of topology. - for a, dq in delta_q.items(): - a = int(a) - if a not in applied and dq != 0.0: - charges.append(dq) - positions.append(pos_of[a]) - scatter.append(((a, 1.0),)) - for v in virtuals: - m1, m2 = v.hosts - w0, w1 = v.weights - positions.append(w0 * pos_of[m1] + w1 * pos_of[m2]) - charges.append(v.charge) - scatter.append(((m1, w0), (m2, w1))) - return (np.asarray(charges, dtype=float), - np.asarray(positions, dtype=float), - scatter) diff --git a/pyoqp/oqp/library/qmmm_driver.py b/pyoqp/oqp/library/qmmm_driver.py index ea2f30ded..55ce8f32a 100644 --- a/pyoqp/oqp/library/qmmm_driver.py +++ b/pyoqp/oqp/library/qmmm_driver.py @@ -12,7 +12,6 @@ from oqp.utils.tb_backends import is_tb_method from oqp.library.qmmm_connectivity import ( detect_link_atoms, link_atom_position, - redistribute_frontier_charges, assemble_embedding_sites, ) import oqp @@ -121,7 +120,6 @@ def __init__( mol=None, Cutoff=app.NoCutoff, Embedding='mechanical', - frontier_scheme='none', ): if oqp_cfg is None and mol is None: raise ValueError("Either 'oqp_cfg' or 'mol' must be provided.") @@ -133,17 +131,7 @@ def __init__( self.positions = positions self.topology = topology self.forcefield = forcefield - # Sort the QM selection into ascending (topology) order. The QM geometry - # handed to the engine is built by iterating self.topology.atoms() filtered - # by membership (i.e. topology order), so gqm / f_qm / pchg come back in - # topology order. The force-scatter loops in _assemble_force / _assemble_ - # force_espf and the link-atom host_row projection index those arrays by - # position in self.qm_atoms; if qm_atoms were given out of order (e.g. - # "5,2,7") those positions would not match topology order and QM gradients, - # coupling forces, and link-atom projections would land on the wrong atoms. - # Sorting makes input order == topology order without changing the QM - # calculation (the engine already sees the atoms in topology order). - self.qm_atoms = np.array(sorted(int(i) for i in qm_atoms), dtype=int) + self.qm_atoms = np.array(qm_atoms, dtype=int) self.Cutoff = Cutoff self.Embedding = Embedding @@ -174,24 +162,17 @@ def __init__( # driver. self.link_atoms = self._detect_link_atoms() - # Frontier (M1) charge treatment across a covalent QM/MM cut. Default - # 'none' = full-field, the validated ESPF baseline: ESPF couples the MM - # potential to QM atomic-charge operators (H += sum_A phi_A Q_A), which - # already suppresses the spill-out that would justify redistribution in a - # density-based embedding. 'rcd' (conserve deleted M1 charge + its dipole - # via virtual midpoint charges), 'rc', and 'z1' are optional refinements. - # A no-op when the QM region is whole molecules (no link atoms). - self.frontier_scheme = str(frontier_scheme or 'none').lower() - if self.link_atoms and self.espf_full: - note = ("full-field ESPF baseline" if self.frontier_scheme == 'none' - else f"'{self.frontier_scheme}' redistribution (optional refinement)") - print( - f"[QM/MM] {len(self.link_atoms)} covalent boundary bond(s); " - f"frontier-charge embedding = {note}." - ) - self.mm_systems = self.prepare_mm() + # Per-QM-centre set of MM atoms to exclude from the QM-MM electrostatics + # (the 1-2/1-3 bonded neighbours of a frontier atom), applied + # consistently to both the embedding potential and the coupling force so + # the analytic gradient stays exact. Link atoms inherit their frontier + # atom's exclusions plus the MM host. + self._qm_mm_excluded = ( + self._build_qm_mm_exclusions() if self.espf_full else None + ) + # --- Internal helpers ------------------------------------------------- def _detect_link_atoms(self): @@ -529,11 +510,7 @@ def compute_force(self, positions, topology, mm_systems, qm_atoms): self.positions = positions self.topology = topology self.mm_systems = mm_systems - # Re-normalize to ascending (topology) order on every call: callers - # (e.g. QMMM_MD) pass their own config-order qm_atoms here, which would - # otherwise overwrite the sorted copy set in __init__ and reintroduce the - # force mis-scatter this fix prevents. See __init__ for the rationale. - self.qm_atoms = np.array(sorted(int(i) for i in qm_atoms), dtype=int) + self.qm_atoms = qm_atoms potmm = potqm = None if self.Embedding in ("electrostatic", "split") or self.espf_full: @@ -757,41 +734,34 @@ def _qm_center_positions_bohr(self): coords.append([c * self._ANG2BOHR for c in pos]) return np.asarray(coords, dtype=float) - def _frontier_hosts(self): - """[(m1_abs_idx, [m2_abs_idx, ...]), ...] : each unique MM host atom - (the MM end of a bond the QM/MM boundary cuts) with its MM neighbours - M2, from the topology. Empty for a whole-molecule QM region.""" - if not self.link_atoms: - return [] + def _build_qm_mm_exclusions(self): + """List (one set per QM centre, in centre order) of MM atom indices to + exclude from the QM-MM electrostatics: the 1-2/1-3 bonded neighbours of + each frontier QM atom (OpenMM full exclusions, chargeProd==0). Link + atoms inherit their frontier atom's set plus the MM host.""" + nb = next(f for f in self.mm_systems["sys0"].getForces() + if isinstance(f, mm.NonbondedForce)) qm_set = set(int(i) for i in self.qm_atoms) - nbrs = {} - for b in self.topology.bonds(): - i, j = int(b[0].index), int(b[1].index) - nbrs.setdefault(i, set()).add(j) - nbrs.setdefault(j, set()).add(i) - hosts = {} - for link in self.link_atoms: - m1 = int(link.mm_index) - if m1 in hosts: + excl = {int(a): set() for a in self.qm_atoms} + for i in range(nb.getNumExceptions()): + p1, p2, cp, _, _ = nb.getExceptionParameters(i) + if cp.value_in_unit(unit.elementary_charge ** 2) != 0.0: continue - hosts[m1] = sorted(n for n in nbrs.get(m1, ()) if n not in qm_set) - return [(m1, hosts[m1]) for m1 in sorted(hosts)] - - def _embedding_sites(self): - """The MM embedding point-charge set the QM density sees, with the - frontier-host (M1) charges redistributed per ``self.frontier_scheme`` - (RCD by default). Returns (charges (S,), positions (S,3) bohr, scatter), - where scatter maps a per-site force back onto the real MM atoms (identity - for a real atom, 0.5/0.5 for a virtual midpoint charge). With no link - atoms this is the identity -> whole-molecule QM/MM is unchanged.""" - mmq, mm_xyz, mm_idx = self._mm_charges_positions_bohr() - q_of = {int(mm_idx[k]): float(mmq[k]) for k in range(len(mm_idx))} - deleted, delta_q, virtuals = redistribute_frontier_charges( - self._frontier_hosts(), lambda a: q_of.get(a, 0.0), - self.frontier_scheme, - ) - return assemble_embedding_sites( - mm_idx, mmq, mm_xyz, deleted, delta_q, virtuals) + if p1 in qm_set and p2 not in qm_set: + excl[p1].add(p2) + elif p2 in qm_set and p1 not in qm_set: + excl[p2].add(p1) + # NOTE: full-field embedding (no exclusions) is used -- excluding the + # bonded MM neighbours was found to worsen the analytic-gradient + # consistency, so the embedding potential and the coupling force both + # see the complete MM charge set (per the ESPF formulation). + rows = [] + for atom in self.topology.atoms(): + if atom.index in qm_set: + rows.append(set()) + for link in self.link_atoms: + rows.append(set()) + return rows def _box_lengths_bohr(self): """Orthorhombic periodic box lengths (bohr), or None when the QM/MM @@ -831,45 +801,46 @@ def _mm_charges_positions_bohr(self): idx.append(i) return np.asarray(q), np.asarray(xyz, dtype=float), np.asarray(idx, dtype=int) + def _center_mm_mask(self, center_row, mm_idx): + """Boolean mask over the MM arrays: False for MM atoms excluded from + this QM centre's electrostatics (bonded 1-2/1-3 neighbours).""" + excluded = self._qm_mm_excluded[center_row] + if not excluded: + return np.ones(len(mm_idx), dtype=bool) + return np.array([int(m) not in excluded for m in mm_idx], dtype=bool) + def _full_field_potmm(self): - """MM electrostatic potential at every QM centre from the (frontier- - redistributed) embedding charge set: phi_A = sum_s Q_s / |r_A - r_s| - (Hartree/e).""" + """MM electrostatic potential at every QM centre, with each frontier + centre's bonded MM neighbours excluded (consistent with the coupling + force), phi_A = sum_{M not excluded} Q_M / |r_A - r_M| (Hartree/e).""" qm_xyz = self._qm_center_positions_bohr() - q_s, xyz_s, _ = self._embedding_sites() + mmq, mm_xyz, mm_idx = self._mm_charges_positions_bohr() box = self._box_lengths_bohr() potmm = np.zeros(len(qm_xyz)) for a in range(len(qm_xyz)): - d = self._min_image(qm_xyz[a] - xyz_s, box) + mask = self._center_mm_mask(a, mm_idx) + d = self._min_image(qm_xyz[a] - mm_xyz[mask], box) r = np.linalg.norm(d, axis=1) - potmm[a] = np.sum(q_s / r) + potmm[a] = np.sum(mmq[mask] / r) return potmm def _coupling_forces(self, pchg): - """Classical Coulomb force (Hartree/bohr) between the QM ESP charges q_A - and the (frontier-redistributed) MM embedding charges Q_s -- the field- - fluctuation term Tr[q dphi/dx]. A force on a virtual midpoint charge is - chain-ruled onto its real hosts (M1, M2) via the site scatter map. - Returns (F on QM centres, F on real MM atoms, real MM idx).""" + """Classical Coulomb force (Hartree/bohr) between the QM ESP charges + q_A and the (unexcluded) MM charges Q_M -- the field-fluctuation term + Tr[q dphi/dx] of eq 7. Uses the same per-centre exclusions as the + embedding potential. Returns (F on QM centres, F on MM atoms, MM idx).""" qm_xyz = self._qm_center_positions_bohr() - _, _, mm_idx = self._mm_charges_positions_bohr() - q_s, xyz_s, scatter = self._embedding_sites() + mmq, mm_xyz, mm_idx = self._mm_charges_positions_bohr() box = self._box_lengths_bohr() f_qm = np.zeros_like(qm_xyz) - f_site = np.zeros_like(xyz_s) + f_mm = np.zeros_like(mm_xyz) for a in range(len(qm_xyz)): - d = self._min_image(qm_xyz[a] - xyz_s, box) # r_A - r_s + mask = self._center_mm_mask(a, mm_idx) + d = self._min_image(qm_xyz[a] - mm_xyz[mask], box) # r_A - r_M r = np.linalg.norm(d, axis=1) - coeff = pchg[a] * q_s / r ** 3 # q_A Q_s / r^3 + coeff = pchg[a] * mmq[mask] / r ** 3 # q_A Q_M / r^3 f_qm[a] = np.sum(coeff[:, None] * d, axis=0) - f_site -= coeff[:, None] * d # Newton's third law - # Scatter each site force onto the real MM atoms it is built from - # (identity for a real atom, chain rule 0.5/0.5 for a midpoint). - row_of = {int(m): j for j, m in enumerate(mm_idx)} - f_mm = np.zeros((len(mm_idx), 3)) - for s, contribs in enumerate(scatter): - for atom, w in contribs: - f_mm[row_of[atom]] += w * f_site[s] + f_mm[mask] -= coeff[:, None] * d # Newton's third law return f_qm, f_mm, mm_idx def electrostatic_potential(self): diff --git a/pyoqp/oqp/library/qmmm_md.py b/pyoqp/oqp/library/qmmm_md.py index 6572eaea1..c2aea49c4 100644 --- a/pyoqp/oqp/library/qmmm_md.py +++ b/pyoqp/oqp/library/qmmm_md.py @@ -224,7 +224,6 @@ def __init__(self, oqp_cfg=None, mol=None): self.cutoff = _resolve_cutoff(qmmm_cfg.get("cutoff", "PME")) self.embedding = str(qmmm_cfg.get("embedding", "electrostatic")) - self.frontier_scheme = str(qmmm_cfg.get("frontier_scheme", "none")) self.n_steps = int(qmmm_cfg.get("n_steps", 1000)) self.timestep = float(qmmm_cfg.get("timestep", 1.0)) * unit.femtoseconds self.temperature = float(qmmm_cfg.get("temperature", 300.0)) * unit.kelvin @@ -272,15 +271,6 @@ def __init__(self, oqp_cfg=None, mol=None): self._apply_xyz_positions(str(qm_atoms_xyz), qm_list) # ------ store QM config / mol for the driver ---------------------- - # The outer runtype=md only selects this ground-state QM/MM MD driver - # (see pyoqp dispatch); the QM subsystem itself runs a single-point - # energy+gradient each step. Force the internal QM runtype to 'energy' - # so the QM engine's input check accepts a config that arrived with - # runtype=md (config mode; harmless when the key is absent). - if isinstance(qm_cfg, dict): - for _k in list(qm_cfg): - if str(_k).split('.')[-1].strip().lower() == 'runtype': - qm_cfg[_k] = 'energy' self.oqp_cfg = qm_cfg self.mol = mol @@ -355,7 +345,6 @@ def _build_oqp_driver(self): mol=self.mol, Cutoff=self.cutoff, Embedding=self.embedding, - frontier_scheme=self.frontier_scheme, ) self.mm_systems = self.oqp_driver.mm_systems diff --git a/pyoqp/oqp/molecule/oqpdata.py b/pyoqp/oqp/molecule/oqpdata.py index 303a9e430..f8efa94c8 100644 --- a/pyoqp/oqp/molecule/oqpdata.py +++ b/pyoqp/oqp/molecule/oqpdata.py @@ -97,11 +97,6 @@ def _perf_parse_float(v): 'cutoff': {'type': str, 'default': 'NoCutoff'}, 'embedding': {'type': str, 'default': 'electrostatic'}, 'temperature': {'type': float, 'default': '300.0'}, - # Frontier (M1) charge treatment across a covalent QM/MM cut for ESPF - # embedding: none (default = full-field, the validated ESPF baseline; - # ESPF's charge-operator coupling already suppresses spill-out) | rcd | - # rc | z1 (optional redistribution refinements). - 'frontier_scheme': {'type': str, 'default': 'none'}, }, 'input': { 'charge': {'type': int, 'default': '0'}, diff --git a/pyoqp/oqp/openqp.py b/pyoqp/oqp/openqp.py index cfc2d8a27..bcf5a2cb8 100644 --- a/pyoqp/oqp/openqp.py +++ b/pyoqp/oqp/openqp.py @@ -457,7 +457,7 @@ def molecule(self, system=None, system2=None, basis=None, charge=None, def qmmm(self, pdb_file=None, forcefield=None, forcefield_files=None, qm_atoms=None, cutoff=None, embedding=None, rigidwater=None, - frontier_scheme=None, **kwargs): + **kwargs): """Enable ESPF QM/MM embedding and configure the ``[qmmm]`` section. Sets ``[input] qmmm_flag=true`` and populates ``[qmmm]``. The QM geometry @@ -470,15 +470,6 @@ def qmmm(self, pdb_file=None, forcefield=None, forcefield_files=None, ``qm_atoms`` accepts a string (``"0-2"`` / ``"0 1 2"``) or a list of indices; a ``forcefield`` list is joined into a comma-separated string. - ``frontier_scheme`` sets the MM frontier-host (M1) charge treatment when - the QM/MM partition cuts a covalent bond in the ESPF path: - ``'none'`` (default, the validated full-field ESPF baseline), or the - optional redistributions ``'rcd'`` / ``'rc'`` / ``'z1'``. It is a no-op - for whole-molecule QM regions. Covalent QM/MM boundaries are handled by - the ground-state QM/MM MD path (``QMMM_MD``); the nonadiabatic - ``runtype=namd`` path does not yet append link atoms to its QM molecule - and raises on a covalent cut. - Any other ``[qmmm]`` keyword can be passed through as a keyword argument. """ if forcefield is not None and forcefield_files is not None: @@ -503,8 +494,6 @@ def qmmm(self, pdb_file=None, forcefield=None, forcefield_files=None, updates["embedding"] = embedding if rigidwater is not None: updates["rigidwater"] = rigidwater - if frontier_scheme is not None: - updates["frontier_scheme"] = frontier_scheme updates.update(kwargs) if updates: self.section("qmmm", **updates) diff --git a/tests/test_openqp_api.py b/tests/test_openqp_api.py index 122daf0fb..b95d9b2fc 100644 --- a/tests/test_openqp_api.py +++ b/tests/test_openqp_api.py @@ -34,7 +34,6 @@ def _string(value): "cutoff": {"type": _string, "default": "NoCutoff"}, "embedding": {"type": _string, "default": "electrostatic"}, "rigidwater": {"type": bool, "default": "False"}, - "frontier_scheme": {"type": _string, "default": "none"}, }, "md": { "nstep": {"type": int, "default": "100"}, @@ -809,18 +808,6 @@ def test_qmmm_enables_flag_and_section(self): self.assertEqual(config["input"]["runtype"], "energy") self.assertEqual(config["qmmm"]["embedding"], "electrostatic") - def test_qmmm_frontier_scheme_sets_section_key(self): - openqp = load_openqp_module() - job = ( - openqp.OpenQP(project="qmmm_rcd") - .molecule("ala.pdb 8 9 16 17 18", basis="6-31g*") - .qmmm(pdb_file="ala.pdb", forcefield="amber14-all.xml", - qm_atoms=[8, 9, 16, 17, 18], embedding="electrostatic", - frontier_scheme="rcd") - ) - config = job.to_input_dict() - self.assertEqual(config["qmmm"]["frontier_scheme"], "rcd") - def test_qmmm_normalizes_lists_and_rejects_duplicate_forcefield(self): openqp = load_openqp_module() job = openqp.OpenQP(project="qmmm_md").molecule("m.pdb 0 1 2", basis="6-31g") diff --git a/tests/test_qmmm_frontier.py b/tests/test_qmmm_frontier.py deleted file mode 100644 index 2d57eca32..000000000 --- a/tests/test_qmmm_frontier.py +++ /dev/null @@ -1,260 +0,0 @@ -"""Unit tests for QM/MM frontier (M1) charge redistribution. - -These exercise the pure-Python RCD/RC/Z1 redistribution used by the QM/MM -driver to treat the electrostatics of a covalent QM/MM boundary: deleting the -MM host (M1) charge and redistributing it to virtual point charges at the M1-M2 -bond midpoints. They load the module by path so they run without the compiled -OpenQP backend or OpenMM. - -The key property tested is gradient consistency: because every virtual charge -sits at a bond midpoint (a linear function of two real atom positions), the -Coulomb force on it must redistribute onto its hosts by the chain rule, so the -analytic coupling force matches a finite-difference of the coupling energy. -""" - -import importlib.util -import sys -import unittest -from pathlib import Path - -import numpy as np - -ROOT = Path(__file__).resolve().parents[1] - - -def load_module(name, relative_path): - spec = importlib.util.spec_from_file_location(name, ROOT / relative_path) - module = importlib.util.module_from_spec(spec) - sys.modules[name] = module - spec.loader.exec_module(module) - return module - - -def _coulomb_energy(qA, rA, qs, rs): - """E = sum_A sum_s qA * qs / |rA - rs| (arbitrary consistent units).""" - E = 0.0 - for a in range(len(qA)): - d = rA[a] - rs - r = np.linalg.norm(d, axis=1) - E += qA[a] * np.sum(qs / r) - return E - - -def _analytic_mm_forces(qA, rA, qs, rs, scatter, n_mm, row_of): - """Force on each real MM atom, f = -dE/dR, via site forces + scatter.""" - f_site = np.zeros_like(rs) - for a in range(len(qA)): - d = rA[a] - rs - r = np.linalg.norm(d, axis=1) - coeff = qA[a] * qs / r ** 3 - f_site -= coeff[:, None] * d - f_mm = np.zeros((n_mm, 3)) - for s, contribs in enumerate(scatter): - for atom, w in contribs: - f_mm[row_of[atom]] += w * f_site[s] - return f_mm - - -class TestFrontierCharge(unittest.TestCase): - def setUp(self): - self.c = load_module( - "qmmm_connectivity_frontier_under_test", - "pyoqp/oqp/library/qmmm_connectivity.py", - ) - - @staticmethod - def _charge_of(charges): - return lambda a: charges[a] - - # -- charge conservation ------------------------------------------- - - def test_rcd_conserves_total_charge(self): - charges = {10: 0.5, 11: -0.3, 12: -0.2, 13: 0.11} # M1=10, M2={11,12} - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of(charges), "rcd") - orig = sum(charges.values()) - new = (sum(charges[a] for a in charges if a not in deleted) - + sum(dq.values()) + sum(v.charge for v in virt)) - self.assertAlmostEqual(new, orig, places=12) - - def test_rc_conserves_total_charge(self): - charges = {10: 0.5, 11: -0.3, 12: -0.2} - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of(charges), "rc") - self.assertEqual(dq, {}) # RC does not touch M2 charges - new = (sum(charges[a] for a in charges if a not in deleted) - + sum(v.charge for v in virt)) - self.assertAlmostEqual(new, sum(charges.values()), places=12) - - def test_z1_deletes_only(self): - charges = {10: 0.5, 11: -0.3, 12: -0.2} - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of(charges), "z1") - self.assertEqual(deleted, {10}) - self.assertEqual(dq, {}) - self.assertEqual(virt, []) - - def test_none_is_noop(self): - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of({10: 1.0}), "none") - self.assertEqual((deleted, dq, virt), (set(), {}, [])) - - def test_isolated_host_falls_back_to_deletion(self): - # M1 with no MM neighbours: cannot redistribute -> delete only. - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [])], self._charge_of({10: 0.4}), "rcd") - self.assertEqual((deleted, dq, virt), ({10}, {}, [])) - - def test_unknown_scheme_raises(self): - with self.assertRaises(ValueError): - self.c.redistribute_frontier_charges( - [(10, [11])], self._charge_of({10: 1.0, 11: 0.0}), "bogus") - - # -- dipole conservation ------------------------------------------- - - def _dipole_about_m1(self, pos, m1, dq, virt): - mu = np.zeros(3) - for v in virt: - r = v.weights[0] * pos[v.hosts[0]] + v.weights[1] * pos[v.hosts[1]] - mu += v.charge * (r - pos[m1]) - for a, dqa in dq.items(): - mu += dqa * (pos[a] - pos[m1]) - return mu - - def test_rcd_conserves_dipole_about_m1_but_rc_does_not(self): - # M1 deleted (dipole of a point charge about itself is zero), so the - # redistributed set must also have zero dipole about M1 for RCD. - pos = {10: np.array([0.0, 0.0, 0.0]), - 11: np.array([1.5, 0.0, 0.0]), - 12: np.array([-0.4, 1.3, 0.2])} - charges = {10: 0.7, 11: -0.2, 12: -0.2} - _, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of(charges), "rcd") - mu_rcd = self._dipole_about_m1(pos, 10, dq, virt) - self.assertTrue(np.allclose(mu_rcd, 0.0, atol=1e-12), mu_rcd) - - _, dq_rc, virt_rc = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of(charges), "rc") - mu_rc = self._dipole_about_m1(pos, 10, dq_rc, virt_rc) - self.assertFalse(np.allclose(mu_rc, 0.0, atol=1e-9), mu_rc) - - def test_virtual_charges_sit_at_bond_midpoints(self): - _, _, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of({10: 0.6, 11: 0.0, 12: 0.0}), "rcd") - for v in virt: - self.assertEqual(v.weights, (0.5, 0.5)) - self.assertEqual(v.hosts[0], 10) - - # -- site assembly -------------------------------------------------- - - def test_assemble_identity_without_frontier(self): - mm_idx = [3, 4, 5] - q = [0.1, -0.2, 0.3] - xyz = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]], float) - qs, rs, scatter = self.c.assemble_embedding_sites( - mm_idx, q, xyz, set(), {}, []) - self.assertTrue(np.allclose(qs, q)) - self.assertTrue(np.allclose(rs, xyz)) - self.assertEqual([s for s in scatter], [((3, 1.0),), ((4, 1.0),), ((5, 1.0),)]) - - def test_assemble_drops_deleted_and_adds_virtuals(self): - mm_idx = [10, 11, 12] - q = [0.6, -0.2, -0.2] - xyz = np.array([[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0]], float) - deleted, dq, virt = self.c.redistribute_frontier_charges( - [(10, [11, 12])], self._charge_of({10: 0.6, 11: -0.2, 12: -0.2}), "rcd") - qs, rs, scatter = self.c.assemble_embedding_sites( - mm_idx, q, xyz, deleted, dq, virt) - # M1 (index 10) has no direct site; 2 real (M2) + 2 virtual midpoints. - self.assertEqual(len(qs), 4) - # a midpoint site is built from M1 and one M2 with 0.5/0.5 weights - mids = [s for s in scatter if len(s) == 2] - self.assertEqual(len(mids), 2) - for s in mids: - self.assertEqual(s[0][1], 0.5) - self.assertEqual(s[1][1], 0.5) - - # -- assembled-charge conservation (through assemble, not redistribute) -- - - def _assembled_total(self, mm_idx, mm_q, mm_xyz, frontier, scheme): - charge_of = lambda a: mm_q[mm_idx.index(a)] - deleted, dq, virt = self.c.redistribute_frontier_charges( - frontier, charge_of, scheme) - qs, _, _ = self.c.assemble_embedding_sites( - mm_idx, mm_q, mm_xyz, deleted, dq, virt) - return float(np.sum(qs)) - - def test_assembled_charge_conserved_isolated_host(self): - mm_idx = [10, 11, 12, 13] - mm_q = [0.5, -0.3, -0.2, 0.11] - xyz = np.array([[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0], [3, 0, 0]], float) - for scheme in ("rc", "rcd"): - self.assertAlmostEqual( - self._assembled_total(mm_idx, mm_q, xyz, [(10, [11, 12])], scheme), - sum(mm_q), places=12, msg=scheme) - - def test_assembled_charge_conserved_adjacent_hosts(self): - # Two frontier hosts bonded to each other (adjacent cuts): each is in the - # other's M2 list, so an RCD delta_q lands on a deleted atom. The residual - # -charge site must keep total charge conserved. - mm_idx = [10, 11, 12, 13] - mm_q = [0.5, 0.3, -0.2, -0.1] - xyz = np.array([[0, 0, 0], [1.5, 0, 0], [-1.4, 0.3, 0], [2.9, 0.2, 0]], float) - frontier = [(10, [11, 12]), (11, [10, 13])] - for scheme in ("rc", "rcd"): - self.assertAlmostEqual( - self._assembled_total(mm_idx, mm_q, xyz, frontier, scheme), - sum(mm_q), places=12, msg=scheme) - - # -- gradient consistency (the decisive test) ----------------------- - - def _scenario(self): - # 2 QM centres (fixed ESP charges/positions); MM = M1(0)+M2a(1)+M2b(2) - # + 2 bulk (3,4). M1 deleted & redistributed by RCD. - qA = np.array([0.3, -0.25]) - rA = np.array([[0.2, 0.1, -0.3], [-0.7, 0.4, 0.5]]) - mm_idx = [0, 1, 2, 3, 4] - mm_q = [0.55, -0.15, -0.18, 0.2, -0.31] - mm_xyz = np.array([ - [1.4, 0.0, 0.0], # M1 - [2.6, 0.3, -0.2], # M2a - [1.1, 1.5, 0.4], # M2b - [3.5, -1.0, 1.2], # bulk - [-2.0, -1.5, 0.8], # bulk - ], float) - return qA, rA, mm_idx, mm_q, mm_xyz - - def test_coupling_force_matches_finite_difference(self): - qA, rA, mm_idx, mm_q, mm_xyz = self._scenario() - row_of = {m: j for j, m in enumerate(mm_idx)} - charge_of = self._charge_of({i: mm_q[j] for j, i in enumerate(mm_idx)}) - - # isolated host (0 -> M2 {1,2}), and two adjacent hosts (0,1 bonded) - for frontier in ([(0, [1, 2])], [(0, [1, 3]), (1, [0, 4])]): - for scheme in ("none", "z1", "rc", "rcd"): - deleted, dq, virt = self.c.redistribute_frontier_charges( - frontier, charge_of, scheme) - - def energy(xyz): - qs, rs, _ = self.c.assemble_embedding_sites( - mm_idx, mm_q, xyz, deleted, dq, virt) - return _coulomb_energy(qA, rA, qs, rs) - - qs, rs, scatter = self.c.assemble_embedding_sites( - mm_idx, mm_q, mm_xyz, deleted, dq, virt) - f_mm = _analytic_mm_forces(qA, rA, qs, rs, scatter, len(mm_idx), row_of) - - h = 1e-6 - fd = np.zeros_like(mm_xyz) - for j in range(len(mm_idx)): - for c in range(3): - xp = mm_xyz.copy(); xp[j, c] += h - xm = mm_xyz.copy(); xm[j, c] -= h - fd[j, c] = -(energy(xp) - energy(xm)) / (2 * h) - self.assertTrue( - np.allclose(f_mm, fd, atol=1e-6), - f"scheme={scheme}: analytic vs FD mismatch\n{f_mm - fd}") - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_qmmm_frontier_openmm.py b/tests/test_qmmm_frontier_openmm.py deleted file mode 100644 index 362afc44d..000000000 --- a/tests/test_qmmm_frontier_openmm.py +++ /dev/null @@ -1,123 +0,0 @@ -"""Integration test: QM/MM frontier (M1) charge redistribution on a real force -field and covalent boundary. - -Builds the alanine dipeptide (ACE-ALA-NH2) with the AMBER-14 force field via -OpenMM, cuts the QM/MM boundary through the ALA C-CA backbone bond (QM = the -C-terminal amide), and checks that: - - * the cut is detected as exactly one hydrogen link atom, and - * RCD / RC frontier-charge redistribution conserves the total MM charge that - the QM density is embedded in (Z1 deletes and so does not). - -This exercises the same pure connectivity core the QM/MM driver uses, on real -AMBER partial charges rather than synthetic ones. It is skipped when OpenMM (or -the bundled amber14-all.xml) is unavailable, so it does not require the compiled -OpenQP backend. -""" - -import importlib.util -import os -import sys -import unittest -from pathlib import Path - -import numpy as np - -ROOT = Path(__file__).resolve().parents[1] - -try: - import openmm as mm - import openmm.app as app - import openmm.unit as u - _HAVE_OPENMM = True -except Exception: # pragma: no cover - optional dependency - _HAVE_OPENMM = False - - -def load_module(name, relative_path): - spec = importlib.util.spec_from_file_location(name, ROOT / relative_path) - module = importlib.util.module_from_spec(spec) - sys.modules[name] = module - spec.loader.exec_module(module) - return module - - -@unittest.skipUnless(_HAVE_OPENMM, "OpenMM not installed") -class TestFrontierChargeOnAmberAlanine(unittest.TestCase): - def setUp(self): - self.c = load_module( - "qmmm_connectivity_amber_under_test", - "pyoqp/oqp/library/qmmm_connectivity.py", - ) - pdb_path = ROOT / "examples" / "QMMM" / "ala.pdb" - if not pdb_path.exists(): - self.skipTest("examples/QMMM/ala.pdb missing") - self.pdb = app.PDBFile(str(pdb_path)) - try: - ff = app.ForceField("amber14-all.xml") - # Keep the System referenced: its Force objects become invalid once - # the owning System is garbage-collected. - self.system = ff.createSystem(self.pdb.topology) - except Exception as err: # pragma: no cover - self.skipTest(f"amber14-all.xml unavailable: {err}") - self.nb = next(f for f in self.system.getForces() - if isinstance(f, mm.NonbondedForce)) - self.charge = { - i: self.nb.getParticleParameters(i)[0].value_in_unit(u.elementary_charge) - for i in range(self.nb.getNumParticles()) - } - self.z = {a.index: a.element.atomic_number for a in self.pdb.topology.atoms()} - self.bonds = [(b[0].index, b[1].index) for b in self.pdb.topology.bonds()] - self.name = {a.index: (a.residue.name, a.name) - for a in self.pdb.topology.atoms()} - # QM = C-terminal amide (ALA C, ALA O + all NH2 atoms): cuts only ALA CA-C. - self.qm = [i for i, n in self.name.items() - if n[0] == "NH2" or n in (("ALA", "C"), ("ALA", "O"))] - - def _frontier(self, links): - qm_set = set(self.qm) - nbrs = {} - for i, j in self.bonds: - nbrs.setdefault(i, set()).add(j) - nbrs.setdefault(j, set()).add(i) - hosts = {l.mm_index: sorted(x for x in nbrs[l.mm_index] if x not in qm_set) - for l in links} - return [(m1, hosts[m1]) for m1 in sorted(hosts)] - - def test_single_link_atom_on_backbone_cut(self): - links = self.c.detect_link_atoms(self.bonds, self.qm, lambda i: self.z[i]) - self.assertEqual(len(links), 1) - la = links[0] - self.assertEqual(self.name[la.qm_index], ("ALA", "C")) - self.assertEqual(self.name[la.mm_index], ("ALA", "CA")) - self.assertTrue(0.0 < la.g < 1.0) - - def test_rcd_and_rc_conserve_amber_charge(self): - links = self.c.detect_link_atoms(self.bonds, self.qm, lambda i: self.z[i]) - frontier = self._frontier(links) - mm_idx = [i for i in range(self.nb.getNumParticles()) if i not in set(self.qm)] - mm_q = [self.charge[i] for i in mm_idx] - mm_xyz = np.zeros((len(mm_idx), 3)) - raw = sum(mm_q) - for scheme in ("rcd", "rc"): - deleted, dq, virt = self.c.redistribute_frontier_charges( - frontier, lambda a: self.charge[a], scheme) - qs, _, _ = self.c.assemble_embedding_sites( - mm_idx, mm_q, mm_xyz, deleted, dq, virt) - self.assertAlmostEqual(float(qs.sum()), raw, places=9, msg=scheme) - - def test_z1_deletes_host_charge(self): - links = self.c.detect_link_atoms(self.bonds, self.qm, lambda i: self.z[i]) - frontier = self._frontier(links) - (m1, _), = frontier - mm_idx = [i for i in range(self.nb.getNumParticles()) if i not in set(self.qm)] - mm_q = [self.charge[i] for i in mm_idx] - deleted, dq, virt = self.c.redistribute_frontier_charges( - frontier, lambda a: self.charge[a], "z1") - qs, _, _ = self.c.assemble_embedding_sites( - mm_idx, mm_q, np.zeros((len(mm_idx), 3)), deleted, dq, virt) - self.assertAlmostEqual(float(qs.sum()), sum(mm_q) - self.charge[m1], places=9) - - -if __name__ == "__main__": - unittest.main()