Quantum-classical hybrid computation of electron transfer in a cryptochrome protein via VQE-PDFT and multiscale modeling

Abstract

Accurate calculation of strongly correlated electronic systems requires proper treatment of both static and dynamic correlations, which remains challenging for conventional methods. To address this, we present VQE-PDFT, a quantum-classical hybrid framework that integrates the variational quantum eigensolver with multiconfiguration pair-density functional theory (MC-PDFT). This framework strategically employs quantum circuits for multiconfigurational wavefunction representation while utilizing density functionals for correlation energy evaluation. The hybrid strategy maintains accurate treatment of static and dynamic correlations while reducing quantum resource requirements compared to highly expressive quantum algorithms. Benchmark validation, performed via a noiseless quantum circuit simulator, on the charge-transfer dataset confirmed that VQE-PDFT achieved results comparable to conventional MC-PDFT. Building upon this, we developed shallow-depth hardware-efficient ansatz circuits and integrated them into a QM/MM multiscale architecture to enable applications in complex biological systems. This extended framework, when applied to electron transfer in the European robin cryptochrome protein ErCRY4 with noiseless simulations, yielded transfer rates that aligned well with experimental measurements. Finally, as a proof-of-concept hardware demonstration, we executed reduced-density-matrix measurements for a single protein conformation on a 13-qubit superconducting device and showed the impact of noise through a comprehensive error analysis.

---

1 Introduction

Evaluating strongly correlated systems is one of the most fundamental challenges in modern science, where conventional theoretical approaches often struggle due to the intricate electronic coupling. The accurate description of these systems has led to major breakthroughs across multiple disciplines: in physics, it has driven the understanding of high-temperature superconductivity;¹ in chemistry, it has enabled the precise modeling of catalytic processes;² in biology, it has provided crucial insights into the mechanisms of nitrogen fixation³ and avian magnetoreception.^4,5

From an electronic structure perspective, the challenge of evaluating strongly correlated systems lies in representing the multiconfigurational wavefunction.⁶ Addressing this problem has spurred the advent of several computational strategies, each with distinct advantages and limitations. Within the multiconfiguration self-consistent field (MCSCF) framework, complete active space configuration interaction (CASCI)⁷ provides a foundation by performing configuration interaction within a predefined active space. Building upon this foundation, the complete active space self-consistent-field (CASSCF) method enhances CASCI by simultaneously optimizing both the CI coefficients and orbital parameters, thereby more accurately capturing static correlation effects—the near-degeneracy effects arising from multiple electronic configurations of comparable energy.^8,9

However, both CASCI and CASSCF often provide insufficient accuracy due to their incomplete treatment of dynamic correlation effects^8,10—the instantaneous electron–electron repulsion effects.¹¹ To remedy this deficiency, researchers have developed various post-SCF approaches that build upon MCSCF reference wavefunctions.¹² Complete active space second-order perturbation theory (CASPT2)^13,14 represents one of the most successful examples, systematically recovering dynamic correlation, though it suffers from computational scaling limitations.

To address these computational limitations, multiconfiguration pair-density functional theory (MC-PDFT)^15,16 offers a promising alternative by combining CASSCF-derived density matrices with an on-top density functional that captures dynamic correlation through opposite-spin electron pair probabilities at identical spatial coordinates. This hybrid approach greatly improves computational scaling relative to traditional post-SCF methods.^6,17

However, MC-PDFT''s efficiency remains fundamentally constrained by its underlying CASSCF calculation. The exponential scaling of CASSCF with active space size renders calculations for large chemical systems computationally prohibitive, even with state-of-the-art hardware and advanced approximation schemes.¹⁸

The emergence of quantum computing presents a potential solution to these exponential scaling challenges.^19,20 Quantum algorithms naturally exploit superposition and entanglement to encode complex many-body wavefunctions, offering theoretical advantages for strongly correlated electronic systems.²¹ Meanwhile, current noisy intermediate-scale quantum (NISQ) devices have motivated the development of hybrid quantum-classical algorithms,²² most notably the variational quantum eigensolver (VQE)²³ and its variants. These include ADAPT-VQE,^24,25 which dynamically constructs quantum circuits through iterative operator selection; cluster-VQE,²⁶ which reduces qubit requirements via system decomposition; DMET-ESVQE,²⁷ which achieves the combined advantages of both circuit optimization and qubit reduction through a distinct embedding-based approach; and ansatz design to maintain physical symmetries including quantum-number-preserving constraints.²⁸ Alternative frameworks such as quantum algorithms for density functional theory (DFT) have also been proposed.²⁹

For large systems, spatial decomposition approaches offer a practical solution by partitioning the system into fragments amenable to quantum hardware. Recent methods combine fragment molecular orbital techniques with quantum embedding³⁰ and utilize localized active space strategies to reduce quantum resource requirements through systematic treatment of inter-fragment correlations.^31,32

Despite these advances, extending VQE-based methods to multiconfigurational systems remains problematic.³³ Existing approaches such as UCCGSD, and k-UpCCGSD struggle with the simultaneous treatment of static and dynamic correlation,^34,35 often requiring prohibitively deep quantum circuits that exceed current NISQ capabilities. Moreover, the qubit counts on near-term devices typically limit such ansatzes to compact active spaces, and dynamic correlation contributions from orbitals outside the active space remain missing. Besides, the inherent noise in NISQ devices raises fundamental questions about computational reliability^36,37 for practical applications, necessitating error analysis in specific calculations. Therefore, the fundamental challenge lies in designing quantum algorithms that can efficiently capture multiconfigurational character while remaining compatible with near-term quantum hardware constraints.

An alternative paradigm leverages quantum devices to compute reduced density matrices (RDMs), delegating correlation energy recovery to classical post-processing. Notable implementations range from coupling self-consistent CASSCF with hardware-efficient sampling to mitigate measurement noise in small active spaces,³⁸ to integrating VQE with adiabatic connection theory for systematic dynamic correlation recovery from orbitals outside the chosen active space.³⁹ Furthermore, Boyn et al. developed a hybrid framework in which a quantum solution of the anti-Hermitian contracted Schrödinger equation directly yields N-representable 2-RDMs, which are then combined with multiconfiguration pair-density functional theory for classical correlation-energy evaluation.⁴⁰

Collectively, the above RDM-based approaches support the viability of separating quantum wavefunction preparation from classical correlation recovery. Building upon this philosophy, we introduce a quantum–classical hybrid framework that integrates VQE with multiconfiguration pair-density functional theory (MC-PDFT), referred to as VQE-PDFT. This framework (described first below) then serves as the basis for two additional developments targeting realistic biological modeling and practical hardware execution.

First, at the methodological level, we employ VQE strictly as a CASCI solver to capture static correlation, while recovering additional dynamic correlation beyond the CASCI description via a classical MC-PDFT on-top functional. In this framework, the total energy is evaluated from a pair-density functional of CASCI-level RDMs rather than directly from the Hamiltonian expectation value, allowing the quantum computation in VQE-PDFT to be reduced to compact active spaces focused on static correlation. By contrast, bare-〈H〉 VQE approaches may access a comparable level of dynamic correlation by enlarging the active space mapped to qubits and/or employing highly expressive ansatzes (e.g., UCCGSD and k-UpCCGSD^34,35), leading to increased qubit requirements and circuit depth.

Second, we embed VQE-PDFT within a multiscale quantum mechanics/molecular mechanics (QM/MM) workflow to study electron transfer in a cryptochrome, computing Marcus parameters and transfer rates from ensembles of protein conformations in a complex biological environment.

Third, we develop shallow, symmetry-preserving hardware-efficient ansatzes tailored to open- and closed-shell tryptophan active spaces of the cryptochrome electron-transfer center, further reducing qubit count and circuit depth to support an initial hardware demonstration.

In this study, we first validated VQE-PDFT on the CT7/04 Charge-Transfer benchmark dataset, confirming accuracy comparable to conventional MC-PDFT. We then presented its QM/MM application to electron transfer in the European robin cryptochrome protein (ErCRY4), where the transfer rates from noiseless simulations agreed well with ultrafast spectroscopy measurements.⁴¹ Finally, as a proof-of-principle hardware validation on a simplified active-space model, a single conformation experiment on a 13-qubit superconducting device illustrated feasibility on current NISQ hardware within well-defined resource limits.

Throughout this work, unless explicitly stated otherwise, most of the quantum circuit calculations were carried out on noiseless classical simulators in limited active spaces, with quantum hardware used only for a single proof-of-principle validation.

2 Results

2.1 The VQE integrated MC-PDFT

