⚛️ Variational Quantum Eigensolver Molecular Energy

From Molecular Geometry to Second-Quantized Electronic Hamiltonian

Every VQE calculation begins entirely classically. Given a molecular geometry and a chosen basis set, a Hartree-Fock self-consistent field calculation produces a set of molecular orbitals and the one- and two-electron integrals that define the electronic structure problem. These integrals are then packaged into the second-quantized fermionic Hamiltonian that the quantum computer will ultimately help solve — the point where quantum resources enter is deliberately delayed as long as possible, since qubits are the most expensive resource in the whole pipeline.

• STO-3G: Basis set (typical NISQ demo) (minimal basis, fewest orbitals)
• 4: H₂ spin-orbitals (2 spatial × 2 spin)
• PySCF / Psi4: Classical package (HF + integral generation)
• frozen core: Active space reduction (drops core electrons from qubits)

Second quantization and the electronic structure Hamiltonian

The non-relativistic, Born-Oppenheimer electronic Hamiltonian in second quantization:

Ĥ = Σ_pq h_pq a†_p a_q + (1/2) Σ_pqrs h_pqrs a†_p a†_q a_r a_s + E_nuc

where: • a†_p, a_q are fermionic creation/annihilation operators for spin-orbital p, q • h_pq = one-electron integrals (kinetic energy + electron-nuclear attraction), computed as ⟨p|−½∇² − Σ_A Z_A/|r−R_A||q⟩ • h_pqrs = two-electron (Coulomb repulsion) integrals ⟨pq|1/r_12|rs⟩ • E_nuc = classical nuclear-nuclear Coulomb repulsion constant, added at the end

Workflow (fully classical, pre-quantum): 1. Specify molecular geometry (bond lengths in Ångström) — e.g., H₂ at experimental equilibrium 0.735 Å 2. Choose basis set: STO-3G (minimal, 1 orbital per atom for H) is standard for NISQ demos since it minimizes qubit count; larger bases (6-31G, cc-pVDZ) give more accurate energies but require far more qubits 3. Run restricted Hartree-Fock (RHF) to obtain molecular orbital coefficients and the reference Slater determinant 4. Transform atomic-orbital integrals into the molecular-orbital basis (h_pq, h_pqrs) 5. Apply active-space reduction: freeze chemically inert core orbitals (e.g., O 1s in H₂O) and truncate high-lying virtual orbitals — this is the single biggest lever for keeping qubit counts NISQ-feasible

Qubit count scaling: • H₂ (STO-3G): 2 spatial orbitals × 2 spins = 4 spin-orbitals → 4 qubits (Jordan-Wigner, no reduction) or 2 qubits with parity/symmetry reduction • LiH (STO-3G, frozen core): 6 active spin-orbitals → 4-6 qubits typically used in demonstrations • BeH₂ (STO-3G, frozen core): 6-8 active spin-orbitals • H₂O (STO-3G, frozen core, active space): 8-10 qubits for a minimally reduced treatment • Full H₂O in cc-pVDZ basis: >50 spin-orbitals → far beyond current NISQ device coherence budgets

This exponential mismatch between chemically interesting basis sets and available qubit counts is precisely why VQE targets small molecules as its proving ground rather than pharmaceutically relevant systems today.

Mapping Fermions to Qubits — the Jordan-Wigner and Bravyi-Kitaev Transforms

Qubits are distinguishable, independently addressable two-level systems; electrons are indistinguishable fermions obeying the Pauli exclusion principle and anticommutation relations {a_p, a†_q} = δ_pq. Before any quantum circuit can represent a molecular wavefunction, the fermionic algebra of creation/annihilation operators must be translated into an isomorphic algebra of Pauli operators (X, Y, Z, I) acting on qubits — a purely mathematical transformation with no approximation, but with real cost in circuit complexity.

• O(N): Jordan-Wigner locality (Pauli string weight per term)
• O(log N): Bravyi-Kitaev locality (binary-tree encoding)
• 15: H₂ Hamiltonian terms (Pauli strings (JW, 4 qubits))
• −2 qubits: Symmetry reduction (parity + Z₂ tapering, typical)

Jordan-Wigner strings, Bravyi-Kitaev trees, and qubit tapering

Jordan-Wigner (JW) transform (Jordan & Wigner, 1928; adapted for quantum simulation by Somma et al. 2002):

