Variational Quantum Eigensolver with PennyLane

The Variational Quantum Eigensolver (VQE) is one of the most important algorithms for near-term quantum hardware. It estimates the ground state energy of a molecule by combining a quantum circuit with a classical optimizer, keeping circuit depth shallow enough to run on noisy devices. In this tutorial you will implement VQE from scratch in PennyLane, compute the ground state energy of molecular hydrogen (H2), and build intuition for every component along the way.

What VQE Actually Computes

The foundation of VQE is the variational principle from quantum mechanics. It states a clean, powerful fact: for any normalized quantum state |psi(theta)>, the expectation value of the Hamiltonian H satisfies

<psi(theta)|H|psi(theta)>  >=  E_0

where E_0 is the true ground state energy. In other words, no matter what state you prepare, the measured energy can never dip below the real answer. It can only be equal to or higher than E_0.

VQE exploits this by parameterizing a quantum circuit with angles theta, measuring the energy expectation value, and then using a classical optimizer to adjust theta downward. Each iteration lowers the energy estimate, and the variational principle guarantees that the result is always a valid upper bound on E_0. When the optimizer converges, you have the best approximation to E_0 that your circuit family (the “ansatz”) can represent.

The quality of the final result depends on two factors:

1. Ansatz expressibility: Can the circuit family represent the true ground state? If the true ground state |psi_0> lies within the set of states reachable by the ansatz, VQE can in principle find the exact answer. If not, you get the closest state the ansatz can produce, which is still a valid upper bound.

2. Optimizer convergence: Does the classical optimizer find the global minimum of the energy landscape, or does it get stuck in a local minimum? The energy landscape for VQE can have local minima, saddle points, and flat regions (barren plateaus), all of which can trap gradient-based optimizers.

For H2 with the UCCSD ansatz (which we use in this tutorial), the true ground state lies within the ansatz family. This means VQE can reach the exact Full Configuration Interaction (FCI) energy, limited only by the optimizer’s ability to converge. H2 is therefore a perfect test case: you know the right answer, so you can verify that your implementation works correctly before moving to harder molecules.

Installation

pip install pennylane pennylane-qchem

PennyLane’s qchem module handles the chemistry-specific pieces: building molecular Hamiltonians from atomic coordinates, computing Hartree-Fock states, and generating excitation operators. The module uses PySCF under the hood for the electronic structure integrals.

Building the Hydrogen Hamiltonian

The hydrogen molecule (H2) is the simplest molecule with a non-trivial electronic structure. At a bond distance of approximately 0.74 angstroms (1.3984 bohr), the Hamiltonian can be expressed as a sum of Pauli operators acting on 4 qubits (representing 4 spin orbitals).

import pennylane as qml
from pennylane import numpy as np

# Define the H2 molecule geometry (atomic units: bohr)
symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.3984])  # ~0.74 angstrom

# Build the molecular Hamiltonian
H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)

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

The molecular_hamiltonian function performs several steps internally. It computes the one-electron and two-electron integrals in the STO-3G basis set, constructs the second-quantized fermionic Hamiltonian, and then applies the Jordan-Wigner transformation to map fermionic creation and annihilation operators onto Pauli operators that act on qubits.

Understanding the Hamiltonian Terms

The resulting qubit Hamiltonian is a weighted sum of Pauli strings. You can inspect the individual terms to understand what each one represents physically:

# Print the full Hamiltonian
print(H)

For H2 in the STO-3G basis with Jordan-Wigner mapping, the Hamiltonian contains several types of terms:

• Constant (Identity) term: This captures the nuclear repulsion energy between the two hydrogen nuclei plus core electronic contributions. It shifts the entire energy spectrum by a fixed amount.

• ZI and IZ terms (single-qubit Z operators): These are single-body terms that represent the energy of individual spin-orbitals. A Z operator on qubit i measures whether spin-orbital i is occupied (eigenvalue -1) or empty (eigenvalue +1), so these terms assign an energy cost or benefit to occupying each orbital.

