Manage a quantum-centric supercomputing workflow with QCSC-Prefect

Usage estimate: 30 seconds of QPU time on an IBM Quantum® Heron r2 processor. This estimate excludes queue time and classical processing. Your runtime might vary.

Learning outcomes

• Trace a quantum-centric supercomputing (QCSC) workflow in Prefect Cloud, with one task per Qiskit pattern stage.
• Configure the workflow by role, with run parameters in the Prefect Cloud run form and the high-performance computing (HPC) and quantum processing unit (QPU) settings in Prefect blocks.
• Resume a failed run without repeating its QPU job, and reuse the samples of one run in a new run.

Prerequisites

• Set up QCSC-Prefect on a Slurm cluster, which prepares the Slurm login node, the Python environments, and the IBM Quantum credentials.
• Sample-based quantum diagonalization of a chemistry Hamiltonian, which explains the algorithm.

Background

A QCSC workflow runs a QPU and an HPC system as one machine. QCSC-Prefect runs QCSC workflows on Prefect and keeps them traceable, configurable, and resumable.

This tutorial demonstrates these features with sample-based quantum diagonalization (SQD) of the N$_2$ molecule. The Quantum Resource Management Interface (QRMI) samples the circuit inside a Slurm job. The Fulqrum eigensolver diagonalizes the sampled subspaces in more Slurm jobs. One Prefect Cloud deployment runs both.

QCSC-Prefect also runs large QCSC workloads, as in the following projects.

• A study of protein-ligand complexes with more than 12,000 atoms used QCSC-Prefect. It sampled quantum-embedding fragments on two 156-qubit IBM Quantum processors, with up to 94 qubits, and diagonalized the fragment subspaces on the Fugaku, Miyabi-G, and ROQUO supercomputers [1].
• A study of iron-sulfur clusters ran an extended QCSC-Prefect SQD workflow with 36 orbitals (72 qubits) and subspace dimensions up to $3 \times 10^{10}$ on the Fugaku, Miyabi, and ROQUO supercomputers [2].

The challenge of QCSC workflows

The SQD workflow is typical of QCSC workloads. QCSC-Prefect answers each of its challenges with one feature, as in the following table.

Each step of Qiskit Pattern in QCSC-Prefect step by step ends with a box that shows how QCSC-Prefect traces, configures, and resumes that stage.

Traceable workflow

A Prefect flow is a Python function made of tasks. QCSC-Prefect makes each stage of a Qiskit pattern (Map, Optimize, Execute, and Postprocess) one task. For every run, Prefect Cloud shows the following information:

• The task graph, with an arrow wherever one task passes its result to another.
• The state, inputs, logs, and duration of each task, filterable by the stage tags map, optimize, execute, and post-process.
• The artifacts, which are tables and links that tasks publish, such as the source of the samples or the energy of each subspace.

The workflow host only submits and monitors work. The compute-intensive work runs as Slurm jobs, as in the following image.

Configurable parameters

QCSC-Prefect keeps each kind of setting in its own place, as in the following table.

The blocks describe each job in the three layers of the following table.

The flow loads blocks by name when it runs, so a change to a block leaves the flow and its run parameters unchanged.

Resumable execution

Prefect stores the results of the first three stages with the flow run. If a later task fails, Retry on the flow run reuses the stored results and runs only what failed, as in the following image.

A new run can also load the samples of an earlier run instead of sampling, because the Execute stage records where its samples are.

SQD at a glance

SQD samples bitstrings from a local unitary cluster Jastrow (LUCJ) circuit built from coupled cluster singles and doubles (CCSD) amplitudes, fixes the bitstrings that noise has corrupted (configuration recovery), and diagonalizes the Hamiltonian in subspaces spanned by the sampled configurations [3]. The SQD tutorial explains the algorithm. In QCSC-Prefect, SQD is the flow sqd, which has the following task graph.

When a run reuses samples, load_qrmi_result (or load_runtime_job for samples from IBM Quantum Platform) replaces both get_lucj_params and quantum_sampling, so the run has no Optimize task. In random mode (Quantum Source set to random), there is no Optimize or Execute task, and Map feeds sqd_loop directly.

