PaperlayerView exact revision on arXiv

Source-linked AI summary

Reliably assessing the electronic structure of cytochrome P450 on today's classical computers and tomorrow's quantum computers

Joshua J. Goings, Alec White, Joonho Lee, Christofer S. Tautermann, Matthias Degroote, Craig Gidney, Toru Shiozaki, Ryan Babbush, Nicholas C. Rubin

arXiv:2202.01244v1quant-phphysics.chem-ph

TL;DR

Reliable quantum-advantage assessment for chemical simulation requires classical benchmarks alongside quantum resource estimates for relevant chemical questions. The paper compares these resources for CYP models using classical electronic-structure methods and qubitized phase estimation compiled to the surface code. Its estimates suggest potential quantum advantage for high-accuracy CYP simulations, while direct full-system simulation remains unlikely to be feasible.

  • Problem

    Reliable quantum-advantage assessment for chemical simulation requires defining the classical frontier alongside resource estimates for relevant chemical questions.

    Source: §I · pp. 1–2

  • Method

    The study benchmarks CYP models with classical electronic-structure methods and estimates qubitized phase-estimation resources compiled to the surface code.

    Source: §IV · p. 11, §V · p. 14

  • Results

    The estimates indicate potential quantum advantage for high-accuracy CYP models, with the largest models predicted to run in under 100 hours under realistic hardware configurations.

    Source: §V · pp. 12, 14

  • Takeaways & Limitations

    Classical calculations help define chemistry's computational frontier and guide quantum-algorithm research.

    Source: §V · p. 14

    Direct simulation of the full Cpd I system is unlikely to be feasible given its very large resource requirements.

    Source: §V · p. 12

Abstract

from arXiv · show

An accurate assessment of how quantum computers can be used for chemical simulation, especially their potential computational advantages, provides important context on how to deploy these future devices. In order to perform this assessment reliably, quantum resource estimates must be coupled with classical simulations attempting to answer relevant chemical questions and to define the classical simulation frontier. Herein, we explore the quantum and classical resources required to assess the electronic structure of cytochrome P450 enzymes (CYPs) and thus define a classical-quantum advantage boundary. This is accomplished by analyzing the convergence of DMRG+NEVPT2 and coupled cluster singles doubles with non-iterative triples (CCSD(T)) calculations for spin-gaps in models of the CYP catalytic cycle that indicate multireference character. The quantum resources required to perform phase estimation using qubitized quantum walks are calculated for the same systems. Compilation into the surface-code provides runtime estimates to compare directly to DMRG runtimes and to evaluate potential quantum advantage. Both classical and quantum resource estimates suggest that simulation of CYP models at scales large enough to balance dynamic and multiconfigurational electron correlation has the potential to be a quantum advantage problem and emphasizes the important interplay between classical simulations and quantum algorithms development for chemical simulation.

I. INTRODUCTION

The paper frames quantum-advantage assessment as a comparison between rigorous classical simulations and realistic quantum resource estimates for biologically relevant CYP models. It uses those models to examine chemical questions including catalytic reactivity and spin-state ordering.

I. INTRODUCTION

  • Quantum-advantage boundaries are difficult to establish because classical methods may solve a problem, quantum resource estimates evolve, and either approach may prove prohibitive.Before scalable fault-tolerant quantum computers exist, estimates are rigorous upper bounds based on prefactors.Source: §I · p. 1
  • The study uses CYP mimics to assess chemical mechanism and reactive-intermediate spin ordering, then estimates phase-estimation runtimes and logical-qubit requirements under surface-code error correction.These questions are framed as representative chemical tasks for comparing classical and quantum algorithms.Source: §I · pp. 1–2
  • Classical benchmarking uses systematically improvable methods because active-space treatments must balance convergence with the cost of accounting for correlation outside the active space.Coupled-cluster methods can fail for significant multireference character, while larger active spaces and balanced external-correlation corrections remain costly challenges.Source: §I · p. 2
  • The quantum estimate uses qubitized phase estimation and compares three Hamiltonian factorizations, with tensor hypercontraction yielding substantially lower CYP runtimes.The estimates also examine tradeoffs between qubit count and Toffoli factories across hardware configurations.Source: §I · p. 2
  • CYPs are a biologically relevant target because they metabolize drugs, and their catalytic cycle includes at least eight intermediates whose spin-state ordering can affect chemical interpretation.The study models water-bound resting, pentacoordinate empty, pyridine-inhibited, and Cpd I states; Cpd I's near-degenerate doublet and quartet ordering remains debated.Source: §I · pp. 2–3, Figure 1 · p. 3
  • The study combines classical electronic-structure analysis, quantum resource estimates, and correlation characterization to assess four CYP models; it identifies Cpd I as showing multiconfigurational character.The authors report that the other three models can be characterized by single-reference methods after dynamic correlation is included.Source: §I · p. 3