• ZZ terms (two-qubit ZZ operators): These represent two-body Coulomb repulsion between electrons. When two spin-orbitals are both occupied, the ZZ term contributes a repulsive energy.

• XX and YY terms (paired X and Y operators): These are the most physically interesting terms. They represent electron hopping and correlation, allowing the quantum state to develop superpositions between different electronic configurations. The correlation energy of H2 comes entirely from these terms. Without them, VQE would just return the Hartree-Fock energy.

You can extract the coefficients and operators to see each term explicitly:

# Inspect individual terms and their coefficients
coeffs, ops = H.terms()
for coeff, op in zip(coeffs, ops):
    print(f"  {coeff:+.6f} * {op}")

The relative magnitudes of these coefficients reveal the physics: the Z-type terms dominate (reflecting the large single-particle energies), while the XX+YY correlation terms are smaller but essential for capturing the electron-electron interaction that Hartree-Fock misses.

The Ansatz: UCCSD

The ansatz is the parameterized circuit that prepares the trial state. For chemistry problems, the Unitary Coupled Cluster Singles and Doubles (UCCSD) ansatz is a strong choice because it is chemically motivated and captures the dominant electron correlation effects.

What Singles and Doubles Mean Physically

UCCSD builds on the Hartree-Fock reference state (the simplest mean-field electron solution) and adds parameterized excitation operators. These operators come in two types:

• Singles excitations: One electron is promoted from an occupied spin-orbital to a virtual (unoccupied) one. For instance, moving an electron from spin-orbital 0 to spin-orbital 2. Singles capture orbital relaxation effects, where the optimal orbitals in the correlated state differ slightly from the Hartree-Fock orbitals.

• Doubles excitations: Two electrons are simultaneously promoted from occupied to virtual spin-orbitals. For instance, moving electrons from spin-orbitals (0, 1) to spin-orbitals (2, 3). Doubles capture the bulk of electron correlation energy, particularly the tendency of electrons to avoid each other (dynamic correlation).

For H2 with 2 electrons and 4 spin-orbitals (indices 0, 1, 2, 3, where 0 and 1 are occupied in the Hartree-Fock state), the possible excitations are constrained by spin selection rules. You can inspect exactly what PennyLane generates:

# Get the active electrons and orbitals
electrons = 2
orbitals = qubits

# Generate singles and doubles excitations
singles, doubles = qml.qchem.excitations(electrons, orbitals)

print(f"Singles excitations: {singles}")
print(f"Doubles excitations: {doubles}")

For H2, the singles list contains excitations like [0, 2] and [1, 3], meaning “move an electron from orbital 0 to orbital 2” and “move an electron from orbital 1 to orbital 3.” These respect spin conservation: orbital 0 and 2 have the same spin (alpha), while orbitals 1 and 3 share the other spin (beta).

The doubles list contains the excitation [0, 1, 2, 3], which simultaneously promotes both electrons from the occupied orbitals (0, 1) to the virtual orbitals (2, 3). This double excitation is the most important one for H2, because it creates the entangled superposition that captures the bulk of the correlation energy.

Mapping Excitations to Circuit Wires

Each excitation gets mapped to a set of qubit wires that define which gates act on which qubits:

# Map excitations to the qubit register
s_wires, d_wires = qml.qchem.excitations_to_wires(singles, doubles, wires=range(qubits))

print(f"Single excitation wires: {s_wires}")
print(f"Double excitation wires: {d_wires}")

The Hartree-Fock State as a Qubit State

The UCCSD circuit starts from the Hartree-Fock reference state, which is the lowest-energy single-determinant state (no electron correlation). For H2 with 2 electrons, the Hartree-Fock state occupies the two lowest spin-orbitals:

# Hartree-Fock initial state: occupy the lowest two spin orbitals
hf_state = qml.qchem.hf_state(electrons, orbitals)

print(f"HF state: {hf_state}")  # [1, 1, 0, 0]
# This corresponds to the computational basis state |1100>
# Qubit i = 1 means spin-orbital i is occupied
# Qubit i = 0 means spin-orbital i is empty

