VQE for H2 Molecule: Complete Worked Example in Qiskit

Why H2?

The hydrogen molecule is the canonical testbed for quantum chemistry on quantum computers. It has the smallest non-trivial electronic structure problem, a known exact ground state energy, and a well-understood mapping to qubits. Mastering VQE on H2 gives you the scaffolding to tackle larger molecules like LiH, BeH2, and water.

The ground state energy of H2 at equilibrium bond length (~0.74 angstrom) is approximately -1.137 Hartree. VQE should reproduce this from a parameterized ansatz optimized by a classical optimizer. Because H2 is so small, you can verify the VQE result against the exact Full Configuration Interaction (FCI) solution, making it the ideal molecule for learning, debugging, and benchmarking.

The Variational Principle

VQE is built on the variational principle from quantum mechanics. For any normalized trial state |psi(theta)>, the expectation value of the Hamiltonian is always greater than or equal to the true ground state energy:

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

The inequality becomes an equality only when |psi(theta)> is the exact ground state. For any other state, including the Hartree-Fock reference, the energy is strictly above E_0. This has three practical consequences for VQE:

1. Upper bound guarantee. Every VQE evaluation gives an upper bound on the true ground state energy. The optimizer walks downhill on the energy landscape, and the lowest energy it finds is still above (or equal to) the exact answer.
2. Ansatz quality metric. Comparing the VQE energy to the FCI energy tells you how expressive your ansatz is. If the gap is large, the ansatz cannot represent the ground state well.
3. Chemical accuracy benchmark. In computational chemistry, “chemical accuracy” means the result is within 1 kcal/mol of the exact answer. In Hartree, that threshold is about 1.6 milliHartree (0.0016 Ha). VQE results within this window are considered reliable for predicting reaction energies and molecular properties.

Later in this tutorial, you will compare the Hartree-Fock energy (the initial point at theta = 0) against the optimized VQE energy to see the variational principle in action.

Setup

pip install qiskit qiskit-nature qiskit-aer qiskit-algorithms pyscf scipy

Prerequisite: All code in this tutorial requires PySCF (pip install pyscf). PySCF runs on Linux and macOS only; on Windows, use WSL2. If you want a VQE example that runs without PySCF, see the Hybrid Quantum-Classical VQE tutorial.

Step 1: Define the Molecular Problem

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper

# H2 at equilibrium bond length (0.735 angstrom)
driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.735",
    basis="sto-3g",
    charge=0,
    spin=0,
)

problem = driver.run()
print(f"Number of spatial orbitals: {problem.num_spatial_orbitals}")  # 2
print(f"Number of particles: {problem.num_particles}")               # (1, 1) = 1 alpha, 1 beta
print(f"Nuclear repulsion energy: {problem.nuclear_repulsion_energy:.6f} Ha")

In the minimal STO-3G basis, H2 has 2 spatial orbitals and 2 electrons (1 alpha, 1 beta). Each spatial orbital splits into 2 spin-orbitals (alpha and beta), yielding 4 spin-orbitals total. After Jordan-Wigner mapping, this becomes a 4-qubit problem.

The nuclear_repulsion_energy is the classical Coulomb repulsion between the two protons. In atomic units, this equals 1/R where R is the internuclear distance in bohr. At 0.735 angstrom = 1.389 bohr, the nuclear repulsion is 1/1.389 = 0.720 Ha. This constant must be added to the electronic energy to get the physically meaningful total energy. The electronic Hamiltonian does not include nuclear repulsion because it only describes electron behavior in the field of fixed nuclei (the Born-Oppenheimer approximation).

Step 2: Map to Qubits and Inspect the Hamiltonian

mapper = JordanWignerMapper()
second_q_op = problem.hamiltonian.second_q_op()
qubit_op = mapper.map(second_q_op)

print(f"Number of qubits: {qubit_op.num_qubits}")    # 4
print(f"Number of Pauli terms: {len(qubit_op)}")      # 15
print("\nFull Pauli Hamiltonian:")
print(qubit_op)

Jordan-Wigner maps fermionic creation/annihilation operators to tensor products of Pauli operators. Each spin-orbital maps to one qubit: occupation = |1>, vacancy = |0>. The resulting Hamiltonian is a sum of Pauli strings whose expectation values give the electronic energy.

Understanding the H2 Pauli Hamiltonian

When you print the qubit Hamiltonian for H2 at the STO-3G level, you see a weighted sum of Pauli strings. The terms group into physically meaningful categories:

• IIII (identity): A constant energy offset that combines parts of the one-electron integrals. This is always present and sets the overall energy scale.
• ZZII, IIZZ, ZIZI, IZIZ, ZIIZ, IZZI (Z-only terms): These encode one-body energies (kinetic energy and nuclear attraction of individual electrons) and two-body diagonal Coulomb repulsion. Z operators measure whether a qubit is occupied (eigenvalue -1 for |1>) or vacant (eigenvalue +1 for |0>), so these terms assign energy contributions based on the occupation pattern.
• XXYY, YYXX, XXXX, YYYY (XX+YY type terms): These are the hopping and correlation terms. They flip pairs of qubits simultaneously, creating superpositions between different occupation configurations. Physically, they represent electron excitations from occupied to virtual orbitals. These terms are what make the problem non-trivial: without them, the Hartree-Fock state would already be exact. The UCCSD ansatz must capture these correlations.

The coefficients of all these terms change with bond length, but the Pauli string structure stays the same for H2 in STO-3G. At equilibrium, the XX+YY hopping terms are relatively small, meaning Hartree-Fock is already a good approximation. At stretched geometries (large R), these terms grow, and correlation becomes essential.

Step 3: Build the UCCSD Ansatz

Unitary Coupled Cluster with Singles and Doubles (UCCSD) is the chemistry-motivated ansatz of choice for VQE. It starts from the Hartree-Fock reference state and applies parameterized single and double excitation operators:

|psi(theta)> = exp(T(theta) - T^dagger(theta)) |HF>

where T includes single excitation operators (promoting one electron) and double excitation operators (promoting two electrons simultaneously).

from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock

hf_state = HartreeFock(
    num_spatial_orbitals=problem.num_spatial_orbitals,
    num_particles=problem.num_particles,
    qubit_mapper=mapper,
)

ansatz = UCCSD(
    num_spatial_orbitals=problem.num_spatial_orbitals,
    num_particles=problem.num_particles,
    qubit_mapper=mapper,
    initial_state=hf_state,
)

print(f"Ansatz parameters: {ansatz.num_parameters}")  # 1

Why Exactly 1 Parameter for H2?

With 2 electrons in 4 spin-orbitals (STO-3G basis), the excitation space is extremely limited. You can inspect the excitation list directly:

# Examine the excitation structure
singles = ansatz.excitation_list
# For UCCSD, access singles and doubles separately:
print(f"Total parameters: {ansatz.num_parameters}")  # 1

There is exactly 1 symmetry-allowed double excitation: promoting both electrons from the occupied spin-orbitals (0, 1) to the virtual spin-orbitals (2, 3). This corresponds to the excitation (0, 1, 2, 3), which takes the Hartree-Fock configuration |1100> to the doubly-excited configuration |0011>. Single excitations are either zero by symmetry or redundant with the double excitation for this two-electron system.

This single parameter theta controls a rotation in the two-dimensional subspace spanned by |1100> and |0011>:

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

For two-electron systems, UCCSD with this one parameter is equivalent to Full Configuration Interaction (FCI). This means the VQE result is exact (up to optimizer convergence) for H2.

Drawing the Ansatz Circuit

To see the actual gates, decompose the high-level ansatz into primitive operations:

decomposed = ansatz.decompose()
print(decomposed.draw(output='text', fold=120))

The circuit has two distinct sections:

1. Hartree-Fock initialization: X gates on qubits 0 and 1 flip them from |0> to |1>, preparing the state |1100>. This is the mean-field reference where both electrons occupy the lowest-energy spin-orbitals.
2. Parameterized excitation operator: A sequence of CNOT gates, single-qubit rotations, and more CNOT gates implement the unitary exp(i * theta * (excitation - de-excitation)). This creates the superposition between |1100> and |0011> controlled by the parameter theta.

If you call ansatz.draw() without decomposing first, you often see a single high-level block labeled “UCCSD” with no internal detail. Always decompose before drawing if you want to inspect the gate-level structure.

Step 4: Run VQE

from qiskit_aer.primitives import Estimator
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import SLSQP
import numpy as np

estimator = Estimator()

optimizer = SLSQP(maxiter=200)

vqe = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=optimizer,
    initial_point=np.zeros(ansatz.num_parameters),
)

vqe_result = vqe.compute_minimum_eigenvalue(qubit_op)