a†_p = (X_p − iY_p)/2 ⊗ Z_1 ⊗ Z_2 ⊗ ... ⊗ Z_{p−1} a_p = (X_p + iY_p)/2 ⊗ Z_1 ⊗ Z_2 ⊗ ... ⊗ Z_{p−1}

• The trailing string of Z operators ("Jordan-Wigner string") enforces fermionic antisymmetry under qubit swap by tracking occupation parity of all lower-indexed orbitals • Simple, deterministic, exact — but a single fermionic operator can produce a Pauli string acting on up to N qubits, so two-electron terms (four fermionic operators) create long, deep strings of CNOTs when implemented as circuits • For H₂ in STO-3G (4 qubits), the qubit Hamiltonian after JW mapping has exactly 15 Pauli terms (identity, single Z, ZZ, and XXYY-type terms from the double-excitation piece) — small enough to write out by hand

Bravyi-Kitaev (BK) transform (Bravyi & Kitaev, 2002; Seeley, Richard, Love 2012 for chemistry): • Encodes partial sums of occupation numbers in a binary-tree structure rather than a fully cumulative string • Reduces the Pauli weight of both creation/annihilation operators and of update/parity operators from O(N) to O(log₂N) • Same physics, same spectrum, different (more circuit-efficient) encoding — becomes advantageous as qubit count grows past ~8-10

Qubit-count reduction techniques (applied after mapping): 1. Parity reduction: exploiting known total electron number and spin symmetry removes 2 qubits automatically (used routinely for H₂: 4 → 2 qubits) 2. Z₂ symmetry tapering (Bravyi et al. 2017): identifies additional Z₂ symmetries of the Hamiltonian (e.g., particle-number parity, spin-z parity) and removes one qubit per independent symmetry via a Clifford change of basis 3. Combined effect: LiH is often demonstrated on IBM hardware with as few as 4 qubits after aggressive tapering, versus 12 in the untapered JW encoding

The choice of mapping is a genuine engineering trade-off: JW gives short-depth single-excitation terms but long two-excitation strings; BK balances term weights more evenly; newer mappings (e.g., generalized superfast encoding) target specific hardware connectivity graphs.

UCCSD and Hardware-Efficient Ansätze — Preparing the Trial Wavefunction

The ansatz is the parameterized quantum circuit U(θ) applied to a reference state (typically the Hartree-Fock determinant) to produce the trial wavefunction |ψ(θ)⟩ = U(θ)|HF⟩. Its design is the central creative and practical challenge of VQE: it must be expressive enough to reach near the true ground state, yet shallow enough to survive decoherence and gate errors on real NISQ hardware — two goals in direct tension.

• O(N⁴): UCCSD gate scaling (excitation operators, N=qubits)
• O(p·N): Hardware-efficient depth (p=repetition layers)
• 1: H₂ UCCSD parameters (single double-excitation θ)
• ~20-40: LiH hardware-efficient CNOTs (typical p=2-3 ansatz)

UCCSD, Trotterization, and hardware-efficient circuit design

Unitary Coupled Cluster Singles and Doubles (UCCSD):

|ψ(θ)⟩ = e^{T(θ) − T†(θ)} |HF⟩

where T(θ) = T₁ + T₂ (singles + doubles excitation operators): T₁ = Σ_ia θ_i^a a†_a a_i (excite electron from occupied i to virtual a) T₂ = Σ_ijab θ_ij^ab a†_a a†_b a_j a_i (double excitations)