The example is N$_2$ at a bond length of 1.1 $\AA$ in the 6-31G basis, with the two 1s core orbitals frozen. This leaves 16 spatial orbitals with five $\alpha$ and five $\beta$ electrons, which is a 32-qubit circuit.

Requirements

To run the code cells of this notebook, install the following packages in the notebook's kernel.

• Qiskit SDK v2.0 or later, with visualization support.
• Qiskit Runtime v0.40 or later (pip install qiskit-ibm-runtime).
• SQD Qiskit addon v0.14 or later (pip install qiskit-addon-sqd).
• ffsim v0.0.75 or later (pip install ffsim).
• Fulqrum v1.0 or later (see the Fulqrum installation guide).

Step 3 of this notebook samples with Qiskit Runtime, so it needs an IBM Quantum account saved on the computer that runs the notebook, as described in Save your login credentials.

Run the following steps one time. Each step names the node that it runs on.

Install the SQD workflow on the workflow host

Set the variables of the setup guide, activate ~/venv/qcsc-host, and install the SQD workflow and the Fulqrum integration.

# Run in a terminal on the workflow host.
# Optional if your ~/.bashrc sets the two variables
export QCSC_PREFECT="$HOME/qcsc-prefect"
export QCSC_WORK_DIR="/shared/work/qcsc-prefect"

source ~/venv/qcsc-host/bin/activate
cd "$QCSC_PREFECT"
uv pip install --no-deps -e packages/qcsc-prefect-fulqrum
uv pip install -e algorithms/qcsc_workflow_utility
uv pip install -e algorithms/open_loop/sqd
which sqd-create-blocks sqd-deploy

Install the Fulqrum solver on a compute node

Open a shell on a compute node. srun passes the variables of the workflow host shell to the compute node.

# Run in a terminal on the workflow host.
srun -N1 --pty bash

The QRMI worker transpiles the circuit on the compute node with the SQD transpilation hook, which imports only Qiskit. Add the hook to ~/venv/qcsc-node and check the imports. The python command prints nothing when the imports succeed.

# Run in the shell on the compute node.
source ~/venv/qcsc-node/bin/activate
cd "$QCSC_PREFECT"
uv pip install --no-deps -e algorithms/open_loop/sqd
python -c "import qrmi, qcsc_prefect.integrations.qiskit.qrmi.worker, sqd_common.transpile"

Build Fulqrum and qcsc-fulqrum-solve into their own environment, packages/qcsc-prefect-fulqrum/native/fulqrum. The build needs a C++ compiler with OpenMP. FULQRUM_WITH_PRIMME=1 also builds the PRIMME eigensolver, which needs a C compiler, Python headers, pkg-config, and a BLAS and LAPACK library with a .pc file.

# Run in the shell on the compute node.
cd "$QCSC_PREFECT"
FULQRUM_WITH_PRIMME=1 bash packages/qcsc-prefect-fulqrum/native/build_fulqrum.sh
realpath packages/qcsc-prefect-fulqrum/native/fulqrum/bin/qcsc-fulqrum-solve
exit

Record the printed path, because the solver blocks point to it. exit returns to the workflow host. For the build options, see Fulqrum integration.

Configure and create the Prefect blocks on the workflow host

Copy the example configuration.

# Run in a terminal on the workflow host.
cd "$QCSC_PREFECT"
cp algorithms/open_loop/sqd/sqd_blocks.example.toml \
   algorithms/open_loop/sqd/sqd_blocks.toml

Edit the settings at the top of sqd_blocks.toml as in the following table.

Keep the following points in mind when you edit the file.

