Our PennyLane VQE guide builds a variational chemistry solver from PennyLane's side. This post covers the same problem, finding a molecule's ground-state energy, from the Qiskit ecosystem, using Qiskit Nature (version 0.8), which handles the chemistry-specific plumbing: driving a classical chemistry package, building the molecular Hamiltonian, and mapping it to qubits.
pip install qiskit-nature pyscf
Driving the classical chemistry
Qiskit Nature doesn't compute molecular orbitals itself. It drives an external classical chemistry package, PySCF here, to do that groundwork, then wraps the result in a form the rest of the pipeline understands:
from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
driver = PySCFDriver(
atom="H 0 0 0; H 0 0 0.735",
basis="sto3g",
charge=0,
spin=0,
unit=DistanceUnit.ANGSTROM,
)
problem = driver.run()
This example builds H₂ at its equilibrium bond length (0.735 Å), the standard first molecule for testing a chemistry pipeline, small enough to solve exactly and check your setup against a known answer. basis="sto3g" sets a minimal basis set, the smallest reasonable choice for orbital representation, which keeps qubit count low at the cost of chemical accuracy, a trade-off worth understanding before scaling to a real molecule of interest.
Mapping fermions to qubits
problem describes the molecule in terms of fermionic operators, the natural language of electronic structure, which then needs mapping onto qubits:
from qiskit_nature.second_q.mappers import JordanWignerMapper
mapper = JordanWignerMapper()
Jordan-Wigner is the most direct mapping: one spin-orbital maps to one qubit, which makes it the easiest to reason about but not the most qubit-efficient. Our companion post on qubit mappers covers the alternatives (ParityMapper, BravyiKitaevMapper) and the qubit-count trade-offs each one makes.
Solving classically first
Before running anything on a quantum circuit, solve the problem with a classical eigensolver. This gives you a known-correct answer to check any quantum approach against, and for a molecule this small, it's also faster:
from qiskit_algorithms import NumPyMinimumEigensolver
from qiskit_nature.second_q.algorithms import GroundStateEigensolver
solver = GroundStateEigensolver(mapper, NumPyMinimumEigensolver())
result = solver.solve(problem)
print(result.total_energies)
For H₂ at equilibrium bond length, this returns a ground-state energy close to -1.137 Hartree, the reference value against which any quantum solution should be checked. If your VQE result later doesn't land near this number, the bug is worth chasing before assuming it's a hardware noise issue.
Swapping in VQE
GroundStateEigensolver takes any compatible minimum eigensolver, which is what makes moving from classical to quantum mostly a one-line change rather than a rewrite:
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Estimator
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
ansatz = UCCSD(
problem.num_spatial_orbitals,
problem.num_particles,
mapper,
initial_state=HartreeFock(problem.num_spatial_orbitals, problem.num_particles, mapper),
)
vqe_solver = VQE(Estimator(), ansatz, COBYLA())
vqe_solver.initial_point = [0.0] * ansatz.num_parameters
solver = GroundStateEigensolver(mapper, vqe_solver)
result = solver.solve(problem)
print(result.total_energies)
UCCSD (Unitary Coupled Cluster Singles and Doubles) is the standard chemistry-motivated ansatz, built to respect the physical structure of electron excitations rather than being a generic parameterized circuit. HartreeFock initializes the circuit in the classical mean-field reference state, a far better starting point than an arbitrary superposition, since the true ground state is usually a comparatively small correction away from it.
Why the two results should (nearly) match
For a small, well-behaved molecule like H₂ in a minimal basis, VQE with UCCSD should converge close to the exact classical answer, since UCCSD is expressive enough to represent the exact ground state for a system this small. That agreement is the actual point of running the comparison: it validates that the mapping, ansatz, and optimizer are all working correctly before you scale to a molecule too large to check against a classical solver at all, which is the actual regime where a quantum approach would need to earn its keep.
Try this next
- Sweep the H-H bond length and plot the ground-state energy curve, the classic dissociation-curve exercise that confirms your solver reproduces real molecular behavior rather than a single lucky data point.
- Swap
JordanWignerMapperforParityMapper(num_particles=problem.num_particles)and compare qubit count, covered in our qubit mappers post. - Read what quantum computers do with molecules in 2026 for where this kind of pipeline stands on real, larger molecules rather than toy H₂.