• Chemically motivated: directly descended from classical coupled-cluster theory (CCSD), the workhorse "gold standard-adjacent" method in quantum chemistry • The unitary version (rather than classical CC's non-unitary e^T) is required because quantum circuits must implement unitary operations • e^{T−T†} is not directly implementable; Trotterization approximates it as a product of exponentials of individual Pauli strings: e^{T−T†} ≈ Π_k e^{iθ_k P_k}, each exponential compiled to a short sequence of CNOT+rotation gates via the standard staircase construction • Trotter error is controllable but nonzero; higher-order Trotter formulas (Suzuki-Trotter) trade circuit depth for accuracy • H₂ minimal case: exactly one double-excitation amplitude θ needed — the famous single-parameter H₂ VQE circuit used in nearly every introductory demonstration (e.g., IBM Qiskit textbook, Google Cirq tutorials) • Gate count grows as O(N⁴) with qubit count — quickly becomes NISQ-prohibitive beyond ~12-16 qubits

Hardware-Efficient Ansatz (HEA) (Kandala et al., Nature 2017 — first multi-qubit VQE on real hardware, BeH₂):

• Alternates layers of single-qubit rotation gates (RY, RZ with free parameters) and a fixed pattern of entangling two-qubit gates (CNOT or CZ) matched to the device's native qubit connectivity • Circuit: U(θ) = Π_{l=1}^{p} [Entangler_l · Π_i RY(θ_{l,i})RZ(θ_{l,i})] • p = number of repeated layers (the "ansatz depth" slider): more layers = more expressive but more accumulated gate error • Not chemically motivated — parameters have no direct physical excitation meaning — but achieves far shallower circuits than UCCSD for a given qubit count, which matters enormously given ~10⁻³ two-qubit gate error rates on today's superconducting devices • Kandala et al. 2017 used exactly this ansatz to compute BeH₂ ground-state energy on 6 qubits on an IBM superconducting processor, the first VQE demonstration beyond textbook H₂

Barren plateau risk: as HEA depth and qubit count grow, gradients of the cost function with respect to θ vanish exponentially (McClean et al. 2018), flattening the optimization landscape and making the classical optimizer's job exponentially hard — a core open problem discussed further under QML.

Measuring ⟨H⟩ — Pauli Grouping, Shot Noise, and Readout on Real Hardware

A quantum computer cannot directly output an energy value; it can only be measured in the computational (Z) basis, collapsing each qubit to 0 or 1. To estimate ⟨ψ(θ)|H|ψ(θ)⟩, every distinct Pauli string in the mapped Hamiltonian must be measured separately (after appropriate basis-rotation gates), each requiring many repeated circuit executions ("shots") to beat down statistical shot noise — turning a single energy evaluation into thousands of physical circuit runs.

• 1,000–8,000: Typical shots per term (per Pauli-string measurement)
• 1/√shots: Shot-noise scaling (standard error of the mean)
• ~5–8×: Qubit-wise commuting groups (reduction in circuits needed)
• ~15: H₂ energy evaluations/run (grouped Pauli measurements)

Pauli grouping, basis rotation, and error mitigation on NISQ devices

Measuring a weighted sum of Pauli strings:

⟨H⟩ = Σ_i c_i ⟨ψ(θ)|P_i|ψ(θ)⟩

Each Pauli string P_i (e.g., X₀Z₁Y₂I₃) must be measured independently because quantum hardware natively measures only in the Z basis: • To measure an X operator on a qubit: apply a Hadamard gate before measurement (rotates X-eigenbasis into Z-eigenbasis) • To measure a Y operator: apply S†·H before measurement • The circuit is executed, sampled, and the parity of the relevant qubit outcomes gives one sample of ⟨P_i⟩; averaging over shots estimates the expectation value with standard error ∝ 1/√(shots)

Measurement grouping (critical efficiency technique): • Naively, each of the ~15 (H₂) to thousands (larger molecules) Pauli terms would require its own circuit — the single largest source of wall-clock cost in VQE • Qubit-wise commuting (QWC) grouping: Pauli strings that commute on every qubit individually can be measured simultaneously in one circuit execution, since they share a common eigenbasis rotation • Typical QWC grouping reduces H₂O-sized problems from ~600 Pauli terms down to ~50-100 distinct measurement circuits • More advanced techniques (general commuting groups, classical shadows — Huang, Kueng, Preskill 2020) push this further, especially valuable as molecule size grows

Error sources and mitigation on real NISQ hardware: • Gate errors: two-qubit CNOT/CZ gate fidelity ~99.0-99.9% on IBM/Google superconducting devices (2023-2024 generation) — still accumulates significantly over ansatz depths of tens of gates • Readout error: ~1-3% probability of misreading qubit state; correctable via measurement calibration matrices • Decoherence: T1/T2 relaxation and dephasing times (~100-300 μs on superconducting qubits) limit total usable circuit depth • Zero-noise extrapolation (ZNE): deliberately amplify noise (via gate folding) at several levels, then extrapolate back to the zero-noise limit • Probabilistic error cancellation and Clifford data regression are further mitigation strategies used in published VQE hardware demonstrations

Even with perfect grouping and mitigation, shot noise sets a fundamental floor: reaching chemical accuracy (1.6 mHa, ~1 kcal/mol) via pure sampling for a molecule with hundreds of Pauli terms can require an estimated 10⁶–10⁸ total shots per single energy evaluation — a major open scaling bottleneck for VQE on pharmaceutically relevant molecules.

Closing the Loop — Classical Optimizers Driving VQE to the Ground State

VQE is a hybrid algorithm precisely because the expensive, exponentially-hard part (representing and sampling a many-body quantum state) runs on the quantum processor, while the comparatively cheap part (updating a modest number of real-valued parameters to reduce a scalar cost function) runs on a classical computer. The variational principle guarantees that any measured energy is an upper bound on the true ground-state energy, so the classical optimizer's entire job is to push that upper bound as low as possible.

• COBYLA, SPSA, L-BFGS-B: Common optimizers (gradient-free & gradient-based)
• exact: Parameter-shift gradient (no finite-difference approximation)
• ~1.137 Ha: H₂ convergence (vs. exact FCI −1.1373 Ha)
• 1.6 mHa: Chemical accuracy target (≈1 kcal/mol)

Optimizer choice, the parameter-shift rule, and convergence versus classical benchmarks

The variational principle underlying VQE (Peruzzo et al., Nature Communications 2014 — the original VQE paper):

⟨ψ(θ)|H|ψ(θ)⟩ ≥ E_ground_state for any trial state |ψ(θ)⟩

So minimizing the measured expectation value over θ, by construction, approaches (but never undershoots) the true ground-state energy — a rigorous guarantee classical heuristic optimization methods lack.

Optimizer choices in practice:

1. COBYLA (Constrained Optimization BY Linear Approximation): gradient-free, builds a local linear model from function evaluations; robust to the shot noise inherent in quantum measurements; the default choice in many Qiskit VQE tutorials for small problems

2. SPSA (Simultaneous Perturbation Stochastic Approximation): estimates a gradient direction using only 2 function evaluations regardless of parameter count (versus 2N for finite differences) — well suited to noisy, expensive quantum cost function evaluations and the standard choice for hardware-efficient ansätze with dozens of parameters

3. Parameter-shift rule (Mitarai et al. 2018; Schuld et al. 2019): for gates of the form e^{−iθP/2}, the exact analytic gradient can be obtained as: ∂⟨H⟩/∂θ = [⟨H⟩(θ+π/2) − ⟨H⟩(θ−π/2)] / 2 • Unlike classical finite-difference gradients, this is exact (not an approximation) even in the presence of shot noise in each term • Enables gradient-based optimizers (L-BFGS-B, Adam) at the cost of 2 extra circuit evaluations per parameter per iteration

Convergence behavior: • Each iteration: prepare ansatz with current θ → measure all Pauli groups → compute ⟨H⟩ → optimizer proposes new θ → repeat • Typical H₂ convergence: 30-100 iterations to reach chemical accuracy (1.6 mHa, ≈1 kcal/mol) from a random parameter start • Barren plateaus (vanishing gradients at scale) and noisy cost landscapes can cause optimizers to stall far from the true minimum — an active research area (natural gradients, layer-wise training, warm-starting from classical CCSD amplitudes)

Benchmark against classical quantum chemistry (why VQE has not yet demonstrated advantage): • Full Configuration Interaction (FCI) — exact within a basis set — gives H₂/STO-3G ground energy of −1.1373 Ha; a well-converged VQE run reproduces this to within mHa on both real and simulated hardware • CCSD(T), the classical "gold standard" for single-reference systems, matches or exceeds FCI-level accuracy for molecules up to hundreds of atoms at polynomial classical cost — VQE offers no accuracy advantage for these small, single-reference molecules • DFT (density functional theory) handles thousands of atoms routinely at a fraction of the cost of any correlated wavefunction method, classical or quantum • VQE's hoped-for advantage lies specifically in strongly correlated, multi-reference systems (e.g., transition-metal catalysts, Fe-S clusters in nitrogenase, bond-breaking regions) where classical methods scale exponentially or fail qualitatively — a regime current NISQ hardware (dozens of noisy qubits) cannot yet reach at useful accuracy

As of 2024, no VQE calculation on real quantum hardware has out-performed classical CCSD(T) or DFT on any chemically or pharmaceutically relevant molecule in either accuracy or cost. The largest hardware VQE demonstrations remain in the 10-20 qubit range (e.g., Google's 2020 Hartree-Fock simulation of diazene on 12 qubits) — genuine quantum advantage for molecular energies is expected to require fault-tolerant, error-corrected qubits, likely thousands of physical qubits per logical qubit, still years away.