II. CLASSICAL CALCULATIONS OF ELECTRON CORRELATION IN P450 MODELS

The classical calculations build a hierarchy of active-space models from CYP3A4 structures and use CCSD(T) and DMRG with NEVPT2 corrections to characterize spin states and support quantum-resource comparisons.

II. CLASSICAL CALCULATIONS OF ELECTRON CORRELATION IN P450 MODELS

  • The study uses CCSD(T) and DMRG with NEVPT2 corrections to characterize electronic structure across experimentally motivated spin states.The calculations' timings and accuracy provide context for quantum resource estimates and the proposed advantage boundary.Source: §II · p. 4
  • The CYP3A4 models derive from experimental X-ray structures after non-iron-coordinating protein entities and solvent are removed.Geometry-construction details are given in Appendix B.Source: §II · p. 4
  • A hierarchy of active spaces is used to examine electronic structure and motivate a potential quantum-advantage boundary.The paper points to Appendix C for computational details.Source: §II · p. 4

A. Active space selection

The active spaces form a nested hierarchy that starts from the high-spin reference's singly occupied orbitals and adds further orbitals to support balanced convergence and cost analysis.

A. Active space selection

  • The hierarchy begins with the five singly occupied orbitals of the high-spin reference and adds occupied and virtual orbitals in successive spaces.The orbitals are localized from the ROHF state using Pipek-Mezey localization.Source: A. Active space selection · p. 4, Table I · p. 4
  • The hierarchy prioritizes balanced active spaces for cost and convergence analysis rather than the smallest possible space at each size.For heme-iron systems, computed spin gaps depend strongly on active-space choice.Source: A. Active space selection · p. 4
  • Each space is a superset of the preceding one, and all spaces have roughly the same filling fraction.The sequence is A, B, C, D, E, F, G, and X.Source: A. Active space selection · p. 4

B. Spin-state ordering from DMRG+NEVPT2

DMRG finds the Cpd I doublet and quartet nearly degenerate, while the doublet shows multiconfigurational character and disagreement with CCSD(T). The near-degeneracy and sensitivity to dynamic-correlation corrections complicate identifying the lowest spin state.

B. Spin-state ordering from DMRG+NEVPT2

C. Cost estimates for classical computation

DMRG cost depends on active-orbital count and bond dimension, which controls calculation quality. The estimated 58-orbital calculation is a major computational investment, and larger spaces quickly become intractable for this model.

C. Cost estimates for classical computation

III. QUANTUM COMPUTING RESOURCE ESTIMATES

The study estimates phase-estimation resources for active-space Hamiltonians using qubitized quantum walks and surface-code error correction. For the largest Cpd I model, the estimates specify substantial physical-qubit and runtime requirements.

III. QUANTUM COMPUTING RESOURCE ESTIMATES

  • At 0.1% physical error, phase estimation for the X-space Cpd I Hamiltonian is estimated to require 4.6 million physical qubits and 73 hours.At 0.001% error, the estimate is 500,000 physical qubits and 25 hours.Source: §III · p. 6
  • Qubitized phase estimation estimates eigenvalues and prepares eigenstates through repeated applications of a quantum-walk operator, scaling as O(λ/ϵ).λ is the L1-norm of the Hamiltonian coefficients and ϵ is the target error.Source: §III · p. 6
  • The resource analysis compares single factorization, double factorization, and tensor hypercontraction for Hamiltonian coefficients.Double factorization and tensor hypercontraction use basis rotations to diagonalize central tensors; THC uses a projection into a larger basis.Source: §III · p. 7
  • Phase-estimation costs depend on factorization: Γ scales as Õ(N^3) for SF, Õ(N^2Ξ) for DF, and Õ(N^2) for THC if its rank grows linearly.Ξ is the average rank of the second factorization and is expected to scale as O(N) in most regimes.Source: §III · p. 7, Table III · p. 7

