QM/MM: full-ESPF covalent-boundary embedding under TB dispatch (fixes DFTB link-atom energy conservation, enables xTB QM/MM)#270
Open
karmachoi wants to merge 84 commits into
Conversation
…OpenQP for computing energy and gradient
…ine and OpenQP for computing energy and gradient" This reverts commit ab7afca.
…OpenQP to compute QM/MM energy and gradient
merge main to qmmm
Bring in the ESPF electrostatic-embedding QM/MM work (Schwinn &
Huix-Rotllant method) as the base for NAMD-QMMM:
- source/modules/qmmm.F90: ESPF QM/MM energy + gradient, ESP charges
- ESPF 1e integrals/gradient hooks in int1.F90/grd1.F90/int1e.F90
- population_analysis.F90: excited-state Mulliken/ESP charges
- Python: qmmm_driver.py, qmmm_md.py, SimulationManager.py,
libopenmm.py, libdlfind.py, libscipy QMMMOpt; oqpdata set_qmmm_*
- examples/QMMM/* decks; qmmm_flag control parameter
Conflict resolutions:
- oqp.h / types.F90: keep both PCM/SD and qmmm_flag fields (C-struct
order kept consistent across both)
- int1e.F90: keep stale-Q-cache removal and add tags_qmmm
- tdhf_mrsf_z_vector.F90: keep main's refactored solver
(CG/MINRES/GMRES + build_mrsf_relaxed_density_and_w); inject the
excited-state mulliken/ESP-charge calls after the relaxed density
is built; drop feat's pre-refactor duplicate post-processing
- population_analysis.F90: drop feat's duplicate mulliken_excited;
keep main's version (with the C-binding wrapper)
- input_checker.py: keep main's modern diagnostic checker; QM/MM
input validation to be re-ported into it (follow-up)
- runfunc.py: keep both optimizer suites; add DLFind/QMMMOpt imports
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
End-to-end MRSF-TDDFT surface-hopping dynamics in OpenQP, ported from the
GAMESS namd.src algorithm. Verified energy conservation (~2e-6 Hartree over
2.5 fs on H2O MRSF) end-to-end through SCF -> MRSF -> gradient -> state overlap
-> Fortran FSSH each step.
Fortran (physics):
- source/modules/namd.F90: FSSH core (TDC from state overlaps, RK4 amplitude
propagation, cumulative Tully probabilities, isotropic energy-conserving
rescaling, frustrated-hop handling) + EDC decoherence (Granucci-Persico)
+ trivial-crossing diabatic following (SC-FSSH style) + C-bound driver
mrsf_namd_hop. Consistent atomic units; amu->a.u. mass conversion.
- tagarray_driver.F90: NAMD state tags (coef/velocity/params/results) as
flat 1-D records; also fix all_tags dimension (was 39, list has 46 after
the QM/MM merge) -> implied-shape (*) so it can't drift again.
- include/oqp.h: declare mrsf_namd_hop.
Python (sequencing):
- library/namd.py: velocity-Verlet FSSH trajectory driver reusing
SinglePoint/BasisOverlap/get_states_overlap/Gradient + oqp.mrsf_namd_hop.
- molecule.py: register NAMD tags. oqpdata.py: [md] config section.
- runfunc.compute_namd + pyoqp dispatch + input_checker namd runtype/_check.
Merge-regression fixes (QM/MM made OpenMM mandatory for ALL runs):
- utils/qmmm.py: make OpenMM import + unit constants optional.
- pyoqp.py main(): detect qmmm_flag via configparser, import qmmm_md
(OpenMM) only when actually in QM/MM-MD mode.
Literature basis recorded in session RESEARCH_ic_isc_methods.md (NPI/EDC/
SC-FSSH for IC; SHARC spin-adiabatic for ISC). NPI TDC and ISC pending.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
For NoCutoff embedding, electrostatic_potential() returned potqm=None, so
OQP::POTQM was never a valid tagarray record and the SCF's
add_potqm_contributions aborted ("Record OQP::POTQM"). The QM-QM Ewald
self-interaction correction is identically zero for non-periodic systems, so
return a zero matrix instead of None. The QM/MM ground-state baseline
(water dimer, electrostatic embedding, NoCutoff) now runs end to end.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…odic) NAMD_QMMM extends the gas-phase FSSH driver to QM/MM: the QM region is the OQP Molecule with ESPF electrostatic embedding; the MM region + QM/MM coupling come from OpenMM via the OpenQpQMMM helpers. Per step: sync positions (QM mol + OpenMM context) -> MM potential at QM atoms -> embedded SCF + MRSF excitation -> active-state embedded gradient (Gradient + grad_esp_qmmm_excited) -> excited ESPF QM charges -> MM forces -> full-system velocity Verlet (QM+MM, atomic units) -> QM-only FSSH hop (rescale only QM velocities, as GAMESS RESCALV). VALIDATED: water dimer (QM=1 water, MM=1 water, electrostatic embedding, NoCutoff). Energy conserved to ~2e-6 Hartree over 8 steps at dt=0.05 fs (the full embedded QM + MM + coupling force assembly is consistent with the total energy). Larger dt drifts only in the strained close-contact region (integrator step size, not a force bug). Wiring: [qmmm] schema keys (pdb_file/forcefield_files/qm_atoms/cutoff/embedding/ temperature); main() routes qmmm_flag+runtype=namd through the Runner to compute_namd (not the ground-state QMMM_MD path); compute_namd branches to NAMD_QMMM when qmmm_flag. Periodic (PME) embedding and hop-on-crossing + GAMESS validation are the remaining QM/MM items. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Generalize NAMD_QMMM beyond NoCutoff: resolve [qmmm] cutoff to the OpenMM method (PME/Ewald/CutoffPeriodic/NoCutoff), and sync the periodic Ewald correction contexts (simew/simor) each step. Periodic embedding now runs end to end (water dimer, 30 A box, PME): embedded SCF with the QM-QM Ewald self-interaction correction (POTQM) + PME MM. KNOWN LIMITATION: periodic forces are not yet energy-conserving — E_tot drifts ~3e-3 Ha even at dt=0.05 fs (where NoCutoff conserves to ~2e-6), a force-energy inconsistency from the MISSING gradient of the QM-QM Ewald correction (POTQM); NoCutoff is unaffected because POTQM=0 there. Completing the periodic-QM/MM gradient (Ewald-correction force) is the remaining task for the NAMD-PBC target. NoCutoff QM/MM NAMD regression-checked: still conserves to ~2e-6 Ha. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…source Implement the QM-QM Ewald self-interaction correction force for periodic QM/MM NAMD (the gradient of the POTQM energy term -1/2 q^T POTQM q), obtained from OpenMM (Ewald - direct) force differences with the ESPF QM charges, the force-space analog of electrostatic_potential's energy construction. DIAGNOSIS: this term is negligible for large boxes (30 A water dimer: images ~30 A away) and does NOT remove the periodic energy drift (~3e-3 Ha/6 steps at dt=0.05 fs). The drift is genuine force-energy inconsistency (excess work: dKE=+0.010 > -dPE=+0.0075), traced to the PME embedding consistency term (the Ewald POTMM coupling / ESPF charge response under PME), NOT the QM-QM correction. NoCutoff remains energy-conserving (~2e-6 Ha). Completing the periodic embedding force (with GAMESS periodic formalism as reference + FD validation) is the remaining NAMD-PBC task. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ings
Add the survey's recommended IC coupling scheme as a selectable option
([md] tdc=npi; default tdc=fd unchanged). The Python driver computes the TDC
from the phase-corrected state overlap and passes it to the Fortran hop via a
flat OQP::namd_tdc tag (Fortran falls back to its finite-difference if a zero
matrix is supplied).
NPI is implemented in its rigorous matrix form (Meek & Levine, JPCL 5, 2351
(2014)): T = logm(s (s^T s)^{-1/2}) / dt, the real antisymmetric logarithm of
the Loewdin-orthonormalised step overlap. Validated:
* reduces EXACTLY to the 2-state identity T*dt = arcsin(s_10);
* antisymmetric; reduces to finite difference in the weak-coupling limit;
* FD-via-tag reproduces the prior all-Fortran FD result bit-for-bit
(ethylene MECI: hop at step 35, energy conserved) -> pipeline layout/sign
verified;
* NPI runs end-to-end, conserves energy, and hops identically in the smooth
regime (NPI and FD agree when couplings don't spike within a step).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
# Conflicts: # include/oqp.h # pyoqp/oqp/molecule/molecule.py # pyoqp/oqp/molecule/oqpdata.py # pyoqp/oqp/pyoqp.py # pyoqp/oqp/utils/input_checker.py # source/modules/int1e.F90 # source/tagarray_driver.F90
…les ISC) Decouple mrsf_namd_hop from tddft%nstate: number of states, absolute state energies, and the state-overlap matrix now arrive via namd_* tags (OQP::namd_eabs, OQP::namd_stas, params slot 13 = n), alongside the existing namd_tdc/coef/velocity/params. The same kernel now serves same-spin MRSF (n = tddft%nstate) and spin-adiabatic SOC NAMD (n = ns + 3*nt). Gas-phase regression: reproduces the prior result bit-for-bit (ethylene MECI, hop at step 35, E_tot conserved) -> the generalization is behaviour-preserving. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
NAMD_SOC ([md] soc=True) drives intersystem-crossing dynamics. Per step: SCF -> singlet MRSF -> triplet MRSF -> oqp.soc_mrsf (builds & diagonalises H = diag(E_MCH) + H_SOC) -> spin-adiabatic energies (soc_eval, cm^-1 rel e0) and eigenvectors U. Nuclei propagate on the active spin-adiabatic surface with the dominant-MCH-component gradient (the pure singlet/triplet MRSF gradient of the largest |U|^2 component -- standard weak/moderate-SOC approximation). VALIDATED (SOC-BOMD): H2O MRSF, active spin-adiabatic state = T1; energy conserved to ~6e-6 Ha over 6 steps (dt=0.1 fs). The spin-adiabatic eigenvalue matches the dominant pure-state energy to ~1e-6 Ha (weak SOC), so the dominant-component gradient is consistent -> conservation holds. Reuses the generalised mrsf_namd_hop kernel (n = ns + 3nt). The spin-adiabatic FSSH hopping layer (amplitude propagation with U-phase tracking + block-diagonal MCH state overlap) is the next increment; a heavy-atom system is needed to exercise actual ISC hops (water SOC is negligible, w=1.000 pure triplet). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
# Conflicts: # pyoqp/oqp/pyoqp.py
- H1 (grd1.F90): grad_elpot now builds the gradient density via prepare_grad_density, so pure-spherical (l>=2) shells are expanded to the Cartesian effective density that comp_coulomb_der1 expects. It previously sliced the raw spherical density with ao_offset, silently corrupting ESPF QM/MM forces for spherical-harmonic bases (cc-pVxZ/def2; 6-31G* Cartesian-d masked it). FD-verified on cc-pVDZ: max |grad-FD| 2.4e-3 -> 2.1e-4 (= grid noise floor); Cartesian unchanged. - M1a (qmmm_driver.py): zero POTQM in forces_qm_openqp so the embedded energy is consistent with the (absent) POTQM force, matching the NAMD-driver fix. - M1b (qmmm_driver.py): correct unit factors 49578.9 -> 49614.75 (Ha/bohr -> kJ/mol/nm) and 2625.5 -> 2625.499639 (Ha -> kJ/mol); drop stray debug print. - M1c (qmmm_md.py): update the QM/MM force at the current positions before the integrator step (was one step stale, breaking energy conservation). - M2 (oqpdata.py): fix set_qmmm_constraints typo (contraints -> constraints) so OpenMM constraints are actually applied. - M3 (population_analysis.F90): restore the log-unit open in mulliken_excited (open was commented out while close(iw) remained). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Resolve the single conflict in pyoqp/oqp/pyoqp.py: keep the PR's NAMD/QMMM initialization path (logfile setup + banner in Runner.__init__, slim run()) while adopting upstream's build-version logging via _openqp_build_label() in the section='start' dump_log. Drops the now-redundant inline git-HEAD build block (duplicated _openqp_build_label) and upstream's logfile/banner block inside run() (already moved into __init__). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Removes docs/pcm_ddx_validation.md and docs/plans/*. Note: these files originate from upstream/main; this drops them on the PR branch per request. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
After removing the top-level docs/ folder, clean up the now-dangling mentions of docs/pcm_ddx_validation.md and docs/plans/2026-06-07-symmetry-reductions-design.md in a code comment, a test docstring, a Fortran comment, and the CI workflow note. No functional change. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Addresses two review findings on PR #205 (iamDongHwan): 1. input_checker._check_system treated the entire first line of [input] system as a filename, so a QM/MM spec like "mol.pdb 9 10 17-19" was checked as the path "mol.pdb 9 10 17-19" and always failed with "Referenced XYZ file does not exist." Now the first whitespace token is taken as the path; a .pdb path is validated for existence and its trailing QM-atom indices are validated (integers/ranges, non-negative, unique) via a new _check_qm_atom_indices helper, mirroring oqpdata.read_system. This removes the need to comment out input validation in pyoqp.py. 2. _openmm_from_oqp compared constraints with `is None`, but the config delivers the string "None" (the default), so it fell through to the error path ("Unknown type of constraint None"). Now None, "", and "none" (case-insensitive) are all treated as no constraint. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
The OpenQpQMMM driver (NAMD / QMMM_MD path) assumed the QM region was a set
of whole molecules, so a QM/MM partition that cuts a covalent bond left an
uncapped dangling bond and gave wrong QM energies/forces. This adds hydrogen
link-atom capping to that driver, matching the scheme already used by the
optimizer path (utils/qmmm.py).
New module pyoqp/oqp/library/qmmm_connectivity.py (OpenMM/OpenQP-free, unit
tested):
- detect_link_atoms(): finds bonds with exactly one endpoint in the QM
region and builds one hydrogen link atom per severed bond; rejects
unphysical partitions (g outside (0,1)).
- link_g_factor()/COVALENT_RADII: scaled-position factor
g = (r_H + r_QM)/(r_QM + r_MM), radii matching utils/qmmm.
- link_atom_position(): R_L = R_QM + g(R_MM - R_QM).
- project_link_gradient(): chain-rule split (1-g)->QM host, g->MM host.
Driver integration (OpenQpQMMM):
- detect link atoms once from the topology at construction;
- append capping H atoms to the QM geometry (config and mol modes);
- size the ESPF embedding arrays to natom = nqm + nlink (link atoms are
not embedded in the MM field and carry no periodic self-image term);
- fold each link atom's ESP charge onto its QM host so the QM->MM charge
is conserved;
- redistribute link-atom gradients onto their QM/MM host atoms.
When the QM region is whole molecules (e.g. a water box) no bond is cut, so
link_atoms is empty and every path is a no-op — behaviour is unchanged.
Tests: tests/test_qmmm_connectivity.py (8 cases: detection, determinism,
multi-link, whole-molecule no-op, unphysical-partition guard, geometry,
gradient conservation). Verified with `python -m unittest`.
Note: the pure-Python connectivity core is unit tested; end-to-end QM/MM
with link atoms still needs validation against a full build (ILP64
BLAS/LAPACK), which is not available in this environment.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
…natom!=3 QM regions Building liboqp and finite-difference-testing the QM/MM forces exposed three latent bugs that were masked whenever the QM region was a single 3-atom water (the only case previously exercised): 1. ESPF gradient orientation. OQP::ESPF_GRAD is declared Fortran (3, natom) but its flat buffer is atom-major, matching the QM gradient. Adding the (3, natom) view directly is only shape-correct when natom == 3 (a square coincidence) and silently wrong otherwise. Reshape the flat buffer to (natom, 3) before adding, in both the HF and TDHF branches. 2. forces_mm charge index. np.where(...)[0] yields a length-1 array, so the QM-MM exception charge product became an array and OpenMM rejected it. This branch only fires for QM-MM exception pairs, which do not exist for separated waters. Extract a scalar index. 3. QM-MM bonded exclusions. A 1-2/1-3 exclusion (chargeProd==0) must stay a full exclusion; turning it into a charged pair changes the set of non-excluded exceptions (OpenMM raises in updateParametersInContext) and double-counts QM->MM electrostatics already carried by ESPF. Skip full-exclusion exceptions. Validation (built with OpenBLAS ILP64): FD force check passes for QM=one water (natom=3, ~8e-4) and QM=two waters (natom=6, no boundary, components match). The pure-QM gradient of a 6-atom link-capped fragment matches FD to ~1e-5. Known limitation: with electrostatic embedding across a *covalent* boundary (hydrogen link atom), the ESPF embedding gradient at the boundary is not yet consistent with the energy (FD mismatch); the link-atom geometry, detection and gradient-projection are correct, but the ESPF boundary term needs further work. Whole-molecule QM regions are unaffected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
OpenQP writes oqp_project.json (default project restart/output JSON) when a calculation runs from the config path, dirtying the working tree during examples/validation runs. Ignore it. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
…(gated)
Adds an opt-in Embedding="espf" path to OpenQpQMMM that routes the entire
QM-MM electrostatic coupling through ESPF instead of splitting it between the
OpenQP QM gradient and OpenMM point charges. The split scheme is only
self-consistent for whole-molecule QM regions; across a covalent QM/MM
boundary its analytic force is wrong (FD rel ~2.3) because the OpenMM 1-2/1-3
exclusions make the charge-fluctuation and field-fluctuation halves of the
coupling gradient inconsistent.
Full-ESPF scheme (non-periodic / NoCutoff):
- embedding potential phi_A computed directly from all MM charges;
- QM-MM coupling energy lives entirely in the embedded SCF (no double-count
subtraction); OpenMM reduced to pure MM-MM (QM charges zeroed);
- the field-fluctuation term Tr[q dphi/dx] (eq 7, Huix-Rotllant & Ferre
JCTC 2021) is added analytically as the classical Coulomb force between the
QM ESP charges and the MM charges, on both QM and MM atoms, with link-atom
projection.
On the alanine C8-CA7 boundary this cuts the analytic-force error from FD
rel ~2.3 to ~0.06 (physical magnitudes; ~40x better). The remaining ~6% is a
convention detail in the Fortran grad_esp_qmmm charge-conserving operator
(the Phi_av (Tr[Q]-Nel) term, eqs 24-27) and the small STO-3G ESPF grid's
spurious charge; closing it needs Fortran-level work in qmmm.F90.
Gated behind Embedding="espf" and OFF by default: the existing electrostatic
and mechanical paths are unchanged (water QM/MM FD still 8e-4, connectivity
unit tests still pass).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
…forces Completes the full-ESPF covalent-boundary electrostatics. The embedded SCF carries only the electronic QM-MM coupling: finite differences of the embedded energy w.r.t. the field give dEqm/dphi_A = -Q_A (the electronic population), i.e. the nuclear-MM interaction sum_A Z_A phi_A is absent from eqm. In the full-ESPF scheme OpenMM holds no QM charge, so that nuclear term was dropped from the total energy, leaving the analytic force ~6% inconsistent with finite differences at the boundary (and the ESP charges' field response mis-split between electronic and net charge). Adding sum_A Z_A phi_A to the embedded energy makes the field response the full net QM charge (Z_A - Q_A = q_A), matching the analytic coupling force already computed. On the alanine C8-CA7 boundary the analytic force now agrees with finite differences to ~1e-3 on the significant forces (down from FD rel ~2.3 in the split scheme, and ~6% in the previous full-ESPF commit). The only residual (~2% on the smallest MM-host force) is the ESPF grid-resolution floor for the closely spaced frontier/link centres -- the same accuracy floor as the already-validated whole-molecule case (two waters: 1.2e-2 on the smallest force). espf_full also validated on whole-molecule QM (no boundary). Still gated behind Embedding="espf" and OFF by default; default electrostatic (water FD 8e-4) and mechanical paths and the connectivity unit tests unchanged. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
(1) Periodic support. The full-ESPF embedding potential and coupling force now
use the minimum-image convention under an orthorhombic periodic box
(_box_lengths_bohr / _min_image), so Embedding="espf" works with PME/periodic
QM/MM, not just NoCutoff. The QM-MM real-space electrostatics stay internally
consistent (potential and force use the same minimum image), so the analytic
force is finite-difference-exact under PBC too. Validated: water dimer in a
30 A box with PME, FD rel 8.6e-4.
(2) Wiring.
- QMMM_MD already consumes OpenQpQMMM.compute_force, so embedding="espf" works
with no code change. Verified end-to-end on the alanine C8-CA7 covalent
boundary (3-step NVE): total energy conserved to ~0.02 kJ/mol.
- NAMD_QMMM._total_force and NAMD_SOC_QMMM._total_force_soc gain a full-ESPF
branch (gated on driver.espf_full) mirroring the validated driver assembly:
pure MM-MM from OpenMM, nuclear-MM energy sum_A Z_A phi_A added to the
embedded active-state energy, and the analytic QM-MM coupling force on QM
and MM atoms (with link projection when present).
Scope: NAMD builds its QM mol from qm_atoms only (no link atoms), so NAMD
full-ESPF is validated for whole-molecule QM regions; a covalent boundary in
NAMD raises NotImplementedError (use QMMM_MD, or add link atoms to the NAMD QM
geometry -- a follow-up). Everything remains gated behind Embedding="espf" and
OFF by default: default electrostatic (water FD 8e-4) and mechanical paths and
the connectivity unit tests are unchanged.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
Embedding="electrostatic" (the default) now uses the full-ESPF scheme, which is finite-difference-exact for both whole-molecule and covalent-boundary QM regions (verified: water FD 8.6e-4; alanine C8-CA7 boundary ~1e-3). The legacy split scheme (QM charges routed through OpenMM point charges, only self- consistent for whole-molecule QM) remains available as Embedding="split" for reference/comparison. "mechanical" is unchanged. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SnqzJHvTchibeD8xd6yf4i
Resolve conflicts between the NAMD/SOC/QM/MM feature branch and main: - source/types.F90, include/oqp.h: union the control_parameters struct — keep main's performance knobs (xc_c2f, xc_phi_cache, xc_incdft, grad_cutoff, mrsf_resp_cutoff, mrsf_fp32, mrsf_zv_warmstart) and the branch's qmmm_flag, in the same field order in both files (cffi reads oqp.h; Fortran bind(c) must match). - source/modules/int1e.F90: take main's alloc_or_die allocation (allocates or reuses in place, so the Python-set OQP_mm_potential tagarray is preserved); drop the legacy remove_records(tags_general). - source/scf.F90: keep both the QM/MM hcore backup (hcore_bk) and main's IncDFT reset. - pyoqp/oqp/openqp.py: take main's rewritten module (canonical OpenQP class + helpers) and reconcile the legacy OPENQP wrapper the ESPF QM/MM driver depends on — accept an optional `silent` arg (OPENQP(cfg, True)) and expose `self.sp` (SinglePoint) for embedded SCF/excitation calls. - README.md: adopt main's table-based capabilities layout; re-add a "Dynamics & QM/MM" subsection for the NAMD/SOC-NAMD/QM/MM features. - .github/workflows/CI.yml: take main's example-validation / feature-coverage / test-suite gates (the branch's stale ddX-OFF note contradicted main's ENABLE_DDX=ON build). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
main's tagarray refactor removed container_t%remove_records and container_t%reserve_data (replaced by %erase and %alloc_or_die). The NAMD/QM/MM modules added by this branch still used the old API and failed to compile after merging main: Error: 'remove_records' is not a member of the 'container_t' structure Error: 'reserve_data' is not a member of the 'container_t' structure Migrate the affected calls (source/modules/namd.F90, source/modules/qmmm.F90) to alloc_or_die (allocates-or-reuses and binds the pointer in one call) and erase, matching the pattern main uses across the rest of the tree. The full tree now builds (liboqp.dylib + the cffi _oqp extension). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add a QM/MM convenience layer to the native `OpenQP` Python wrapper so ESPF
QM/MM and (SOC-)NAMD jobs can be built in the same section-style API as
`job.theory.mrsf(...)` / `job.workflow.soc(...)`:
- `job.qmmm(pdb_file=, forcefield=, qm_atoms=, cutoff=, embedding=, rigidwater=,
**[qmmm] keys)` sets `[input] qmmm_flag=true` and the `[qmmm]` section in one
call. `forcefield` is an alias for `forcefield_files`; a `forcefield` list is
joined to a comma string and a `qm_atoms` list to a space string.
- `job.workflow.namd(...)` selects `runtype=namd` and the `[md]` section
(MRSF-TDDFT only, like the other workflows). Pass `soc=True`/`soc_basis=...`
for SOC-NAMD; it composes with `job.qmmm(...)` for NAMD-QMMM.
Example:
job = OpenQP("chromophore")
job.molecule("system.pdb 0 1 2 3 4", basis="6-31g*")
job.theory.mrsf(functional="bhhlyp", nstate=3)
job.qmmm(pdb_file="system.pdb",
forcefield=["amber14-all.xml", "amber14/tip3p.xml"],
qm_atoms="0-4", cutoff="PME")
job.workflow.namd(soc=True, soc_basis="mch", nstep=200, dt=0.5, init_state="S1")
mol = job.run()
Also make the legacy `OPENQP.sp` (SinglePoint, used by the ESPF QM/MM driver) a
lazy property so importing `oqp.openqp` / constructing `OPENQP` no longer
requires the compiled library at that point. Adds 4 unit tests and the `[qmmm]`
/`[md]` sections to the API-test schema stub; full test_openqp_api.py passes
(39 tests).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…full suite Merging main enabled the example-coverage and full-suite gates, which the branch's new NAMD/QM/MM surface did not satisfy: - check_feature_coverage failed on 7 new opt-in bool flags with no example. Classify them (audited against their code usage): md.restart, md.econs, md.soc_du_dt_corr, md.soc_tdc_grad_corr are EXEMPT (IO restart / numerical stabilizer / diagnostic gradient sub-knobs); md.soc, md.dt_adaptive, and qmmm.rigidwater are KNOWN_UNCOVERED tracked gaps (real capabilities that need a runtype=namd / QM/MM example, PR #205). - run_tests all would then error on examples/QMMM/*.inp: those require the optional OpenMM backend plus external PDB/force-field files that the isolated per-example runner does not stage (PDB-not-found), so they are not self-contained regression tests. Exclude qmmm_flag=true examples from the full run (like numerical Hessians / IRC); they still run when invoked explicitly. Feature coverage now passes; HF/DFT/MRSF/SOC/SF example suites still pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add two small (nstep=1) NAMD-QMMM examples on formaldehyde solvated by 5 TIP3P waters, verified to run end-to-end: - H2CO-water_BHHLYP-MRSF-NAMD-QMMM.inp -- MRSF-TDDFT FSSH (internal conversion), ESPF QM/MM, NoCutoff. - H2CO-water_BHHLYP-SOC-NAMD-QMMM.inp -- SOC-NAMD (intersystem crossing) on the spin-adiabatic manifold, ESPF QM/MM. with the companion topology/force-field files (formaldehyde_water.pdb, formaldehyde.xml, tip3p.xml) and a README. These require the optional OpenMM backend plus the staged auxiliary files, and NAMD writes a trajectory log rather than a regression .json, so they are runnable demonstrations, not self-contained regression tests: they are skipped by `openqp --run_tests all` (qmmm_flag=true) and, having no committed .json, ignored by `--validate_examples`. The SOC example now exercises the md.soc feature flag, so md.soc is dropped from regression.KNOWN_UNCOVERED. Single-point QM/MM energy (runtype=energy + qmmm_flag) is not yet functional in this branch (the SCF's add_potqm_contributions needs OQP::ESPF_CORR, which is only created by the NAMD/QM-MM driver path), so no single-point QM/MM example is added here. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…no OpenMM)
Make the NAMD-QMMM examples first-class members of `openqp --run_tests all`
instead of excluding them:
- oqp_tester.run_single_test: when a run raises ModuleNotFoundError('openmm'),
report the example SKIPPED (like a build without an optional feature) rather
than ERROR, so a build without the optional OpenMM backend still yields a green
suite. SKIPPED does not affect the suite exit code.
- oqp_tester._skip_in_full_run: stop skipping runtype=namd QM/MM examples; keep
skipping the single-point / OpenMM-MD QM/MM decks (qmmm_flag without
runtype=namd), whose single-point ESPF path is not yet functional.
- namd.NAMD_QMMM: resolve the [qmmm] pdb_file and local force-field files
relative to the input file's directory when a bare name does not resolve
against the CWD, so the examples run from any working directory (e.g. the
per-example test runner). OpenMM built-in force fields are left untouched.
Verified: with OpenMM the two H2CO-water NAMD-QMMM examples PASS through the
tester; without OpenMM they report SKIPPED (suite status stays 0). Feature
coverage and test_openqp_api.py still pass.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…e build The ESPF QM/MM code (source/modules/qmmm.F90 espf_weights) called LAPACK dgels directly, bypassing the oqp_<routine>_i64 wrapper layer that every other BLAS/LAPACK call goes through. Consequences: - macOS Accelerate ILP64 build FAILED the check_accelerate_aliases gate: the raw dgels (and dgetrf, called from the oqp_dgetrf_i64 wrapper) were not interposed onto Accelerate's $NEWLAPACK$ILP64 interface -- an un-aliased LP64 symbol would silently run the 4-byte-integer BLAS and corrupt results. - On OQP_BLAS_INT=4 builds the raw call passed default (4-byte) integers to the ILP64 BLAS. Fix: - Add oqp_dgels_i64 to source/mathlib/lapack_wrap.F90 (mirrors oqp_dgeqrf_i64: converts the integer arguments to blas_int and calls dgels). - Rename dgels => oqp_dgels_i64 in source/mathlib/oqp_linalg.F90, so qmmm.F90's existing `call dgels(...)` transparently uses the wrapper (like dgesv/dgeqrf); no change needed in qmmm.F90. - Add dgels and dgetrf to OQP_ACCELERATE_ILP64_SYMS in cmake/oqp_functions.cmake so the raw symbols (still linked from inside the wrappers) are interposed onto Accelerate ILP64. Verified locally on macOS arm64 with LINALG_LIB=auto -> Apple Accelerate ILP64: liboqp.dylib links and "Accelerate ILP64 alias check OK: no LP64 BLAS/LAPACK leaks"; the H2CO-water NAMD-QMMM example runs to E_tot=-114.26343417, identical to the netlib build. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Lead the Dynamics & QM/MM functionality subsection with a description of SOC-NAMD-QMMM (excited-state MRSF-TDDFT surface hopping + spin-orbit intersystem crossing + ESPF QM/MM embedding in an explicit MM solvent), and link the SOC-NAMD-QMMM manual guide plus the compact job.qmmm/job.workflow.namd Python API. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… TB dispatch (dftb + xtb)
Reconcile the two diverged QM/MM lineages so method=dftb AND method=xtb get the
correct covalent (link-atom) QM/MM energy conservation.
- Take soc-namd-options's full-ESPF QM/MM driver internals (full-field ESPF with
per-centre MM exclusions, nuclear-MM -> FD-exact covalent-boundary forces, PBC;
frontier-charge redistribution removed).
- Re-apply the TB dispatch (oqp.utils.tb_backends: is_tb_method/tb_section_name/
tb_config/make_tb_adapter) on top of the full-ESPF driver:
* qmmm_driver.OpenQpQMMM._forces_qm mol/config dispatch + _forces_qm_dftb
(sets mol.dftb_external_potential=POTMM, reads relaxed OQP::partial_charges,
skips the ESPF hcore/grad mutation; no nuclear term for TB).
* namd NAMD_QMMM electronic/gradient/energy TB arms + SOC(-QMMM) TB arms and
the _dftb_soc_tags/_dftb_spatial_overlap helpers.
* runfunc.compute_soc, single_point (energy/reference/excitation/gradient/
state-overlap), input_checker NAMD+QM/MM gates, oqpdata [dftb]/[xtb] schema,
openqp [dftb]/[xtb] API, oqp_tester skip handling.
- Keep the DFTB/XTB adapters (openqp_dftb.py, openqp_xtb.py, tb_backends.py).
- Remove stale frontier_scheme artifacts superseded by full-ESPF
(test_qmmm_frontier[_openmm].py, ala-dipeptide RCD example).
Static verification: all changed .py parse; TB schema/dispatch/interface unit
tests pass (54 passed, 20 skipped); full suite shows the same 51 pre-existing
(liboqp-gated) failures as the pre-merge baseline (no new failures). Covalent-
boundary energy-vs-force FD for dftb/xtb remains to be run on the integrated
runtime; see docs/qmmm-tb-reconciliation.md.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: b217f46d57
ℹ️ 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".
| for _k in list(qm_cfg): | ||
| if str(_k).split('.')[-1].strip().lower() == 'runtype': | ||
| qm_cfg[_k] = 'energy' | ||
| self.oqp_cfg = qm_cfg |
There was a problem hiding this comment.
When users run the ground-state QM/MM MD path via an input file with [input] runtype=md, pyoqp.py dispatches directly to QMMM_MD, and this assignment now preserves input.runtype=md in qm_cfg. The first compute_force() then constructs OPENQP(self.oqp_cfg_base), whose check_input_values() rejects md as a recognized-but-not-implemented runtype before any MD step can run. The deleted normalization to energy was needed because each QM subcall is only an energy/gradient evaluation inside the outer MD driver.
Useful? React with 👍 / 👎.
| # 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) |
There was a problem hiding this comment.
For an out-of-order qm_atoms selection, the QM geometry is still built by iterating self.topology.atoms() in topology order, but gradients/charges are scattered back by enumerating self.qm_atoms and link-atom host_row is computed from this list. Preserving the caller order here (and again in compute_force) therefore sends QM forces and link projections to the wrong atoms for inputs such as qm_atoms = 16,8,18,9,17; the existing tests/test_qmmm_qm_atoms_order.py documents that the driver must sort this list.
Useful? React with 👍 / 👎.
Comment on lines
+7
to
+9
| 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. |
There was a problem hiding this comment.
This change explicitly introduces method=dftb/method=xtb QM/MM at covalent link-atom boundaries, but the diff does not add any examples/** input that exercises that new user-facing path (and it deletes the only dedicated covalent-boundary QMMM-MD deck). AGENTS.md rule 2 requires new capabilities to ship with a small, fast example, so users and CI have no example-level coverage for the feature this reconciliation enables.
Useful? React with 👍 / 👎.
Comment on lines
+7
to
+9
| 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. |
There was a problem hiding this comment.
Because this change advertises new user-facing method=dftb/method=xtb QM/MM support at covalent boundaries, AGENTS.md rule 4 requires the PR description to link the companion openqp-docs manual update. No such link is present in the submitted description, so the keyword/workflow manual can drift from the behavior enabled here.
Useful? React with 👍 / 👎.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Reconciles the two diverged OpenQP lineages so that both
method=dftbandmethod=xtbget energy-conserving covalent-boundary (link-atom) QM/MM. Merges the full-ESPF QM/MM work fromcodex/soc-namd-options(FD-exact covalent-boundary forces, PBC, full-ESPF default) onto the TB-adapter branch (#269), re-applying the TB method dispatch on top of the newer full-ESPF driver.Stacked on #269 (base =
feat/openqp-xtb-adapter); the net delta here is the full-ESPF QM/MM driver + its TB (dftb/xtb) dispatch. Retarget tomainafter #269 merges.The problem this fixes
upstream/mainshipped the DFTB backend (#266) but its released ESPF QM/MM does not conserve energy at a covalent QM/MM boundary (force ≠ −dE/dR at link atoms). The fix — full-ESPF with FD-exact covalent-boundary forces — lived only oncodex/soc-namd-options, which never had the DFTB/xTB adapters. Neither branch had both. This PR is the reconciliation.Runtime validation (integrated OpenQP + TB libraries)
Built liboqp from this branch +
libopenqp_xtb_c(openqp-xtba00d5db) +libopenqp_dftb_c(openqp-dftbmain) and ran the definitive checks:Covalent link-atom QM/MM energy-vs-force FD (central FD of the total ESPF QM/MM energy vs projected analytic force; ethane CH3/CH3 boundary, OpenMM):
The repo's own authoritative arbiter
tests/test_qmmm_dftb.pypasses all 3 tests (ethane link-atom + H2/water ground + H2/water MRSF) with the stable DFTB library — direct proof the full-ESPF reconciliation fixes the DFTB covalent-boundary bug.Gas-phase MRSF-xTB FSSH NVE (ethylene, velocity-Verlet, no thermostat): no secular drift; conservation error scales O(dt²) (peak-peak ×4.8 for dt 1→2 fs); ~0.001 eV/ps at dt=0.5 fs — the analytic MRSF-xTB force is F = −dE/dR.
Scope / caveats
_forces_qm_dftb(serves both dftb and xtb via the shared adapter +mol.dftb_external_potentialhandoff), gated byespf_full.docs/qmmm-tb-reconciliation.md.🤖 Generated with Claude Code