Our quantum-classical hybrid approach VQE-PDFT replaces the computationally expensive CASSCF optimization in MC-PDFT with a VQE solver for the CASCI active-space Hamiltonian. In all calculations reported here, this CASCI Hamiltonian is constructed using Hartree–Fock canonical molecular orbitals for the chosen active spaces, without further orbital optimization included. Consequently, the VQE state should be viewed as a CASCI-level wavefunction rather than a fully self-consistent CASSCF state.

The modified workflow proceeds as follows: VQE optimizes a parameterized quantum circuit ansatz (UCCSD for the CT7/04 benchmark, and ROUCCSD or the empirical HEA circuits for the ErCRY4 application) to approximate the CASCI ground state in the chosen active space, from which one-particle and two-particle reduced density matrices (1-RDM and 2-RDM) are extracted by measuring the corresponding matrix element expectation values on the quantum circuit as shown in Fig. 1. These reduced density matrices are then utilized to compute the total energy following the MC-PDFT formalism.

Fig. 1 Workflow of VQE-PDFT. Quantum computing is used as a solver for the CASCI active-space Hamiltonian and for the subsequent evaluation of 1-RDM and 2-RDM elements, where “CASCI” denotes the complete active space configuration interaction active-space Hamiltonian and the associated one- and two-electron integrals. The VQE loop and the quantum circuit U(θ) below it represent the numerical solution of this CASCI Hamiltonian, yielding the state |ψ(θ)〉 from which the 1- and 2-RDMs are sampled and passed to further calculate the PDFT energy. Here, “CASCI” refers to the level of theory of the active-space Hamiltonian, while the particular VQE ansatz (UCCSD/ROUCCSD/HEA) specifies how the CASCI wavefunction is parameterized on the quantum circuit.

The 1-RDM provides the kinetic energy T, nuclear–electron interaction V_ne, total electron density ρ, and classical Coulomb repulsion V_ee(ρ). The 2-RDM yields the on-top pair density Π, which combines with ρ to determine the on-top density functional energy E_ot(ρ, Π) that captures exchange–correlation effects. Finally, together with the nuclear–nuclear repulsion V_nn, the total energy follows the standard MC-PDFT expression:

E = T + V_ne + V_nn + V_ee(ρ) + E_ot(ρ, Π). (1)

This hybrid strategy is more compatible with NISQ constraints because the quantum computation is confined to a CASCI active-space solver that captures static correlation, while dynamic correlation is recovered classically through the on-top functional.¹⁵ It therefore avoids treating both correlation regimes variationally on the quantum circuit, which would generally require a highly expressive yet prohibitively deep ansatz (e.g., UCCGSD³⁴). Moreover, capturing a large part of the dynamic correlation beyond a compact active space would require an expanded orbital space, further increasing the qubit counts.

2.1.1 Benchmarks on dissociation energy. To validate our VQE-PDFT approach, we evaluated its performance on the CT7/04 dataset⁴² with a noiseless simulated quantum circuit. To ensure computational feasibility, we employed a reduced active space up to (10e, 10o), requiring up to 20 qubits (see Table 1). Detailed quantum resource information for each dimer and monomer is summarized in SI Table S1.

Table 1 Dissociation energy of charge-transfer dimer dataset CT7/04. Energies are given in kcal mol⁻¹. The geometries of dimers and their reference values were fetched from ref. 42 and 45. Here, the column labeled ‘VQE’ corresponds to the CASCI energy 〈ψ(θ)|Ĥ_CASCI|ψ(θ)〉 evaluated from the same VQE-optimized active-space wavefunction that is used to construct the RDMs in VQE-PDFT. All ‘VQE’ and ‘VQE-PDFT’ results were obtained from noiseless classical simulations of UCCSD ansatz in limited active spaces (see SI Table S1). The column “N_q (dimer)” reports the number of qubits used for each dimer calculation (20 qubits for six systems and 16 qubits for one system). The detailed active space specifications and quantum resource summary (for both dimers and monomers) are provided in SI Table S1. The results of MC-PDFT were obtained from Ghosh et al.'s research.¹⁶ The jul-cc-pVTZ⁴⁶ basis set was used in all calculations. The tPBE functional was adopted in VQE-PDFT, which aligned with ref. 16. MUE stands for the mean unsigned error. More explicit absolute energies of all systems are listed in SI Table S3

Dimers	N q (dimer)	VQE	VQE-PDFT	MC-PDFT	W1-reference
NH3⋯FCl	20	5.38	11.56	12.42	10.62
NH3⋯Cl2	20	1.57	5.09	4.55	4.88
NH3⋯F2	20	−0.08	1.88	1.11	1.81
HCN⋯FCl	16	2.19	3.13	3.67	4.86
H2O⋯FCl	20	10.01	3.30	4.38	5.36
C2H2⋯FCl	20	4.13	3.07	4.01	3.81
C2H4⋯F2	20	−1.14	0.85	0.33	1.06
MUE		2.896	0.853	0.847	—

---

This dataset comprises 7 dimers that exhibit significant charge transfer effects upon dissociation. These systems are notoriously difficult to calculate accurately due to two key factors: their inherent multiconfigurational character and the dramatic electronic reorganization that occurs during dissociation.^43,44 Such characteristics make these dimers ideal test cases for evaluating multiconfigurational methods in further applications.¹⁶

Table 1 summarizes the CT7/04 dissociation energies obtained with CASCI-level VQE (‘VQE’), VQE-PDFT, and classical MC-PDFT, together with high-level Weizmann-1 theory reference values (‘W1-reference’).⁴² Utilizing the Unitary Coupled-Cluster Singles and Doubles (UCCSD) ansatz, our VQE-PDFT method achieved a mean unsigned error (MUE) of 0.853 kcal mol⁻¹ relative to W1 theory, closely matching the classical MC-PDFT result of 0.847 kcal mol⁻¹.¹⁶

We emphasize that this close agreement between VQE-PDFT and classical MC-PDFT for CT7/04 was achieved within the specific limited active spaces and numerical settings adopted here (SI Table S1), and different active-space choices or numerical approximations could modify the quantitative agreement even though both protocols share the same on-top functional.

Notably, for two dimers (NH₃⋯F₂ and C₂H₄⋯F₂), the CASCI-level VQE (column ‘VQE’ in Table 1) yielded negative dissociation energies. This indicated that, within these active spaces, the bare CASCI-level energies could be insufficiently correlated for these charge-transfer interactions. When the identical VQE-derived 1- and 2-RDMs were instead combined with the MC-PDFT functional (VQE-PDFT column), the dissociation energies restored the correct sign and moved closer to the reference values, and the MUE was reduced to 0.853 kcal mol⁻¹. This comparison implied that the improvement of VQE-PDFT over CASCI-level VQE arose primarily from the functional treatment of dynamic correlation, instead of a change in the underlying VQE state or additional quantum resources.

This interpretation is consistent with the analysis of Ghosh et al. for CT7/04, where classical CASSCF (an orbital-optimizing version of CASCI) was found to have the largest MUE (3.92 kcal mol⁻¹) among the methods considered and yielded negative dissociation energies on C₂H₂⋯FCl and C₂H₄⋯F₂ dimers in CT7/04. Ghosh et al. attributed this to the absence of an explicit approximation to the full dynamic correlation energy.¹⁶ Accordingly, the CASCI-level energies in the ‘VQE’ column may mainly reflect the known limitations of the CASCI framework with respect to dynamic correlation, rather than an intrinsic deficiency of the VQE algorithm itself. Conversely, the improved VQE-PDFT dissociation energies result from supplementing this missing dynamic correlation through the MC-PDFT functional applied to the same VQE-derived RDMs (even though the resulting energy is no longer a strict variational 〈Ψ|Ĥ|Ψ〉 for that state).

From this perspective, the ‘VQE’ column in Table 1 may be viewed as the CASCI-level analogue of the classical CASSCF reference step in conventional MC-PDFT, whereas the ‘VQE-PDFT’ and ‘MC-PDFT’ columns correspond to applying the same on-top pair-density functional to quantum- and classically obtained multiconfigurational references, respectively.

To better understand the error patterns, we further analyzed the deviations besides the MUE metric. The root-mean-square errors (RMSE) relative to the W1 reference were 1.129 kcal mol⁻¹ for VQE-PDFT and 0.985 kcal mol⁻¹ for MC-PDFT, showing that VQE-PDFT exhibited a slightly broader error distribution yet averaged out in the MUE metric. Additionally, inspection of the absolute total energies (SI Table S3) revealed a method-dependent systematic shift of VQE-PDFT compared to MC-PDFT results. However, since the dissociation energy is defined as an energy difference (E_dimer − ∑E_monomer), these systematic shifts partially canceled out, leaving residual errors that were comparable to the classical MC-PDFT reference. This suggested that the agreement in dissociation energies reflects a combination of similar dynamic correlation treatments via the tPBE functional and the beneficial cancellation of systematic offset.

