Quantum Chemistry Simulation with Q#

Quantum chemistry simulation is one of the most promising near-term applications of quantum computing. Classical methods for computing molecular electronic structure, such as Full Configuration Interaction (FCI), scale exponentially with the number of electrons. Quantum computers can represent electronic wavefunctions naturally using qubits and, in principle, find ground state energies with polynomial resources. Microsoft’s Q# chemistry library provides a complete pipeline from molecular data to quantum circuits.

The Quantum Chemistry Problem

The electronic structure problem asks: given a molecule’s geometry (nuclear positions), what is the ground state energy of its electrons? This energy determines bond lengths, reaction rates, and material properties.

The molecular Hamiltonian in second quantization is:

H = Σ_{pq} h_{pq} a†p a_q + (1/2) Σ{pqrs} h_{pqrs} a†_p a†_q a_r a_s

Here, a†p and a_p are fermionic creation and annihilation operators for molecular orbital p, and the coefficients h{pq} (one-electron integrals) and h_{pqrs} (two-electron integrals) come from classical Hartree-Fock calculations. The challenge is that the Hilbert space grows as 2^N where N is the number of spin-orbitals.

The Jordan-Wigner Transformation

Qubits obey bosonic commutation relations, but electrons are fermions. The Jordan-Wigner transformation maps fermionic operators to qubit operators while preserving the correct antisymmetric statistics:

a†p → (1/2)(X_p - iY_p) ⊗ Z{p-1} ⊗ Z_{p-2} ⊗ … ⊗ Z_0

Each qubit represents one spin-orbital: |1⟩ means occupied, |0⟩ means unoccupied. The string of Z operators enforces the fermionic sign rule (antisymmetry under particle exchange). After this mapping, the molecular Hamiltonian becomes a sum of Pauli tensor products that can be measured on a quantum computer.

Other mappings exist (Bravyi-Kitaev, parity encoding) that can reduce the Z-string length, but Jordan-Wigner is the most intuitive and is the default in the Q# chemistry library.

The Broombridge Format

Microsoft’s chemistry stack uses the Broombridge YAML schema (named after the bridge in Dublin where Hamilton conceived quaternions). A Broombridge file contains all the data needed to construct a molecular Hamiltonian:

• Molecular geometry (atom types and positions)
• Basis set used (e.g., STO-6G, cc-pVDZ)
• One-electron and two-electron integrals from a Hartree-Fock calculation
• Suggested initial states (Hartree-Fock reference, UCCSD trial states)

Classical quantum chemistry packages like NWChem can export Broombridge files. Microsoft provides several precomputed example files for small molecules including H2, LiH, and BeH2.

Here is a simplified example of Broombridge structure for H2:

schema: broombridge-0.2
problem_description:
  - metadata:
      molecule: H2
    basis_set:
      type: sto-6g
    geometry:
      atoms:
        - {name: H, coords: [0.0, 0.0, 0.0]}
        - {name: H, coords: [0.0, 0.0, 1.4]}
    hamiltonian:
      one_electron_integrals:
        units: hartree
        format: sparse
        values: [[1, 1, -1.2528], [2, 2, -0.4760]]
      two_electron_integrals:
        units: hartree
        format: sparse
        values: [[1, 1, 1, 1, 0.6746], ...]
    initial_state_suggestions:
      - state:
          label: "|HF>"
          method: hartree_fock
          superposition:
            - [1.0, "(1a)↑", "(1a)↓"]

Loading Molecular Data in Q#

The Microsoft.Quantum.Chemistry package provides operations for loading Broombridge data and constructing the qubit Hamiltonian.

ChemistrySimulation.qs

namespace ChemistryDemo {
    open Microsoft.Quantum.Intrinsic;

    /// Prepares the Hartree-Fock initial state for a given
    /// number of electrons in the Jordan-Wigner encoding.
    /// In JW encoding, occupied orbitals are |1⟩.
    operation PrepareHartreeFock(
        nElectrons : Int,
        qubits : Qubit[]
    ) : Unit is Adj + Ctl {
        // Occupy the lowest nElectrons spin-orbitals
        for i in 0..nElectrons - 1 {
            X(qubits[i]);
        }
    }
}

Running Chemistry with the Python Host

The typical workflow uses Python to load the Broombridge file, generate the Jordan-Wigner Hamiltonian, and pass the resulting data to Q# operations for quantum simulation.

run_chemistry.py

Note: The following code uses the qsharp.chemistry API from QDK 0.x. The qsharp 1.x package does not include this module. To run this code, install the legacy qsharp==0.28 package or use the Azure Quantum Resource Estimator for chemistry workflows.

import qsharp
from qsharp.chemistry import load_broombridge, load_fermion_hamiltonian
from qsharp.chemistry import load_input_state, encode

# Step 1: Load molecular data from Broombridge file
broombridge_data = load_broombridge("h2_sto6g.yaml")

# Step 2: Extract the fermion Hamiltonian (one- and two-electron integrals)
fermion_hamiltonian = load_fermion_hamiltonian(broombridge_data)

# Step 3: Load the suggested initial state (Hartree-Fock)
input_state = load_input_state(broombridge_data, "|HF>")

# Step 4: Encode into Jordan-Wigner qubit Hamiltonian
# This produces the JWEncodedData structure consumed by Q# operations
jw_encoded = encode(fermion_hamiltonian, input_state)

print(f"Number of qubits required: {jw_encoded.n_qubits}")
print(f"Number of Hamiltonian terms: {jw_encoded.n_terms}")

