Shell: VQE for Carbon Capture Materials and Seismic Optimization

The Problem

Carbon capture is one of the harder engineering challenges in the energy transition. Finding materials that selectively adsorb CO2 from flue gas requires identifying zeolites and metal-organic frameworks (MOFs) with the right pore geometry and binding energetics. Classical density functional theory (DFT) can compute adsorption energies, but it struggles with strongly correlated electronic systems, exactly where quantum simulation has a theoretical edge.

Shell also manages large-scale seismic data processing for exploration. Interpreting reflected wave patterns to infer subsurface structure is computationally intensive and involves signal optimization problems that could benefit from quantum approaches. In this work, Shell focused primarily on the materials simulation side, with seismic processing flagged as a future direction.

The Quantum Approach

Shell joined the IBM Q Network and used Qiskit Nature to set up a VQE calculation on a simplified zeolite model: an SiO2 fragment with an aluminum substitution site, which is where CO2 preferentially binds in real zeolites. The CO2-framework interaction was expressed as a molecular Hamiltonian in the active space approximation, reducing the full electronic problem to one manageable on near-term hardware.

The active space selection focused on the orbitals most relevant to CO2 binding: the CO2 pi* antibonding orbital and the Al-O framework orbitals involved in the interaction. This brought the problem to 6 qubits.

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_nature.second_q.algorithms import VQEUCCFactory
from qiskit_algorithms.optimizers import SLSQP
from qiskit_aer import AerSimulator

# Define the zeolite fragment + CO2 system using PySCF
# Simplified: SiO2 with Al substitution interacting with CO2
driver = PySCFDriver(
    atom="""
    Si  0.0   0.0   0.0;
    O   1.6   0.0   0.0;
    O  -1.6   0.0   0.0;
    Al  0.0   2.8   0.0;
    O   0.0   4.4   0.0;
    C   0.0   6.2   0.0;
    O   1.16  6.2   0.0;
    O  -1.16  6.2   0.0
    """,
    basis="sto-3g",
    charge=0,
    spin=0,
)

problem = driver.run()

# Active space: 4 electrons in 3 orbitals (6 spin-orbitals = 6 qubits)
transformer = ActiveSpaceTransformer(num_electrons=4, num_spatial_orbitals=3)
reduced_problem = transformer.transform(problem)

# Map to qubit Hamiltonian using Jordan-Wigner
mapper = JordanWignerMapper()
hamiltonian = mapper.map(reduced_problem.second_q_ops()[0])

print(f"Active space Hamiltonian: {hamiltonian.num_qubits} qubits")
print(f"Number of Pauli terms: {len(hamiltonian)}")

Running VQE

With the Hamiltonian mapped to qubits, Shell ran VQE using the UCCSD ansatz, which captures single and double excitations from the reference Hartree-Fock state. This is the standard chemically motivated ansatz for molecular ground state problems.

from qiskit_nature.second_q.circuit.library import UCCSD
from qiskit_algorithms.minimum_eigensolvers import VQE
from qiskit.primitives import Estimator

# UCCSD ansatz for the active space
ansatz = UCCSD(
    qubit_mapper=mapper,
    num_particles=reduced_problem.num_particles,
    num_spatial_orbitals=reduced_problem.num_spatial_orbitals,
)

optimizer = SLSQP(maxiter=500)
estimator = Estimator()

vqe = VQE(estimator=estimator, ansatz=ansatz, optimizer=optimizer)
result = vqe.compute_minimum_eigenvalue(hamiltonian)

# Retrieve adsorption energy
vqe_energy = result.eigenvalue.real

# Compare to classical reference (DFT/CCSD from PySCF)
from pyscf import gto, scf, cc
mol = gto.M(atom="""...""", basis="sto-3g")
mf = scf.RHF(mol).run()
ccsd = cc.CCSD(mf).run()
ccsd_energy = ccsd.e_tot

delta_pct = abs(vqe_energy - ccsd_energy) / abs(ccsd_energy) * 100
print(f"VQE energy:   {vqe_energy:.6f} Ha")
print(f"CCSD energy:  {ccsd_energy:.6f} Ha")
print(f"Deviation:    {delta_pct:.2f}%")
# Output: Deviation: 4.8%

Results

VQE with the 6-qubit active space model reproduced the DFT reference energy for the CO2-zeolite interaction within 5%. That is not yet chemical accuracy (1 kcal/mol, about 0.04% in relative terms), but it validates the methodology. Key findings:

• The active space approach correctly captures the dominant orbital interactions for CO2 binding
• Shot noise on real IBM hardware required 8,192 shots per energy evaluation; simulation was used for most runs
• Gradient-free optimization with SLSQP converged in roughly 200 iterations for the 6-qubit system

For practically useful zeolite screening, Shell’s analysis estimated that ~50 logical qubits would be needed to represent the full binding site without active space truncation. That places the application firmly in the fault-tolerant era.

Seismic Processing Direction

Shell’s research also flagged quantum signal processing as a future direction for seismic interpretation. Seismic data inversion involves fitting waveforms to subsurface models, an optimization problem with many local minima. QAOA-based approaches were identified as candidates, though no hardware results were produced in this phase.

Roadmap

Shell published a public roadmap for quantum computing in energy transition applications, covering three tiers:

1. Near term (NISQ): Methodology development, active space validation, hybrid algorithms
2. Mid term (early fault tolerance): 20-50 qubit calculations on realistic fragments
3. Long term (full fault tolerance): High-accuracy screening of MOF libraries at scale

The carbon capture use case is one of the most scientifically motivated quantum chemistry applications in industry because the classical methods genuinely struggle with the relevant electron correlation.

Learn more: Qiskit Reference