print(f"VQE energy (electronic): {vqe_result.eigenvalue.real:.6f} Ha")

Setting initial_point=np.zeros(...) starts the optimization at the Hartree-Fock state, since all excitation amplitudes are zero. This is standard practice in quantum chemistry VQE because Hartree-Fock is typically a reasonable starting point. Starting from random parameters risks landing in local minima or requiring many more iterations.

Comparing Hartree-Fock to VQE

To see the variational principle at work, evaluate the energy at the initial point (theta = 0, which is the Hartree-Fock state) and compare it to the optimized result:

# Hartree-Fock energy = energy at theta = 0 (the initial point)
from qiskit.primitives import StatevectorEstimator

sv_estimator = StatevectorEstimator()
hf_circuit = ansatz.assign_parameters(np.zeros(ansatz.num_parameters))
hf_job = sv_estimator.run([(hf_circuit, qubit_op)])
hf_electronic_energy = hf_job.result()[0].data.evs.item()

nuclear_repulsion = problem.nuclear_repulsion_energy

print(f"Hartree-Fock electronic energy: {hf_electronic_energy:.6f} Ha")
print(f"Hartree-Fock total energy:      {hf_electronic_energy + nuclear_repulsion:.6f} Ha")
print(f"VQE electronic energy:          {vqe_result.eigenvalue.real:.6f} Ha")
print(f"VQE total energy:               {vqe_result.eigenvalue.real + nuclear_repulsion:.6f} Ha")
print(f"Correlation energy captured:    {hf_electronic_energy - vqe_result.eigenvalue.real:.6f} Ha")

The difference between the Hartree-Fock energy and the VQE energy is the correlation energy. For H2 at equilibrium, this is roughly 0.015 Ha, which is small but chemically significant. The variational principle guarantees that the VQE result is lower (better) than Hartree-Fock, and for this system, it reaches the exact FCI energy.

Monitoring Convergence with a Callback

VQE accepts a callback function that receives intermediate results at each optimizer step. This is invaluable for diagnosing convergence issues:

iteration_energies = []

def vqe_callback(eval_count, params, energy, std):
    iteration_energies.append(energy)
    if eval_count % 10 == 0:
        print(f"Iteration {eval_count}: E = {energy:.6f} Ha")

vqe_with_callback = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=SLSQP(maxiter=200),
    initial_point=np.zeros(ansatz.num_parameters),
    callback=vqe_callback,
)

result_cb = vqe_with_callback.compute_minimum_eigenvalue(qubit_op)

print(f"\nFinal energy: {result_cb.eigenvalue.real:.6f} Ha")
print(f"Total iterations: {len(iteration_energies)}")
print(f"Energy drop from first to last: {iteration_energies[0] - iteration_energies[-1]:.6f} Ha")

For H2 with SLSQP, convergence is fast: the optimizer typically finds the minimum within 10-20 iterations. For larger molecules with more parameters, convergence can take hundreds of iterations, and plotting the callback energies helps identify whether the optimizer is stuck in a local minimum or still making progress.

Step 5: Convert to Total Energy

The VQE solver returns the electronic energy only. To get the total ground-state energy including nuclear repulsion, use GroundStateEigensolver, which handles this addition automatically:

from qiskit_nature.second_q.algorithms import GroundStateEigensolver

solver = GroundStateEigensolver(mapper, vqe)
ground_state_result = solver.solve(problem)

print(f"Total ground state energy: {ground_state_result.total_energies[0]:.6f} Ha")
print(f"Reference (FCI):           -1.137270 Ha")

The VQE result should agree with the FCI value to within numerical precision on the statevector simulator, because UCCSD is exact for two-electron systems.

A common mistake is to report the electronic eigenvalue directly as “the ground state energy.” The electronic energy alone is not physically meaningful for comparing to experiment. You always need to add nuclear repulsion:

Total energy = Electronic energy + Nuclear repulsion energy

The GroundStateEigensolver wrapper handles this, but if you use the raw VQE result, remember to add problem.nuclear_repulsion_energy yourself.

Step 6: Optimizer Comparison

The choice of classical optimizer significantly affects VQE performance, especially on real hardware. Here are three commonly used optimizers and when to choose each:

SLSQP (Sequential Least-Squares Programming)

Gradient-based optimizer that uses exact gradients or finite differences. It converges quickly on smooth energy landscapes and is the best choice for statevector simulation where gradient evaluations are cheap and noise-free.