A. Quantum resource scaling

Resource estimates depend on the accuracy cutoff used to represent the two-electron integrals and on the selected Hamiltonian factorization. The study uses CCSD(T) convergence to set cutoffs and extrapolates resource scaling from the modeled active spaces.

A. Quantum resource scaling

  • The factorization cutoff is chosen so reconstructed-integral CCSD(T) energies for high-spin geometries are within one milliHartree.The authors note that additive errors could require a lower cutoff.Source: A. Quantum resource scaling · p. 7
  • The observed THC-rank-versus-orbital-number slope is 4.7; Table IV tracks high-spin CCSD(T) convergence by THC rank in the 58-orbital Cpd I space.The table highlights the largest THC rank used for resource estimates.Source: A. Quantum resource scaling · p. 8, Table IV · p. 8
  • Resource scaling is extrapolated from active-space estimates; Figure 5 plots Toffoli counts and logical-qubit counts against orbital number for four compounds and three factorizations.The plotted empirical scaling comes from least-squares fits on log-log plots.Source: Figure 5 · p. 8, A. Quantum resource scaling · p. 8

B. Compilation into Surface Code

Surface-code compilation yields resource estimates for phase estimation under specified hardware assumptions, with runtime scaling favoring tensor hypercontraction at small sizes and lower-error forecasts reducing requirements.

B. Compilation into Surface Code

  • For the largest X Hamiltonians, dynamic-correlation corrections may be small, making their spin-gap estimates more likely to be accurate.The runtime analysis begins with these systems.Source: B. Compilation into Surface Code · p. 8
  • The estimates assume a 0.1% per-gate error rate, a 1 µs code cycle, a 10 µs reaction time, and four Toffoli factories.Configurations minimize spacetime volume while correcting all errors at least 90% of the time; factory waiting time bottlenecks execution.Source: B. Compilation into Surface Code · p. 9
  • Tensor-hypercontraction qubitization shows an asymptotic runtime-scaling advantage even at small sizes, despite higher prefactors.At 0.1% physical error, four factories, and the stated timings, the largest-system estimate is 4,624,440 physical qubits and 73 hours; the resource estimates are upper bounds.Source: B. Compilation into Surface Code · p. 9, Figure 6 · p. 9
  • At 0.001% physical error, the largest system is estimated to require 500,000 qubits and 25 hours.The forecast assumes 1 µs surface-code cycles, 10 µs reaction times, and four Toffoli factories.Source: B. Compilation into Surface Code · p. 9, Figure 7 · p. 10

C. Demarcating the quantum advantage boundary

The estimated phase-estimation runtimes suggest a computational advantage over DMRG for large active spaces when computing ground-state energies, but the comparison depends on accuracy assumptions and does not directly provide other observables.

C. Demarcating the quantum advantage boundary

  • For large active spaces and bond dimensions M≥1000, phase estimation already has a computational advantage for ground-state energy simulation.The comparison estimates QPU runtimes with two Toffoli factories, corresponding to 4.9 million qubits and 135 hours for the largest systems.Source: C. Demarcating the quantum advantage… · p. 10, Figure 8 · p. 10
  • The DMRG accuracy comparison assumes that extrapolated energies are highly accurate, because neither DMRG nor coupled-cluster spin-gap calculations converged.The comparison is made against extrapolated DMRG energies for each active space.Source: C. Demarcating the quantum advantage… · p. 10
  • Phase estimation provides energy estimates only; dynamic-correlation corrections outside the active space usually require reduced density matrices.Extrapolating to 500 orbitals, the authors estimate approximately 9,000 logical qubits and 1.5 trillion Toffoli gates for phase estimation on the entire space.Source: C. Demarcating the quantum advantage… · pp. 10–11
  • The largest computational-basis-state overlap falls by a factor of two from 10 to 84 qubits, but is never considered problematic for phase estimation.The plotted overlaps use the DMRG M = 1500 wavefunction for three spin states.Source: C. Demarcating the quantum advantage… · p. 11, Figure 9 · p. 11

