Lithium Hydride VQE with PennyLane: From Molecule to Energy

Why Lithium Hydride?

LiH (lithium hydride) is a standard benchmark in quantum chemistry simulations. It is simple enough to fit on near-term quantum hardware (the STO-3G basis gives 12 spin orbitals and 4 electrons, requiring 12 qubits with Jordan-Wigner encoding) yet complex enough to exhibit genuine electron correlation effects that Hartree-Fock misses. The exact ground state energy is well characterized, so we can compare our VQE result against Full Configuration Interaction (FCI) to quantify the accuracy of the approximation.

This tutorial builds a complete VQE pipeline for LiH using PennyLane’s qchem module.

Background: UCCSD Ansatz

The Unitary Coupled Cluster Singles and Doubles (UCCSD) ansatz is the standard chemistry-inspired ansatz for VQE. It applies a product of single and double excitation operators to the Hartree-Fock reference state:

|psi(theta)> = exp(T1(theta) - T1(theta)^dagger) * exp(T2(theta) - T2(theta)^dagger) * |HF>

where T1 contains single excitations (one electron moved from occupied to virtual orbital) and T2 contains double excitations (two electrons moved). The parameters theta are optimized classically to minimize the energy expectation value.

UCCSD is chemically motivated and usually converges close to FCI for small molecules, but T-count and circuit depth scale as O(N^5) in the number of orbitals, which limits scalability to larger systems.

Setup

pip install pennylane pennylane-qchem scipy

import pennylane as qml
from pennylane import numpy as np
import pennylane.qchem as qchem

Step 1: Build the Molecular Hamiltonian

PennyLane’s qchem.molecular_hamiltonian generates the qubit Hamiltonian for a molecule in a given basis set using OpenFermion and PySCF under the hood:

# LiH at equilibrium bond length (1.57 Angstrom)
symbols  = ["Li", "H"]
geometry = np.array([[0.0, 0.0, 0.0],
                     [0.0, 0.0, 1.57]])   # in Angstroms

hamiltonian, num_qubits = qchem.molecular_hamiltonian(
    symbols,
    geometry,
    basis="sto-3g",
    charge=0,
    mult=1,                # singlet ground state
    active_electrons=2,    # freeze core, use 2 active electrons
    active_orbitals=5,     # 5 active orbitals = 10 spin orbitals = 10 qubits
)

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

Using the full STO-3G active space of LiH requires 12 qubits. Restricting to 2 active electrons in 5 orbitals (the standard frozen-core approximation) reduces this to 10 qubits while retaining the most important correlation effects.

Step 2: Hartree-Fock Initial State and UCCSD Excitations

The Hartree-Fock state fills the lowest-energy spin orbitals. PennyLane provides a helper to compute the HF state and generate all single and double excitations:

electrons = 2
hf_state = qchem.hf_state(electrons, num_qubits)
print(f"HF state: {hf_state}")  # binary occupation vector

singles, doubles = qchem.excitations(electrons, num_qubits)
print(f"Single excitations: {len(singles)}")
print(f"Double excitations: {len(doubles)}")
print(f"Total UCCSD parameters: {len(singles) + len(doubles)}")

Step 3: Define the VQE Circuit

Build a QNode that prepares the UCCSD ansatz and returns the Hamiltonian expectation value:

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

@qml.qnode(dev)
def circuit(params, excitations):
    # Hartree-Fock initial state
    qml.BasisState(hf_state, wires=range(num_qubits))
    
    # Apply all UCCSD excitation operators
    for i, excitation in enumerate(excitations):
        if len(excitation) == 2:
            # Single excitation
            qml.SingleExcitation(params[i], wires=excitation)
        elif len(excitation) == 4:
            # Double excitation
            qml.DoubleExcitation(params[i], wires=excitation)
    
    return qml.expval(hamiltonian)

# Combine singles and doubles into one list
all_excitations = singles + doubles
num_params = len(all_excitations)

# Initialize parameters at zero (HF starting point)
params = np.zeros(num_params, requires_grad=True)

print(f"Initial VQE energy (HF): {circuit(params, all_excitations):.6f} Ha")

Step 4: VQE Optimization with Adam

Use PennyLane’s Adam optimizer, which handles the parameter-shift gradient computation automatically:

opt = qml.AdamOptimizer(stepsize=0.02)

max_iterations = 20
energy_history = []
prev_energy = circuit(params, all_excitations)

print("Starting VQE optimization...")
print(f"{'Iteration':>10} {'Energy (Ha)':>15} {'Delta E':>12}")
print("-" * 40)

for iteration in range(max_iterations):
    params, energy = opt.step_and_cost(
        lambda p: circuit(p, all_excitations),
        params,
    )
    energy_history.append(float(energy))
    delta = energy - prev_energy
    prev_energy = energy
    
    if iteration % 20 == 0:
        print(f"{iteration:>10} {float(energy):>15.8f} {float(delta):>12.2e}")
    
    if abs(delta) < 1e-8 and iteration > 20:
        print(f"\nConverged at iteration {iteration}.")
        break

final_energy = float(energy_history[-1])
print(f"\nVQE ground state energy: {final_energy:.8f} Ha")

Step 5: Compare to Hartree-Fock and FCI

# Hartree-Fock energy (parameters = 0, no correlation)
hf_energy = float(circuit(np.zeros(num_params), all_excitations))

# FCI reference energy for LiH/STO-3G with 2 active electrons, 5 orbitals
# (this value comes from a classical FCI calculation)
fci_energy = -7.882377  # Ha (approximate, varies with active space)