In the circuit, qml.BasisState(hf_state, wires=wires) prepares the state |1100>, meaning spin-orbitals 0 and 1 are occupied and spin-orbitals 2 and 3 are empty.

When the UCCSD circuit acts on this state, it creates a superposition. For H2, the dominant effect is the double excitation, which mixes |1100> with |0011>:

|psi(theta)> = cos(theta/2)|1100> + sin(theta/2)|0011>

At the optimal value of theta, this superposition gives the exact ground state of H2. The parameter theta controls how much weight the correlated configuration |0011> receives relative to the Hartree-Fock configuration |1100>.

# Total number of variational parameters
n_params = len(singles) + len(doubles)
print(f"Variational parameters: {n_params}")

The Cost Function

The cost function is the expectation value of the molecular Hamiltonian with respect to the circuit output. Minimizing it drives the circuit toward the ground state.

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

@qml.qnode(dev)
def circuit(params, wires, s_wires, d_wires, hf_state):
    # Prepare the Hartree-Fock reference state
    qml.BasisState(hf_state, wires=wires)

    # Apply UCCSD ansatz
    qml.UCCSD(params, wires=wires, s_wires=s_wires, d_wires=d_wires, init_state=hf_state)

    # Measure the expectation value of H
    return qml.expval(H)

Understanding the Gradient Computation

PennyLane computes exact gradients of this circuit using the parameter-shift rule. For a parameter theta_i, the gradient is:

d/d(theta_i) <H> = [ f(theta_i + pi/2) - f(theta_i - pi/2) ] / 2

where f evaluates the full circuit expectation value. This formula is exact (not a finite-difference approximation) and works on real quantum hardware, not just simulators.

The cost of computing gradients scales linearly with the number of parameters. For each parameter, you need two circuit evaluations (one shifted up, one shifted down). You can calculate the total cost:

print(f"Parameters: {n_params}")
print(f"Circuit evaluations per gradient step: {2 * n_params}")
# On hardware with 1000 shots per evaluation:
shots_per_step = 1000 * 2 * n_params
print(f"Shots per gradient step: {shots_per_step}")

For H2 with 3 UCCSD parameters, each gradient step requires 6 circuit evaluations. On a quantum device running 1000 measurement shots per evaluation, that translates to 6000 shots per optimization step. For 100 optimization steps, you need 600,000 total shots. This shot budget is a key consideration when planning hardware experiments.

Running the Optimization

We wrap the circuit call in a cost function and use the Adam optimizer to minimize it.

Why Adam Works Well for VQE

Plain gradient descent updates each parameter by a fixed step proportional to the gradient:

theta_i = theta_i - lr * grad_i

This approach can oscillate in narrow valleys, overshoot past minima, and move slowly through flat regions. The Adam optimizer (Adaptive Moment Estimation) addresses all of these issues by maintaining running estimates of both the gradient mean and the gradient variance for each parameter. It adapts the effective step size per-parameter based on the gradient history.

Key hyperparameters for Adam:

• stepsize (learning rate): The base step size, typically 0.01 to 0.1 for VQE.
• beta1 (default 0.9): Controls momentum on the gradient mean. Higher values give more smoothing.
• beta2 (default 0.99): Controls momentum on the squared gradient. Higher values give more stable step sizes.

The Optimization Loop

import matplotlib.pyplot as plt

# Define the cost function
def cost_fn(params):
    return circuit(params, wires=range(qubits), s_wires=s_wires,
                   d_wires=d_wires, hf_state=hf_state)

# Initialize parameters to zero (starts at Hartree-Fock state)
params = np.zeros(n_params, requires_grad=True)

# Adam optimizer handles the saddle-point landscape better than plain gradient descent
opt = qml.AdamOptimizer(stepsize=0.02)

energies = []
max_iterations = 100
conv_threshold = 1e-6  # Convergence threshold in Hartree