IV. NUMERICAL SUPPORT FOR A MULTIREFERENCE TREATMENT OF P450

The diagnostics give mixed evidence of multireference character: improved Kohn–Sham references reduce several single-reference warning signs, but Cpd I’s doublet retains a triradical signature. Reliable spin energetics therefore require multireference methods that balance dynamic and static correlation.

IV. NUMERICAL SUPPORT FOR A MULTIREFERENCE TREATMENT OF P450

  • Kohn–Sham orbitals reduce max(|t1|) for three compounds, but not for Cpd I, where a single-reference approach may struggle.With Hartree–Fock orbitals, every compound has max(|t1|) > 0.19; the diagnostic characterizes reference-orbital quality.Source: §IV · pp. 11–12
  • Kohn–Sham references nearly eliminate spin contamination in most cases, while the Cpd I doublet remains an exception.Hartree–Fock references show large spin contamination across all compounds and spin states.Source: §IV · p. 12
  • max(|t2|) stays near 0.05 across compounds and spin states, below the authors’ stronger-correlation marker of 0.1.The authors therefore find no evidence that the P450 models are strongly correlated by this diagnostic.Source: §IV · p. 12
  • Three natural orbital occupations near one identify triradical character in the Cpd I doublet, which a single determinant cannot faithfully represent.Spin-dependent properties such as the spin gap are often difficult to compute accurately without spin-pure states.Source: §IV · p. 12
  • Accurate spin energetics for the nearly degenerate Cpd I quartet and doublet require multireference methods that approach the exact solution.The authors emphasize that Cpd I is not strongly correlated in the sense of having many strongly entangled orbitals.Source: §IV · p. 12

V. CONCLUSION

The resource estimates indicate potential quantum advantage for CYP models at scales that balance dynamic and static correlation, while direct simulation of the full Cpd I system remains infeasible. The authors conclude that classical benchmarks and further fault-tolerant algorithm development are both needed to define and advance this computational frontier.

V. CONCLUSION

  • Classical and quantum resource assessments suggest potential quantum advantage even for problems that are not strongly correlated, when dynamic and static correlation must be balanced.Directly simulating the entire system to avoid partitioning those correlations is unlikely to be feasible with the estimated phase-estimation resources.Source: §V · p. 12
  • Cpd I’s doublet has three open shells, yet its ground-state spin and spin gaps remain unresolved because dynamic and static correlation must be treated together.The open shells involve Fe–O π* xz/yz orbitals and a sulfur-ligand lone pair; this triradical character agrees with prior calculations [26, 50, 69].Source: §V · p. 13, Figure 11 · p. 13
  • Qubitized phase estimation has more favorable scaling than DMRG with bond dimension, and the largest CYP models are estimated to run in under 100 hours.The full Cpd I extrapolation instead requires 1.5 trillion Toffoli gates, deemed unfeasible.Source: §V · p. 14
  • Tensor hypercontraction consistently outperformed the SF and DF factorizations in qubitized phase estimation.Surface-code compilation provided an upper-bound runtime estimate.Source: §V · pp. 13–14
  • Classical calculations help define chemistry’s computational frontier and guide quantum-algorithm research, while transformative simulations require further fault-tolerant algorithm development.The authors identify observable extraction and scaling as areas for development.Source: §V · pp. 12, 14

DATA AVAILABILITY

The paper makes its data, code, geometries, active-space Hamiltonians, and calculation inputs available through a public repository. Appendix A explains why correlation type and several complementary diagnostics matter when deciding how to model electronic structure.

DATA AVAILABILITY

  • A public Zenodo repository provides the study’s data and code, including geometries, active-space Hamiltonians, and calculation inputs.The passage also identifies the software and repository locations used for tensor factorizations and phase-estimation estimates.Source: DATA AVAILABILITY · p. 14

