Crossing the 12,000-atom barrier with heterogeneous quantum-classical supercomputing: quantum chemistry of protein-ligand complexes

Motivation

• Chemical-accuracy predictions for biomolecules require scalable methods that can treat thousands of atoms while retaining accurate energy differences.This is especially relevant to protein-ligand binding energies and structure-based drug discovery.Source: §II · p. 1
• Wavefunction methods offer systematically improvable accuracy, but solving the many-electron Schrödinger equation remains computationally costly and limits accessible system sizes.DFT is more scalable but may lack the accuracy required for drug-discovery applications.Source: §II · p. 1, §III · p. 2

Approach

• Quantum embedding partitions a molecule into independently treated fragments, allowing quantum processors to act as specialized co-processors within a classical workflow.The approach exploits spatial locality and coarse-grained parallelism.Source: §II · p. 2

Challenges

• Embedding at biomolecular scale requires compact accurate fragment bases and systematically improvable fragment solutions without overwhelming memory bandwidth, accelerators, or inter-node communication.These constraints motivate algorithmic and software advances in the HQC workflow.Source: §II · p. 2

Limitations

• O(M^5) fragment-construction cost and O(M) fragment growth make biochemical-scale EWF calculations computationally demanding.Conventional EWF also faces inaccurate SQD and ExtSQD fragment solutions and requires GPU-accelerated distributed diagonalization.Source: §IV · pp. 3–4

Innovations

• The workflow combines lower-scaling EWF with higher-accuracy TrimSQD fragment solutions to improve scalability and accuracy.These innovations target the principal limitations of conventional EWF with HQC fragment solvers.Source: §IV · p. 4

Impact

• 12,635 atoms represents more than 40× the previous HQC size of 303 atoms.The calculations use quantum embedding and TrimSQD fragment solutions.Source: §IV · p. 4
• 0.1 mHa deviation from DMRG for the largest trypsin fragment corresponds to a 210× accuracy improvement.The largest Trp-cage fragment has a 21 mHa deviation.Source: §IV · p. 4

Input and decomposition

• O(M^3) Hartree–Fock generates molecular orbitals that serve as quantum-embedding inputs.Here, M is the number of orbitals in the full system.Source: A. Overarching algorithm · p. 4
• EWF decomposes the molecule into fragments, with construction reduced from O(M^5) to O(1) per fragment and parallelized across fragments.The reduced-scaling construction is described as embarrassingly parallel.Source: A. Overarching algorithm · p. 4

Fragment solvers

• Fragments use FCI for small hydrogen-centered cases, TrimSQD for heavier-atom cases, and CCSD as a classical comparison.TrimSQD is used for fragments with 13 ≤ M_c ≤ 45.Source: A. Overarching algorithm · p. 4

HQC fragment calculation

• Quantum sampling, subspace diagonalization, and subspace extension generate fragment solutions before parallel energy reconstruction.Quantum circuits produce electronic configurations that guide subspace diagonalization.Source: A. Overarching algorithm · p. 4

Construction strategy

• The construction innovations eliminate data dependencies at full-system scale and make fragment processing embarrassingly parallel.They target the cost and memory bottlenecks of conventional EWF construction.Source: B. Innovations in EWF fragment… · p. 4

Localized MP2

• O(M^5) MP2 cost is reduced by restricting each fragment calculation to orbitals within a 7 Å radius.The restriction is centered on the atom defining the fragment.Source: 1) Spatially localized orbital… · p. 5

Localized ERIs

• O(M^4) ERI construction and its large memory and I/O requirements are addressed by restricting ERIs to orbitals within Rcut + 3 Å.The buffer is intended to avoid artificial boundary effects while controlling computational cost.Source: 2) Localized electron repulsion… · p. 5

Impact

• Localized orbital truncation and ERI construction reduce fragment-construction cost to O(1) per fragment.The resulting operation is embarrassingly parallel.Source: 2) Localized electron repulsion… · p. 5

TrimSQD workflow