for step in range(max_iterations):
    params, energy = opt.step_and_cost(cost_fn, params)
    energies.append(float(energy))

    # Check convergence
    if step > 0 and abs(energies[-1] - energies[-2]) < conv_threshold:
        print(f"Converged at step {step}: E = {energy:.6f} Ha")
        break

    if step % 20 == 0:
        print(f"Step {step:3d}: E = {energy:.6f} Ha")

# Step   0: E = -1.117349 Ha
# Step  20: E = -1.136101 Ha
# Step  40: E = -1.136779 Ha
# Step  60: E = -1.136890 Ha
# Step  80: E = -1.136923 Ha

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

Effect of Learning Rate on Convergence

The choice of learning rate significantly affects how quickly VQE converges and whether it reaches the global minimum. You can experiment with different values:

for lr in [0.01, 0.05, 0.1, 0.3]:
    opt_test = qml.AdamOptimizer(stepsize=lr)
    params_test = np.zeros(n_params, requires_grad=True)
    for _ in range(50):
        params_test, _ = opt_test.step_and_cost(cost_fn, params_test)
    print(f"lr={lr}: final energy = {cost_fn(params_test):.6f} Ha")

Too small a learning rate (0.01) converges slowly and may not reach chemical accuracy within a reasonable number of steps. Too large a learning rate (0.3) can overshoot the minimum and oscillate. For H2 with UCCSD, values between 0.02 and 0.1 typically work well.

Plotting and Interpreting Energy Convergence

plt.figure(figsize=(8, 4))
plt.plot(energies, color="steelblue", linewidth=2, label="VQE energy")
plt.axhline(y=-1.136190, color="crimson", linestyle="--", label="Exact FCI")
plt.axhline(y=-1.117349, color="gray", linestyle=":", alpha=0.7, label="Hartree-Fock")
plt.xlabel("Optimization step")
plt.ylabel("Energy (Hartree)")
plt.title("VQE Energy Convergence for H2")
plt.legend()
plt.tight_layout()
plt.savefig("vqe_convergence.png", dpi=150)
plt.show()

What to Look For in the Convergence Plot

The plot reveals several important features:

1. Starting point: The energy begins at the Hartree-Fock value (approximately -1.117 Ha for H2 at equilibrium). This happens because the initial parameters are all zero, and UCCSD with zero parameters produces exactly the Hartree-Fock state. The excitation operators reduce to identity when their angles are zero.

2. Rapid initial descent: During the first 10-20 steps, the energy drops quickly. The optimizer finds the dominant direction in parameter space that lowers the energy, which for H2 is primarily the double excitation parameter.

3. Slow final convergence: After the initial rapid decrease, the curve flattens. The optimizer is now in a shallow region near the minimum, making small adjustments. This is normal behavior for gradient-based optimization near a minimum.

4. The correlation energy gap: The difference between Hartree-Fock and FCI energies is the correlation energy. For H2:

e_hf = -1.117349    # Hartree-Fock energy
e_fci = -1.136190   # Exact FCI energy
e_corr = e_fci - e_hf
print(f"Correlation energy: {e_corr:.6f} Ha")
print(f"Correlation energy: {e_corr * 627.509:.2f} kcal/mol")
# Correlation energy: -0.018841 Ha
# Correlation energy: -11.82 kcal/mol

This 12 kcal/mol gap is significant in chemistry. Bond energies and reaction barriers often differ by only a few kcal/mol, so capturing correlation energy accurately is essential for useful predictions.

5. Convergence criterion: Convergence is typically defined as |E(step) - E(step-1)| < 1e-6 Ha (sub-millihartree accuracy). At that level, the result is converged to well within “chemical accuracy” (1 kcal/mol = 0.0016 Ha).

Interpreting the Result

Energies are reported in Hartree (Ha). One Hartree is about 627.5 kcal/mol, so even small differences matter in chemistry. The ground state energy of H2 at equilibrium is approximately -1.1363 Ha. Values above this indicate the circuit has not fully converged or the ansatz cannot represent the true ground state.

