- Demos/
- Algorithms/
Quantifying Resource Requirements for Vibronic Dynamics Simulation
Quantifying Resource Requirements for Vibronic Dynamics Simulation
Published: May 05, 2024. Last updated: September 22, 2025.
Standard quantum chemistry often relies on the adiabatic Born-Oppenheimer approximation, assuming nuclei move on a single potential energy surface (PES) well-separated from others. But in the world of photochemistry, where molecules absorb light, and break bonds, this approximation breaks down. The timescales of electronic and nuclear motion become comparable, giving rise to vibronic coupling. This coupling, drives critical processes like photosynthesis, and solar cell efficiency and accurately simulating this requires a Hamiltonian that treats the dynamics of electrons and nuclei simultaneously.
However, modeling such dynamics is computationally demanding. Before we attempt to run such complex algorithms on future fault-tolerant hardware, we need to answer a fundamental question: Is this algorithm actually feasible?
In this demo, we assume the role of an algorithm architect. We will construct a simulation workflow
for the vibronic Hamiltonian from the ground up, utilizing the modular building blocks within PennyLane’s estimator.
By constructing this pipeline ourselves, we can perform a feasibility analysis, determining exactly what resources are required to simulate
these complex non-adiabatic processes on future hardware.
The Vibronic Hamiltonian
To perform this feasibility check, we must first define the system. We use the vibronic coupling Hamiltonian, which describes a set of \(N\) electronic states interacting with \(M\) vibrational modes.
Unlike standard electronic structure problems, this model mixes discrete electronic levels with continuous vibrational motion, and the Hamiltonian takes the form:
where \(T\) represents the kinetic energy of the nuclei, and is diagonal in the electronic subspace while \(V\) contains off-diagonal electronic couplings.
To simulate the time evolution \(U(t) = e^{-iHt}\), we employ a product formula (Trotterization) approach as outlined in D. Motlagh et al. [1]. A key challenge in Trotterization is decomposing the Hamiltonian efficiently. Here, we use the fragmentation scheme proposed in the reference, grouping the terms \(\ket{j} \bra{i} \otimes V_{ji}\) such that they differ by a fixed bitstring \(m\). The original grouping method results in N different fragments, which can be viewed as blocks of the potential energy matrix as shown in Figure:1. Each of these fragments can then be block-diagonalized by using only Clifford gates and implemented as a sequence of evolutions controlled by the corresponding electronic states.
Figure 1: Fragmentation of the Vibronic Hamiltonian into N blocks based on bitstring.¶
Defining the Hamiltonian
With the algorithmic approach defined, the next step is to instantiate the system we wish to benchmark.
While access to the full Hamiltonian coefficients would allow us to further optimize costs by leveraging the commutativity of electronic parts, we can still derive a reliable baseline using only structural parameters. By defining the number of modes, electronic states, and grid size, we can map out the cost topology of a real system without needing to generate and store the full Hamiltonian.
As a case study, we select (NO):math:_4-Anth, a molecule created by introducing four N-oxyl radical fragments to anthracene. Proposed for its theoretically record-breaking singlet fission speed [2], we use this system to define our simulation parameters:
num_modes = 19 # Number of vibrational modes
num_states = 5 # Number of electronic states
k_grid = 4 # Number of qubits per mode (discretization)
taylor_degree = 2 # Truncate to Quadratic Vibronic Coupling
In our model, we truncate the interaction terms for potential energy fragment to linear and quadratic terms only. This approximation, known as the Quadratic Vibronic Coupling (QVC) model, captures the dominant physical effects while simplifying the potential energy circuits significantly.
Constructing Circuits for Single Time-Step
Next step is to define what the circuit will look like for evolving a single time step. Based on the term-based fragmentation scheme, the single Trotter step is composed of two distinct types of quantum circuits interleaved together:
Kinetic Energy Fragment: This implements the nuclear kinetic energy. Since the kinetic operator depends on momentum \(P\), this circuit uses the Quantum Fourier Transform (QFT) to switch to the momentum basis, applies a phase rotation, and then switches back.
Potential Energy Fragments: These implement the interaction terms. For each fragment, the algorithm loads coefficients using a QROM (Quantum Read-Only Memory), computes the vibrational monomial product using quantum arithmetic, and applies a phase gradient.
Let’s first see what the circuit will look for the kinetic energy fragment:
import pennylane.estimator as qre
def kinetic_circuit(mode_wires, phase_wires, scratch_wires, coeff_wires):
ops = []
num_phase_wires = len(phase_wires)
k_grid = len(mode_wires[0])
for mode in range(num_modes):
ops.append(qre.AQFT(order=1, num_wires=k_grid, wires=mode_wires[mode]))
for i in range(num_modes):
ops.append(
qre.OutOfPlaceSquare(
register_size=k_grid, wires=scratch_wires + mode_wires[i]
)
)
for j in range(2 * k_grid):
ctrl_wire = [scratch_wires[j]]
target_wires = (
coeff_wires[: len(phase_wires) - j]
+ phase_wires[: len(phase_wires) - j]
)
ops.append(
qre.Controlled(
qre.SemiAdder(max_register_size=num_phase_wires - j),
num_ctrl_wires=1,
num_zero_ctrl=0,
wires=target_wires + ctrl_wire,
)
)
ops.append(
qre.Adjoint(
qre.OutOfPlaceSquare(
register_size=k_grid, wires=scratch_wires + mode_wires[i]
)
)
)
for mode in range(num_modes):
ops.append(
qre.Adjoint(
qre.AQFT(order=1, num_wires=k_grid, wires=mode_wires[mode])
)
)
return qre.Prod(ops)
Similarly, we can define the structure for the potential energy fragments. For a QVC model truncated to quadratic terms, each fragment will implement a monomial of degree up to 2
where \(c_r\) and \(\tilde{c}_{r}\) are the linear and quadratic coupling coefficients respectively. The exponential of each term is implemented by the modular circuit shown in Figure 2.
Figure 2: Circuit for implementing the exponential of a monomial term in the potential energy fragment.¶
The circuit executes the following steps:
Loading the coefficients from QROM controlled by the electronic state.
Computing the monomial product by using a square operation for quadratic terms.
Accumulating the phase by decomposing the multiplication into a series of controlled adders that add the product directly onto the phase gradient register.
Uncomputing the intermediate steps to clean up ancilla wires.
We can define this circuit using two different segments, one for linear terms and one for quadratic terms:
def linear_circuit(num_states, elec_wires, phase_wires, coeff_wires, scratch_wires):
ops = []
ops.append(
qre.QROM(
num_bitstrings=num_states,
size_bitstring=len(phase_wires),
restored=False,
wires=elec_wires + coeff_wires,
)
)
for i in range(k_grid):
ctrl_wire = [scratch_wires[i]]
target_wires = (
coeff_wires[: len(phase_wires) - i] + phase_wires[: len(phase_wires) - i]
)
ops.append(
qre.Controlled(
qre.SemiAdder(max_register_size=len(phase_wires) - i),
num_ctrl_wires=1,
num_zero_ctrl=0,
wires=target_wires + ctrl_wire,
)
)
ops.append(
qre.Adjoint(
qre.QROM(
num_bitstrings=num_states,
size_bitstring=len(phase_wires),
restored=False,
wires=elec_wires + coeff_wires,
)
)
)
return qre.Prod(ops)
def quadratic_circuit(
num_states, elec_wires, phase_wires, coeff_wires, mode_wires, scratch_wires
):
ops = []
k_grid = len(mode_wires)
ops.append(
qre.QROM(
num_bitstrings=num_states,
size_bitstring=len(phase_wires),
restored=False,
wires=elec_wires + coeff_wires,
)
)
qre.OutOfPlaceSquare(register_size=k_grid, wires=mode_wires + scratch_wires)
for i in range(2 * k_grid):
ctrl_wire = [scratch_wires[i]]
target_wires = (
coeff_wires[: len(phase_wires) - i] + phase_wires[: len(phase_wires) - i]
)
ops.append(
qre.Controlled(
qre.SemiAdder(max_register_size=len(phase_wires) - i),
num_ctrl_wires=1,
num_zero_ctrl=0,
wires=target_wires + ctrl_wire,
)
)
ops.append(
qre.Adjoint(
qre.OutOfPlaceSquare(
register_size=k_grid, wires=mode_wires + scratch_wires
)
)
)
ops.append(
qre.Adjoint(
qre.QROM(
num_bitstrings=num_states,
size_bitstring=len(phase_wires),
restored=False,
wires=elec_wires + coeff_wires,
)
)
)
return qre.Prod(ops)
Finally, to ensure precise resource tracking, we explicitly label our wire registers. This avoids ambiguity about which circuit operations are mapped to which quantum registers and ensures the resource estimator captures the full width of the circuit.
def get_wire_labels(num_modes, num_states, k_grid, phase_prec):
"""Generates the wire map for the full system."""
num_elec_qubits = (num_states - 1).bit_length()
elec_wires = [
f"e_{i}" for i in range(num_elec_qubits)
] # Electronic State Register
phase_wires = [
f"pg_{i}" for i in range(phase_prec)
] # Resource State For Phase Gradients
coeff_wires = [f"c_{i}" for i in range(phase_prec)] # Coefficient Register
mode_wires = []
for m in range(num_modes):
mode_wires.append([f"m{m}_{w}" for w in range(k_grid)]) # Mode m Register
scratch_wires = [
f"s_{i}" for i in range(2 * k_grid)
] # Scratch Space for Arithmetic
return elec_wires, phase_wires, coeff_wires, mode_wires, scratch_wires
We now define the full circuit by assembling our fragments into a second-order Suzuki-Trotter expansion
To construct the expansion, we loop over all fragments: the kinetic fragment is added once, while the number of potential fragments is determined by the size of the electronic register. This ensures that every block of the Hamiltonian identified in our fragmentation scheme (Figure 1) is accounted for in the simulation.
Additionally, an important part of the algorithm is the preparation of phase gradient register. The size of this register is dictated by the desired precision for phase rotations in both the potential and kinetic steps, directly influencing the simulation’s accuracy.
import math
def circuit(num_modes, num_states, k_grid, taylor_degree, num_steps, phase_grad_prec=1e-6):
fragments = []
phase_grad_wires = int(math.ceil(math.log2(1 / phase_grad_prec)))
elec_wires, phase_wires, coeff_wires, mode_wires, scratch_wires = (
get_wire_labels(num_modes, num_states, k_grid, phase_grad_wires)
)
qre.PhaseGradient(num_wires=len(phase_wires), wires=phase_wires) # Prepare Phase Gradient State
kinetic_fragment = kinetic_circuit(
mode_wires, phase_wires, scratch_wires, coeff_wires
)
fragments.append(kinetic_fragment) # Add Kinetic Fragment
num_fragments = 2**len(elec_wires)
for i in range(num_fragments):
frag_op = []
for mode in range(num_modes):
if taylor_degree >= 1:
frag_op.append(
linear_circuit(
num_states, elec_wires, phase_wires, coeff_wires, scratch_wires
)
)
if taylor_degree >= 2:
frag_op.append(
quadratic_circuit(
num_states,
elec_wires,
phase_wires,
coeff_wires,
mode_wires[mode],
scratch_wires,
)
)
fragments.append(qre.Prod(frag_op))
qre.TrotterProduct(first_order_expansion=fragments, num_steps=num_steps, order=2)
Finally, we can estimate the resource requirements for this full circuit using the estimate function. To ensure the simulation reaches a total time of 100 fs with sufficient accuracy, we set the number of Trotter steps to 500. This value is chosen based on benchmark data from Motlagh et al. [1], where this step count was shown to maintain reliable accuracy for a 100 fs simulation of the (NO)$_4$-Anth system.
print(qre.estimate(circuit)(num_modes, num_states, k_grid, taylor_degree, num_steps=500))
--- Resources: ---
Total wires: 146
algorithmic wires: 127
allocated wires: 19
zero state: 19
any state: 0
Total gates : 5.210E+8
'Toffoli': 6.082E+7,
'T': 749,
'CNOT': 2.712E+8,
'X': 4.560E+6,
'Z': 1,
'S': 1,
'Hadamard': 1.844E+8
We observe that while our results follow the same scaling and orders of magnitude reported in Motlagh et al. [1], the Toffoli gate count is slightly higher. The higher gate counts in this estimation occur because we assume a dense Hamiltonian, whereas the reference work leverages system-specific sparsity by only implementing non-zero coupling terms. These numbers can therefore be viewed as a reliable upper bound for the cost of simulating (NO):math:_4-Anth dynamics.
Conclusions
By constructing the full simulation workflow for vibronic dynamics from the ground up, we have gained critical insights into the resource requirements for simulating these complex processes on future quantum hardware. This approach allows us to break down complex algorithms into manageable building blocks, providing a transparent view of how algorithmic choices directly impact gate counts and qubit overhead. This demonstration highlights the power of modular quantum algorithm design, enabling researchers to build, analyze, and optimize complex workflows for next-generation quantum simulations. Ultimately, this framework is designed to be an open sandbox for researchers. We encourage you to use the building blocks provided here to explore alternative algorithms tailored to your specific research interests.
References
About the author
Diksha Dhawan
Developing Tools to Simulate Chemistry Using Quantum Computers
Total running time of the script: (0 minutes 0.192 seconds)
