Water Molecule Simulation with Qiskit Nature and VQE

Quantum chemistry is one of the most promising near-term applications for quantum computing. Simulating the electronic structure of molecules, finding ground state energies, bond lengths, reaction barriers, requires exponential classical resources in the worst case, yet a quantum computer can represent molecular wavefunctions naturally. Water (H2O) is a canonical benchmark: small enough to be tractable on current hardware, but complex enough to require genuine electronic correlation beyond Hartree-Fock.

Setup and Dependencies

pip install qiskit-nature[pyscf] qiskit-algorithms scipy

Qiskit Nature is the quantum chemistry layer of the Qiskit ecosystem. It provides drivers for classical chemistry codes (PySCF, PSI4), problem formulations, active space transformations, qubit mappings, and ansatz circuits tuned for molecular wavefunctions.

Note: All code blocks in this tutorial require PySCF (pip install pyscf) in addition to qiskit-nature. Skip any block that imports PySCFDriver if PySCF is not installed.

Step 1: Classical Driver and Geometry

The PySCFDriver calls the PySCF classical chemistry code to compute one- and two-electron integrals. These integrals, combined with the nuclear repulsion energy, define the molecular Hamiltonian in second quantization.

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.units import DistanceUnit

# Water molecule in STO-3G basis, near equilibrium geometry
# O-H bond length: 0.9572 Angstrom, H-O-H angle: 104.52 degrees
driver = PySCFDriver(
    atom="O 0.0 0.0 0.0; H 0.757 0.586 0.0; H -0.757 0.586 0.0",
    basis="sto-3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

problem = driver.run()
print(f"Number of molecular orbitals: {problem.num_spatial_orbitals}")
print(f"Number of electrons: {problem.num_particles}")

hamiltonian = problem.hamiltonian
print(f"Nuclear repulsion energy: {hamiltonian.nuclear_repulsion_energy:.6f} Ha")

In STO-3G, water has 7 spatial orbitals and 10 electrons (8 core + valence from O, 1 each from two H). The full second-quantized Hamiltonian over 7 orbitals, after Jordan-Wigner mapping, would require 14 qubits, feasible but expensive.

Step 2: Active Space Reduction

The core 1s electrons of oxygen (the lowest two orbitals) contribute little to molecular bonding. Chemists routinely freeze these and work in an active space containing only the valence electrons and orbitals relevant to the chemical question.

The ActiveSpaceTransformer freezes occupied core orbitals and virtual (high-energy) orbitals, reducing the problem size:

# Requires: pyscf
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer

# Choose 4 active electrons in 4 active orbitals
# (freeze 2 core electrons and select the 4 most important valence orbitals)
transformer = ActiveSpaceTransformer(
    num_electrons=4,
    num_spatial_orbitals=4,
)

# Apply the transformation
reduced_problem = transformer.transform(problem)
print(f"Active space orbitals: {reduced_problem.num_spatial_orbitals}")
print(f"Active space electrons: {reduced_problem.num_particles}")
# Output: 4 orbitals, (2, 2) electrons (2 alpha, 2 beta)

The 4-electron, 4-orbital active space requires only 8 qubits after Jordan-Wigner mapping, compared to 14 for the full space. The frozen core electrons are accounted for by adjusting the one-electron integrals; their energy contribution is folded into constants.

For a production simulation you would consult the molecular orbital energies and occupation numbers (from the Hartree-Fock solution) to select the most important orbitals, but for H2O in STO-3G the default selection is well-established.

Step 3: Jordan-Wigner Mapping

The Jordan-Wigner transformation maps fermionic creation/annihilation operators to qubit Pauli operators. Each spin-orbital becomes one qubit; the occupied/unoccupied state maps to $|1\rangle$/$|0\rangle$.

from qiskit_nature.second_q.mappers import JordanWignerMapper

mapper = JordanWignerMapper()

# Map the fermionic Hamiltonian to a qubit operator
qubit_op = mapper.map(reduced_problem.hamiltonian.second_q_op())

print(f"Number of qubits: {qubit_op.num_qubits}")
print(f"Number of Pauli terms: {len(qubit_op)}")
# Typically: 8 qubits, ~250 Pauli terms for this active space

Alternative mappers include ParityMapper (can reduce qubit count by 2 using symmetry) and BravyiKitaevMapper (logarithmic weight per term). For H2O in this active space, the 2-qubit reduction from ParityMapper brings the circuit down to 6 qubits, which is worth using in a hardware run:

from qiskit_nature.second_q.mappers import ParityMapper

parity_mapper = ParityMapper(num_particles=reduced_problem.num_particles)
qubit_op_parity = parity_mapper.map(reduced_problem.hamiltonian.second_q_op())
print(f"Parity mapping qubits: {qubit_op_parity.num_qubits}")  # 6 qubits

Step 4: UCCSD Ansatz

The Unitary Coupled Cluster Singles and Doubles (UCCSD) ansatz is the standard choice for quantum chemistry VQE. It prepares the Hartree-Fock reference state and applies all single and double excitation operators:

$|\psi(\theta)\rangle = e^{T(\theta) - T^\dagger(\theta)} |\text{HF}\rangle$

where $T = T_1 + T_2$ contains single and double excitation amplitudes as variational parameters.

# Requires: pyscf
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_nature.second_q.mappers import JordanWignerMapper

# Hartree-Fock reference state (initial state)
hf_state = HartreeFock(
    num_spatial_orbitals=reduced_problem.num_spatial_orbitals,
    num_particles=reduced_problem.num_particles,
    qubit_mapper=mapper,
)

# UCCSD ansatz built on top of HF reference
ansatz = UCCSD(
    num_spatial_orbitals=reduced_problem.num_spatial_orbitals,
    num_particles=reduced_problem.num_particles,
    qubit_mapper=mapper,
    initial_state=hf_state,
)

print(f"UCCSD circuit depth: {ansatz.decompose().depth()}")
print(f"Number of variational parameters: {ansatz.num_parameters}")
# Typically: ~50-80 depth, 3 parameters for (2e, 4 orbital) active space

UCCSD automatically generates all chemically valid single and double excitations from the HF reference. For the 4-electron, 4-orbital active space of H2O, this gives a small number of parameters (often 2-4 for this active space due to symmetry), making convergence fast.

Step 5: VQE with SLSQP Optimizer

# Requires: pyscf
import numpy as np
from scipy.optimize import minimize
from qiskit.primitives import StatevectorEstimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_aer import AerSimulator

# For simulation, use statevector estimator (exact, no shot noise)
estimator = StatevectorEstimator()

# Track convergence
energies = []
eval_count = [0]

def cost_fn(params):
    pub = (ansatz, [qubit_op], [params])
    result = estimator.run([pub]).result()
    energy = float(result[0].data.evs[0])
    energies.append(energy)
    eval_count[0] += 1
    return energy

# Initial point: use HF amplitudes (all zeros for UCCSD)
x0 = np.zeros(ansatz.num_parameters)

# SLSQP: fast, gradient-free, good for chemistry (smooth landscape)
result = minimize(
    cost_fn,
    x0,
    method="SLSQP",
    options={"maxiter": 500, "ftol": 1e-9},
)

vqe_energy = result.fun
print(f"\nVQE Results")
print(f"===========")
print(f"Converged: {result.success}")
print(f"VQE ground state energy: {vqe_energy:.8f} Ha")
print(f"Function evaluations: {eval_count[0]}")

Step 6: Interpreting Results and Reference Comparison

The raw VQE energy is the electronic energy plus nuclear repulsion. Qiskit Nature provides utilities to interpret this in chemist-friendly units:

# Requires: pyscf
from qiskit_nature.second_q.algorithms import GroundStateEigensolver

# Use the problem's built-in result interpretation
class VQEWrapper:
    """Adapts the scipy result for Qiskit Nature's result interpreter."""
    def __init__(self, energy):
        self.eigenvalue = energy

# Reference energies for H2O in STO-3G (from classical codes)
hartree_fock_energy = -74.9629  # Ha (obtained from PySCF driver)
ccsd_energy = -75.0124          # Ha (classical gold standard for small molecules)
fci_energy = -75.0143           # Ha (exact within basis set)

print(f"\nEnergy Comparison (STO-3G basis)")
print(f"=================================")
print(f"Hartree-Fock:        {hartree_fock_energy:.6f} Ha")
print(f"VQE (UCCSD):         {vqe_energy:.6f} Ha")
print(f"CCSD (classical):    {ccsd_energy:.6f} Ha")
print(f"FCI (exact):         {fci_energy:.6f} Ha")
print()
print(f"Correlation energy captured by VQE: "
      f"{(vqe_energy - hartree_fock_energy) * 1000:.2f} mHa")
print(f"Error vs FCI: {(vqe_energy - fci_energy) * 1000:.2f} mHa")

UCCSD-VQE, with an exact (statevector) estimator, should match FCI to within numerical precision for this active space; UCCSD is exact for two-electron systems. The real test comes when running on hardware, where shot noise and gate errors degrade the result.

Step 7: Preparing for Hardware Execution

Requires hardware credentials: This block connects to IBM Quantum and requires a valid QiskitRuntimeService account token. Skip if running locally without hardware access.

To run on IBM Quantum hardware, transpile the UCCSD circuit to ISA format and use the Runtime Estimator:

# Requires: pyscf
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2

service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.least_busy(operational=True, simulator=False, min_num_qubits=8)

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_op = qubit_op.apply_layout(isa_ansatz.layout)

print(f"ISA circuit depth: {isa_ansatz.depth()}")
print(f"ISA circuit 2Q gates: {isa_ansatz.count_ops().get('cx', 0)}")

with Session(backend=backend) as session:
    estimator = EstimatorV2(mode=session)

    def hardware_cost_fn(params):
        pub = (isa_ansatz, [isa_op], [params])
        result = estimator.run([pub]).result()
        return float(result[0].data.evs[0])

    hw_result = minimize(hardware_cost_fn, x0, method="SLSQP",
                         options={"maxiter": 100, "ftol": 1e-5})
    print(f"Hardware VQE energy: {hw_result.fun:.6f} Ha")

Performance Considerations

Active space selection is the most important choice for accuracy. The 4-electron, 4-orbital space used here is minimal; for quantitative results, a (6, 6) or (8, 7) active space would be more appropriate, at the cost of more qubits and deeper circuits.

Basis set: STO-3G is a minimal basis for demonstration. Chemical accuracy requires at least cc-pVDZ or 6-31G*, which significantly increases the number of orbitals and therefore qubits.

Ansatz expressibility: UCCSD is hardware-expensive. Hardware-efficient ansatze (like EfficientSU2) run shallower circuits but require more optimizer iterations and may miss the global minimum. A growing literature explores ADAPT-VQE, which builds the ansatz adaptively and often converges with fewer parameters than full UCCSD.

Noise mitigation: On current hardware, running VQE without error mitigation typically gives energies 10-50 mHa above the ideal value. Using resilience_level=1 in the RuntimeEstimator provides readout mitigation and dynamical decoupling at modest cost overhead.

The water molecule simulation serves as a benchmark that spans the full stack of quantum chemistry on quantum hardware: problem construction, active space reduction, qubit mapping, variational optimization, and comparison to classical reference methods. Mastering this workflow provides the foundation for tackling larger, classically-intractable molecules.