• Write the paths out, because the file does not expand environment variables. The ~ character works. Keep every top-level key, such as python_executable, before the first table header.
• The [qrmi] table needs no change, because the qrmi-runtime block supplies the instance and the QPU. Set resource only to override the backend. The credentials stay in qrmi_config.json, not in this file.
• The slurm profile holds the scheduler and resource settings of a generic Slurm cluster, and any key in your file overrides it. The commented [sqd_advanced] table sets the sqd-advanced-settings block, the commented [advanced] table sets block names and timeouts, and the table of each solver holds its tuning. To choose the nodes, ranks, and threads of the solver jobs, see Size the SQD diagonalization for your HPC resources.
• The QRMI job and each solver job have a walltime of 12 hours (#SBATCH --time), and the flow waits up to 13 hours (timeout_seconds = 46800) for each job. To change the walltime, set walltime at the top of the file for the solver jobs, or in the [qrmi] table for the QRMI job, for example walltime = "02:00:00". A shorter walltime can start sooner on a busy cluster.

Create the blocks.

# Run in a terminal on the workflow host.
sqd-create-blocks --config algorithms/open_loop/sqd/sqd_blocks.toml

The command saves the blocks in the following table.

Add --dry-run to print the resolved settings without saving anything. prefect block ls and the Blocks page of Prefect Cloud list the blocks, as in the following image.

Deploy the workflow on the workflow host

sqd-deploy serves the flow from the workflow host and picks up the runs that you start in Prefect Cloud, so keep the terminal open while you run the workflow. The following command also writes the output to a new log file under $QCSC_WORK_DIR/logs/deploy at each start, named after the start time in UTC, for example sqd-deploy-20261109T051530Z.log.

# Run in a terminal on the workflow host.
mkdir -p "$QCSC_WORK_DIR/logs/deploy"
sqd-deploy 2>&1 | tee "$QCSC_WORK_DIR/logs/deploy/sqd-deploy-$(date -u +%Y%m%dT%H%M%SZ).log"

When the output shows that the flow sqd is served, the deployment appears in Prefect Cloud under Deployments. Press Ctrl+C to stop serving.

Setup

Import the dependencies.

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
from fulqrum.addon_sqd_support import sci_solver_from_integrals
from pyscf import ao2mo, cc, mcscf, scf
from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy
from qiskit_ibm_runtime import QiskitRuntimeService, Sampler

Qiskit Pattern in QCSC-Prefect step by step

The following steps compute the four stages in this notebook with PySCF, ffsim, Qiskit, Qiskit Runtime, and qiskit-addon-sqd, without Prefect. See the SQD tutorial for the details of each computation.

Step 1: Map classical inputs to a quantum problem

The Map stage, the task compute_molecular_integrals in the workflow, turns the molecule into the inputs of all later stages. Restricted Hartree-Fock (RHF) gives the molecular orbitals, and the active-space integrals leave out the two frozen 1s cores. CCSD gives the $t_2$ amplitudes of the ansatz, and its density matrix gives the orbital occupancies that start configuration recovery. CCSD iterates CCSD_MAXITER times, each time from the previous amplitudes, so the result does not depend on the DIIS history of the first iteration. This notebook uses 4 iterations, and the workflow uses 10.

ATOM = "N 0 0 0; N 0 0 1.1"
BASIS = "6-31g"
N_FROZEN = 2
CCSD_MAXITER = 4

mol = pyscf.gto.M(atom=ATOM, basis=BASIS, symmetry=False, verbose=0)
mf = scf.RHF(mol).run()
norb = mol.nao_nr() - N_FROZEN
nelec = (mol.nelec[0] - N_FROZEN, mol.nelec[1] - N_FROZEN)

# Active-space Hamiltonian with the 1s cores frozen
cas = mcscf.CASCI(mf, ncas=norb, nelecas=nelec)
one_body_tensor, core_energy = cas.get_h1eff()
two_body_tensor = ao2mo.restore(1, cas.get_h2eff(), norb)

# CCSD, iterated from the previous amplitudes
ccsd = cc.CCSD(mf, frozen=N_FROZEN)
eris = ccsd.ao2mo()
ccsd.kernel(eris=eris)
for _ in range(CCSD_MAXITER - 1):
    ccsd.kernel(ccsd.t1, ccsd.t2, eris=eris)
t2 = ccsd.t2

occupancies = np.diag(ccsd.make_rdm1())[N_FROZEN:] / 2
initial_occupancies = (occupancies, occupancies)

print(f"norb = {norb}, nelec = {nelec}")
print(f"E(RHF)  = {mf.e_tot:.8f} Ha")
print(f"E(CCSD) = {ccsd.e_tot:.8f} Ha")

Output:

norb = 16, nelec = (5, 5)
E(RHF)  = -108.86761837 Ha
E(CCSD) = -109.09392115 Ha

Step 2: Optimize problem for quantum hardware execution

The Optimize stage, the task get_lucj_params in the workflow, fits the ansatz to the QPU before any circuit exists. ffsim builds the LUCJ operator from $t_2$ with only the interactions that a heavy-hex lattice runs without SWAP gates. These are the same-spin pairs of neighboring orbitals and the opposite-spin pairs on every fourth orbital, at most five on a heavy-hex (Heron) processor. ffsim then compresses the operator (optimize=True, on by default in the workflow) so that one layer reproduces $t_2$ more closely.

N_REPS = 1

pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = [(p, p) for p in range(0, norb, 4)][:5]

lucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2,
    n_reps=N_REPS,
    interaction_pairs=(pairs_aa, pairs_ab),
    optimize=True,
    options={"maxiter": 50},
)
n_params = ffsim.UCJOpSpinBalanced.n_params(
    norb, N_REPS, interaction_pairs=(pairs_aa, pairs_ab)
)