In summary, this benchmark suggests that VQE-PDFT may effectively describe the electronic correlation in these charge-transfer systems. This indicates its potential for handling similar complex interactions found in bio-molecules, where multiconfigurational character is ubiquitous. Moreover, the observation that method-dependent shifts partially canceled in energy differences provides the rationale for our subsequent application to the electron transfer process in ErCRY4 protein. Since the key parameters (reorganization energy and driving force in Marcus theory) are also defined as differences between single-point energies, we anticipate a similar degree of error cancellation, enabling reliable predictions of electron transfer rates even within simplified active-space models.

2.2 Quantum-classical hybrid framework for biological electron transfer

Building upon the validated accuracy of our VQE-PDFT method, we developed a comprehensive framework to apply quantum–classical hybrid calculations to a representative biological system, focusing on the electron transfer process in ErCRY4. In Marcus Theory, the reorganization energy and driving force are constructed from differences of closely related single-point energies; therefore, the previously observed partial cancellation of method-dependent shifts in energy differences now provides a direct motivation for the following transfer rate calculations.

Overall, this framework addresses two key challenges: adapting a multiscale QM/MM architecture, based on standard electrostatic embedding and link-atom techniques, for complex biological environments; and designing hardware-efficient quantum circuits optimized for NISQ constraints in electron transfer calculations.

Regarding the multiscale coupling, the QM/MM machinery itself follows well-established practice. In our workflow, the quantum circuit is used specifically as a CASCI solver for the multiconfigurational QM region, replacing the classical CASSCF step within an otherwise MC-PDFT/QM/MM treatment. Because the present ErCRY4 proof-of-concept employs deliberately compact active spaces that are also classically tractable, we do not claim a demonstrated quantum speed-up for this application. Rather, the motivation is that, as qubit numbers and gate fidelities improve, a quantum CASCI solver could help alleviate the main bottleneck of classical multiconfigurational solvers, and thereby enable larger active spaces in the QM region within the same multiscale framework. This long-term perspective is consistent with quantum resource-estimate studies of other strongly correlated bioinorganic centers, such as the nitrogenase FeMo cofactor, where Fe–S clusters pose challenges for conventional multiconfigurational treatments and have been proposed as natural targets for future quantum simulations.^3,48

2.2.1 QM/MM multiscale architecture. Building on the hybrid quantum computing pipeline previously developed by Li et al.,⁴⁹ the multiscale architecture employs a QM/MM partitioning scheme where the quantum region, containing the multiconfigurational active sites, is treated by our VQE-PDFT method, while the surrounding protein environment and solvent are described using classical molecular mechanics force fields, as illustrated in Fig. 2.

Fig. 2 Multiscale calculation framework of ErCRY4 protein. The quantum mechanical (QM) region is set to be two indole rings on the adjacent tryptophan residues (TrpB and TrpC), which are treated in separate QM calculations, while the rest of the systems are considered at the molecular mechanical (MM) level. The standard electrostatic embedding and link-atom methodology are adopted for evaluating interactions between the QM and MM regions. The figures were drawn with the help of ChimeraX.⁴⁷

Extending this architecture for the subsequent study, the standard electrostatic embedding and link-atom methodology are implemented to handle the interface between QM and MM regions. In particular, it involves capping the covalent ‘C–C’ bond cut at the QM/MM boundary, ensuring proper treatment of boundary effects.⁵⁰

To provide a proof-of-concept application of this multiscale framework, we selected the electron transfer process between adjacent tryptophan residues (TrpB and TrpC in Fig. 2) in the European robin cryptochrome protein (ErCRY4) as our target study. This system presents an ideal test case due to its well-characterized experimental properties^5,41 and the multiconfigurational nature. The electron transfer process involves photo-excitation of the FAD cofactor followed by sequential electron transfer among adjacent tryptophan residues, during which tryptophan residues alternate between neutral and cationic states.

In this QM/MM setup, the QM region is represented by two separated fragments corresponding to the indole rings of TrpB and TrpC. Each single-point QM calculation is carried out on one fragment (TrpB or TrpC), while the remainder of the amino-acid side chain is treated at the MM level. For each QM fragment, the VQE-PDFT calculations use a compact CASCI active space consisting of three frontier orbitals localized on the tryptophan indole ring (SI Fig. S1), with 4 or 3 active electrons for the neutral and cationic states, respectively. The corresponding CASCI Hamiltonian and active orbitals are generated by the CASCI module of PySCF⁵¹ as interfaced in TenCirChem,⁵² starting from Hartree–Fock canonical orbitals and using the default active-space construction based on the specified (N_e, N_o) pattern.

However, the cationic states exhibit open-shell electronic character with unpaired electrons, making the previously used standard UCCSD ansatz inadequate, as it assumes paired α and β electrons in each spatial orbital. This limitation, combined with NISQ device constraints, necessitates the development of specialized quantum circuits.

2.3 Shallow-depth empirical ansatz design

Conventional UCCSD ansatz performs well on evaluating closed-shell systems, and can be extended to restricted open-shell UCCSD (ROUCCSD) by explicitly considering single-electron occupations and the corresponding excitations to accommodate the unpaired electron scenarios in cationic tryptophan. However, ROUCCSD's prohibitively deep circuit structure and extensive non-local gate operations render it impractical for NISQ device implementation. To address this limitation, we developed shallow depth empirical hardware-efficient ansatz (HEA) circuits specifically tailored for the electronic structure and symmetry of the tryptophan system.

Our approach employs distinct quantum circuits for the two electronic states encountered in the transfer process: a closed-shell HEA (CHEA, Fig. 3A) for neutral tryptophan and an open-shell HEA (OHEA, Fig. 3B) for cationic tryptophan. Both circuits are empirically designed to maintain the particle number conservation for each spin (i.e., within the target N_α/N_β symmetry sector after tapering), while achieving minimal circuit depth compared to the ROUCCSD ansatz. The design methodology is detailed in Section 4.2.

Fig. 3 The shallow-depth empirical ansatzes for electron transfer. (A) Closed-shell HEA circuit (CHEA) with a circuit depth of 4. (B) Open-shell HEA circuit (OHEA) with a circuit depth of 6. These two circuits are dramatically shallower than ROUCCSD with a circuit depth of 35. The circuit diagrams were drawn with the PennyLane library.⁵³

To validate the accuracy of the empirical HEA circuits, we randomly extracted a protein conformation from ErCRY4 molecular dynamics simulations reported in ref. 5. We then applied our multiscale framework by defining two adjacent tryptophan residues as separated QM regions for computational analysis. We employed a systematic approach that generated eight single-point energy calculations for two tryptophan residues, encompassing four open-shell cationic states and four closed-shell neutral states. These states are typically encountered in electron transfer processes, and the specific procedure of constructing these states follows the four-point scheme⁵⁴ (see also Section 2.4.2).

For this validation study, we performed calculations using three different methods: the empirical HEA, ROUCCSD, and full configuration interaction (FCI), each applied to the active space within the CASCI framework. The comparative results are presented in Table 2.

Table 2 Energy calculations on the active space of tryptophan. The 1st, 3rd, 5th, and 7th single-point calculations are open-shell with a (3e, 3o) active space, while the others are closed-shell with a (4e, 3o) active space. The 6-31G basis set was used in all calculations. All energies are given in Hartree

	E FCI	E ROUCCSD	E HEA	ΔEHEA-FCI
1	−2.1968	−2.1968	−2.1947	0.0021
2	−2.5535	−2.5535	−2.5500	0.0035
3	−2.1902	−2.1902	−2.1899	0.0004
4	−2.5424	−2.5424	−2.5391	0.0033
5	−2.1912	−2.1912	−2.1899	0.0013
6	−2.5621	−2.5621	−2.5584	0.0036
7	−2.1963	−2.1963	−2.1958	0.0005
8	−2.5310	−2.5310	−2.5279	0.0032

---

The results demonstrate that the proposed empirical HEA circuits achieve accuracy of energy differences typically ranging from 10⁻⁴ to 10⁻³ Hartree with respect to the FCI results. Moreover, ROUCCSD energies match FCI values within excellent precision (at least <10⁻⁴ Hartree). We note that this near identity between ROUCCSD and FCI may arise from the compact (3e, 3o) and (4e, 3o) active spaces adopted in this validation, and should not be interpreted as a general property of ROUCCSD for larger multiconfigurational active spaces, where it remains an approximate high-level reference to FCI.

Also in Table 2, one may observe the slightly smaller HEA-FCI deviations for the open-shell cases (rows 1, 3, 5, and 7), compared to closed-shell results. This is consistent with our deliberate circuit design: to accommodate the more complex open-shell configurations, OHEA is constructed with a modestly larger depth and more variational parameters than CHEA (depth 6 vs. 4, and 5 vs. 2 parameters; Fig. 3). At the same time, the empirical HEA circuits are intentionally constrained to a fixed shallow depth and a small number of variational parameters to prioritize NISQ feasibility. Therefore, the residual HEA-FCI discrepancies in Table 2 mainly reflect this resource-oriented ansatz choice rather than an intrinsic limitation of the general HEA framework.

