Quantum Chemistry with PennyLane | QuantumComputingCourses.com

Quantum chemistry is one of the most credible near-term applications for quantum computers. Simulating the electronic structure of molecules scales exponentially with system size on classical hardware, but polynomially on a quantum computer. The Variational Quantum Eigensolver (VQE) is the leading near-term approach: use a parameterized quantum circuit as an ansatz, optimize it classically to minimize the expectation value of the molecular Hamiltonian.

PennyLane’s qchem module wraps PySCF and OpenFermion to give a high-level Python API that goes from molecular geometry to ground-state energy with minimal boilerplate.

Installation

pip install pennylane pennylane-qchem pyscf openfermion-pyscf scipy

PySCF handles the classical Hartree-Fock reference calculation and the integral generation. OpenFermion handles the fermion-to-qubit mapping (Jordan-Wigner or Bravyi-Kitaev). PennyLane ties it together and provides the VQE optimization loop.

Building the H2 Hamiltonian

The hydrogen molecule H2 at its equilibrium bond length is the standard benchmark for quantum chemistry algorithms. It is small enough to solve exactly classically but captures the essential structure of more complex molecules.

import pennylane as qml
from pennylane import numpy as np

# Define molecular geometry: list of (symbol, [x, y, z]) in Angstroms
symbols = ["H", "H"]
coordinates = np.array([
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 0.74],   # equilibrium bond length ~0.74 Angstrom
])

# Build the molecular Hamiltonian
H, qubits = qml.qchem.molecular_hamiltonian(
    symbols,
    coordinates,
    charge=0,
    mult=1,           # spin multiplicity: 2S+1 = 1 (singlet)
    basis="sto-3g",   # minimal basis set
    mapping="jordan_wigner",
    unit="angstrom",  # coordinates are in Angstrom
)

print(f"Number of qubits: {qubits}")
print(f"Number of Hamiltonian terms: {len(H.operands)}")
print(f"Hamiltonian:\n{H}")

For H2 in the STO-3G basis with Jordan-Wigner mapping, you get 4 qubits and a Hamiltonian with around 15 Pauli terms. The exact FCI ground state energy is approximately -1.136 Hartree.

Getting Electrons and Excitations

The UCCSD ansatz is built from single and double excitations of electrons out of occupied orbitals into virtual orbitals. PennyLane can enumerate these for you:

electrons = 2   # H2 has 2 electrons
orbitals  = 4   # STO-3G gives 2 spatial orbitals = 4 spin-orbitals

# Singles: list of [occupied, virtual] orbital index pairs
# Doubles: list of [occ1, occ2, virt1, virt2] pairs
singles, doubles = qml.qchem.excitations(electrons, orbitals)

print(f"Single excitations: {singles}")
print(f"Double excitations: {doubles}")

# Map excitations to qubit operations (Jordan-Wigner)
s_wires, d_wires = qml.qchem.excitations_to_wires(singles, doubles, wires=range(orbitals))

For H2/STO-3G, there is 1 double excitation and 0 relevant singles (by symmetry). This means the optimal UCCSD circuit has just 1 parameter.

The UCCSD Ansatz

UCCSD (Unitary Coupled Cluster Singles and Doubles) is the quantum chemistry standard for near-exact ground states. The ansatz prepares the Hartree-Fock reference state and applies a parameterized unitary that generates all single and double excitations:

# Reference state: Hartree-Fock for 2 electrons in 4 spin-orbitals
# Occupies the two lowest-energy orbitals
hf_state = qml.qchem.hf_state(electrons, orbitals)
print(f"HF reference state: {hf_state}")  # [1, 1, 0, 0]

dev = qml.device("default.qubit", wires=qubits)

@qml.qnode(dev)
def vqe_circuit(params):
    # UCCSD prepares the HF reference state internally via init_state,
    # then applies the parameterized excitation operators on top.
    qml.UCCSD(params, wires=range(qubits), s_wires=s_wires, d_wires=d_wires, init_state=hf_state)
    # Return the expectation value of the Hamiltonian
    return qml.expval(H)

Running VQE

Initialize the parameters near zero (the HF reference is the zero-excitation limit) and optimize:

# Initialize parameters
num_params = len(singles) + len(doubles)
params = np.zeros(num_params, requires_grad=True)

print(f"Number of variational parameters: {num_params}")
print(f"HF energy (initial): {vqe_circuit(params):.6f} Ha")

# Optimize with Adam
opt = qml.AdamOptimizer(stepsize=0.05)
energies = []

for step in range(30):
    params, energy = opt.step_and_cost(vqe_circuit, params)
    energies.append(float(energy))
    if step % 10 == 0:
        print(f"Step {step:3d}: E = {energy:.6f} Ha")

print(f"\nFinal VQE energy: {energies[-1]:.6f} Ha")

Comparing to the Exact FCI Energy

