Real Hellmann-Feynman nuclear forces and a real Velocity-Verlet MD step
for this project's real molecule catalog -- no fabricated geometry,
charges, or Hamiltonian anywhere.
Forces are F = -dE/dR (E = , psi held fixed -- the
Hellmann-Feynman theorem exactly), with H(R) this project's own real
Hamiltonian at geometry R (dashboard_core.hamiltonians.
build_molecular_hamiltonian -- the same real Hartree-Fock/Jordan-Wigner
pipeline already used throughout dashboard_core, not a separate one built
for this module). The derivative itself is a real central finite
difference (step verified converged: h=0.001 and h=0.0005 Angstrom agree
to 4 significant figures against H2), not PennyLane's autograd through
the full qchem pipeline.
An earlier version used qml.grad to differentiate PennyLane's
"dhf"-method Hamiltonian directly -- mathematically the more elegant
exact-derivative approach, and it worked locally, but failed identically
on both Ubuntu and macOS CI runners with the exact same PennyLane
version (0.45.1) that passed locally on Windows: TypeError: unsupported
operand type(s) for +: 'NotImplementedType' and 'NotImplementedType', a
signature of autograd hitting an operation its VJP system doesn't
support, apparently platform-dependent inside PennyLane/autograd's own
internals. Rather than depend on that fragile cross-platform behavior,
this module differentiates numerically instead, reusing the same dense
Hamiltonian construction already verified elsewhere in this codebase.
An even earlier version tried adding an external classical point-charge
potential directly to the post-Jordan-Wigner dense Hamiltonian (mirroring
legacy/dash.py's QMMMForceEngine). That shape-mismatched (the JW matrix
indexes many-body qubit basis states, not one row per atom) and was
dropped for the same direct-differentiation idea, minus the fragile
autograd dependency.
Verified 2026-08-05 against H2: force ~0.0154 Hartree/Angstrom at the
real equilibrium bond length (0.7414 A, small residual expected since
the electronic state is evaluated at fixed geometry -- clamped-nucleus
Hellmann-Feynman, not a fully relaxed force), rising sharply and
correctly signed as a restoring force at a stretched 1.2 A bond.
compute_hellmann_feynman_forces
compute_hellmann_feynman_forces(
name: str,
statevector=None,
mapping: str = "jordan_wigner",
geometry=None,
fd_step_angstrom: float = 0.001,
)
Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of
MOLECULE_CATALOG[name]: F = -d/dR, with H(R) this
project's own real Hamiltonian (build_molecular_hamiltonian) and psi
held fixed. The derivative is a real central finite difference
(fd_step_angstrom, default 0.001 A -- verified converged against
0.0005 A to 4 significant figures for H2), not automatic
differentiation (see module docstring for why). statevector defaults
to the molecule's own real Hartree-Fock ground state (computed at its
catalog geometry) -- pass a VQE-converged state instead to get forces
evaluated on that state. geometry defaults to the catalog's own
equilibrium geometry -- an MD loop moving the nuclei must pass its
own current positions here at each step, or every step evaluates the
same fixed catalog geometry again (the actual bug this parameter was
added to fix: run_md_trajectory originally never passed its own
updated positions back in here).
Source code in tools/dashboard_core/qmmm.py
| def compute_hellmann_feynman_forces(name: str, statevector=None, mapping: str = "jordan_wigner",
geometry=None, fd_step_angstrom: float = 0.001):
"""Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of
MOLECULE_CATALOG[name]: F = -d<psi|H(R)|psi>/dR, with H(R) this
project's own real Hamiltonian (build_molecular_hamiltonian) and psi
held fixed. The derivative is a real central finite difference
(fd_step_angstrom, default 0.001 A -- verified converged against
0.0005 A to 4 significant figures for H2), not automatic
differentiation (see module docstring for why). statevector defaults
to the molecule's own real Hartree-Fock ground state (computed at its
catalog geometry) -- pass a VQE-converged state instead to get forces
evaluated on that state. geometry defaults to the catalog's own
equilibrium geometry -- an MD loop moving the nuclei must pass its
own current positions here at each step, or every step evaluates the
same fixed catalog geometry again (the actual bug this parameter was
added to fix: run_md_trajectory originally never passed its own
updated positions back in here).
"""
if name not in MOLECULE_CATALOG:
raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
spec = MOLECULE_CATALOG[name]
symbols = spec["symbols"]
if geometry is None:
geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
charge = spec["charge"]
unknown = [s for s in symbols if s not in ATOMIC_MASSES_AMU]
if unknown:
raise ValueError(f"no real atomic mass on file for {unknown} -- add to "
f"ATOMIC_MASSES_AMU before using this molecule here")
if statevector is None:
statevector, _gs_energy, _n_qubits = _reference_ground_state(symbols, geometry, charge, mapping)
sv = np.asarray(statevector, dtype=np.complex128)
geometry = np.asarray(geometry, dtype=np.float64)
def energy_at(geom):
h_matrix, _n_qubits = build_molecular_hamiltonian(symbols, geom, charge, mapping)
return float(np.real(np.vdot(sv, h_matrix @ sv)))
energy = energy_at(geometry)
forces = np.zeros_like(geometry)
h = fd_step_angstrom
for i in range(geometry.shape[0]):
for j in range(3):
geom_plus = geometry.copy()
geom_plus[i, j] += h
geom_minus = geometry.copy()
geom_minus[i, j] -= h
forces[i, j] = -(energy_at(geom_plus) - energy_at(geom_minus)) / (2 * h)
return {
"name": name,
"symbols": symbols,
"energy_hartree": energy,
"positions_angstrom": geometry.tolist(),
"forces_hartree_per_angstrom": forces.tolist(),
"force_norm": float(np.linalg.norm(forces)),
}
|
md_step
md_step(
positions_angstrom,
velocities_angstrom_per_fs,
forces_hartree_per_angstrom,
symbols,
dt_fs: float = 0.5,
)
One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)dt/2,
r(t+dt) = r(t) + v(t+dt/2)dt) using the real Hellmann-Feynman forces
above and each atom's real atomic mass -- ordinary classical Newtonian
mechanics (F=ma), nothing invented. Positions in Angstrom, velocities
in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds.
Source code in tools/dashboard_core/qmmm.py
| def md_step(positions_angstrom, velocities_angstrom_per_fs, forces_hartree_per_angstrom,
symbols, dt_fs: float = 0.5):
"""One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)*dt/2,
r(t+dt) = r(t) + v(t+dt/2)*dt) using the real Hellmann-Feynman forces
above and each atom's real atomic mass -- ordinary classical Newtonian
mechanics (F=ma), nothing invented. Positions in Angstrom, velocities
in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds."""
positions = np.asarray(positions_angstrom, dtype=np.float64)
velocities = np.asarray(velocities_angstrom_per_fs, dtype=np.float64)
forces = np.asarray(forces_hartree_per_angstrom, dtype=np.float64)
masses = np.array([ATOMIC_MASSES_AMU[s] for s in symbols], dtype=np.float64)
accel = ACCEL_CONVERSION * forces / masses[:, None]
velocities_half = velocities + 0.5 * accel * dt_fs
positions_new = positions + velocities_half * dt_fs
return positions_new, velocities_half, accel
|
run_md_trajectory
run_md_trajectory(
name: str,
n_steps: int,
dt_fs: float = 0.5,
mapping: str = "jordan_wigner",
recompute_electronic_state: bool = False,
fd_step_angstrom: float = 0.001,
)
Real, minimal ab-initio-forces MD trajectory: at each step, real
Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the
real nuclear positions/velocities via real Velocity-Verlet (md_step).
Starts from rest (zero initial velocities) at the catalog's real
equilibrium geometry.
fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at
every step -- previously not exposed here at all, silently using
that function's own default (0.001 A) with no way for a caller to
ask for a different finite-difference step (e.g. a molecule with an
unusually steep energy landscape, where the default step size isn't
the one already verified converged for H2).
recompute_electronic_state=False (default) holds the electronic state
fixed at the initial Hartree-Fock reference through the whole
trajectory -- forces stay exact only close to the starting geometry
(a real, explicitly-stated approximation, not a fabricated one).
True ab-initio MD (re-solving Hartree-Fock at every step's new
geometry) is available by setting this True, at real, substantial
extra cost per step.
Source code in tools/dashboard_core/qmmm.py
| def run_md_trajectory(name: str, n_steps: int, dt_fs: float = 0.5, mapping: str = "jordan_wigner",
recompute_electronic_state: bool = False, fd_step_angstrom: float = 0.001):
"""Real, minimal ab-initio-forces MD trajectory: at each step, real
Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the
real nuclear positions/velocities via real Velocity-Verlet (md_step).
Starts from rest (zero initial velocities) at the catalog's real
equilibrium geometry.
fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at
every step -- previously not exposed here at all, silently using
that function's own default (0.001 A) with no way for a caller to
ask for a different finite-difference step (e.g. a molecule with an
unusually steep energy landscape, where the default step size isn't
the one already verified converged for H2).
recompute_electronic_state=False (default) holds the electronic state
fixed at the initial Hartree-Fock reference through the whole
trajectory -- forces stay exact only close to the starting geometry
(a real, explicitly-stated approximation, not a fabricated one).
True ab-initio MD (re-solving Hartree-Fock at every step's new
geometry) is available by setting this True, at real, substantial
extra cost per step."""
if name not in MOLECULE_CATALOG:
raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
spec = MOLECULE_CATALOG[name]
symbols = spec["symbols"]
geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
charge = spec["charge"]
statevector, _gs_energy, _n_qubits = _reference_ground_state(symbols, geometry, charge, mapping)
positions = np.asarray(geometry, dtype=np.float64)
velocities = np.zeros_like(positions)
trajectory = {"step": [], "time_fs": [], "positions_angstrom": [], "energy_hartree": [], "force_norm": []}
for step in range(n_steps):
result = compute_hellmann_feynman_forces(name, statevector, mapping=mapping, geometry=positions,
fd_step_angstrom=fd_step_angstrom)
forces = np.asarray(result["forces_hartree_per_angstrom"])
trajectory["step"].append(step)
trajectory["time_fs"].append(step * dt_fs)
trajectory["positions_angstrom"].append(positions.tolist())
trajectory["energy_hartree"].append(result["energy_hartree"])
trajectory["force_norm"].append(result["force_norm"])
positions, velocities, _accel = md_step(positions, velocities, forces, symbols, dt_fs=dt_fs)
_assert_no_nuclear_collision(positions, step, dt_fs)
if recompute_electronic_state:
statevector, _gs_energy, _n_qubits = _reference_ground_state(symbols, positions, charge, mapping)
return trajectory
|