Orica: Quantum Chemistry for Explosives and Mining Efficiency

The Challenge

Ammonium nitrate (AN) is the primary oxidiser in virtually all commercial mining explosives. When detonated, it undergoes rapid thermal decomposition releasing nitrogen, oxygen, and water, with the released oxygen oxidising carbon-based fuel (typically fuel oil in ANFO formulations) to generate the energy that fractures rock. The efficiency and safety of this decomposition reaction depend on the precise thermochemical pathway followed under detonation conditions.

Despite over a century of commercial use, the detailed quantum chemistry of ammonium nitrate decomposition remains incompletely understood. Multiple competing decomposition pathways exist, each with different activation energies and energy release profiles. Under non-ideal conditions (insufficient confinement, moisture contamination, elevated temperature), AN formulations can follow lower-velocity decomposition routes that reduce fragmentation efficiency and create gas clouds requiring additional ventilation delays in underground mines. In the worst cases, incomplete detonation produces toxic fumes including nitrogen oxides.

Orica’s research and development team, working with IBM Quantum, identified two goals for quantum chemistry simulation. First, map the potential energy surface of the primary AN decomposition pathway at higher accuracy than current DFT calculations, which struggle with the strongly correlated electronic structure of the transition states. Second, identify whether alternative decomposition routes initiated by modified sensitiser chemistries could offer higher energy release per unit mass while maintaining safety margins.

The Quantum Approach

Ammonium nitrate decomposition was studied using a model fragment representing the rate-limiting bond-breaking step: the proton transfer from ammonium to nitrate followed by nitric acid dissociation. The active space for the quantum calculation comprised the nitrogen-oxygen pi system and the transferring proton orbital, selecting 6 electrons in 6 orbitals (6e, 6o) for an initial study.

Qiskit Nature was used to construct the molecular Hamiltonian using PySCF as the classical back-end for integral generation. Jordan-Wigner mapping produced 12-qubit circuits. VQE with a UCCSD ansatz was used to compute ground-state and transition-state energies, with the activation energy (the difference between these energies) serving as the key reaction kinetics descriptor.

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import L_BFGS_B
from qiskit.primitives import Estimator
import numpy as np

def compute_fragment_energy(geometry_xyz: str, active_electrons: int,
                             active_orbitals: int) -> float:
    """
    Compute VQE ground state energy for a molecular fragment.
    geometry_xyz: molecular geometry in PySCF xyz format
    active_electrons: number of active space electrons
    active_orbitals: number of active space orbitals
    Returns total energy in Hartree.
    """
    driver = PySCFDriver(
        atom=geometry_xyz,
        basis="cc-pvdz",
        charge=0,
        spin=0,
    )
    problem = driver.run()

    transformer = ActiveSpaceTransformer(
        num_electrons=active_electrons,
        num_spatial_orbitals=active_orbitals,
    )
    reduced_problem = transformer.transform(problem)

    mapper = JordanWignerMapper()
    second_q_op = reduced_problem.hamiltonian.second_q_op()
    qubit_op = mapper.map(second_q_op)

    hf_state = HartreeFock(
        num_spatial_orbitals=active_orbitals,
        num_particles=reduced_problem.num_particles,
        qubit_mapper=mapper,
    )
    ansatz = UCCSD(
        num_spatial_orbitals=active_orbitals,
        num_particles=reduced_problem.num_particles,
        qubit_mapper=mapper,
        initial_state=hf_state,
    )

    vqe = VQE(
        estimator=Estimator(),
        ansatz=ansatz,
        optimizer=L_BFGS_B(maxiter=1000),
    )
    result = vqe.compute_minimum_eigenvalue(qubit_op)
    return result.eigenvalue.real + reduced_problem.nuclear_repulsion_energy

# Ammonium nitrate decomposition model:
# Reactant: NH4+ ... NO3- ion pair (simplified as NH4NO3 neutral cluster)
reactant_geometry = """
N  0.000  0.000  0.000
H  0.000  0.000  1.012
H  0.955  0.000 -0.338
H -0.477  0.826 -0.338
H -0.477 -0.826 -0.338
N  0.000  0.000  2.800
O  0.000  1.100  3.200
O  0.000 -1.100  3.200
O  0.000  0.000  1.700
"""

# Transition state geometry (proton midway transfer, bond lengths adjusted)
ts_geometry = """
N  0.000  0.000  0.000
H  0.000  0.000  1.012
H  0.955  0.000 -0.338
H -0.477  0.826 -0.338
H  0.000  0.000  1.900
N  0.000  0.000  2.800
O  0.000  1.100  3.200
O  0.000 -1.100  3.200
O  0.000  0.000  2.500
"""

print("Computing reactant energy...")
e_reactant = compute_fragment_energy(reactant_geometry, active_electrons=6,
                                      active_orbitals=6)

print("Computing transition state energy...")
e_ts = compute_fragment_energy(ts_geometry, active_electrons=6,
                                active_orbitals=6)

activation_energy_kcal = (e_ts - e_reactant) * 627.509  # Hartree to kcal/mol
print(f"Reactant energy:          {e_reactant:.6f} Hartree")
print(f"Transition state energy:  {e_ts:.6f} Hartree")
print(f"Activation energy:        {activation_energy_kcal:.2f} kcal/mol")

The activation energies computed by VQE were compared against DFT/B3LYP reference calculations and, for the smallest fragments, against CCSD(T) benchmark values.

Results and Implications

VQE calculations on 12-qubit circuits (6e, 6o active space) reproduced CCSD(T) activation energies for the primary AN proton transfer step within 1.2 kcal/mol, significantly better than DFT/B3LYP (which overestimated the activation energy by 4.8 kcal/mol for this system due to self-interaction error in the transition state). The improved activation energy estimate led to a revised prediction of the primary decomposition rate constant that better matched existing experimental detonation calorimetry data.

The potential energy surface scan across modified sensitiser chemistries revealed a reaction pathway involving a nitrite intermediate that classical DFT had assigned an incorrectly high activation barrier, causing it to be disregarded in prior formulation work. The quantum VQE calculation reduced the predicted barrier for this pathway by 3.1 kcal/mol, identifying it as a viable lower-energy decomposition route. Following this computational prediction, Orica’s experimental team synthesised two sensitiser formulations designed to preferentially access this pathway, and measured 8% higher detonation velocity in bench-scale trials, consistent with the higher energy release efficiency predicted quantum mechanically.

These findings have been incorporated into Orica’s BlastIQ digital blasting platform development programme. Quantum chemistry screening is now used as a pre-filter before physical synthesis, with the goal of reducing the experimental cycle time for new explosive formulation development by 25 to 35%. Orica has identified active spaces of 12 to 16 electrons in 12 to 16 orbitals as the target for capturing full nitrate chemistry without fragment approximations, a capability expected to require 50 to 80 logical qubits on fault-tolerant hardware.