• TrimSQD introduces determinant trimming based on diagonalization results into ExtSQD to improve configuration selection.The approach draws on closed-loop SQD and TrimCI.Source: 1) TrimSQD · p. 5

Configuration recovery

• The configuration-recovery loop merges recovered configurations, partitions them across subgroups, and diagonalizes each subgroup.The loop carries selected configurations between iterations.Source: 1) TrimSQD · p. 5

Configuration trimming

• Top k1-percent configurations by contribution to |Ψ_x| are retained from subgroup diagonalizations, followed by global retention of the top k2 percent.The global subspace diagonalization spans all MPI ranks.Source: 1) TrimSQD · p. 6

Subspace extension

• The subspace is extended with configurations coupled to retained configurations by H_xy ≠ 0 and Hamming distance at most two.Final diagonalization retains configurations satisfying |Ψ_x| > ε.Source: 1) TrimSQD · p. 6
• The final expanded-subspace diagonalization typically dominates execution time because its dimension becomes very large.It is performed on all MPI ranks using the user-defined threshold ε.Source: 1) TrimSQD · p. 6

• Subspace diagonalization dominates ExtSQD and TrimSQD computational cost, motivating an implementation tailored for large-scale parallel execution.Source: 2) Selected-basis diagonalization · p. 6
• SBD partitions each vector across computing nodes, reducing per-node memory requirements and enabling scalable large-problem computations.Source: 2) Selected-basis diagonalization · p. 6
• GPU acceleration redesigns the dominant matrix–vector multiplication kernel with fine-grained data-parallel primitives and GPU-resident layouts.Source: 2) Selected-basis diagonalization · p. 6
• TrimSQD introduces interactions between subgroups of otherwise independent configurations, improving configuration selection and fragment-energy accuracy.Source: 2) Selected-basis diagonalization · p. 6
• The experiments benchmark two contrasting protein–ligand complexes to assess the workflow’s portability.Source: A. Benchmark molecular systems · p. 6
• The trypsin–benzamidine benchmark contains 12,635 atoms and 31,795 molecular orbitals, while T4-Lysozyme–n-butyl-benzene contains 11,608 atoms and 28,844 molecular orbitals.Source: 1) trypsin, benzamidine · p. 6, 2) T4-Lysozyme, n-butyl-benzene · p. 6

• Quantum sampling is distributed across ibm_cleveland and ibm_kobe, two 156-qubit Heron r2 QPUs.Source: 1) Quantum sampling · p. 7
• ExtSQD and TrimSQD subspace diagonalization run on both Fugaku and Miyabi-G.Source: 2) Subspace diagonalization · p. 7
• Fugaku provides 158,976 Arm-based nodes with distributed-memory interconnects for large-scale parallel computation.Source: 2) Subspace diagonalization · p. 7
• Miyabi-G provides 1,120 GH200 nodes, each combining a 72-core Grace CPU with an NVIDIA Hopper H100 GPU.Source: 2) Subspace diagonalization · p. 7
• Classical timing uses synchronized MPI wall-clock measurements, while QPU runtime includes circuit execution and inter-circuit intervals but excludes compilation and network transmission.Source: C. Measurement methodology · p. 7
• QPU performance is characterized by qubit coherence times and 1- and 2-qubit gate and measurement error rates.Source: C. Measurement methodology · p. 7

• The software stack combines ORCA for Hartree–Fock, in-house EWF and diagonalization codes, and PySCF, DICE, and Block2 for classical comparisons.Source: D. Setting of experiments · p. 7
• Fugaku uses ExtSQD for fragments with Mc ≤26 and TrimSQD for fragments with Mc ≥27, whereas Miyabi-G uses TrimSQD for all fragments.Source: D. Setting of experiments · pp. 7–8
• Detailed performance and accuracy analyses focus on three representative trypsin fragments.Source: Table VI · p. 7, D. Setting of experiments · p. 8
• EWF-CCSD and EWF-TrimSQD produce comparable ground-state and binding-energy calculations for both benchmark complexes.Source: A. Results · p. 8
• The binding energies are too positive, with the minimal STO-3G basis and selected η values identified as likely causes.Source: A. Results · p. 8, Table VII · p. 8
• Future improvements include more realistic compact basis sets and a more flexible η choice to increase computational cost and accuracy selectively in the binding region.Source: A. Results · p. 8
• Reported timing distinguishes total QPU time in hours from GPU-hours on Miyabi-G.Source: A. Results · p. 8, Table VIII · p. 8