Overall, this validation confirms that our HEA circuits maintain sufficient expressibility for the active space of the current electron transfer system while offering notable quantum resource savings.

2.4 Evaluation of electron transfer in the ErCRY4

2.4.1 Marcus theory. Having demonstrated the accuracy of the empirical HEA circuits, we applied them to compute electron transfer rates in the ErCRY4 protein system. Electron transfer kinetics in biological systems may be depicted by Marcus theory, with the transfer rate expressed as:which incorporates three key parameters: 〈|H_DA|〉, λ, and ΔG^o. The electronic coupling 〈|H_DA|〉 represents the interaction between two tryptophan residues, reflecting the orbital overlap and electronic communication pathways. Then the reorganization energy λ quantifies the energetic cost of nuclear rearrangement accompanying electron transfer. Lastly, the driving force ΔG^o corresponds to the thermodynamic bias for the transfer process, determined by the energy difference between reactant and product states.

Our VQE-PDFT framework within the multiscale QM/MM architecture enables direct evaluation of the reorganization energy and driving force.

2.4.2 Four-point scheme for λ and ΔG^o. With the Marcus theory framework outlined above, we now detail our computational approach. We focused specifically on the electron transfer between tryptophan residues TrpB (W372) and TrpC (W318) in the ErCRY4 protein. The initial state corresponded to cationic TrpB and neutral TrpC, while the final state involved neutral TrpB and cationic TrpC following electron transfer.

To quantify both the reorganization energy λ and driving force ΔG^o, we employed the four-point scheme,⁵⁴ a well-established thermodynamic approach that captures the essential physics of Marcus theory. This method is based on the fundamental assumption that the system responds instantaneously to changes in electronic charge distribution.⁵⁵ This approximation, although it neglects potential non-Markovian effects arising from slower environmental responses,⁵⁶ enables us to bypass costly dynamical simulations and focus directly on energetic differences between key electronic configurations.

As illustrated in Fig. 4, the four-point scheme evaluates single-point energies at four distinct state-geometry combinations: each electronic state (initial and final) calculated at both its own optimized nuclear conformation and at the geometry optimized for the other state. This approach captures the energetic penalty arising from the mismatch between optimal nuclear arrangements for different charge distributions. Thereafter, the reorganization energy and the driving force of a given protein conformation can be extracted through:

ΔG^o = |Eⁱ_i − E^f_f|, (3)

λ = |Eⁱ_f − Eⁱ_i| + |E^f_i − E^f_f|, (4)

where E^geometry_state denotes the energy of a given electronic state calculated at a specific nuclear geometry, while “i” and “f” stand for “initial” and “final”, respectively.

Fig. 4 The four-point scheme for evaluating reorganization energy and driving force. (a) TrpB⁺/TrpC⁰ state on the initial geometry; (b) TrpB⁰/TrpC⁺ state on the initial geometry; (c) TrpB⁰/TrpC⁺ state on the final geometry; (d) TrpB⁺/TrpC⁰ state on the final geometry.

To capture the dynamic behavior of the protein in biological environments, we sampled 20 distinct protein conformations from molecular dynamics simulations of ErCRY4 and applied the four-point scheme to each conformation using our VQE-PDFT/MM framework. We calculated λ and ΔG^o using both ROUCCSD and the proposed empirical HEA circuits, obtaining averaged value of λ = 0.5701 eV and ΔG^o = 0.0689 eV for ROUCCSD, and λ = 0.4356 eV and ΔG^o = 0.0724 eV for HEA, as listed in Table 3 and detailed in the Supplementary Information. The good agreement between HEA and ROUCCSD results, considering the great reduction in circuit depth in HEA, validates the feasibility of the empirical HEA for the current biological application, achieving a balance between accuracy and quantum resource efficiency, although with a compact active space.

Table 3 The averaged λ, ΔG^o, and 〈|H_DA|〉. The reorganization energy and driving force were evaluated by the four-point scheme, utilizing ROUCCSD as well as the empirical HEA, while the electronic coupling was determined by the direct coupling scheme. The detailed results for individual configurations are provided in the SI. λ, ΔG^o, and 〈|H_DA|〉 are given in eV, while 〈|H_DA|²〉 is in (eV)².The 6-31G basis set and the tPBE functional were used in all calculations

Ansatz	λ	ΔG^o	〈|HDA|〉	〈|HDA|^2〉
ROUCCSD	0.5701	0.0689	6.352 × 10^−3	1.1431 × 10^−4
HEA	0.4356	0.0724

---

2.4.3 Electronic coupling. With the reorganization energy and driving force determined, we now evaluate the electronic coupling 〈|H_DA|〉.

We computed the 〈|H_DA|〉 through transfer integral calculations between the TrpB and TrpC indole rings. The approach employed Boys localization⁵⁷ to transform canonical molecular orbitals into spatially localized orbitals centered on individual tryptophan sites, enabling clear identification of the relevant frontier orbitals for each fragment. This localization procedure was essential for defining well-separated donor and acceptor orbitals in weakly coupled systems where standard canonical orbitals are typically delocalized across multiple sites.

The electronic coupling is extracted using the direct coupling (DC) scheme,⁵⁸ which accounts for orbital non-orthogonality through:

H_DA = [T_DA − (e_D + e_A)S_DA/2]/(1 − S_DA²), (5)

where T_DA is the Hamiltonian coupling element, S_DA is the overlap integral, and e_D and e_A are localized orbital energies for donor and acceptor sites, respectively.

Applying this methodology to the same 20 protein configurations, we obtained an averaged electronic coupling of 〈|H_DA|〉 = 6.35 × 10⁻³ eV (6.35 meV) that aligned well with the 5 meV reported by Timmer et al.,⁴¹ also listed in Table 3.

2.4.4 Electron transfer rate. With all Marcus theory parameters determined from the 20 sampled protein configurations, we proceeded to evaluate the electron transfer rates using eqn (2).

The computed rates were derived by substituting the averaged reorganization energies, driving forces, and electronic couplings into the Marcus expression (eqn (2)). As a result, our empirical HEA approach yielded k_ET = 0.944 × 10¹⁰ s⁻¹, aligning well with both ROUCCSD calculations (0.864 × 10¹⁰ s⁻¹) and the experimental ultrafast transient absorption measurements⁴¹ (0.709 × 10¹⁰ s⁻¹), as shown in Fig. 5. Notably, these results were approximately one order of magnitude slower than previous DFT predictions⁵ ((5.0 ± 1.8) × 10¹⁰ s⁻¹), as also pointed out in ref. 41.

Fig. 5 Electron transfer rate comparison. The red marked transfer rate was evaluated on quantum hardware. The single conformation results show rates one order of magnitude lower than the averaged values of 20 conformations, reflecting typical conformational fluctuations in MD simulations where tryptophan positioning variations can alter transfer rates by several orders of magnitude.⁵⁶ The DFT prediction was obtained from ref. 5. Explicit parameters calculated for each conformation are listed in the SI.

2.4.4.1 Experimental validation on quantum hardware. Moreover, to validate our approach on actual NISQ hardware, we executed the empirical HEA circuits on a superconducting quantum device for a randomly selected protein conformation. With VQE's optimization performed on a noiseless simulator and density matrix evaluations on quantum hardware with readout error mitigation,^37,59 we obtained the single-point energies for the four-point scheme, eventually leading to the transfer rate k_ET ≈ 2.62 × 10⁸ s⁻¹, as listed in Table 4.

Table 4 Single-point energies, Marcus parameters, and the transfer rate of the single conformation. The conformation was randomly selected as the 2071 frame from the MD simulation and calculated in the QM/MM framework, from which the QM energies were extracted and listed here. The calculations were conducted by empirical HEA on a 13-qubit superconducting quantum hardware. The single-point energies are given in Hartree. The 6-31G basis set and the tPBE functional were used in all calculations. Explicit single-point calculations and Marcus parameters of ROUCCSD and noiseless HEA can be found in the SI

Index	System	Empirical HEA on quantum hardware
1	point-a TrpB^+	−362.92966
2	point-a TrpC^0	−363.17866
3	point-b TrpB^0	−363.20554
4	point-b TrpC^+	−362.91475
5	point-c TrpB^0	−363.19954
6	point-c TrpC^+	−362.91126
7	point-d TrpB^+	−362.92451
8	point-d TrpC^0	−363.18130

---

λ	0.47021 eV
ΔG^o	0.06764 eV
k ET	2.61742 × 10^8 s^−1

---