The full configuration interaction (FCI) energy is the exact ground-state energy within the chosen basis set. For verification, diagonalize the Hamiltonian matrix directly:

import numpy as np

# Get the Hamiltonian as a dense matrix
H_matrix = qml.matrix(H, wire_order=range(qubits))

# Diagonalize to find all eigenvalues
eigenvalues = np.linalg.eigvalsh(H_matrix)
fci_energy = eigenvalues[0]

print(f"FCI ground state energy: {fci_energy:.6f} Ha")
print(f"VQE result:              {energies[-1]:.6f} Ha")
print(f"Difference:              {abs(energies[-1] - fci_energy)*1000:.3f} mHa")

For H2/STO-3G, a well-converged VQE should match FCI to within 0.1 millihartree (chemical accuracy is 1.6 mHa). With a single double-excitation parameter, UCCSD is exact for H2 in this basis.

Active Space Approximation

For larger molecules, the full orbital space generates too many qubits. The active space approximation freezes core electrons and removes high-energy virtual orbitals, retaining only the chemically relevant frontier orbitals.

# LiH example: 4 electrons, 6 spatial orbitals in STO-3G -> 12 qubits
# Active space: 2 active electrons, 2 active orbitals -> 4 qubits

symbols_lih = ["Li", "H"]
coordinates_lih = np.array([
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 1.594],   # LiH equilibrium ~1.594 Angstrom
])

# Full calculation: 12 qubits
H_lih_full, qubits_full = qml.qchem.molecular_hamiltonian(
    symbols_lih,
    coordinates_lih,
    charge=0, mult=1, basis="sto-3g",
    unit="angstrom",
)
print(f"LiH full space: {qubits_full} qubits")

# Active space: select 2 active electrons, 4 active orbitals
core, active = qml.qchem.active_space(
    electrons=4,
    orbitals=6,
    mult=1,
    active_electrons=2,
    active_orbitals=2,
)

H_lih_as, qubits_as = qml.qchem.molecular_hamiltonian(
    symbols_lih,
    coordinates_lih,
    charge=0, mult=1, basis="sto-3g",
    active_electrons=2,
    active_orbitals=2,
    unit="angstrom",
)
print(f"LiH active space: {qubits_as} qubits")

The active space reduces LiH from 12 to 4 qubits, at the cost of some accuracy. The frozen-core approximation (keeping 1s electrons of Li frozen) accounts for most of the error, which is typically small for ground-state energetics.

Full VQE for LiH Active Space

electrons_as = 2
orbitals_as  = 4

hf_lih = qml.qchem.hf_state(electrons_as, orbitals_as)
singles_lih, doubles_lih = qml.qchem.excitations(electrons_as, orbitals_as)
s_wires_lih, d_wires_lih = qml.qchem.excitations_to_wires(
    singles_lih, doubles_lih, wires=range(orbitals_as)
)

dev_lih = qml.device("default.qubit", wires=qubits_as)

@qml.qnode(dev_lih)
def vqe_lih(params):
    qml.UCCSD(params, wires=range(qubits_as),
              s_wires=s_wires_lih, d_wires=d_wires_lih, init_state=hf_lih)
    return qml.expval(H_lih_as)

num_params_lih = len(singles_lih) + len(doubles_lih)
params_lih = qml.numpy.zeros(num_params_lih, requires_grad=True)

opt_lih = qml.AdamOptimizer(stepsize=0.05)
for step in range(20):
    params_lih, energy_lih = opt_lih.step_and_cost(vqe_lih, params_lih)
    if step % 10 == 0:
        print(f"Step {step:3d}: E(LiH) = {energy_lih:.6f} Ha")

print(f"\nFinal LiH VQE energy (active space): {energy_lih:.6f} Ha")

Practical Notes

Gradient computation: PennyLane uses the parameter-shift rule by default for hardware-compatible analytic gradients. On a real QPU, each gradient evaluation requires 2 circuit evaluations per parameter. For UCCSD with many parameters, this becomes expensive.

Ansatz choice: UCCSD is accurate but circuit-depth heavy. Hardware-efficient ansatze (HEA) use shallower circuits with less physical motivation. For NISQ devices, HEA often performs better in practice because it generates less noise.

Convergence: VQE can get stuck in local minima. Strategies include: running multiple random initializations, using gradient-free optimizers (COBYLA, Nelder-Mead) for rough convergence followed by gradient-based refinement, and using natural gradient (QNGO).

What to Try Next

Use @qml.qnode(dev, interface="torch") with a "default.qubit" device to route gradients through PyTorch’s autograd engine, enabling GPU acceleration of the classical gradient computation
Implement a bond dissociation curve by scanning the H-H distance from 0.5 to 3.0 Angstrom and plotting VQE energy vs. bond length
Try the qml.AdaptiveOptimizer for ADAPT-VQE, which builds the ansatz adaptively by adding the most important excitation operators one at a time
Read the PennyLane qchem documentation for the full API reference