Appendix A: Electron correlation and the computational cost of electronic structure calculations

  • Dynamic correlation comprises many small contributions to the wave function and can be systematically improved by increasing the order of coupled-cluster theory.The number of excited determinants may still make these methods costly.Source: Appendix A · p. 17
  • Multireference systems need several comparably important determinants; balancing static and dynamic correlation matters because errors in either can distort spin-state ordering.The passage describes single-reference and multireference behavior as useful classification extremes.Source: Appendix A · p. 17
  • Because strong correlation cannot be identified reliably a priori and has no single clear measure, several probes should be considered together.Determining whether multireference methods are necessary itself poses a scoping challenge.Source: Appendix A · p. 17
  • Spin contamination can reflect essential multireference character or artificial symmetry breaking, while max(|t1|) primarily diagnoses reference-orbital quality.Spin symmetry may be restored by changing the SCF reference or using orbital-optimized methods with dynamic correlation.Source: Appendix A · pp. 17–18
  • Large max(|t2|) indicates multireference character, while correlated natural orbital occupations can reveal additional open shells.The paper uses max(|t2|) > 0.1 as an example threshold and κ-OOMP2 NOONs as an affordable diagnostic.Source: Appendix A · p. 18

Appendix B: Cytochrome P450

The appendix sets out CYP context and model construction alongside classical calculations and quantum-oriented tensor factorization. It describes how active spaces are selected and notes that NEVPT2 corrections can alter spin ordering in smaller Cpd I spaces.

Appendix B: Cytochrome P450

  • CYPs catalyze oxygen insertion into substrate C-H bonds, supporting drug detoxification; the catalytic cycle proceeds through reactive Compound I.The cycle returns to the resting state after alcohol dissociation and water coordination.Source: Appendix B · p. 18, Figure 12 · p. 19

a. Geometries

  • CYP3A4 model geometries came from water- and pyridine-bound crystal structures, with non-iron-coordinating protein and solvent removed before optimization.Cysteine was truncated to methyl-thiolate; the pentacoordinate model omits water, and Cpd I models omit both water-bound hydrogens.Source: a. Geometries · p. 20

a. Full valence space calculations

  • Full-valence κ-OOMP2 and CCSD calculations probe multireference character using diagnostics including CCSD amplitudes, spin contamination, and natural-orbital occupations.The calculations used the B3LYP high-spin geometry; CCSD tested UHF, UB3LYP, and UBP86 reference orbitals.Source: a. Full valence space calculations · p. 21

1. Active space calculations

  • DMRG energies for doublet, quartet, and sextet states were extrapolated, with NEVPT2 estimating dynamic correlation; active-space CCSD(T) calculations also compared spin states.For CCSD(T), orbitals were re-optimized within each active space for the different spin multiplicities.Source: §1 · p. 21

Appendix D: Details of DMRG calculations on Cpd I

  • DMRG cost empirically follows O(k^3) for fixed bond dimension over the studied sizes; active spaces A-C were exact, while X was too large for practical Cpd I calculations.The fixed-bond-dimension scaling is an estimate for the system sizes considered.Source: Table V · p. 21, Appendix D · p. 22

Appendix E: DMRG+NEVPT2 calculations on Cpd I

  • In active spaces D and E, NEVPT2 corrections can change the spin ordering, though the shift may reflect approximation in NEVPT2 or the DMRG wavefunction.The figure caption specifies that the corrected ordering places the high-spin state lowest.Source: Appendix E · p. 22, Figure 15 · p. 22

Appendix F: Tensor Hypercontraction by Optimized CP3

  • Optimized THC factorization targets both integral accuracy and a smaller central-tensor L1 norm; the rank meeting a CCSD(T) error below 1 mHartree was about five times the basis size.L1 regularization was used because unregularized optimization could produce erratically large central-tensor elements and fault-tolerant cost estimates.Source: Appendix F · p. 23

Appendix G: Active space sizes

  • The active-space hierarchy starts with five open-shell orbitals, then expands through iron and axial ligand orbitals to heme and more distal orbitals; the same orbitals were used across spin states.The hierarchy sought relevant ligand-binding and spin-character orbitals while maintaining roughly 50% filling; Table VI identifies the modeled resting, empty, Cpd1, and inhibited systems.Source: Appendix G · p. 24, Table VI · p. 24
Loading 2202.01244v1…