This result compared favorably with noiseless calculations: 1.89 × 10⁸ s⁻¹ (HEA) and 2.49 × 10⁸ s⁻¹ (ROUCCSD), demonstrating that quantum hardware may conduct electron transfer evaluations within a compact active space despite inherent noise. More explicit single-point energies and derived parameters are detailed in the Supplementary Information.

2.5 Error analysis

The successful quantum hardware validation demonstrated current feasibility, but might raise questions about the reliability of noisy quantum devices for the electron transfer process. To address these concerns, we conducted comprehensive error analysis comparing single-point energies and derived Marcus parameters across different computational approaches.

2.5.1 Systematic error in single-point energies. To assess the computational accuracy of quantum hardware implementation, we examined single-point energy errors of the protein conformation (frame 2071) presented in Table 4. Our analysis employed ROUCCSD calculations as the reference standard, building upon our earlier validation that demonstrated excellent agreement between ROUCCSD and FCI results (Table 2). The quantum hardware calculations revealed systematic deviations from this ROUCCSD reference, as illustrated in Fig. 6A and B, where all eight energy calculations exhibit consistent positive shifts.

Fig. 6 Error analysis of single-point energies and Marcus parameters. ROUCCSD results are set as the comparison baseline. (A) Error comparison on 8 single-point energies of quantum hardware and noiseless empirical HEA, with quantum hardware results showing systematic errors (65.3 ± 14.2 mHartree positive deviation) compared to noiseless HEA relative to the ROUCCSD reference. The system indices follow the index column in Table 4. (B) Error distributions of quantum hardware and noiseless empirical HEA. (C) Error magnitude comparison between single-point energies and Marcus parameters. Explicit single-point calculations can be found in the Supplementary Information.

This systematic bias indicated a non-random error source. The errors likely originated from correlated uncertainties in noisy matrix element measurements or from VQE-PDFT's systematic response to noisy RDMs, rather than from purely random quantum fluctuations.

Moreover, the systematic nature of these hardware errors is noteworthy for biological electron transfer applications. Unlike random noise that would propagate unpredictably and through energy difference calculations, systematic errors might enable partial cancellation when computing the energy differences that underlie Marcus theory parameters. This observation motivated examining how these correlated single-point energy errors influence the derived reorganization energy and driving force that determine electron transfer rates.

2.5.2 Error cancellation in Marcus parameters. The feasibility of Marcus parameter calculations on quantum hardware can be understood through their mathematical structure. Both the driving force and reorganization energy depend on energy differences, as defined in eqn (3) and (4). Consider the driving force calculation on quantum hardware:where systematic errors δE relative to the ROUCCSD reference partially cancel in the energy difference calculation, reducing the overall uncertainty to the net difference (δE_i − δE_f). A similar partial cancellation also functions in the reorganization energy calculation.

This error cancellation effect was demonstrated in our calculations, as illustrated in Fig. 6B and C. While individual single-point energies exhibited errors of 61–86 mHartree on quantum hardware, the derived Marcus parameters showed notably reduced absolute uncertainties: λ differed from the ROUCCSD reference by 2.0 mHartree, and ΔG^o by 0.9 mHartree. However, it should be noted that the relative errors of Marcus parameters actually increased due to their smaller absolute magnitudes (∼12% for λ, ∼36% for ΔG^o).

Despite these relative error increases, the systematic nature of quantum hardware errors enables partial cancellation when computing energy differences. This error cancellation allowed the derived electron transfer rates (2.62 × 10⁸ s⁻¹) to achieve consistency with classical predictions (1.89 ∼ 2.49 × 10⁸ s⁻¹), demonstrating that the approach remains viable for the current electron transfer application. Additionally, together with the similar effect shown in the CT7/04 benchmark test, the observed error cancellation in energy differences may thus provide a strategic pathway for NISQ devices to contribute meaningfully to computational biochemistry, by focusing on difference-based observables rather than absolute electronic energies.

3 Conclusion

In this work, we developed and validated VQE-PDFT, a quantum-classical hybrid framework that integrates the variational quantum eigensolver with multiconfiguration pair-density functional theory for strongly correlated electronic systems. In VQE-PDFT, the quantum circuit is used as a CASCI active-space solver to optimize the multiconfigurational wavefunction and to sample RDMs, while the total correlation energy is evaluated classically via the MC-PDFT on-top functional as a post-processing step. This strategy delivers accuracy comparable to conventional MC-PDFT on our benchmarks, while reducing quantum resource requirements (qubit count and circuit depth) relative to highly expressive VQE approaches that seek to capture full dynamic correlation variationally by enlarging the orbital space.

Validation on the CT7/04 dataset demonstrated excellent agreement with reference calculations (MUE = 0.853 kcal mol⁻¹ with respect to W-1 theory), confirming the method's reliability for multiconfigurational systems. We also extended this framework to a biological electron transfer application by developing specialized hardware-efficient ansatz circuits optimized for NISQ device constraints, achieving circuit depths of 4–6, compared to 35 for ROUCCSD, while maintaining high accuracy.

Our comprehensive application to electron transfer in the European robin cryptochrome ErCRY4 protein yielded transfer rates (0.944 × 10¹⁰ s⁻¹ for HEA and 0.864 × 10¹⁰ s⁻¹ for ROUCCSD) that aligned well with experimental ultrafast spectroscopy measurements (0.709 × 10¹⁰ s⁻¹), validating our quantum-classical approach for this complex biological environment. Finally, as a proof-of-concept hardware validation within this quantum-classical multiscale framework, we executed the reduced-density-matrix measurements for a single ErCRY4 protein conformation on a 13-qubit superconducting device and obtained an estimated electron transfer rate (2.62 × 10⁸ s⁻¹) that aligned well with noiseless simulations, even in the presence of hardware noise.

Moreover, this work also mitigates three challenges for NISQ-era quantum biology applications: explicit treatment of multiconfigurational electronic correlations through purpose-designed shallow-depth quantum circuits, enhancing the ability for calculating strongly correlated system; systematic error cancellation through focus on energy differences rather than absolute energies, where we discovered that errors partially cancel in Marcus parameter calculations with quantum hardware, alleviating the quantum noise's impact on derived observables; and integration within a scalable multiscale framework adaptable to future hardware improvements, minimizing software and algorithmic modification while maintaining utility via the flexible VQE calculation.

Several limitations warrant future investigation. Firstly, current quantum hardware constraints may restrict active space sizes and computational accuracy. Secondly, efficient optimization and measurement grouping methods should be developed to lower the costs of quantum computation, enabling sampling across more protein conformations on superconducting quantum hardware, which was the main obstacle while conducting experimental validation. Thirdly, our separate treatment of indole rings on two tryptophan residues, while justified by weak electronic coupling, represents an approximation that could be improved by simultaneous evaluation of multiple rings to capture longer-range interactions and intermediate transfer states. Besides, in this proof-of-principle study we did not perform a systematic sensitivity analysis with respect to the precise QM/MM partitioning or to enlarging the tryptophan active spaces, as such an investigation would require a substantially larger set of high-level quantum calculations; exploring these alternatives will be an important direction for future work.

Additionally, the empirical HEA circuits adopted here are constructed for the present compact tryptophan active spaces under fixed (N_α, N_β) constraints and are not expected to be directly transferable. Therefore, developing more systematic and automated strategies for symmetry- and configuration-guided HEA constructions, and assessing their performance in larger active spaces and broader scenarios (e.g., bond-dissociation/potential-energy curves), would also be a valuable direction for future work.

Future developments will further focus on constraining electronic populations during VQE optimization through penalty functions or Lagrangian multipliers with information from Mulliken population analysis, enabling treatment of multiple simultaneous transfer pathways. As quantum hardware develops, our scalable framework can accommodate increased qubit counts and improved fidelities for more comprehensive biological system calculations.

In conclusion, VQE-PDFT illustrates a practical pathway for quantum computing applications in biochemistry, suggesting that despite current limitations, quantum-classical hybrid approaches may yield reasonably accurate predictions for simplified active-space models of particular biological processes and can serve as an initial step toward exploring quantum utility in molecular sciences.

4 Methods

4.1 Quantum computing reduced density matrices

In VQE-PDFT, VQE is used as a CASCI solver to optimize the CASCI active-space Hamiltonian, which is performed before the PDFT energy evaluation. Each VQE run is initialized in a product state consistent with the chosen ansatz (Hartree–Fock determinant for UCCSD/ROUCCSD and the all-zero computational state for the HEA circuits). Once the optimization has converged, the parameters of circuits are frozen and passed for subsequently evaluating 1-RDM and 2-RDM. The matrix elements of 1-RDM and 2-RDM may be written as the expectation value in terms of creation/annihilation operators under the second quantization, such as for 1-RDM,

γ_pq = 〈ψ|a^†_pa_q|ψ〉,