from qiskit_algorithms.optimizers import SLSQP
# Already used above
optimizer_slsqp = SLSQP(maxiter=200)

COBYLA (Constrained Optimization by Linear Approximation)

Gradient-free optimizer that uses only function evaluations (no gradients). It builds a linear approximation of the objective at each step. Because it never computes gradients, it handles noisy energy landscapes better than SLSQP. This makes it a strong choice for shot-based simulation or near-term hardware where each energy evaluation has statistical noise.

from qiskit_algorithms.optimizers import COBYLA

vqe_cobyla = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=COBYLA(maxiter=200),
    initial_point=np.zeros(ansatz.num_parameters),
)
result_cobyla = vqe_cobyla.compute_minimum_eigenvalue(qubit_op)
print(f"COBYLA energy: {result_cobyla.eigenvalue.real:.6f} Ha")

SPSA (Simultaneous Perturbation Stochastic Approximation)

Gradient-free stochastic optimizer that estimates the gradient using only 2 function evaluations per step, regardless of the number of parameters. This property makes SPSA scale well for circuits with many parameters. It is the standard choice for VQE on real quantum hardware.

from qiskit_algorithms.optimizers import SPSA

vqe_spsa = VQE(
    estimator=estimator,
    ansatz=ansatz,
    optimizer=SPSA(maxiter=200),
    initial_point=np.zeros(ansatz.num_parameters),
)
result_spsa = vqe_spsa.compute_minimum_eigenvalue(qubit_op)
print(f"SPSA energy: {result_spsa.eigenvalue.real:.6f} Ha")

For the 1-parameter H2 problem, all three optimizers converge to the same answer. The differences become pronounced for larger molecules: SLSQP is fastest in simulation, COBYLA is more robust to noise, and SPSA is most scalable on hardware.

Step 7: Parity Mapping for Qubit Reduction

Jordan-Wigner is the simplest qubit mapping, but it uses one qubit per spin-orbital (4 qubits for H2). The Parity mapping can exploit particle number and spin symmetries to reduce the qubit count:

from qiskit_nature.second_q.mappers import ParityMapper

parity_mapper = ParityMapper(num_particles=problem.num_particles)
qubit_op_parity = parity_mapper.map(second_q_op)

print(f"Jordan-Wigner qubits: {qubit_op.num_qubits}")        # 4
print(f"Parity (reduced) qubits: {qubit_op_parity.num_qubits}")  # 2

When you pass num_particles to ParityMapper, it automatically identifies and removes 2 qubits corresponding to the conserved total particle number and total spin symmetry. For H2, this reduces the problem from 4 qubits to 2 qubits, cutting the circuit depth roughly in half.

The 2-qubit Hamiltonian contains the same physics as the 4-qubit version but operates in the symmetry-restricted subspace. You can build a corresponding reduced ansatz and run VQE on it:

hf_parity = HartreeFock(
    num_spatial_orbitals=problem.num_spatial_orbitals,
    num_particles=problem.num_particles,
    qubit_mapper=parity_mapper,
)

ansatz_parity = UCCSD(
    num_spatial_orbitals=problem.num_spatial_orbitals,
    num_particles=problem.num_particles,
    qubit_mapper=parity_mapper,
    initial_state=hf_parity,
)

vqe_parity = VQE(
    estimator=estimator,
    ansatz=ansatz_parity,
    optimizer=SLSQP(maxiter=200),
    initial_point=np.zeros(ansatz_parity.num_parameters),
)

result_parity = vqe_parity.compute_minimum_eigenvalue(qubit_op_parity)
total_parity = result_parity.eigenvalue.real + problem.nuclear_repulsion_energy
print(f"Parity-mapped VQE total energy: {total_parity:.6f} Ha")

On real hardware, fewer qubits means fewer CNOT gates, less decoherence, and better results. Always use symmetry reduction when possible.

Step 8: Potential Energy Surface

Running VQE at multiple bond lengths traces out the potential energy surface (PES), which reveals the equilibrium geometry and dissociation behavior:

import numpy as np

bond_lengths = np.linspace(0.3, 3.0, 20)
energies = []