print(f"Opposite-spin pairs: {pairs_ab}")
print(f"Ansatz parameters: {n_params}")

Output:

Opposite-spin pairs: [(0, 0), (4, 4), (8, 8), (12, 12)]
Ansatz parameters: 275

Step 3: Execute using Qiskit primitives

The Execute stage, the task quantum_sampling in the workflow, builds the circuit (Hartree-Fock state, LUCJ operator, and measurement), transpiles it for the QPU, and samples it. The ffsim PRE_INIT passes run first in the transpiler and turn the ffsim gates into efficient native-gate sequences. This notebook samples with Qiskit Runtime on ibm_fez. Replace ibm_fez with a backend that your instance can use. The workflow samples with QRMI instead, in a Slurm job that holds the QPU.

service = QiskitRuntimeService()
backend = service.backend("ibm_fez")

qreg = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qreg, ClassicalRegister(2 * norb, name="meas"))
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qreg)
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(lucj_op), qreg)
circuit.measure(qreg, circuit.cregs[0])

pass_manager = generate_preset_pass_manager(
    optimization_level=3, target=backend.target, seed_transpiler=24
)
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Backend: {backend.name}, ISA circuit depth: {isa_circuit.depth()}")

Output:

Backend: ibm_fez, ISA circuit depth: 640

SHOTS = 100_000

sampler = Sampler(mode=backend)
job = sampler.run([isa_circuit], shots=SHOTS)
quantum_bool_array = job.result()[0].data.meas.to_bool_array()
print(f"Job {job.job_id()}: {len(quantum_bool_array)} shots")

Output:

Job dbpl975t7fqc7389koog: 100000 shots

As a baseline, draw the same number of uniform random bitstrings, as the workflow does with Quantum Source set to random. The following cell also prints the fraction of each sample set with the correct electron numbers. Step 4 compares SQD on the two sample sets.

random_bool_array = generate_bit_array_uniform(
    SHOTS, 2 * norb, rand_seed=24
).to_bool_array()


def correct_electron_fraction(bool_array):
    """Fraction of the bitstrings with the correct numbers of alpha and beta electrons."""
    # The right half of each bitstring is alpha, the left half is beta.
    n_alpha, n_beta = nelec
    valid = (bool_array[:, norb:].sum(axis=1) == n_alpha) & (
        bool_array[:, :norb].sum(axis=1) == n_beta
    )
    return valid.mean()


print("Correct electron numbers:")
print(
    f"  Quantum samples: {correct_electron_fraction(quantum_bool_array):.4f}"
)
print(
    f"  Random samples:  {correct_electron_fraction(random_bool_array):.4f}"
)

Output:

Correct electron numbers:
  Quantum samples: 0.0336
  Random samples:  0.0040

Step 4: Postprocess and return result in desired classical format

The Postprocess stage, the task sqd_loop in the workflow, runs self-consistent configuration recovery. Each round corrects the sampled bitstrings with the current orbital occupancies, diagonalizes the Hamiltonian in several batches of configurations, and updates the occupancies from the best result. The following cell runs diagonalize_fermionic_hamiltonian on both sample sets with the same settings. The subspace dimension is the workflow's 1,000,000 determinants (1,000 strings per spin). The notebook uses two batches per round and ten rounds, and the workflow uses four batches and five rounds. As in the workflow, Fulqrum diagonalizes each subspace. sci_solver_from_integrals builds the Fulqrum Hamiltonian one time and returns the sci_solver that the function calls for each batch. A TrimPolicy, as in the SQD tutorial, merges the heaviest strings of the batches into one subspace and diagonalizes it. Of the strings above carryover_threshold, the policy keeps the heaviest half (carryover_ratio=0.5) for the next round, like the workflow's Carryover Fraction.