where the electronic state |ψ〉 can be efficiently reconstructed with the previously frozen circuit parameters. Therefore,

γ_pq = 〈ψ(θ)|a^†_pa_q|ψ(θ)〉,

where

|ψ(θ)〉 = U(θ)|ψ₀〉.

This enables us to compute each element of 1-RDM by executing the circuit with frozen parameters and sampling the expectation value of the corresponding operators after fermion-to-qubit mapping such as the parity transformation.

The same procedure is also valid for computing 2-RDM Γ_pq,rs.

In this element-wise measurement protocol, the additional post-optimization measurement overhead is governed by the number of distinct RDM elements to be sampled. For an active space comprising N_so spin orbitals, the 1-RDM γ_pq and 2-RDM Γ_pq,rs would formally contain O(N_so²) and O(N_so⁴) matrix elements, respectively (up to constant prefactors from Hermiticity and index symmetries). After the fermion-to-qubit transformation, each fermionic operator is mapped to a constant number of Pauli words to be measured, so the number of distinct measurement operators required for sampling the full 1- and 2-RDM scales as O(N_so²) for the 1-RDM and O(N_so²) for the 2-RDM, respectively.

Given this measurement overhead, we note that alternative strategies have been proposed to reduce the quantum measurement cost of RDM evaluations. For instance, natural orbital functional (NOF) methods reconstruct an approximate 2-RDM from the measured 1-RDM, reducing formal sampling requirements from O(N_so⁴) to O(N_so²).⁶⁰ However, because MC-PDFT depends explicitly on 2-RDM information through the on-top pair density Π, introducing an additional 2-RDM reconstruction approximation may bias the functional evaluation. Our approach therefore directly measures the 1-RDM and 2-RDM elements in the present work. Separately, for evaluating the CASCI energy one may further reduce measurement cost by truncating Pauli terms with small Hamiltonian coefficients. These term-selection techniques are complementary to VQE-PDFT and could be combined with our framework.

4.2 Designing the empirical HEA

The hardware-efficient ansatz (HEA) circuits are designed based on the electronic structure characteristics of tryptophan indole systems, optimizing quantum resources while maintaining computational accuracy for NISQ devices. For clarity, the empirical HEA construction follows a simple workflow:

• Identify the active space and fermion-to-qubit mapping (parity transformation here), and determine the fixed N_α/N_β symmetry sector;

• Reduce the qubit register by tapering qubits associated with these symmetries;

• Determine the set of computational-basis configurations that satisfy particle conservation;

• Construct a shallow symmetry-preserving circuit that couples these configurations with minimal depth and parameter count.

4.2.1 Qubit reduction through symmetry constraints. The (4e, 3o) closed-shell and (3e, 3o) open-shell active spaces require 6 qubits under parity transformation. Exploiting α and β electron number conservation reduces this to 4 qubits, as qubits 3 and 6 (representing total α-electron and total electron occupation parity) remain fixed throughout calculations. This symmetry-based tapering reduces quantum resource requirements while preserving all physically accessible states.

4.2.2 Determination of the symmetry-allowed configuration subspace. After tapering, particle-number conservation further restricts the computational basis accessible in the reduced qubit register. For example, in open-shell (3e, 3o) systems with 2α and 1β electrons, the first two qubits must not access |00〉 (violating α-electron conservation) while qubits 3 and 4 must not access |10〉 (violating β-electron conservation), restricting the total accessible space to 3 × 3 = 9 states. Consistent with this symmetry restriction, the converged ROUCCSD wavefunctions have support only within this nine-state sector for the present active spaces.

4.2.3 Circuit construction based on the accessible subspace. With the symmetry constraints enforced and the accessible configuration subspace identified, we construct the CHEA/OHEA circuits in a pragmatic, empirical manner tailored to the present compact active spaces for the proof-of-concept demonstration. In addition, we may even select the highest-weight configurations within the symmetry-allowed subspace as the dominant configurations to be captured, via inspecting the classical reference in the same active space (such as ROUCCSD), and further reduce the target subspace. Then, we iteratively assemble a shallow circuit from a restricted gate set {X, R_y, CNOT, CR_y}: at each step, candidate gate additions are evaluated classically on the reduced register to ensure that the resulting state maintains particle number conservation and to increase the circuit's ability to populate the target subspace.

Owing to the strong particle-number constraints and the small, factorized structure of the α- and β-spin qubit registers in the present (4e, 3o)/(3e, 3o) cases, this iterative construction remains tractable and yields depth-4 (CHEA, Fig. 3A) and depth-6 (OHEA, Fig. 3B) circuits that are sufficiently flexible for the subsequent validation and hardware execution, compared to ROUCCSD's 35-layer requirement. A quantitative comparison of the logical circuit depths, total numbers of single- and two-qubit gates for ROUCCSD and the HEA circuits is provided in Supplementary Table S2(a).

This empirical approach represents a compromise necessitated by current NISQ limitations rather than fundamental algorithmic constraints that can be directly transferred for arbitrary systems. The system-specific ansatz optimized for tryptophan electronic structure patterns requires extensive reanalysis for different chemical environments or larger active spaces. For other systems, the same design workflow may be followed, but the resulting circuit structure will generally change with (a) the active-space size and occupation pattern (hence the number of qubits after symmetry tapering) and (b) the number and connectivity of dominant configurations within the target symmetry sector, which may require additional entangling layers and variational parameters to maintain accuracy. As quantum hardware develops to support deeper circuits (hundreds of circuit-depths) with high fidelity, more general approaches like ROUCCSD or highly expressive multi-layer HEAs may become preferable over these tailored solutions.⁶¹

4.3 Computational details

4.3.1 Quantum hardware implementation. VQE optimizations were performed on classical simulators to avoid noise accumulation during iterative parameter updates. Once parameter convergence was achieved, the optimized parameters were saved and transferred to quantum hardware for reduced density matrix evaluations.

4.3.2 Hardware specification. Quantum computation employed a customized 13-qubit superconducting quantum device with the following performance metrics: single-qubit gate fidelity of 99.93%, two-qubit CZ gate fidelity of 99.13%, coherence times of T₁ = 83.8 µs and T₂ = 45.6 µs, and readout fidelities of F₀ = 98.35% and F₁ = 95.88%. The connecting topology is shown in Fig. 7.

Fig. 7 Topology of a 13-qubit superconducting quantum processor. During our validations, we only used 4 qubits indexing {0, 3, 11, 12} and their couplings.

4.3.3 RDM matrix element evaluation. The computation of 1-RDM and 2-RDM matrix elements followed a systematic six-step protocol:

(1) Pauli operator preparation. Pauli measurement operators were obtained through parity transformation of the corresponding matrix elements.

(2) Measurement grouping. Compatible measurement operators were grouped using a simple grouping strategy to reduce quantum measurement overhead and noise impact. For instance, operators Z₁Z₂I₃I₄ and I₁I₂Z₃Z₄ can be simultaneously measured using the single operation Z₁Z₂Z₃Z₄.

(3) Circuit compilation. The empirical HEA circuits and grouped measurements were uploaded to the quantum hardware platform. These circuits then underwent qubit mapping and gate decomposition to match the chip topology and native gate set. As shown in Fig. 7, the resulting compiled circuits utilized only 4 physical qubits, benefiting from our restrained use of two-qubit gates that eliminated the need for additional connectivity qubits. Quantitative information for the compiled HEA circuits on the 13-qubit device is summarized in SI Table S2(b).

(4) Circuit execution. Each compiled circuit was executed with 2048 measurement shots, returning raw bit-string distribution data.

(5) Readout error mitigation. Raw measurement data underwent readout error correction using the correlated Markovian noise model approach. This correction required calibration circuits executed once per single-point energy calculation with 8192 shots. The resulting calibration data were then stored for all subsequent measurements within that calculation. Finally, the error mitigation module processed raw bit-string distribution, calibration data, and target measurement operators to yield the corrected expectation value as the matrix elements.

(6) RDM assembly. Steps 1–5 were repeated for all required matrix elements to construct the entire 1-RDM and 2-RDM matrices.

This protocol enabled reliable extraction of reduced density matrices from noisy quantum hardware while maintaining computational efficiency through strategic measurement grouping and error mitigation.

4.3.4 Molecular dynamics sampling and statistical analysis. Twenty protein conformations were randomly sampled from molecular dynamics simulations of ErCRY4.⁵ All indices of these conformations are detailed in the SI. Marcus theory parameters (reorganization energy λ, driving force ΔG^o, and electronic coupling 〈|H_DA|〉) were calculated independently for each conformation, with statistical averages of these parameters input into the Marcus rate expression. This approach ensures thermodynamically consistent averaging of Marcus theory parameters, representing the actual biological environment where conformational fluctuations modulate the underlying electronic structure properties.