Estimating Ground State Energy with Phase Estimation

Once the Hamiltonian is encoded, QPE can estimate the ground state energy. The Q# chemistry library provides TrotterStepOracle, which constructs the time-evolution operator e^(-iHt) via Trotterization, suitable for use as the unitary in phase estimation.

EnergyEstimation.qs

namespace ChemistryDemo {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Chemistry.JordanWigner;
    open Microsoft.Quantum.Simulation;
    open Microsoft.Quantum.Characterization;

    /// Estimates the ground state energy of a molecule
    /// using quantum phase estimation with Trotterized
    /// time evolution.
    operation EstimateGroundStateEnergy(
        jwEncoded : JordanWignerEncodingData,
        nBitsPrecision : Int,
        trotterStepSize : Double,
        trotterOrder : Int
    ) : Double {
        // Construct the Trotter step oracle from the JW Hamiltonian.
        // This creates a unitary U = e^(-iHt) via product formula.
        let (nQubits, (rescaleFactor, oracle)) =
            TrotterStepOracle(jwEncoded, trotterStepSize, trotterOrder);

        // Run phase estimation to find the eigenvalue.
        // The oracle encodes the energy as a phase.
        let phase = RobustPhaseEstimation(
            nBitsPrecision,
            oracle,
            PrepareTrialState(jwEncoded, _)
        );

        // Convert phase back to energy (in Hartree)
        let energy = phase * rescaleFactor + IdentityCoefficient(jwEncoded);
        return energy;
    }
}

Python host for energy estimation

Note: The following code uses the qsharp.chemistry API and the from ChemistryDemo import ... / .simulate() patterns from QDK 0.x. Neither qsharp.chemistry nor .simulate() exist in the qsharp 1.x package. To run this code, install the legacy qsharp==0.28 package or use the Azure Quantum Resource Estimator for chemistry workflows.

import qsharp
from qsharp.chemistry import load_broombridge, load_fermion_hamiltonian
from qsharp.chemistry import load_input_state, encode
from ChemistryDemo import EstimateGroundStateEnergy

# Load H2 molecular data
broombridge_data = load_broombridge("h2_sto6g.yaml")
fermion_hamiltonian = load_fermion_hamiltonian(broombridge_data)
input_state = load_input_state(broombridge_data, "|HF>")
jw_encoded = encode(fermion_hamiltonian, input_state)

# Estimate ground state energy
# 10 bits of precision, Trotter step size 0.5, first-order Trotter
energy = EstimateGroundStateEnergy.simulate(
    jwEncoded=jw_encoded,
    nBitsPrecision=10,
    trotterStepSize=0.5,
    trotterOrder=1
)

print(f"Estimated ground state energy: {energy:.6f} Hartree")
print(f"Expected (exact FCI): -1.137260 Hartree")

Expected Output

Estimated ground state energy: -1.137200 Hartree
Expected (exact FCI): -1.137260 Hartree

The small discrepancy comes from the finite Trotter step size and limited precision bits. Increasing either parameter improves accuracy at the cost of deeper circuits.

Understanding the Trotterization

The time evolution operator e^(-iHt) cannot be implemented directly as a quantum gate. Trotterization approximates it by splitting H into a sum of individually simulable terms:

H = H_1 + H_2 + … + H_m

The first-order Trotter formula gives:

e^(-iHt) ≈ (e^(-iH_1 t/n) · e^(-iH_2 t/n) · … · e^(-iH_m t/n))^n

Each term e^(-iH_k t/n) is a product of Pauli rotations, which map directly to quantum gates. Higher-order formulas (second-order Suzuki-Trotter and beyond) achieve better accuracy for the same number of Trotter steps by interleaving the terms in a symmetric pattern.

For H2 in an STO-6G basis, the Jordan-Wigner Hamiltonian requires 4 qubits (4 spin-orbitals) and contains roughly 15 Pauli terms. Each Trotter step applies rotations for all 15 terms. A larger molecule like LiH in the same basis needs 12 qubits and hundreds of terms per step.

Connection to Classical Quantum Chemistry

The quantum simulation pipeline builds on decades of classical quantum chemistry. Classical packages handle the computationally affordable steps:

1. NWChem, PySCF, or Gaussian compute the molecular orbital integrals (the h_{pq} and h_{pqrs} coefficients). This step scales polynomially and is well-suited to classical computers.
2. The integrals are exported in Broombridge format.
3. Q# chemistry library constructs the Jordan-Wigner Hamiltonian and Trotter circuits.
4. Quantum phase estimation finds the ground state energy eigenvalue.

The quantum advantage lies specifically in step 4: solving the eigenvalue problem for strongly correlated systems where classical methods like CCSD(T) break down. Molecules with transition metals, open-shell configurations, or stretched bonds are prime candidates where quantum simulation could provide value beyond classical capabilities.

Variational Quantum Eigensolver (VQE) Alternative

For near-term hardware without error correction, VQE provides a shallower-circuit alternative to QPE. VQE uses a parameterized quantum circuit (ansatz) to prepare trial states and a classical optimizer to minimize the energy expectation value:

E(θ) = ⟨ψ(θ)|H|ψ(θ)⟩

The Q# chemistry library supports UCCSD (Unitary Coupled Cluster Singles and Doubles) ansatz preparation, which is the most common chemistry-motivated ansatz. While VQE requires fewer qubits and shorter circuits, it trades these benefits for many more circuit executions and dependence on the classical optimizer converging to the global minimum. For production-scale quantum chemistry, fault-tolerant QPE remains the target algorithm.