SUBSPACE_DIM = 1_000_000
NUM_BATCHES = 2
MAX_ITERATIONS = 10

policy = TrimPolicy(
    trim_ratio=0.8,
    carryover_threshold=1e-4,
    carryover_ratio=0.5,
)
sci_solver = sci_solver_from_integrals(one_body_tensor, two_body_tensor)


def run_sqd(bool_array, label):
    """Run SQD on the samples, and return the final energy and the energy of each round."""
    round_energies = []

    def callback(results):
        round_energies.append(
            min(result.energy for result in results) + core_energy
        )
        print(
            f"{label}, round {len(round_energies)}: {round_energies[-1]:.8f} Ha"
        )

    result = diagonalize_fermionic_hamiltonian(
        one_body_tensor,
        two_body_tensor,
        bool_array,
        samples_per_batch=1000,
        norb=norb,
        nelec=nelec,
        num_batches=NUM_BATCHES,
        max_iterations=MAX_ITERATIONS,
        policy=policy,
        sci_solver=sci_solver,
        symmetrize_spin=True,
        max_dim=int(np.sqrt(SUBSPACE_DIM)),
        initial_occupancies=initial_occupancies,
        callback=callback,
        seed=np.random.default_rng(24),
    )
    return result.energy + core_energy, round_energies


quantum_energy, quantum_round_energies = run_sqd(
    quantum_bool_array, "Quantum samples"
)
random_energy, random_round_energies = run_sqd(
    random_bool_array, "Random samples"
)
print(f"SQD energy with quantum samples: {quantum_energy:.8f} Ha")
print(f"SQD energy with random samples:  {random_energy:.8f} Ha")

Output:

Quantum samples, round 1: -109.09698391 Ha
Quantum samples, round 2: -109.10046517 Ha
Quantum samples, round 3: -109.10193633 Ha
Quantum samples, round 4: -109.10186225 Ha
Quantum samples, round 5: -109.10196116 Ha
Quantum samples, round 6: -109.10190268 Ha
Quantum samples, round 7: -109.10199765 Ha
Quantum samples, round 8: -109.10196037 Ha
Quantum samples, round 9: -109.10199167 Ha
Quantum samples, round 10: -109.10205287 Ha
Random samples, round 1: -109.09582359 Ha
Random samples, round 2: -109.09962562 Ha
Random samples, round 3: -109.10059383 Ha
Random samples, round 4: -109.10115380 Ha
Random samples, round 5: -109.10131183 Ha
Random samples, round 6: -109.10148002 Ha
Random samples, round 7: -109.10154865 Ha
Random samples, round 8: -109.10159227 Ha
Random samples, round 9: -109.10161699 Ha
Random samples, round 10: -109.10158185 Ha
SQD energy with quantum samples: -109.10205287 Ha
SQD energy with random samples:  -109.10161699 Ha

The plot compares both sample sets, round by round, with the full configuration interaction (FCI) energy of this active space, $-109.10336546$ Ha, which PySCF computes in about a minute (cas.run() on the preceding CASCI object). The dashed line is chemical accuracy, 1.6 mHa.

E_FCI = -109.10336546

# SQD stops early if it converges, so each run can have fewer rounds.
for label, round_energies in [
    ("Quantum samples", quantum_round_energies),
    ("Random samples", random_round_energies),
]:
    rounds = range(1, len(round_energies) + 1)
    plt.semilogy(
        rounds, np.abs(np.array(round_energies) - E_FCI), "o-", label=label
    )
plt.axhline(
    abs(ccsd.e_tot - E_FCI), color="tab:gray", linestyle=":", label="CCSD"
)
plt.axhline(
    1.6e-3, color="tab:red", linestyle="--", label="Chemical accuracy"
)
plt.xlabel("Recovery round")
plt.ylabel("|E - E(FCI)| (Ha)")
plt.xticks(range(1, MAX_ITERATIONS + 1))
plt.legend()
plt.show()

Output:

Summary of the stages

The following image summarizes the four stages.

Run the full workflow