For H2, UCCSD is essentially exact: the ansatz can represent the true ground state, and the optimizer converges to it. For larger molecules, more expressive ansatze or active space approximations become necessary to manage the qubit count.

Potential Energy Surface Scan

A single VQE calculation gives the energy at one molecular geometry. To understand the full bonding behavior of H2, you can scan across multiple bond lengths and plot the potential energy surface (PES). This reveals the equilibrium bond length (the energy minimum) and dissociation behavior (what happens as the atoms move apart).

bond_lengths_bohr = np.linspace(0.5, 5.0, 20)  # in bohr
pes_energies = []

for r in bond_lengths_bohr:
    # Place two hydrogen atoms along the z-axis separated by distance r
    coords = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r])
    H_r, qubits_r = qml.qchem.molecular_hamiltonian(symbols, coords)

    # Set up a fresh device and circuit for this geometry
    dev_r = qml.device("default.qubit", wires=qubits_r)

    @qml.qnode(dev_r)
    def circuit_r(params):
        qml.BasisState(hf_state, wires=range(qubits_r))
        qml.UCCSD(params, wires=range(qubits_r), s_wires=s_wires,
                  d_wires=d_wires, init_state=hf_state)
        return qml.expval(H_r)

    # Optimize for this bond length
    opt_r = qml.AdamOptimizer(stepsize=0.05)
    params_r = np.zeros(n_params, requires_grad=True)

    for _ in range(60):
        params_r, energy_r = opt_r.step_and_cost(circuit_r, params_r)

    pes_energies.append(float(energy_r))
    print(f"r = {r:.2f} bohr ({r * 0.529177:.3f} A): E = {energy_r:.6f} Ha")

# Find equilibrium bond length
min_idx = np.argmin(pes_energies)
r_eq_bohr = bond_lengths_bohr[min_idx]
r_eq_angstrom = r_eq_bohr * 0.529177
print(f"\nEquilibrium bond length: {r_eq_bohr:.2f} bohr = {r_eq_angstrom:.3f} angstrom")
print(f"Energy at equilibrium: {pes_energies[min_idx]:.6f} Ha")

The conversion factor is 1 bohr = 0.529177 angstrom. For H2, the PES minimum falls near 1.4 bohr (0.74 angstrom), matching the experimentally known bond length. At short distances, strong nuclear repulsion drives the energy up. At large distances, the energy plateaus at the dissociation limit (two isolated hydrogen atoms).

# Plot the potential energy surface
plt.figure(figsize=(8, 4))
plt.plot(bond_lengths_bohr * 0.529177, pes_energies, "o-", color="steelblue",
         linewidth=2, markersize=5)
plt.xlabel("Bond length (angstrom)")
plt.ylabel("Energy (Hartree)")
plt.title("H2 Potential Energy Surface from VQE")
plt.axvline(x=r_eq_angstrom, color="gray", linestyle=":", alpha=0.5, label="Equilibrium")
plt.legend()
plt.tight_layout()
plt.savefig("h2_pes.png", dpi=150)
plt.show()

Note that molecular_hamiltonian expects coordinates in bohr, not angstrom. Using angstrom without converting (multiplying by 1.8897) gives a Hamiltonian for the wrong geometry and produces incorrect energies. This is one of the most common mistakes in quantum chemistry calculations.

Hardware-Efficient Ansatz Comparison

UCCSD is motivated by chemistry and guarantees that physically relevant excitations are included. However, it produces relatively deep circuits with many CNOT gates, which makes it challenging to run on noisy hardware. An alternative is a hardware-efficient ansatz (HEA), which uses shallow, regular layers of single-qubit rotations and entangling gates.

PennyLane provides qml.StronglyEntanglingLayers as a built-in hardware-efficient ansatz:

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

@qml.qnode(dev_hea)
def hea_circuit(params):
    # No HF state preparation -- HEA explores the full Hilbert space
    qml.StronglyEntanglingLayers(params, wires=range(qubits))
    return qml.expval(H)