for r in bond_lengths:
    driver_r = PySCFDriver(atom=f"H 0 0 0; H 0 0 {r}", basis="sto-3g")
    problem_r = driver_r.run()
    qubit_op_r = mapper.map(problem_r.hamiltonian.second_q_op())

    vqe_r = VQE(
        estimator=Estimator(),
        ansatz=UCCSD(
            num_spatial_orbitals=problem_r.num_spatial_orbitals,
            num_particles=problem_r.num_particles,
            qubit_mapper=mapper,
            initial_state=HartreeFock(
                problem_r.num_spatial_orbitals,
                problem_r.num_particles,
                mapper,
            ),
        ),
        optimizer=SLSQP(maxiter=200),
        initial_point=[0.0],
    )

    result_r = vqe_r.compute_minimum_eigenvalue(qubit_op_r)
    # Add nuclear repulsion for the total energy
    nuclear_repulsion = problem_r.nuclear_repulsion_energy
    energies.append(result_r.eigenvalue.real + nuclear_repulsion)

for r, e in zip(bond_lengths, energies):
    print(f"r = {r:.2f} A  E = {e:.6f} Ha")

Reading the Potential Energy Surface

The PES encodes all the key molecular properties:

Equilibrium bond length. The minimum of the PES curve gives the equilibrium geometry. Find it programmatically:

min_idx = np.argmin(energies)
r_eq = bond_lengths[min_idx]
e_eq = energies[min_idx]
print(f"Equilibrium bond length: {r_eq:.2f} angstrom")
print(f"Equilibrium energy: {e_eq:.6f} Ha")

For H2 in STO-3G, the equilibrium bond length is approximately 0.74 angstrom, matching the experimental value of 0.741 angstrom.

Dissociation energy. The well depth, from the minimum to the energy at large separation, gives the bond dissociation energy. At large R, the molecule dissociates into two isolated hydrogen atoms. In the STO-3G basis, each isolated H atom has energy -0.4668 Ha, so the dissociation limit is approximately -0.9335 Ha. The dissociation energy is:

e_dissociation = energies[-1] - e_eq  # Energy at largest R minus equilibrium
print(f"Dissociation energy: {e_dissociation:.4f} Ha")
print(f"Dissociation energy: {e_dissociation * 627.509:.1f} kcal/mol")

Vibrational frequency. Near the minimum, the PES is approximately parabolic. Fitting a harmonic potential E = E_0 + (1/2)k(r - r_0)^2 gives the force constant k, from which the vibrational frequency follows. This harmonic approximation breaks down far from the minimum.

Nuclear repulsion at each point. At each bond length r, the nuclear repulsion energy is 1/r (in atomic units, with r in bohr). At short bond lengths, nuclear repulsion dominates and the total energy shoots upward. At large bond lengths, nuclear repulsion goes to zero and the total energy approaches the sum of isolated atom energies.

Step 9: FreezeCoreTransformer for Larger Molecules

For molecules beyond H2, some electrons occupy tightly bound core orbitals (inner-shell) that barely participate in chemistry. The FreezeCoreTransformer removes these orbitals to reduce the qubit count:

from qiskit_nature.second_q.transformers import FreezeCoreTransformer

# Example: LiH with frozen Li 1s core orbital
driver_lih = PySCFDriver(atom="Li 0 0 0; H 0 0 1.6", basis="sto-3g")
problem_lih = driver_lih.run()

print(f"LiH full spatial orbitals: {problem_lih.num_spatial_orbitals}")      # 6
print(f"LiH full particles: {problem_lih.num_particles}")                    # (2, 2)

freeze_core = FreezeCoreTransformer(freeze_core=True)
reduced_problem = freeze_core.transform(problem_lih)

print(f"LiH active spatial orbitals: {reduced_problem.num_spatial_orbitals}")  # 5
print(f"LiH active particles: {reduced_problem.num_particles}")               # (1, 1)

Freezing the Li 1s core orbital removes 1 spatial orbital (2 spin-orbitals = 2 qubits) and accounts for 2 electrons. The frozen electrons contribute a constant energy that gets added back automatically. This is the same core approximation used routinely in classical quantum chemistry.

For H2 itself, there are no core orbitals to freeze, so this transformer has no effect. But for any molecule with second-row atoms (Li, Be, B, C, N, O, F) or heavier, freezing core orbitals is standard practice and can save 2 or more qubits per heavy atom.

Step 10: Running on Real IBM Hardware

To run VQE on actual quantum hardware through IBM Quantum, replace the local Estimator with IBM Runtime’s estimator. The key steps are: transpile the ansatz for the target backend, then map the Hamiltonian to the transpiled circuit’s qubit layout.

from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2, Session
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_algorithms.optimizers import COBYLA