4.3.5 Software implementation. The VQE-PDFT framework was implemented in Python through local modifications of PySCF⁵¹ and TenCirChem⁵² libraries. QM/MM calculations employed the ASH⁶² framework for the link-atom and electrostatic embedding between the QM regions and the MM regions. The quantum hardware validation process was performed via the TensorCircuit⁶³ library for measurement grouping, quantum device communication, and readout error mitigation.

Author contributions

Conceptualization: Y. C. and W. L.; methodology: Y. C., Z. S., and W. L.; software: Y. C., Z. S., and W. L.; validation: Y. C., Z. S., and W. L.; formal analysis: Y. C. and W. L.; investigation: Y. C.; resources: J.-H. H. and Y. L.; data curation: Y. C.; writing – original draft: Y. C. and Z. S.; writing – review & editing: Y. C., Z. S., W. L., Y. Z., X. X., J.-H. H., and Y. L.; visualization: Y. C.; supervision: W. L., Y. Z., X. X., J.-H. H., and Y. L.; project administration: J.-H. H.; funding acquisition: W. L., J.-H. H., Y. Z., and Y. L.

Conflicts of interest

The authors declare no competing interests.

Data availability

The CT7/04 dimers reported in ref. 42 and 45 can be fetched from https://comp.chem.umn.edu/db/dbs/ncce31.html. The MD simulation file of ErCRY4 reported in ref. 5 can be obtained from https://cloud.uol.de/s/NrTYpoEzL6RbPq7. Other relevant data supporting this study are available from the corresponding author upon reasonable request. The custom code and local modifications for implementation are available in the GitHub repository: https://github.com/yiboch/VQE-PDFT.

Supplementary information (SI): visualizations of active space orbitals and detailed summaries of quantum resource requirements, including qubit counts and circuit depths for the computational schemes; tabulated data on absolute energies for the charge-transfer dataset, alongside calculated Gibbs free energies, reorganization energies, and electronic coupling values for the sampled protein conformations. See DOI: https://doi.org/10.1039/d5sc07528a.

Acknowledgements

This work was supported by the State Key Laboratory of Genome and Multi-omics Technologies, the Guangdong Bigdata Engineering Technology Research Center for Life Sciences and the Shenzhen City Peacock Team Project (grant no. KQTD20240729102028011 to W. L.).

References