• Quantum sampling uses up to 94 qubits on ibm_cleveland and sustains performance for 56 h on trypsin and 44 h on T4-Lysozyme across ibm_cleveland and ibm_kobe.Source: B. QPU capacity · p. 8
• Fugaku and Miyabi-G sustain 95.7% and 98.6% of their available nodes, respectively, during EWF fragment calculations.Source: C. HPC capacity · p. 8
• Using two QPUs and two leadership-class supercomputers demonstrates workflow portability and flexibility across CPU- and GPU-based platforms.Source: C. HPC capacity · p. 8
• Parallelization targets subspace diagonalization and its dominant matrix–vector multiplication routine, based on strong-scaling analysis.Source: D. Scaling · p. 8, Figure 4 · p. 9
• 72.5% parallel efficiency is achieved with Nnode = 64 in the practically relevant production regime.Source: 1) Strong scaling of TrimSQD · pp. 8–9, Figure 3 · p. 9
• Increasing node counts reduces time-to-solution but eventually deviates from ideal linear scaling, making moderate node counts preferable for these problems.Source: 1) Strong scaling of TrimSQD · pp. 8–9

Strong scaling of matrix-vector multiplication:

• ≈11 s or less per matrix-vector operation is achieved for a subspace dimension of 4.29 · 10^9.Increasing node counts eventually yields diminishing time-reduction returns because scaling becomes sublinear.Source: 2) Strong scaling of matrix-vector… · p. 9
• 72.5% parallel efficiency enables scalable distributed execution of EWF fragment calculations across Fugaku and Miyabi-G.Strong scaling determines the optimal parallelization of subspace diagonalizations across these large-scale platforms.Source: 2) Strong scaling of matrix-vector… · p. 9
• SBD-G supports accurate fragment treatment alongside strong scaling, motivating the assessment of EWF fragment energies after performance evaluation.The workflow treats performance gains as meaningful only when the underlying fragments are accurately solved.Source: E. Accuracy · p. 9
• TrimSQD’s accuracy advantage over ExtSQD becomes more pronounced as fragment size increases, confirming benefits from its subspace-diagonalization innovations.The comparison evaluates energy differences for 403 non-water fragments with large molecular-orbital spaces.Source: 1) Comparison between TrimSQD and ExtSQD · pp. 9–10

Accuracy-time tradeoff in TrimSQD and ExtSQD:

• TrimSQD achieves high accuracy with substantially reduced time-to-solution versus ExtSQD and selected CPU-based SCI and DMRG implementations.SBD-G avoids storing the Hxy matrix and distributed Ψx vectors, reducing memory demands for large-scale diagonalization.Source: 2) Accuracy-time tradeoff in TrimSQD… · p. 10
• 12,000 atoms are exceeded in two protein-ligand complex simulations using quantum embedding and heterogeneous quantum-classical fragment solvers.The approach combines lower-scaling fragment construction with algorithmic and implementation advances in fragment solutions.Source: §VII · p. 10
• 1.3·10^9 measurement outcomes from 9,200 quantum circuits support the reported large-scale fragment calculations.The circuits include 1,016 runs using 58 to 94 qubits.Source: §VII · p. 10
• The study establishes a practical baseline for future HQC applications in biochemistry and other systems with comparable size and chemical intricacy.Proposed applications include covalent-binding drugs, catalytic reactions, transition states, and phase diagrams in quantum liquids and solids.Source: §VII · p. 10