# Connect to IBM Quantum (requires saved credentials)
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(f"Using backend: {backend.name}")

# Transpile the ansatz for the target hardware
pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_op = qubit_op.apply_layout(isa_ansatz.layout)

# Run VQE on hardware
with Session(backend=backend) as session:
    estimator = EstimatorV2(mode=session)
    vqe_hw = VQE(
        estimator=estimator,
        ansatz=isa_ansatz,
        optimizer=COBYLA(maxiter=50),
        initial_point=np.zeros(ansatz.num_parameters),
    )
    result_hw = vqe_hw.compute_minimum_eigenvalue(isa_op)

total_hw = result_hw.eigenvalue.real + problem.nuclear_repulsion_energy
print(f"Hardware VQE total energy: {total_hw:.6f} Ha")

Important notes for hardware runs:

• Transpile first. The ansatz must be compiled into the native gate set and connectivity of the target device. The generate_preset_pass_manager handles this.
• Apply layout to the operator. After transpilation, qubits may be reordered. apply_layout maps the Hamiltonian to match the physical qubit assignment.
• Use COBYLA or SPSA. Gradient-based optimizers like SLSQP perform poorly on hardware because each energy evaluation has shot noise. COBYLA is the recommended starting point for hardware VQE.
• Reduce iterations. Hardware jobs are expensive in both time and cost. Start with 50 iterations and increase only if needed.
• Expect noise. Hardware results will be above the exact energy (the variational principle still holds for the ideal circuit, but noise can push the measured energy in either direction). Error mitigation techniques, available through the EstimatorV2 options, can improve accuracy.

Common Mistakes

PySCF not available on Windows

PySCF does not support Windows natively. Use WSL2 (Windows Subsystem for Linux) or a Docker container. Alternatively, construct the Hamiltonian manually using qiskit_nature.second_q.hamiltonians.ElectronicEnergy with pre-computed one- and two-electron integrals.

Reporting electronic energy as total energy

The most common beginner mistake is reporting the VQE eigenvalue directly without adding nuclear repulsion. The electronic energy at equilibrium H2 is around -1.86 Ha, while the total energy is around -1.14 Ha. If your “ground state energy” is below -1.5 Ha, you are likely reporting the electronic component only.

Starting from random parameters

VQE starts at random parameters by default if you do not set initial_point. Random starting points can land in local minima, especially for multi-parameter ansatze. For chemistry problems, always start at zeros:

initial_point = np.zeros(ansatz.num_parameters)

This corresponds to the Hartree-Fock state, which is a physically motivated and usually reasonable starting point.

Too few optimizer iterations

For H2 with 1 parameter, SLSQP converges in under 20 iterations. But for larger molecules with 10-100 parameters, the optimizer may need 500 or more iterations. If your VQE energy is noticeably above the expected value, try increasing maxiter before concluding the ansatz is inadequate.

Not decomposing the ansatz before drawing

Calling ansatz.draw() directly often displays a single opaque block labeled “UCCSD.” To see the actual gate structure, always decompose first:

# This shows a black box:
print(ansatz.draw(output='text'))

# This shows the actual gates:
print(ansatz.decompose().draw(output='text', fold=120))

You may need to call .decompose() multiple times to reach the primitive gate level.

Key Takeaways

1. UCCSD is exact for H2. With 2 electrons in 4 spin-orbitals, UCCSD captures the full correlation energy through a single parameter. This makes H2 an ideal validation target: if your VQE implementation does not match FCI for H2, something is wrong.

2. The variational principle guarantees upper bounds. Every VQE energy is above or equal to the true ground state energy. Comparing to FCI quantifies ansatz quality.

3. Nuclear repulsion is not optional. Always add problem.nuclear_repulsion_energy to the electronic eigenvalue when reporting total energies.

4. Symmetry reduction saves qubits. ParityMapper with particle number specified reduces H2 from 4 to 2 qubits. For larger molecules, this savings is critical for fitting the problem onto real hardware.

5. Optimizer choice depends on the execution environment. Use SLSQP for noiseless simulation, COBYLA for noisy or shot-based runs, and SPSA for many-parameter circuits on hardware.

6. The PES reveals molecular properties. Bond length, dissociation energy, and vibrational frequency all come from scanning VQE across different geometries.

7. Scale up carefully. FreezeCoreTransformer, symmetry reduction, and hardware-efficient ansatze are the tools for moving beyond H2 to chemically interesting molecules.