The blocks from Configure and create the Prefect blocks on the workflow host already describe the cluster and the QPU, so a run needs only the run parameters.

Submit the run

In Prefect Cloud, go to Deployments, select sqd, then select Run -> Custom run. Set FCIDump File, the only field without a default, as in the following table. The other rows show the defaults. The form marks the three SQD fields as required, although they have defaults. Select Submit.

A path relative to the ~/qcsc-prefect checkout also works. algorithms/open_loop/sqd/default_params/n2.json holds the same parameters for scripted runs. Checkpoint Recovery adds a folder of less than 1 MB per run under work_dir/sqd_checkpoints/, even for 1000 orbitals, and nothing removes it.

Read the results

The flow run page shows the task graph with the state of each task, the logs, and, on the Artifacts tab, the five artifacts in the following table.

The energy that the flow reports is the lowest energy in the sqd-configuration-recovery table. Each solver job leaves a job_<timestamp>_<id> directory under work_dir. The directory holds the rendered Slurm script, input.npz, solver.json, and the Slurm output. It also holds output.npz and result.json on success, or error.json on failure. Use sacct to see the jobs of the run.

# Run in a terminal on the workflow host.
sacct --format=JobID,JobName,State,ExitCode,Elapsed,WorkDir%60

Optional: resume and reuse the workflow

Retry a failed run, and reuse its samples

Retry: Suppose that, in a run with Max Iteration set to 10, a solver job hits its walltime in round 5 and sqd_loop fails. Raise the walltime in the exec-sqd-fulqrum-slurm block under Blocks, then select Retry on the failed flow run. The retried run loads the stored results of the first three stages and submits no QPU job. It runs only sqd_loop again, with the new walltime, and starts at round 5 from the checkpoint of round 4. Its log shows [ckpt] RESUME from …: completed 4 step(s), then SQD iteration 5/10, as in the following image.

Reuse: A new run computes everything again unless you point it to earlier samples. To try another subspace dimension with the same samples, start a custom run with the same parameters, set Load from QRMI Job Directory (qrmi_result_dir) to the work directory in the sqd-sampling-source artifact of the first run, and change Subspace Dimension. The new run replaces get_lucj_params and quantum_sampling with load_qrmi_result, which reads the samples and submits nothing to the QPU. Reset Mitigation in sqd-advanced-settings must match the first run, because the reload applies it. For samples taken through IQP, set Load from Qiskit Runtime Job ID (qiskit_runtime_job_id) instead, and the stage is load_runtime_job. It reads the job with the credentials of the prefect-qiskit QuantumRuntime block named in Quantum Runtime Block (reload) (ibm-runner by default). Create that block one time with prefect block create quantum-runtime.

Continue a recovery with more rounds: Start a custom run with the same parameters, set Resume from Checkpoint (resume_checkpoint) to the path in the sqd-checkpoint artifact of the first run, and raise Max Iteration. The new run continues after the completed rounds, from the stored occupancies and best energy. For example, a checkpoint with 32 of 100 rounds done and Max Iteration set to 200 continues at round 33/200. The new run reads that checkpoint and writes its own, and the following rules apply.

• The samples, the molecule, Subspace Dimension, Batch Number, and the solver type (the fulqrum_solver_job part of the solver block) must match the first run, so reuse its samples as in Reuse.
• The name and the job settings of the solver block, such as its walltime, can change.
• If a setting that must match differs, the run fails before any solver job and names the difference.

Next steps

References

[1] Merz Jr., K. M., Shajan, A., Kaliakin, D., Liang, F., Otsuka, Y., Shirakawa, T., et al. (2026). Crossing the 12,000-atom barrier with heterogeneous quantum-classical supercomputing: quantum chemistry of protein-ligand complexes. arXiv:2605.01138

[2] Kumar, K., Shirakawa, T., Kim, B., Wakizaka, R., Hamamura, I., Sato, M., et al. (2026). Broken-symmetry reference states for sample-based quantum diagonalization with near-DMRG accuracy. arXiv:2611.04496

[3] Robledo-Moreno, J., Motta, M., Haas, H., Javadi-Abhari, A., Jurcevic, P., Kirby, W., et al. (2025). Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer. Science Advances, 11(25), eadu9991. DOI, arXiv:2405.05068