VQE for H2 Molecule Ground State

What You'll Learn:

• How the Variational Quantum Eigensolver (VQE) combines quantum circuits with classical optimization
• How a molecular Hamiltonian is mapped to qubit operators via Jordan-Wigner transformation
• How to measure energy expectation values using Pauli basis rotations
• What "chemical accuracy" means and why it matters for quantum chemistry

Level: Intermediate | Time: 25 minutes | Qubits: 2 | Framework: Qiskit

Prerequisites

• Bell State — entanglement with CX gates
• Hadamard Test — expectation value estimation

The Idea

Imagine you want to know the lowest energy state of a hydrogen molecule — the energy at which its two electrons settle into the most stable configuration. Classically, solving this exactly requires diagonalizing a matrix that grows exponentially with system size. For H2 it's still tractable (4x4 matrix), but for larger molecules like caffeine (24 atoms, 160 electrons) the matrix is astronomically large.

VQE takes a different approach: instead of solving the full eigenvalue problem, it guesses a quantum state, measures its energy, and adjusts the guess to make the energy lower. The quantum computer prepares the trial state (something a classical computer struggles to represent efficiently), while a classical optimizer handles the parameter updates. This hybrid loop converges to the ground state energy — the lowest eigenvalue of the molecular Hamiltonian.

The variational principle guarantees this works: any trial state gives an energy at least as high as the true ground state. So we can only go down, never accidentally overshoot.

How It Works

The VQE Loop

CODECopy
Classical Computer                    Quantum Computer
┌─────────────────┐                  ┌──────────────────┐
│                  │  parameters θ   │                   │
│  COBYLA          │ ──────────────> │  Prepare |ψ(θ)>  │
│  Optimizer       │                 │  Measure in Z,X,Y │
│                  │  energy E(θ)    │  bases             │
│  Update θ to     │ <────────────── │                   │
│  minimize E(θ)   │                 └──────────────────┘
└─────────────────┘
        │
        ▼
  Converged? → E(θ*) ≈ E₀

Step 1: Build the Ansatz

The ansatz is a parameterized quantum circuit that prepares trial wavefunctions. For H2, we use a hardware-efficient ansatz with 4 parameters:

CODECopy
     ┌──────────┐          ┌──────────┐
q_0: ┤ RY(θ₀)  ├────■─────┤ RY(θ₂)  ├
     ├──────────┤  ┌─┴─┐   ├──────────┤
q_1: ┤ RY(θ₁)  ├──┤ X ├───┤ RY(θ₃)  ├
     └──────────┘  └───┘   └──────────┘

• RY gates rotate each qubit's state on the Bloch sphere
• CX gate creates entanglement — essential for capturing electron correlation
• Without the CX, we could only represent product states (independent electrons), missing the correlation energy entirely

PYTHONCopy
from circuit import create_h2_ansatz

# Create ansatz with specific angles
ansatz = create_h2_ansatz([0.0, 3.14159, 0.0, 0.0])
print(ansatz.draw())

Step 2: Map the Molecule to Qubits

The H2 Hamiltonian in the STO-3G basis has 6 terms. Each term is a Pauli string — a tensor product of Pauli operators (I, X, Y, Z) acting on the 2 qubits:

These coefficients come from the parity mapping with two-qubit reduction of the 4-qubit Jordan-Wigner Hamiltonian. Nuclear repulsion (V_nn = 0.7199 Ha) is added separately to get the total molecular energy.

Step 3: Measure in Pauli Bases

To compute the energy, we need the expectation value of each Pauli term. The trick: different terms require different measurement bases.

ZZ basis — just measure in the computational basis. One circuit gives us three expectation values (⟨Z₀⟩, ⟨Z₁⟩, ⟨Z₀Z₁⟩) simultaneously:

PYTHONCopy
# Computational basis: no rotation needed
qc.measure_all()

XX basis — apply Hadamard (H) before measurement to rotate from X eigenbasis to Z eigenbasis:

PYTHONCopy
qc.h(0)
qc.h(1)
qc.measure_all()

Only 2 circuits are needed per energy evaluation — the parity mapping with two-qubit reduction eliminates the Y₀Y₁ term from the Hamiltonian.

Step 4: Compute the Energy

From measurement counts, compute the expectation value of each Pauli string using the parity rule:

For a 2-qubit Pauli string, bitstring b has eigenvalue (-1)^parity(b), where parity = number of 1-bits. So:

⟨P⟩ = Σ_b (-1)^parity(b) × count(b) / shots

The total molecular energy combines all terms plus nuclear repulsion:

CODECopy
E_elec(θ) = c₀ + c₁⟨Z₀⟩ + c₂⟨Z₁⟩ + c₃⟨Z₀Z₁⟩ + c₄⟨X₀X₁⟩
E_total(θ) = E_elec(θ) + V_nn

Step 5: Optimize

The classical optimizer (COBYLA) adjusts θ to minimize E(θ). After ~50-100 iterations, the energy converges to the ground state:

PYTHONCopy
from circuit import optimize_vqe

result = optimize_vqe(max_iterations=100, shots=2048, seed=42)
print(f"Optimized energy: {result['optimal_energy']:.4f} Ha")
print(f"Exact:            {result['exact_ground_state']:.4f} Ha")
print(f"Error:            {result['error']:.4f} Ha")
print(f"Chemical accuracy: {result['chemical_accuracy']}")

The Math

Variational Principle

For any normalized trial state |ψ(θ)⟩:

⟨ψ(θ)|H|ψ(θ)⟩ ≥ E₀

where E₀ is the true ground state energy. Equality holds when |ψ(θ)⟩ = |ψ₀⟩ (the ground state). This guarantees VQE can only overestimate the energy, so minimizing always moves toward the true answer.

H2 Hamiltonian Construction

Starting from the second-quantized electronic Hamiltonian:

H = Σ_{pq} h_{pq} a†_p a_q + ½ Σ_{pqrs} h_{pqrs} a†_p a†_q a_s a_r

where a†_p, a_p are fermionic creation/annihilation operators.

The parity mapping converts these to qubit operators. For H2 in STO-3G (4 spin-orbitals → 4 qubits), two qubits encode conserved symmetries (particle number parity) and can be removed, yielding a 5-term Hamiltonian on 2 qubits. Nuclear repulsion V_nn = 1/R is added classically.

Ansatz State

With default parameters θ = [π, π, 0, 0] (Hartree-Fock):

CODECopy
|ψ⟩ = RY(π)⊗RY(π) → CX → |00⟩
     = CX |1⟩|1⟩           (since RY(π)|0⟩ = |1⟩)
     = |1⟩|0⟩              (CX flips target when control is |1⟩)

This gives the Hartree-Fock state |10⟩ (in Qiskit notation: bitstring "01"), where one electron occupies each spin-orbital. The energy is -1.117 Ha — missing ~20 mHa of correlation energy. VQE optimization finds the superposition (α|01⟩ + β|10⟩ in the Hamiltonian eigenbasis) that captures this correlation.

Expected Output

With 2048 shots and COBYLA optimization, VQE typically reaches within 5-15 mHa of the exact value — close to but not always within chemical accuracy, due to shot noise.

Running the Circuit

PYTHONCopy
from circuit import run_circuit, optimize_vqe, verify_vqe

# Quick evaluation with near-optimal parameters
result = run_circuit()
print(f"Energy: {result['energy']:.4f} Ha (error: {result['error']:.4f})")

# Full optimization
opt = optimize_vqe(max_iterations=100, shots=2048, seed=42)
print(f"Optimized: {opt['optimal_energy']:.4f} Ha")
print(f"Chemical accuracy: {opt['chemical_accuracy']}")

# Verification suite
v = verify_vqe()
for check in v["checks"]:
    status = "PASS" if check["passed"] else "FAIL"
    print(f"[{status}] {check['name']}: {check['detail']}")

Try It Yourself

1. Change the bond length: The Hamiltonian coefficients above are for equilibrium (0.735 A). Look up H2 coefficients at 1.0 A and 2.0 A — what happens to the ground state energy? (Hint: it approaches -1.0 Ha as the bond breaks.)

2. Remove the entangling gate: Delete the CX from the ansatz and re-optimize. You'll get the Hartree-Fock energy (-1.1168 Ha) but never reach the exact ground state — the missing ~20 mHa is the correlation energy.

3. Add more layers: Duplicate the RY-CX-RY block to create a 2-layer ansatz (8 parameters). Does the extra expressibility improve convergence speed?

4. Try different optimizers: Replace COBYLA with method="Nelder-Mead" or method="Powell" in optimize_vqe(). Which converges faster? Which is more stable with finite shots?

5. Increase shots: Run with 256, 1024, 4096, and 16384 shots. Plot the energy variance — you should see it decrease as 1/sqrt(shots).

What's Next

• LiH Ground State — Scale to 4 qubits with a 2-layer ansatz for lithium hydride
• HeH+ Ground State — The first molecule in the universe, with an asymmetric Hamiltonian
• QAOA MaxCut — Another variational algorithm, applied to combinatorial optimization instead of chemistry

Applications

References

• Peruzzo, A. et al. (2014). "A variational eigenvalue solver on a photonic quantum processor." Nature Communications 5, 4213. DOI: 10.1038/ncomms5213
• O'Malley, P.J.J. et al. (2016). "Scalable quantum simulation of molecular energies on a superconducting qubit processor." Physical Review X 6, 031007. DOI: 10.1103/PhysRevX.6.031007
• Kandala, A. et al. (2017). "Hardware-efficient variational quantum eigensolver for small molecules." Nature 549, 242-246. DOI: 10.1038/nature23879
• Nielsen, M.A. & Chuang, I.L. Quantum Computation and Quantum Information, Chapter 4.