Noneq-MQC-QComp
Nonequilibrium Quantum-Classical Dynamics on Quantum Computers
Research Proposal for Unitary Foundation Microgrant
Motivation
Simulating nonadiabatic molecular dynamics at electrochemical interfaces requires treating coupled electron-nuclear motion where the Born-Oppenheimer approximation breaks down. Classical approaches combining the Mapping Approach to Surface Hopping (MASH) [1] with Nonequilibrium Green’s Function Density Functional Theory (NEGF-DFT) [2] can capture this physics, but electronic structure calculations at every timestep scale as \(O(N^3)\) to \(O(N^6)\), making realistic simulations (100+ atoms, microsecond timescales) intractable.
Recent quantum algorithms for nonadiabatic dynamics@ollitrault2020nonadiabatic, [3] show polynomial scaling, but existing work focuses on isolated molecules in Born-Oppenheimer or simple nonadiabatic regimes. Work by Li et al.@li2025multiset on multi-set VQD (MS-VQD) and Walters et al.@walters2025variational on ensemble-averaged Ehrenfest trajectories provides building blocks, but no implementation exists for electrochemical systems with applied voltage and continuous electron exchange with reservoirs.
This 6-month project delivers a proof-of-concept implementation of MS-VQD for MASH-NEGF dynamics, validated on the spin-boson model and two-site Hubbard model, with working code released as open-source contribution to Qiskit Nature.
Research Program
Core Algorithm: MS-VQD for MASH Dynamics
The wavefunction ansatz employs separate parameterized circuits for each electronic state:
\[|\Psi(t)\rangle = \sum_{i=1}^{N_{\text{states}}} c_i(t) U_i(\boldsymbol{\theta}_i(t)) |\phi_i^{\text{nuc}}\rangle \otimes |i\rangle_e\]
where \(U_i(\boldsymbol{\theta}_i(t))\) is a hardware-efficient ansatz with depth \(L_i\) optimized separately for each state \(i\), \(|\phi_i^{\text{nuc}}\rangle\) encodes nuclear configuration via basis function representation on \(n_{\text{nuc}} = \log_2(N_{\text{grid}})\) qubits, and \(c_i(t)\) are classical electronic amplitudes evolved via MASH equations. For nuclear DOFs, we use discrete variable representation (DVR) where position eigenstates \(|x_k\rangle\) map to computational basis states \(|k\rangle\), enabling momentum operators via discrete Fourier transform and kinetic energy via \(\hat{T} = -\frac{\hbar^2}{2m}\hat{p}^2\).
MASH Surface Hopping on Quantum Circuits
MASH maps electronic populations to action variables \((n_i, q_i)\) satisfying
\[\dot{n}_i = -\sum_{j \neq i} \frac{2}{\hbar}\text{Im}[H_{ij}\sqrt{n_i n_j}e^{i(q_i - q_j)}]\]
\[\dot{q}_i = \frac{H_{ii}}{\hbar} - \sum_{j \neq i} \frac{\text{Re}[H_{ij}e^{i(q_i - q_j)}]}{\hbar}\sqrt{\frac{n_j}{n_i}}\]
where \(H_{ij}\) are matrix elements in the adiabatic basis. At each timestep \(\Delta t\): (1) evolve nuclear DOFs via \(U_i(\boldsymbol{\theta}_i(t))\) under forces \(\boldsymbol{F}_i = -\nabla_{\boldsymbol{R}} H_{ii}\), (2) measure electronic qubits to obtain populations \(\{n_i\}\), (3) compute hopping probabilities \(g_{i \to j} = \max(0, -2\dot{n}_i \Delta t / n_i)\), (4) accept hops stochastically and rescale nuclear momenta to conserve energy. Circuit depth per timestep scales as \(O(L_{\text{max}} \cdot N_{\text{states}})\) where \(L_{\text{max}} = \max_i L_i\), compared to \(O(L \cdot 2^{N_{\text{states}}})\) for single-circuit VQD.
NEGF Extension for Applied Voltage
For electrochemical systems, the molecular Hamiltonian couples to electrode reservoirs characterized by self-energies \(\Sigma^r_L(E)\) and \(\Sigma^r_R(E)\) at chemical potentials \(\mu_L = E_F + eV_{\text{bias}}/2\) and \(\mu_R = E_F - eV_{\text{bias}}/2\). The retarded Green’s function satisfies
\[G^r(E) = [E - H_{\text{mol}} - \Sigma^r_L(E) - \Sigma^r_R(E)]^{-1}\]
yielding bias-dependent forces \(\boldsymbol{F}_i(\boldsymbol{R}, V_{\text{bias}}) = -\nabla_{\boldsymbol{R}}\text{Tr}[G^r H_{\text{mol}}]\) used in step (1). For the microgrant scope, we implement NEGF with wide-band limit approximation \(\Sigma^r_{L/R}(E) \approx -i\Gamma_{L/R}/2\) (constant coupling), sufficient for validation on molecular junction benchmarks. Electronic structure inputs \((H_{\text{mol}}, \Gamma_{L/R})\) are precomputed classically via GPAW or PySCF.
Concrete Implementation: Spin-Boson Model
The spin-boson Hamiltonian reads
\[H = -\frac{\Delta}{2}\sigma_x - \frac{\epsilon}{2}\sigma_z + \sum_j \left[\frac{p_j^2}{2m_j} + \frac{1}{2}m_j\omega_j^2\left(x_j - \frac{c_j}{m_j\omega_j^2}\sigma_z\right)^2\right]\]
with system-bath coupling characterized by spectral density \(J(\omega) = \sum_j \frac{c_j^2}{2m_j\omega_j}\delta(\omega - \omega_j)\). For MS-VQD implementation: (1) map \(\sigma_z = \pm 1\) states to 1 electronic qubit, (2) represent each bath mode \(x_j\) on \(n_{\text{bath}} = 6\) qubits (64-point DVR grid spanning \(\pm 5\sigma\) where \(\sigma^2 = \hbar/(2m_j\omega_j)\)), (3) use \(N_{\text{bath}} = 5\) modes with frequencies \(\{\omega_j\}\) and couplings \(\{c_j\}\) chosen to discretize Ohmic spectral density \(J(\omega) = \eta\omega e^{-\omega/\omega_c}\). Total: 31 qubits (1 electronic + 30 bath). Circuit depth per timestep: \(L = 20\) layers (10 rotation layers + 10 entangling layers) using RealAmplitudes ansatz in Qiskit.
Validation compares against exact path integral (QuAPI) for parameters \(\Delta = 1.0\), \(\epsilon = 0\), \(\eta = 0.1\), \(\omega_c = 2.5\), \(\beta = 5.0\) (units where \(\hbar = 1\)). Observables: population dynamics \(P_{\uparrow}(t) = \langle\frac{1 + \sigma_z}{2}\rangle(t)\) and decoherence rate \(\Gamma = -\lim_{t\to\infty}\frac{1}{t}\ln|P_{\uparrow}(t) - 0.5|\).
Validation: Two-Site Hubbard Model
\[H = -t\sum_{\sigma}(c^\dagger_{1\sigma}c_{2\sigma} + \text{h.c.}) + U\sum_i n_{i\uparrow}n_{i\downarrow} + \sum_{i,\sigma}\epsilon_i n_{i\sigma}\]
coupled to harmonic bath at each site. Implement as: 4 electronic qubits (2 sites × 2 spins), 2 nuclear DOFs (site separation \(R\) on 6 qubits, relative angle \(\theta\) on 6 qubits). Parameters: \(t = 1.0\), \(U = 4.0\), \(\epsilon_1 = -\epsilon_2 = 1.0\), bath coupling \(\alpha = 0.2\). Compare electron transfer rates against numerical renormalization group (NRG) results from Ref. [4]. Metrics: pair transfer probability \(P_{\text{pair}}(t) = \langle n_{2\uparrow}n_{2\downarrow}\rangle(t)\) vs sequential transfer.
Timeline
Month 1: Infrastructure Implement DVR encoding for nuclear DOFs in Qiskit (300 lines). Implement MASH equations for action-angle variables with stochastic hopping (200 lines). Unit tests against analytic solutions for 2-level system.
Month 2: Spin-Boson Construct MS-VQD circuits for spin-boson with 5 bath modes. Optimize circuit depth via COBYLA on IBM Quantum simulator. Generate QuAPI reference data using qutip or quapi package. Validate \(P_{\uparrow}(t)\) with target error \(<10\%\) for \(t \leq 10/\omega_c\).
Month 3: Hubbard Model Extend to 4-state Hubbard model with 2 nuclear DOFs. Implement energy-conserving momentum rescaling for surface hops. Compare against NRG data for pair vs sequential transfer. Target circuit depth \(L \leq 30\) layers.
Month 4: NEGF Integration Implement wide-band NEGF with precomputed \(H_{\text{mol}}\) from PySCF. Test on H₂ molecular junction with Au electrodes (\(V_{\text{bias}} = 0.1\)-\(1.0\) V). Validate conductance \(G(V) = \frac{2e^2}{h}\text{Tr}[\Gamma_L G^r \Gamma_R G^a]\) against literature.
Month 5: Hardware Testing Run optimized circuits on IBM Quantum hardware (127-qubit Eagle processor via Unitary Foundation QPU access). Apply error mitigation (zero-noise extrapolation + readout correction). Document noise impact on hopping probabilities.
Month 6: Release & Documentation Package code as qiskit-namd module with tutorials. Write paper draft comparing circuit depth scaling (MS-VQD vs single-circuit VQD) and accuracy (VQD vs exact). Submit to Quantum journal (open access). Present at Qiskit Advocate workshop.
Funding Breakdown ($15,000 total):
- IBM Quantum Premium access (3 months): $6,000
- AWS Braket for larger simulations: $2,000
- Classical HPC for QuAPI/NRG reference (50k CPU-hrs @ $0.04/hr): $2,000
- Open access publication (Quantum journal): $1,500
- Conference registration + travel (Qiskit Advocate workshop): $2,500
- Hardware (workstation GPU for circuit optimization): $1,000
Outcome
Primary Deliverable: qiskit-namd Package
Open-source Python package on GitHub providing:
DVREncodingclass: Maps nuclear DOFs to qubit basis via discrete variable representationMASHCircuitclass: Constructs MS-VQD ansatz with separate circuits per electronic state
SurfaceHoppingclass: Implements stochastic hopping with energy-conserving momentum rescalingNEGFIntegrationclass: Computes bias-dependent forces from wide-band NEGF- Tutorial notebooks: Spin-boson model, Hubbard model, H₂ junction
Target: 2000 lines of documented code with 80% test coverage. Modular design allows users to swap ansatzes (RealAmplitudes, EfficientSU2, custom) or NEGF implementations (wide-band → full frequency-dependent via future extension).
Benchmark Dataset
Zenodo repository containing:
- Spin-boson: Exact QuAPI dynamics for 20 parameter sets varying \((\eta, \omega_c, \beta)\)
- Hubbard: NRG transfer rates for 10 parameter sets varying \((U/t, \alpha)\)
- H₂ junction: Experimental conductance vs \(V_{\text{bias}}\) from literature + our VQD results
- Circuit specifications: Optimal depths \(\{L_i\}\) achieving \(<5\%\) error, gate counts, CNOT counts
Scientific Validation
Paper draft (12 pages + appendices) demonstrating:
- Circuit depth reduction: MS-VQD requires \(L_i \sim 20\) per state vs single-circuit VQD requiring \(L \sim 60\) for spin-boson, achieving 66% reduction consistent with Li et al.’s@li2025multiset 50% target
- Accuracy: MS-VQD with MASH hopping reproduces exact \(P_{\uparrow}(t)\) within 5% error vs 15% error for Ehrenfest (no hopping)
- Scalability: Hardware demonstration on 31-qubit spin-boson shows feasible noise levels with error mitigation
- NEGF integration: H₂ conductance matches experimental \(G(0.5\text{ V}) \approx 0.02 \times 2e^2/h\) within factor of 2
Ecosystem Impact
Contribution to Qiskit Nature enables computational chemists to explore nonadiabatic dynamics without quantum computing expertise, quantum algorithm developers to test new ansatzes on realistic chemical problems, and the electrochemistry community to benchmark against first-principles dynamics.
Differentiation from existing Qiskit Nature modules: Current modules focus on ground-state VQE or time-independent dynamics. Our contribution adds time-dependent nonadiabatic dynamics with surface hopping, filling gap identified in Qiskit Nature roadmap discussions (GitHub issues #150, #287).
Why This Project Needs Microgrant Funding
IBM Quantum Premium access ($6k) is essential—free tier limits runtime to 10 minutes per job, insufficient for converged MS-VQD optimization requiring 50+ iterations × 5 minutes each. Classical HPC for exact reference data ensures rigorous validation but requires dedicated allocation. Open access publication ensures broader community adoption critical for ecosystem impact. Without this funding, project reduces to toy demonstrations on simulators without hardware validation or proper benchmarking.
Follow-on Work (Beyond Grant Scope)
Success of this 6-month project enables future extensions: full frequency-dependent NEGF (removes wide-band approximation), non-Markovian baths via ensemble-averaged trajectories, larger systems (3-site Hubbard with 8 qubits, benzene-Au junction with 15 qubits), and integration with battery modeling (PyBAMM) for Li⁺ intercalation dynamics. Each extension builds on core qiskit-namd infrastructure delivered by microgrant.