# Shape: (n_layers, n_wires, 3 rotations per wire)
n_layers = 2
shape = qml.StronglyEntanglingLayers.shape(n_layers=n_layers, n_wires=qubits)
print(f"HEA parameter shape: {shape}")
print(f"Total HEA parameters: {np.prod(shape)}")

# Initialize with small random values
np.random.seed(42)
params_hea = np.random.uniform(-0.1, 0.1, shape, requires_grad=True)

opt_hea = qml.AdamOptimizer(stepsize=0.05)
hea_energies = []

for step in range(200):
    params_hea, energy_hea = opt_hea.step_and_cost(hea_circuit, params_hea)
    hea_energies.append(float(energy_hea))

print(f"HEA final energy:  {hea_energies[-1]:.6f} Ha")
print(f"UCCSD final energy: {energies[-1]:.6f} Ha")
print(f"Exact FCI energy:  -1.136190 Ha")

The trade-offs between the two approaches are clear:

• UCCSD uses fewer parameters (3 for H2) but deeper circuits. It converges reliably because the ansatz is physically motivated and starts from the Hartree-Fock state.
• HEA uses more parameters (24 for 2 layers on 4 qubits) but shallower, more regular circuits. It may require more optimization steps and can get stuck in local minima because it lacks chemistry-specific structure.

For H2 (a small, easy molecule), both approaches typically reach near-exact energies. The differences become more pronounced for larger molecules, where UCCSD’s chemical structure gives it a significant advantage in convergence, while HEA’s shallow depth gives it an advantage on noisy hardware.

Shot Noise and Finite-Shot Evaluation

In simulation mode, default.qubit computes exact expectation values analytically. On real quantum hardware (or when simulating finite-shot measurements), the expectation value is estimated from a finite number of measurement samples, introducing statistical uncertainty.

You can observe this effect by comparing exact simulation with finite-shot evaluation:

optimal_params = params  # Use the converged UCCSD parameters

for shots in [None, 10000, 1000, 100]:
    dev_shots = qml.device("default.qubit", wires=qubits, shots=shots)

    @qml.qnode(dev_shots)
    def circuit_shots(params):
        qml.BasisState(hf_state, wires=range(qubits))
        qml.UCCSD(params, wires=range(qubits), s_wires=s_wires,
                  d_wires=d_wires, init_state=hf_state)
        return qml.expval(H)

    # Evaluate 5 times to see variance
    evals = [float(circuit_shots(optimal_params)) for _ in range(5)]
    label = "exact" if shots is None else str(shots)
    print(f"shots={label:>5s}: mean={np.mean(evals):.4f}, std={np.std(evals):.4f}")

With exact simulation (shots=None), every evaluation returns the identical value and the standard deviation is zero. With 10,000 shots, the standard deviation is on the order of 0.001 Ha. With 100 shots, it grows to roughly 0.01 Ha or more, which is larger than the correlation energy we are trying to capture.

This has practical implications for VQE on hardware:

• Optimization: Shot noise makes gradient estimates noisy, which slows convergence and can prevent the optimizer from reaching chemical accuracy. Use exact simulation (shots=None) for prototyping and algorithm development.
• Final evaluation: Once you have optimized parameters, evaluate the final energy with a large shot count to reduce statistical error bars.
• Shot budget: The total number of shots required scales as O(1/epsilon^2) to achieve precision epsilon in the energy estimate.

Comparing PennyLane and Qiskit Approaches

If you have used Qiskit for quantum chemistry, you may find it helpful to see how the two frameworks compare for the same H2 VQE problem:

Aspect	PennyLane	Qiskit Nature
Hamiltonian	qml.qchem.molecular_hamiltonian	PySCFDriver + JordanWignerMapper
Ansatz	qml.UCCSD	UCCSD from circuit library
Optimizer	qml.AdamOptimizer	scipy.optimize or SLSQP
Gradient	Parameter-shift rule (automatic)	Parameter-shift via Estimator
Hardware	PennyLane device plugins	Qiskit Runtime
Best for	Quick prototyping, autodiff ecosystem	IBM hardware runs, production chemistry