1. T. Yanagisawa, Mechanism of High-Temperature Superconductivity in Correlated-Electron Systems, Condens. Matter, 2019, 4(2), 57.
2. C. Biz, M. Fianchini and J. Gracia, Strongly Correlated Electrons in Catalysis: Focus on Quantum Exchange, ACS Catal., 2021, 11(22), 14249–14261.
3. M. Reiher, N. Wiebe, K. M. Svore, D. Wecker and M. Troyer, Elucidating Reaction Mechanisms on Quantum Computers, Proc. Natl. Acad. Sci. U. S. A., 2017, 114(29), 7555–7560.
4. Y. Zhang, G. P. Berman and S. Kais, The Radical Pair Mechanism and the Avian Chemical Compass: Quantum Coherence and Entanglement, Int. J. Quantum Chem., 2015, 115(19), 1327–1341.
5. J. Xu, L. E. Jarocha, T. Zollitsch, M. Konowalczyk, K. B. Henbest and S. Richert, et al., Magnetic Sensitivity of Cryptochrome 4 from a Migratory Songbird, Nature, 2021, 594(7864), 535–540.
6. C. Zhou, M. R. Hermes, D. Wu, J. J. Bao, R. Pandharkar and D. S. King, et al., Electronic Structure of Strongly Correlated Systems: Recent Developments in Multiconfiguration Pair-Density Functional Theory and Multiconfiguration Nonclassical-Energy Functional Theory, Chem. Sci., 2022, 13(26), 7685–7706.
7. B. G. Levine, A. S. Durden, M. P. Esch, F. Liang and Y. Shu, CAS without SCF-Why to Use CASCI and Where to Get the Orbitals, J. Chem. Phys., 2021, 154(9), 090902.
8. J. Olsen, The CASSCF Method: A Perspective and Commentary, Int. J. Quantum Chem., 2011, 111(13), 3267–3272.
9. A. J. Wallace and D. L. Crittenden, Optimal Composition of Atomic Orbital Basis Sets for Recovering Static Correlation Energies, J. Phys. Chem. A, 2014, 118(11), 2138–2148.
10. S. Pathak, L. Lang and F. Neese, A Dynamic Correlation Dressed Complete Active Space Method: Theory, Implementation, and Preliminary Applications, J. Chem. Phys., 2017, 147(23), 234109.
11. C. L. Benavides-Riveros, N. N. Lathiotakis and M. A. L. Marques, Towards a Formal Definition of Static and Dynamic Electronic Correlations, Phys. Chem. Chem. Phys., 2017, 19(20), 12655–12664.
12. T. Helgaker, P. Jørgensen and J. Olsen, in Multiconfigurational Self-Consistent Field Theory, 1st edn, 2000, pp. 598–647.
13. B. O. Roos, K. Andersson, M. P. Fülscher, P. Malmqvist, L. Serrano-Andrés, K. Pierloot, et al., in Multiconfigurational Perturbation Theory: Applications in Electronic Spectroscopy, 1st edn, 1996, pp. 219–331.
14. S. Battaglia, I. Fdez Galván and R. Lindh, in Multiconfigurational Quantum Chemistry: The CASPT2 Method, 2023, pp, 135–162.
15. M. G. Li, R. K. Carlson, S. Luo, D. Ma, J. Olsen and D. G. Truhlar, et al., Multiconfiguration Pair-Density Functional Theory, J. Chem. Theory Comput., 2014, 10(9), 3669–3680.
16. S. Ghosh, A. L. Sonnenberger, C. E. Hoyer, D. G. Truhlar and L. Gagliardi, Multiconfiguration Pair-Density Functional Theory Outperforms Kohn-Sham Density Functional Theory and Multireference Perturbation Theory for Ground-State and Excited-State Charge Transfer, J. Chem. Theory Comput., 2015, 11(8), 3643–3649.
17. L. Gagliardi, D. G. Truhlar, G. Li Manni, R. K. Carlson, C. E. Hoyer and J. L. Bao, Multiconfiguration Pair-Density Functional Theory: A New Way To Treat Strongly Correlated Systems, Acc. Chem. Res., 2017, 50(1), 66–73.
18. D. S. Levine, D. Hait, N. M. Tubman, S. Lehtola, K. B. Whaley and M. Head-Gordon, CASSCF with Extremely Large Active Spaces Using the Adaptive Sampling Configuration Interaction Method, J. Chem. Theory Comput., 2020, 16(4), 2340–2354.
19. A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink and J. M. Chow, et al., Hardware-Efficient Variational Quantum Eigensolver for Small Molecules and Quantum Magnets, Nature, 2017, 549(7671), 242–246.
20. Google AI Quantum and Collaborators, F. Arute, K. Arya, R. Babbush, D. Bacon and J. C. Bardin, et al., Hartree-Fock on a Superconducting Qubit Quantum Computer, Science, 2020, 369(6507), 1084–1089.
21. R. N. Tazhigulov, S. N. Sun, R. Haghshenas, H. Zhai, A. T. K. Tan and N. C. Rubin, et al., Simulating Models of Challenging Correlated Molecules and Materials on the Sycamore Quantum Processor, PRX Quant., 2022, 3(4), 040318.
22. J. Preskill, Quantum Computing in the NISQ Era and Beyond, Quantum, 2018, 2, 79.
23. A. Peruzzo, J. McClean, P. Shadbolt, M. H. Yung, X. Q. Zhou and P. J. Love, et al., A Variational Eigenvalue Solver on a Photonic Quantum Processor, Nat. Commun., 2014, 5(1), 4213.
24. H. R. Grimsley, S. E. Economou, E. Barnes and N. J. Mayhall, An Adaptive Variational Algorithm for Exact Molecular Simulations on a Quantum Computer, Nat. Commun., 2019, 10(1), 3007.
25. H. L. Tang, V. O. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall and E. Barnes, et al., Qubit-ADAPT-VQE: An Adaptive Algorithm for Constructing Hardware-Efficient Ansatze on a Quantum Processor, PRX Quant., 2021, 2(2), 020310.
26. Y. Zhang, L. Cincio, C. F. A. Negre, P. Czarnik, P. J. Coles and P. M. Anisimov, et al., Variational Quantum Eigensolver with Reduced Circuit Complexity, npj Quantum Inf., 2022, 8(1), 1–10.
27. W. Li, Z. Huang, C. Cao, Y. Huang, Z. Shuai and X. Sun, et al., Toward Practical Quantum Embedding Simulation of Realistic Chemical Systems on Near-Term Quantum Computers, Chem. Sci., 2022, 13(31), 8953–8962.
28. G. L. R. Anselmetti, D. Wierichs, C. Gogolin and R. M. Parrish, Local, Expressive, Quantum-Number-Preserving VQE Ansätze for Fermionic Systems, New J. Phys., 2021, 23(11), 113010.
29. B. Senjean, S. Yalouz and M. Saubanère, Toward Density Functional Theory on Quantum Computers? SciPost, Physics, 2023, 14(3), 055.
30. T. S. Hardikar, K. Heitritter, J. Brown, R. D'Cunha, A. Mitra, S. Weatherly, et al., Quanta-Bind: A Quantum Computing Pipeline for Modeling Strongly Correlated Metal-Protein Interactions, in 2024 IEEE International Conference on Quantum Computing and Engineering (QCE), 2024, pp. 538–544.
31. M. Otten, M. R. Hermes, R. Pandharkar, Y. Alexeev, S. K. Gray and L. Gagliardi, Localized Quantum Chemistry on Quantum Computers, J. Chem. Theory Comput., 2022, 18(12), 7205–7217.
32. A. Mitra, R. D'Cunha, Q. Wang, M. R. Hermes, Y. Alexeev and S. K. Gray, et al., The Localized Active Space Method with Unitary Selective Coupled Cluster, J. Chem. Theory Comput., 2024, 4c00528.
33. K. Sugisaki, T. Kato, Y. Minato, K. Okuwaki and Y. Mochizuki, Variational Quantum Eigensolver Simulations with the Multireference Unitary Coupled Cluster Ansatz: A Case Study of the C_2v Quasi-Reaction Pathway of Beryllium Insertion into a H₂ Molecule, Phys. Chem. Chem. Phys., 2022, 24(14), 8439–8452.
34. D. Wecker, M. B. Hastings and M. Troyer, Progress towards Practical Quantum Variational Algorithms, Phys. Rev. A, 2015, 92(4), 042303.
35. J. Lee, W. J. Huggins, M. Head-Gordon and K. B. Whaley, Generalized Unitary Coupled Cluster Wave Functions for Quantum Computation, J. Chem. Theory Comput., 2019, 15(1), 311–324.
36. S. Dasgupta and T. Humble, Impact of Unreliable Devices on Stability of Quantum Computations, ACM Trans. Quantum Comput., 2024, 5(4), 1–23.
37. P. Lolur, M. Skogh, W. Dobrautz, C. Warren, J. Biznárová and A. Osman, et al., Reference-State Error Mitigation: A Strategy for High Accuracy Quantum Computation of Chemistry, J. Chem. Theory Comput., 2023, 19(3), 783–789.
38. J. Tilly, P. V. Sriluckshmy, A. Patel, E. Fontana, I. Rungger and E. Grant, et al., Reduced Density Matrix Sampling: Self-Consistent Embedding and Multiscale Electronic Structure on Current Generation Quantum Computers, Phys. Rev. Res., 2021, 3(3), 033230.
39. M. Matoušek, K. Pernal, F. Pavošević and L. Veis, Variational Quantum Eigensolver Boosted by Adiabatic Connection, J. Phys. Chem. A, 2024, 128(3), 687–698.
40. J. N. Boyn, A. O. Lykhin, S. E. Smart, L. Gagliardi and D. A. Mazziotti, Quantum-Classical Hybrid Algorithm for the Simulation of All-Electron Correlation, J. Chem. Phys., 2021, 155(24), 244106.
41. D. Timmer, A. Frederiksen, D. C. Lünemann, A. R. Thomas, J. Xu and R. Bartölke, et al., Tracking the Electron Transfer Cascade in European Robin Cryptochrome 4 Mutants, J. Am. Chem. Soc., 2023, 145(21), 11566–11578.
42. Y. Zhao and D. G. Truhlar, Benchmark Databases for Nonbonded Interactions and Their Use To Test Density Functional Theory, J. Chem. Theory Comput., 2005, 1(3), 415–432.
43. J. G. Vitillo, C. J. Cramer and L. Gagliardi, Multireference Methods Are Realistic and Useful Tools for Modeling Catalysis, Isr. J. Chem., 2022, 62(1–2), e202100136.
44. H. Lischka, D. Nachtigallová, A. J. A. Aquino, P. G. Szalay, F. Plasser and F. B. C. Machado, et al., Multireference Approaches for Excited States of Molecules, Chem. Rev., 2018, 118(15), 7293–7361.
45. Y. Zhao and D. G. Truhlar, Design of Density Functionals That Are Broadly Accurate for Thermochemistry, Thermochemical Kinetics, and Nonbonded Interactions, J. Phys. Chem. A, 2005, 109(25), 5656–5667.
46. E. Papajak, J. Zheng, X. Xu, H. R. Leverentz and D. G. Truhlar, Perspectives on Basis Sets Beautiful: Seasonal Plantings of Diffuse Basis Functions, J. Chem. Theory Comput., 2011, 7(10), 3027–3034.
47. E. C. Meng, T. D. Goddard, E. F. Pettersen, G. S. Couch, Z. J. Pearson and J. H. Morris, et al., ChimeraX, Protein Sci., 2023, 32(11), e4792.
48. Z. Li, J. Li, N. S. Dattani, C. J. Umrigar and G. K. L. Chan, The Electronic Complexity of the Ground-State of the FeMo Cofactor of Nitrogenase as Relevant to Quantum Simulations, J. Chem. Phys., 2019, 150(2), 024302.
49. W. Li, Z. Yin, X. Li, D. Ma, S. Yi and Z. Zhang, et al., A Hybrid Quantum Computing Pipeline for Real World Drug Discovery, Sci. Rep., 2024, 14(1), 16942.
50. H. M. Senn and W. Thiel, QM/MM Methods for Biomolecular Systems, Angew. Chem., Int. Ed., 2009, 48(7), 1198–1229.
51. Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry and N. S. Blunt, et al., Recent Developments in the PySCF Program Package, J. Chem. Phys., 2020, 153(2), 024109.
52. W. Li, J. Allcock, L. Cheng, S. X. Zhang, Y. Q. Chen and J. P. Mailoa, et al., TenCirChem: An Efficient Quantum Computational Chemistry Package for the NISQ Era, J. Chem. Theory Comput., 2023, 19(13), 3966–3981.
53. V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, S. Ahmed and V. Ajith, et al., PennyLane: Automatic Differentiation of Hybrid Quantum-Classical Computations, arXiv, 2022, preprint, arXiv:1811.04968 [quant-ph],  DOI:10.48550/arXiv.1811.04968.
54. O. López-Estrada, H. G. Laguna, C. Barrueta-Flores and C. Amador-Bedolla, Reassessment of the Four-Point Approach to the Electron-Transfer Marcus-Hush Theory, ACS Omega, 2018, 3(2), 2130–2140.
55. R. A. Marcus and N. Sutin, Electron Transfers in Chemistry and Biology, Biochim. Biophys. Acta Rev. Bioenerg., 1985, 811(3), 265–322.
56. T. Firmino, E. Mangaud, F. Cailliez, A. Devolder, D. Mendive-Tapia and F. Gatti, et al., Quantum Effects in Ultrafast Electron Transfers within Cryptochromes, Phys. Chem. Chem. Phys., 2016, 18(31), 21442–21457.
57. J. E. Subotnik, S. Yeganeh, R. J. Cave and M. A. Ratner, Constructing Diabatic States from Adiabatic States: Extending Generalized Mulliken-Hush to Multiple Charge Centers with Boys Localization, J. Chem. Phys., 2008, 129(24), 244101.
58. C. P. Hsu, The Electronic Couplings in Electron Transfer and Excitation Energy Transfer, Acc. Chem. Res., 2009, 42(4), 509–518.
59. S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay and J. M. Gambetta, Mitigating Measurement Errors in Multiqubit Experiments, Phys. Rev. A, 2021, 103(4), 042605.
60. J. F. H. Lew-Yee and M. Piris, Efficient Energy Measurement of Chemical Systems via One-Particle Reduced Density Matrix: A NOF-VQE Approach for Optimized Sampling, J. Chem. Theory Comput., 2025, 21(5), 2402–2413.
61. S. Guo, J. Sun, H. Qian, M. Gong, Y. Zhang and F. Chen, et al., Experimental Quantum Computational Chemistry with Optimized Unitary Coupled Cluster Ansatz, Nat. Phys., 2024, 20(8), 1240–1246.
62. R. Björnsson, ASH – a Multiscale Modelling Program, https://ash.readthedocs.io/.
63. S. X. Zhang, J. Allcock, Z. Q. Wan, S. Liu, J. Sun and H. Yu, et al., TensorCircuit: A Quantum Software Framework for the NISQ Era, Quantum, 2023, 7, 912.

---