print(f"\nEnergy comparison:")
print(f"  Hartree-Fock:   {hf_energy:.6f} Ha")
print(f"  UCCSD-VQE:      {final_energy:.6f} Ha")
print(f"  FCI (ref):      {fci_energy:.6f} Ha")
print(f"\nCorrelation energy captured by VQE: {hf_energy - final_energy:.6f} Ha")
print(f"Error vs FCI:                        {abs(final_energy - fci_energy):.6f} Ha")

# Chemical accuracy is 1 kcal/mol = 0.0016 Ha
chemical_accuracy = 0.0016
print(f"\nWithin chemical accuracy: {abs(final_energy - fci_energy) < chemical_accuracy}")

Step 6: Potential Energy Surface Scan

A PES scan computes the ground state energy as a function of bond length. This reveals the equilibrium geometry, dissociation energy, and how well VQE handles bond breaking:

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

bond_lengths = np.linspace(1.0, 3.5, 5)  # Angstroms
pes_energies = []

print("\nScanning potential energy surface...")
for r in bond_lengths:
    geo = np.array([[0.0, 0.0, 0.0],
                    [0.0, 0.0, float(r)]])
    
    h, nq = qchem.molecular_hamiltonian(
        symbols, geo,
        basis="sto-3g",
        active_electrons=2,
        active_orbitals=5,
    )
    
    dev_r = qml.device("default.qubit", wires=nq)
    
    @qml.qnode(dev_r)
    def pes_circuit(p):
        qml.BasisState(hf_state, wires=range(nq))
        for i, exc in enumerate(all_excitations):
            if len(exc) == 2:
                qml.SingleExcitation(p[i], wires=exc)
            elif len(exc) == 4:
                qml.DoubleExcitation(p[i], wires=exc)
        return qml.expval(h)
    
    # Quick optimization at each geometry
    p = np.zeros(num_params, requires_grad=True)
    opt_r = qml.AdamOptimizer(stepsize=0.05)
    for _ in range(20):
        p, e = opt_r.step_and_cost(pes_circuit, p)
    
    pes_energies.append(float(e))
    print(f"  r = {r:.2f} A: E = {float(e):.6f} Ha")

# Plot the PES
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(bond_lengths, pes_energies, "o-", label="UCCSD-VQE")
ax.set_xlabel("Li-H Bond Length (Angstrom)")
ax.set_ylabel("Energy (Hartree)")
ax.set_title("LiH Potential Energy Surface (STO-3G, 2e/5o)")
ax.legend()
fig.tight_layout()
fig.savefig("lih_pes.png")
print("\nSaved lih_pes.png")

Step 7: Analyzing Optimal Parameters

The optimal parameters reveal which excitations contribute most to electron correlation:

# Re-run final optimization to get optimal params
params_opt = np.zeros(num_params, requires_grad=True)
opt_final = qml.AdamOptimizer(stepsize=0.02)
for _ in range(max_iterations):
    params_opt, _ = opt_final.step_and_cost(
        lambda p: circuit(p, all_excitations), params_opt
    )

print("\nLargest UCCSD amplitudes:")
sorted_indices = np.argsort(np.abs(params_opt))[::-1]
for rank, idx in enumerate(sorted_indices[:5]):
    exc = all_excitations[idx]
    exc_type = "single" if len(exc) == 2 else "double"
    print(f"  {rank+1}. {exc_type} {exc}: theta = {float(params_opt[idx]):.6f}")

Large amplitudes indicate strongly correlated electron pairs. For LiH, the dominant double excitation typically involves the HOMO-LUMO pair.

Advanced: Adaptive VQE (ADAPT-VQE)

Instead of including all excitations from the start, ADAPT-VQE greedily adds only the excitation with the largest gradient at each step. This can dramatically reduce circuit depth:

def adapt_vqe_step(current_params, current_excitations, candidate_excitations, hamiltonian):
    """Add the excitation with the largest gradient magnitude."""
    best_gradient = 0
    best_excitation = None
    
    for exc in candidate_excitations:
        if exc in current_excitations:
            continue
        
        trial_excitations = current_excitations + [exc]
        trial_params = np.append(current_params, 0.0)
        trial_params = np.array(trial_params, requires_grad=True)
        
        grad_fn = qml.grad(lambda p: circuit(p, trial_excitations))
        grad = grad_fn(trial_params)
        grad_norm = abs(float(grad[-1]))
        
        if grad_norm > best_gradient:
            best_gradient = grad_norm
            best_excitation = exc
    
    return best_excitation, best_gradient

# Run one step of ADAPT-VQE selection
initial_excitations = []
initial_params = np.array([], requires_grad=True)
next_exc, grad_norm = adapt_vqe_step(
    initial_params, initial_excitations, all_excitations, hamiltonian
)
print(f"\nADAPT-VQE selects: {next_exc} (gradient norm = {grad_norm:.6f})")

Key Takeaways

• qchem.molecular_hamiltonian handles the entire pipeline from atomic coordinates to qubit Hamiltonian via Jordan-Wigner mapping.
• The frozen-core approximation reduces qubit count significantly while retaining the key physics.
• UCCSD with PennyLane’s SingleExcitation and DoubleExcitation gates provides hardware-efficient gradient computation via the parameter-shift rule.
• The PES scan is a direct output of the VQE that classical Hartree-Fock can also produce, making it a useful comparison target.
• ADAPT-VQE selects only the most important excitations, making it more suitable for near-term hardware with circuit depth constraints.