PennyLane’s strength is its tight integration with automatic differentiation. You write a quantum circuit, decorate it with @qml.qnode, and PennyLane handles gradient computation transparently. Qiskit Nature offers deeper integration with IBM’s quantum hardware and a broader set of classical chemistry solvers for benchmarking.

Both frameworks use the same underlying physics (Jordan-Wigner mapping, UCCSD excitations, parameter-shift gradients), so the results should match to numerical precision for the same molecular geometry and basis set.

Common Mistakes

Using regular NumPy instead of PennyLane’s NumPy

# WRONG: breaks gradient tracking
import numpy as np
params = np.zeros(n_params)

# CORRECT: PennyLane's numpy wraps arrays to track gradients
from pennylane import numpy as np
params = np.zeros(n_params, requires_grad=True)

PennyLane’s numpy module wraps standard NumPy arrays with gradient-tracking metadata. If you use regular NumPy, PennyLane cannot compute gradients through parameter-shift, and the optimizer will not update the parameters. This is the single most common source of “my VQE isn’t converging” bugs.

Forgetting requires_grad=True

# WRONG: PennyLane cannot differentiate with respect to this array
params = np.zeros(n_params)

# CORRECT: explicitly mark the array as differentiable
params = np.zeros(n_params, requires_grad=True)

Without requires_grad=True, PennyLane treats the parameters as fixed constants. The optimizer’s step will have zero gradient, and the energy will never change from its initial value.

Using finite shots during optimization

# WRONG for prototyping: shot noise makes gradients noisy
dev = qml.device("default.qubit", wires=qubits, shots=1000)

# CORRECT for prototyping: exact expectation values
dev = qml.device("default.qubit", wires=qubits)

Shot noise adds variance to gradient estimates, which slows convergence dramatically. Use exact simulation (shots=None, which is the default) when developing and debugging your VQE workflow. Only switch to finite shots when you specifically want to study the effect of sampling noise or when running on hardware.

Setting coordinates in angstrom instead of bohr

# WRONG: 0.74 is in angstrom, but the function expects bohr
coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.74])

# CORRECT: convert angstrom to bohr (1 angstrom = 1.8897 bohr)
coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.74 * 1.8897])  # = 1.3984 bohr

The molecular_hamiltonian function expects atomic coordinates in bohr (atomic units). Passing angstrom values without conversion places the atoms at the wrong distance, producing a Hamiltonian for a different molecular geometry and giving incorrect energies. The resulting energy may still look reasonable (it is a valid energy for some geometry), making this bug hard to catch without careful validation.

Starting from random initial parameters

# SUBOPTIMAL: random start has no physical motivation
params = np.random.uniform(-np.pi, np.pi, n_params, requires_grad=True)

# BETTER: zero parameters give the Hartree-Fock state
params = np.zeros(n_params, requires_grad=True)

For UCCSD, setting all parameters to zero produces exactly the Hartree-Fock state, which is already a reasonable approximation (capturing about 99% of the total energy for H2). Starting from this known good state, the optimizer only needs to find the relatively small corrections that capture correlation energy. Random initial parameters place you at an arbitrary point in parameter space, potentially far from the minimum, and the optimizer may converge to a local minimum or take many more steps to reach the global one.

What to Try Next

• Run the potential energy surface scan above and compare with the Hartree-Fock PES (set all UCCSD params to zero) to visualize where correlation energy matters most
• Increase the number of StronglyEntanglingLayers and observe how HEA accuracy improves with circuit depth
• Try a different molecule like LiH (6 qubits in minimal basis) to see how VQE scales
• Explore excited states using the folded spectrum method or the subspace search VQE (SSVQE) approach
• Add noise to the simulation using qml.device("default.mixed", wires=qubits) and study how decoherence affects energy estimates
• See the PennyLane Reference for the full qchem API