Sample IBM Qiskit Code for the Potential Use Case
Sample IBM Qiskit Code for the Potential Use Case
Embedded Files
This module implements a complete VQE baseline for molecular
energy calculations on IBM quantum hardware.
Requirements:
pip install qiskit qiskit-ibm-runtime qiskit-nature
pyscf scipy numpy matplotlib
Usage:
python milestone2_vqe.py --molecule water --backend ibm_boston
# ── Standard Library ────────────────────────────────────────
import os
import json
import logging
import argparse
from datetime import datetime
from dataclasses import dataclass, field
from typing import Optional
# ── Numerical / Scientific ───────────────────────────────────
import numpy as np
from scipy.optimize import minimize
# ── Qiskit Core ──────────────────────────────────────────────
from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.circuit.library import TwoLocal, EfficientSU2
# ── Qiskit IBM Runtime ───────────────────────────────────────
from qiskit_ibm_runtime import (
QiskitRuntimeService,
Session,
Estimator,
Options,
)
# ── Qiskit Nature (Electronic Structure) ────────────────────
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import (
JordanWignerMapper,
ParityMapper,
BravyiKitaevMapper,
)
from qiskit_nature.second_q.circuit.library import (
HartreeFock,
UCCSD,
)
from qiskit_nature.second_q.transformers import (
FreezeCoreTransformer,
ActiveSpaceTransformer,
)
from qiskit_nature.second_q.algorithms import (
GroundStateEigensolver,
VQEUCCFactory,
)
# ── PySCF (Classical Reference) ──────────────────────────────
from pyscf import gto, scf, fci
# ── Visualization ────────────────────────────────────────────
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# ── Logging Setup ────────────────────────────────────────────
logging.basicConfig(
level=logging.INFO,
format="%(asctime)s | %(levelname)s | %(message)s",
handlers=[
logging.StreamHandler(),
logging.FileHandler(
f"vqe_run_{datetime.now().strftime('%Y%m%d_%H%M%S')}.log"
),
],
)
log = logging.getLogger(__name__)
# ════════════════════════════════════════════════════════════
# 1. CONFIGURATION
# ════════════════════════════════════════════════════════════
@dataclass
class MoleculeConfig:
"""
Molecular system configuration.
Defines the molecule, basis set, and active space.
"""
name: str
geometry: str # XYZ string or PySCF format
basis: str = "6-31g*"
charge: int = 0
spin: int = 0 # 2S (number of unpaired electrons)
freeze_core: bool = True
active_electrons: Optional[int] = None
active_orbitals: Optional[int] = None
def __post_init__(self):
log.info(f"MoleculeConfig created: {self.name} | "
f"basis={self.basis} | charge={self.charge}")
@dataclass
class VQEConfig:
"""
VQE hyperparameters and runtime settings.
"""
mapper: str = "jordan_wigner" # jordan_wigner | parity | bravyi_kitaev
ansatz_type: str = "uccsd" # uccsd | two_local | efficient_su2
optimizer: str = "COBYLA"
max_iterations: int = 1000
convergence_tol: float = 1e-6
shots: int = 8192
# TwoLocal / EfficientSU2 specific
reps: int = 3
entanglement: str = "linear" # linear | full | circular
use_simulator: bool = True # False = real IBM hardware
backend_name: str = "ibm_boston"
ibm_token: str = "" # Set via env var IBM_TOKEN
results_dir: str = "results/"
def __post_init__(self):
os.makedirs(self.results_dir, exist_ok=True)
if not self.ibm_token:
self.ibm_token = os.environ.get("IBM_TOKEN", "")
log.info(f"VQEConfig: ansatz={self.ansatz_type} | "
f"optimizer={self.optimizer} | shots={self.shots}")
# ── Predefined Molecules ─────────────────────────────────────
MOLECULES = {
"water": MoleculeConfig(
name = "water",
geometry = "O 0.0 0.0 0.0; H 0.757 0.586 0.0; H -0.757 0.586 0.0",
basis = "6-31g", # 14 qubits after freeze-core + JW
charge = 0,
spin = 0,
freeze_core = True,
),
"hydrogen": MoleculeConfig(
name = "hydrogen",
geometry = "H 0.0 0.0 0.0; H 0.0 0.0 0.735",
basis = "sto-3g", # 4 qubits — good sanity check
charge = 0,
spin = 0,
freeze_core = False,
),
"lih": MoleculeConfig(
name = "lih",
geometry = "Li 0.0 0.0 0.0; H 0.0 0.0 1.6",
basis = "sto-3g",
charge = 0,
spin = 0,
freeze_core = True,
),
}
# ════════════════════════════════════════════════════════════
# 2. IBM BACKEND CONNECTION
# ════════════════════════════════════════════════════════════
class IBMBackendManager:
"""
Manages connection to IBM Quantum hardware or simulator.
Falls back gracefully to local Aer simulator if no token.
"""
def __init__(self, cfg: VQEConfig):
self.cfg = cfg
self.backend = None
self.service = None
self._connect()
def _connect(self):
if self.cfg.use_simulator:
self._use_local_simulator()
else:
self._use_ibm_hardware()
def _use_local_simulator(self):
"""Use Qiskit Aer statevector simulator (no token needed)."""
try:
from qiskit_aer import AerSimulator
self.backend = AerSimulator(method="statevector")
log.info("Connected to local Aer statevector simulator.")
except ImportError:
log.warning(
"qiskit-aer not installed. "
"Run: pip install qiskit-aer"
)
raise
def _use_ibm_hardware(self):
"""Connect to real IBM Quantum hardware."""
if not self.cfg.ibm_token:
raise ValueError(
"IBM_TOKEN environment variable not set. "
"Export your IBM Quantum token:\n"
" export IBM_TOKEN='your_token_here'"
)
try:
self.service = QiskitRuntimeService(
channel = "ibm_quantum",
token = self.cfg.ibm_token,
)
self.backend = self.service.backend(
self.cfg.backend_name
)
f"Connected to IBM backend: {self.cfg.backend_name} | "
f"Qubits: {self.backend.num_qubits}"
)
except Exception as e:
log.error(f"IBM connection failed: {e}")
log.info("Falling back to local simulator.")
self._use_local_simulator()
def get_backend(self):
return self.backend
def get_backend_info(self) -> dict:
info = {"backend_name": str(self.backend)}
if hasattr(self.backend, "num_qubits"):
info["num_qubits"] = self.backend.num_qubits
if hasattr(self.backend, "properties"):
try:
props = self.backend.properties()
if props:
gates = props.gates
cx_errors = [
g.parameters[0].value
for g in gates
if g.gate == "cx" and g.parameters
]
if cx_errors:
info["avg_cx_error"] = float(np.mean(cx_errors))
except Exception:
pass
return info
# ════════════════════════════════════════════════════════════
# 3. HAMILTONIAN GENERATION
# ════════════════════════════════════════════════════════════
class HamiltonianBuilder:
"""
Generates the qubit Hamiltonian from molecular geometry
using PySCF + Qiskit Nature.
Pipeline:
Molecule geometry
→ PySCF RHF
→ Qiskit Nature ElectronicStructureProblem
→ Freeze core / active space reduction
→ Fermion → Qubit mapping (Jordan-Wigner etc.)
→ SparsePauliOp Hamiltonian
"""
MAPPERS = {
"jordan_wigner": JordanWignerMapper,
"parity": ParityMapper,
"bravyi_kitaev": BravyiKitaevMapper,
}
def __init__(self, mol_cfg: MoleculeConfig, vqe_cfg: VQEConfig):
self.mol_cfg = mol_cfg
self.vqe_cfg = vqe_cfg
self.problem = None
self.mapper = None
self.hamiltonian = None
self.num_qubits = None
self.num_particles = None
self.num_spatial_orbs = None
self.nuclear_repulsion = None
def build(self) -> SparsePauliOp:
log.info(f"Building Hamiltonian for {self.mol_cfg.name} ...")
# ── PySCF Driver ─────────────────────────────────────
driver = PySCFDriver(
atom = self.mol_cfg.geometry,
basis = self.mol_cfg.basis,
charge = self.mol_cfg.charge,
spin = self.mol_cfg.spin,
)
self.problem = driver.run()
# ── Transformers (freeze core / active space) ────────
transformers = []
if self.mol_cfg.freeze_core:
transformers.append(FreezeCoreTransformer())
if (self.mol_cfg.active_electrons is not None
and self.mol_cfg.active_orbitals is not None):
transformers.append(
ActiveSpaceTransformer(
self.mol_cfg.active_electrons,
self.mol_cfg.active_orbitals,
)
)
for t in transformers:
self.problem = t.transform(self.problem)
# ── Mapper ───────────────────────────────────────────
mapper_cls = self.MAPPERS[self.vqe_cfg.mapper]
self.mapper = mapper_cls()
# ── Second-quantized → Qubit Hamiltonian ─────────────
second_q_ops = self.problem.second_q_ops()
fermionic_op = second_q_ops[0]
self.hamiltonian = self.mapper.map(fermionic_op)
# ── Metadata ─────────────────────────────────────────
self.num_qubits = self.hamiltonian.num_qubits
ep = self.problem.properties.electronic_integrals
self.num_particles = self.problem.num_particles
self.num_spatial_orbs = self.problem.num_spatial_orbitals
self.nuclear_repulsion = (
self.problem.nuclear_repulsion_energy
)
f"Hamiltonian built | "
f"Qubits: {self.num_qubits} | "
f"Particles: {self.num_particles} | "
f"Spatial orbitals: {self.num_spatial_orbs} | "
f"Pauli terms: {len(self.hamiltonian)}"
)
return self.hamiltonian
def get_hf_energy(self) -> float:
"""Returns classical Hartree-Fock reference energy."""
if self.problem is None:
self.build()
return self.problem.reference_energy
def get_fci_energy(self) -> Optional[float]:
"""
Computes exact FCI energy classically via PySCF.
Only feasible for small systems (H2, LiH, water/sto-3g).
"""
try:
mol = gto.Mole()
mol.atom = self.mol_cfg.geometry
mol.basis = self.mol_cfg.basis
mol.charge = self.mol_cfg.charge
mol.spin = self.mol_cfg.spin
mol.build()
mf = scf.RHF(mol).run(verbose=0)
fci_solver = fci.FCI(mf)
fci_energy, _ = fci_solver.kernel()
log.info(f"FCI energy: {fci_energy:.8f} Ha")
return fci_energy
except Exception as e:
log.warning(f"FCI calculation failed: {e}")
return None
# ════════════════════════════════════════════════════════════
# 4. ANSATZ CONSTRUCTION
# ════════════════════════════════════════════════════════════
class AnsatzBuilder:
"""
Constructs the parameterized quantum circuit (ansatz).
Options:
uccsd — Physically motivated, compact, hardware-intensive
two_local — Hardware-efficient, less physically motivated
efficient_su2 — Hardware-efficient with SU(2) rotations
"""
def __init__(
self,
hamiltonian_builder: HamiltonianBuilder,
vqe_cfg: VQEConfig,
):
self.hb = hamiltonian_builder
self.cfg = vqe_cfg
self.ansatz = None
self.hf_circuit = None
def build(self) -> QuantumCircuit:
log.info(f"Building ansatz: {self.cfg.ansatz_type}")
if self.cfg.ansatz_type == "uccsd":
return self._build_uccsd()
elif self.cfg.ansatz_type == "two_local":
return self._build_two_local()
elif self.cfg.ansatz_type == "efficient_su2":
return self._build_efficient_su2()
else:
raise ValueError(
f"Unknown ansatz type: {self.cfg.ansatz_type}. "
"Choose: uccsd | two_local | efficient_su2"
)
def _build_uccsd(self) -> QuantumCircuit:
"""
Unitary Coupled Cluster Singles and Doubles ansatz.
Physically motivated — best accuracy for small molecules.
More circuit depth than hardware-efficient alternatives.
"""
hf_state = HartreeFock(
num_spatial_orbitals = self.hb.num_spatial_orbs,
num_particles = self.hb.num_particles,
qubit_mapper = self.hb.mapper,
)
uccsd = UCCSD(
num_spatial_orbitals = self.hb.num_spatial_orbs,
num_particles = self.hb.num_particles,
qubit_mapper = self.hb.mapper,
initial_state = hf_state,
)
self.hf_circuit = hf_state
self.ansatz = uccsd
f"UCCSD ansatz | "
f"Parameters: {uccsd.num_parameters} | "
f"Depth: {uccsd.decompose().depth()}"
)
return uccsd
def _build_two_local(self) -> QuantumCircuit:
"""
Hardware-efficient TwoLocal ansatz.
Lower circuit depth — better for NISQ hardware.
"""
n = self.hb.num_qubits
ansatz = TwoLocal(
num_qubits = n,
rotation_blocks = ["ry", "rz"],
entanglement_blocks = "cx",
entanglement = self.cfg.entanglement,
reps = self.cfg.reps,
insert_barriers = True,
)
# Prepend Hartree-Fock initial state
hf = HartreeFock(
num_spatial_orbitals = self.hb.num_spatial_orbs,
num_particles = self.hb.num_particles,
qubit_mapper = self.hb.mapper,
)
self.hf_circuit = hf
self.ansatz = hf.compose(ansatz)
f"TwoLocal ansatz | "
f"Parameters: {ansatz.num_parameters} | "
f"Depth: {ansatz.decompose().depth()}"
)
return self.ansatz
def _build_efficient_su2(self) -> QuantumCircuit:
"""
EfficientSU2 — hardware-efficient with SU(2) rotations.
Good balance between expressibility and circuit depth.
"""
n = self.hb.num_qubits
ansatz = EfficientSU2(
num_qubits = n,
entanglement = self.cfg.entanglement,
reps = self.cfg.reps,
)
hf = HartreeFock(
num_spatial_orbitals = self.hb.num_spatial_orbs,
num_particles = self.hb.num_particles,
qubit_mapper = self.hb.mapper,
)
self.hf_circuit = hf
self.ansatz = hf.compose(ansatz)
f"EfficientSU2 ansatz | "
f"Parameters: {ansatz.num_parameters} | "
f"Depth: {ansatz.decompose().depth()}"
)
return self.ansatz
def get_initial_parameters(self) -> np.ndarray:
"""
Smart initialization:
- UCCSD: zeros (known good starting point)
- Others: small random perturbation around zero
to break symmetry and avoid barren plateaus
"""
n_params = self.ansatz.num_parameters
if self.cfg.ansatz_type == "uccsd":
params = np.zeros(n_params)
else:
rng = np.random.default_rng(seed=42)
params = rng.uniform(-0.1, 0.1, n_params)
log.info(f"Initialized {n_params} parameters.")
return params
# ════════════════════════════════════════════════════════════
# 5. VQE RUNNER
# ════════════════════════════════════════════════════════════
class VQERunner:
"""
Executes the VQE optimization loop.
Supports:
- Local Aer simulator (statevector — noiseless)
- Real IBM quantum hardware via Qiskit Runtime
- Multiple classical optimizers (COBYLA, L-BFGS-B, SLSQP)
- Energy history tracking and convergence monitoring
"""
def __init__(
self,
hamiltonian: SparsePauliOp,
ansatz: QuantumCircuit,
backend,
cfg: VQEConfig,
nuclear_repulsion: float = 0.0,
):
self.hamiltonian = hamiltonian
self.ansatz = ansatz
self.backend = backend
self.cfg = cfg
self.nuclear_repulsion = nuclear_repulsion
self.energy_history: list = []
self.param_history: list = []
self.iteration_count: int = 0
self.optimal_params: Optional[np.ndarray] = None
self.optimal_energy: Optional[float] = None
def _evaluate_energy(self, params: np.ndarray) -> float:
"""
Core energy evaluation function.
Called by the classical optimizer at each iteration.
"""
self.iteration_count += 1
try:
if self.cfg.use_simulator:
energy = self._statevector_energy(params)
else:
energy = self._hardware_energy(params)
# Add nuclear repulsion for total energy
total_energy = energy + self.nuclear_repulsion
self.energy_history.append(total_energy)
self.param_history.append(params.copy())
if self.iteration_count % 10 == 0:
f"Iteration {self.iteration_count:4d} | "
f"Energy: {total_energy:.8f} Ha"
)
return total_energy
except Exception as e:
log.error(f"Energy evaluation failed: {e}")
return 1e10 # Return large value to signal failure
def _statevector_energy(self, params: np.ndarray) -> float:
"""
Fast exact energy via statevector simulation.
No shot noise — used for development and validation.
"""
bound_circuit = self.ansatz.assign_parameters(params)
sv = Statevector(bound_circuit)
energy = sv.expectation_value(self.hamiltonian).real
return float(energy)
def _hardware_energy(self, params: np.ndarray) -> float:
"""
Energy from real IBM hardware via Qiskit Runtime Estimator.
Includes transpilation to backend native gate set.
"""
options = Options()
options.execution.shots = self.cfg.shots
options.optimization_level = 1
options.resilience_level = 1 # Basic error mitigation
transpiled = transpile(
self.ansatz,
backend = self.backend,
optimization_level = 1,
)
with Session(backend=self.backend) as session:
estimator = Estimator(session=session, options=options)
job = estimator.run(
[(transpiled, self.hamiltonian, params)]
)
result = job.result()
energy = result[0].data.evs
return float(energy)
def run(
self,
initial_params: np.ndarray,
) -> dict:
"""
Main VQE optimization loop.
Returns full results dictionary.
"""
f"Starting VQE | "
f"Optimizer: {self.cfg.optimizer} | "
f"Max iterations: {self.cfg.max_iterations}"
)
start_time = datetime.now()
# ── Run optimization ─────────────────────────────────
result = minimize(
fun = self._evaluate_energy,
x0 = initial_params,
method = self.cfg.optimizer,
options = {
"maxiter": self.cfg.max_iterations,
"xatol": self.cfg.convergence_tol,
"fatol": self.cfg.convergence_tol,
"rhobeg": 0.1, # COBYLA initial step size
},
)
elapsed = (datetime.now() - start_time).total_seconds()
self.optimal_params = result.x
self.optimal_energy = result.fun
vqe_result = {
"optimal_energy": self.optimal_energy,
"optimal_params": self.optimal_params.tolist(),
"energy_history": self.energy_history,
"iterations": self.iteration_count,
"converged": result.success,
"optimizer_message": result.message,
"elapsed_seconds": elapsed,
"nuclear_repulsion": self.nuclear_repulsion,
}
f"VQE complete | "
f"Energy: {self.optimal_energy:.8f} Ha | "
f"Iterations: {self.iteration_count} | "
f"Converged: {result.success} | "
f"Time: {elapsed:.1f}s"
)
return vqe_result
# ════════════════════════════════════════════════════════════
# 6. RESULTS ANALYZER
# ════════════════════════════════════════════════════════════
class ResultsAnalyzer:
"""
Analyzes and benchmarks VQE results against
classical reference methods (HF, FCI).
Generates plots and saves results to disk.
"""
def __init__(
self,
vqe_result: dict,
mol_cfg: MoleculeConfig,
vqe_cfg: VQEConfig,
hf_energy: Optional[float] = None,
fci_energy: Optional[float] = None,
):
self.result = vqe_result
self.mol_cfg = mol_cfg
self.vqe_cfg = vqe_cfg
self.hf_energy = hf_energy
self.fci_energy = fci_energy
def print_summary(self):
"""Prints a formatted benchmark summary table."""
vqe_e = self.result["optimal_energy"]
print("\n" + "═" * 55)
print(f" VQE RESULTS — {self.mol_cfg.name.upper()}")
print("═" * 55)
print(f" Basis set : {self.mol_cfg.basis}")
print(f" Ansatz : {self.vqe_cfg.ansatz_type}")
print(f" Mapper : {self.vqe_cfg.mapper}")
print(f" Optimizer : {self.vqe_cfg.optimizer}")
print(f" Iterations : {self.result['iterations']}")
print(f" Converged : {self.result['converged']}")
print(f" Wall time : {self.result['elapsed_seconds']:.1f}s")
print("─" * 55)
print(f" VQE energy : {vqe_e:>16.8f} Ha")
if self.hf_energy is not None:
hf_err = abs(vqe_e - self.hf_energy) * 1000
corr = (self.hf_energy - vqe_e) * 1000
print(f" HF energy : {self.hf_energy:>16.8f} Ha")
print(f" Correlation E : {corr:>+14.4f} mHa")
if self.fci_energy is not None:
fci_err = abs(vqe_e - self.fci_energy) * 1000
chem = 1.594 # mHa = 1 kcal/mol
print(f" FCI energy : {self.fci_energy:>16.8f} Ha")
print(f" VQE error vs FCI : {fci_err:>14.4f} mHa")
status = "✓ Chemical accuracy" if fci_err < chem else "✗ Above chemical accuracy"
print(f" Status : {status}")
print("═" * 55 + "\n")
def validate_convergence(self) -> bool:
"""
Checks whether VQE has converged to chemical accuracy
(< 1.594 mHa = 1 kcal/mol) relative to FCI.
"""
if self.fci_energy is None:
log.warning("FCI energy not available for validation.")
return False
err = abs(self.result["optimal_energy"] - self.fci_energy)
chem_accuracy = 1.594e-3 # Hartree
passed = err < chem_accuracy
f"Validation: error={err*1000:.4f} mHa | "
f"Chemical accuracy: {'PASSED' if passed else 'FAILED'}"
)
return passed
def plot_convergence(self, save: bool = True):
"""
Plots VQE energy convergence history with
HF and FCI reference lines.
"""
history = self.result["energy_history"]
if not history:
log.warning("No energy history to plot.")
return
fig = plt.figure(figsize=(12, 8))
gs = gridspec.GridSpec(2, 2, figure=fig)
ax1 = fig.add_subplot(gs[0, :])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[1, 1])
iterations = range(1, len(history) + 1)
# ── Main convergence plot ─────────────────────────────
ax1.plot(
iterations, history,
"b-o", ms=3, lw=1.5,
label="VQE energy",
)
if self.hf_energy:
ax1.axhline(
self.hf_energy, color="orange",
ls="--", lw=1.5, label=f"HF ({self.hf_energy:.5f} Ha)"
)
if self.fci_energy:
ax1.axhline(
self.fci_energy, color="green",
ls="--", lw=1.5, label=f"FCI ({self.fci_energy:.5f} Ha)"
)
ax1.set_xlabel("VQE Iteration")
ax1.set_ylabel("Energy (Hartree)")
ax1.set_title(
f"VQE Convergence — {self.mol_cfg.name.capitalize()} "
f"({self.vqe_cfg.ansatz_type.upper()})"
)
ax1.legend()
ax1.grid(True, alpha=0.3)
# ── Error vs FCI ─────────────────────────────────────
if self.fci_energy:
errors = [abs(e - self.fci_energy) * 1000 for e in history]
ax2.semilogy(iterations, errors, "r-", lw=1.5)
ax2.axhline(
1.594, color="green", ls="--",
label="Chemical accuracy (1 kcal/mol)"
)
ax2.set_xlabel("Iteration")
ax2.set_ylabel("|VQE − FCI| (mHa)")
ax2.set_title("Error vs FCI (log scale)")
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)
# ── Energy distribution (last 20%) ───────────────────
tail = history[int(0.8 * len(history)):]
ax3.hist(tail, bins=15, color="steelblue", edgecolor="white")
ax3.axvline(
self.result["optimal_energy"],
color="red", ls="--",
label=f"Optimal: {self.result['optimal_energy']:.6f} Ha"
)
ax3.set_xlabel("Energy (Hartree)")
ax3.set_ylabel("Count")
ax3.set_title("Energy Distribution (final 20% of iterations)")
ax3.legend(fontsize=8)
plt.tight_layout()
if save:
fname = os.path.join(
self.vqe_cfg.results_dir,
f"vqe_convergence_{self.mol_cfg.name}_"
f"{self.vqe_cfg.ansatz_type}.png"
)
plt.savefig(fname, dpi=150, bbox_inches="tight")
log.info(f"Convergence plot saved: {fname}")
plt.show()
def save_results(self):
"""Saves full results to JSON for reproducibility."""
output = {
"timestamp": datetime.now().isoformat(),
"molecule": self.mol_cfg.name,
"basis": self.mol_cfg.basis,
"ansatz": self.vqe_cfg.ansatz_type,
"mapper": self.vqe_cfg.mapper,
"optimizer": self.vqe_cfg.optimizer,
"vqe_energy": self.result["optimal_energy"],
"hf_energy": self.hf_energy,
"fci_energy": self.fci_energy,
"iterations": self.result["iterations"],
"converged": self.result["converged"],
"elapsed_s": self.result["elapsed_seconds"],
"energy_history": self.result["energy_history"],
}
if self.fci_energy:
output["error_vs_fci_mHa"] = (
abs(self.result["optimal_energy"] - self.fci_energy) * 1000
)
fname = os.path.join(
self.vqe_cfg.results_dir,
f"vqe_results_{self.mol_cfg.name}_"
f"{self.vqe_cfg.ansatz_type}.json"
)
with open(fname, "w") as f:
json.dump(output, f, indent=2)
log.info(f"Results saved: {fname}")
return fname
# ════════════════════════════════════════════════════════════
# 7. BOND LENGTH SCAN
# ════════════════════════════════════════════════════════════
class BondLengthScanner:
"""
Scans VQE energy across a range of bond lengths to
generate a potential energy surface (PES).
Replicates the water PES shown in Figure 6 of the
QUICK proposal.
"""
def __init__(
self,
mol_cfg: MoleculeConfig,
vqe_cfg: VQEConfig,
backend,
):
self.mol_cfg = mol_cfg
self.vqe_cfg = vqe_cfg
self.backend = backend
def scan_oh_bond(
self,
bond_lengths: np.ndarray,
) -> dict:
"""
Scans O-H bond length in water and computes
VQE energy at each geometry.
"""
f"Starting PES scan over {len(bond_lengths)} "
f"geometries ..."
)
vqe_energies = []
hf_energies = []
fci_energies = []
# Reuse parameters from previous geometry (warm start)
prev_params = None
for idx, r in enumerate(bond_lengths):
log.info(f"Geometry {idx+1}/{len(bond_lengths)} | O-H = {r:.3f} Å")
# Update geometry
cfg = MoleculeConfig(
name = f"water_r{r:.3f}",
geometry = (
f"O 0.0 0.0 0.0; "
f"H {r:.4f} 0.0 0.0; "
f"H {-r*0.5:.4f} {r*0.866:.4f} 0.0"
),
basis = self.mol_cfg.basis,
freeze_core = self.mol_cfg.freeze_core,
)
try:
hb = HamiltonianBuilder(cfg, self.vqe_cfg)
H = hb.build()
hf_e = hb.get_hf_energy()
fci_e = hb.get_fci_energy()
ab = AnsatzBuilder(hb, self.vqe_cfg)
ansatz = ab.build()
params = (prev_params
if prev_params is not None
and len(prev_params) == ansatz.num_parameters
else ab.get_initial_parameters())
runner = VQERunner(
H, ansatz, self.backend,
self.vqe_cfg, hb.nuclear_repulsion
)
result = runner.run(params)
prev_params = result["optimal_params"]
vqe_energies.append(result["optimal_energy"])
hf_energies.append(hf_e)
if fci_e:
fci_energies.append(fci_e)
except Exception as e:
log.error(f"Scan failed at r={r:.3f}: {e}")
vqe_energies.append(None)
hf_energies.append(None)
fci_energies.append(None)
return {
"bond_lengths": bond_lengths.tolist(),
"vqe_energies": vqe_energies,
"hf_energies": hf_energies,
"fci_energies": fci_energies,
}
def plot_pes(self, scan_result: dict, save: bool = True):
"""Plots the potential energy surface."""
r = np.array(scan_result["bond_lengths"])
vqe = np.array([e for e in scan_result["vqe_energies"] if e])
hf = np.array([e for e in scan_result["hf_energies"] if e])
fci = np.array([e for e in scan_result["fci_energies"] if e])
plt.figure(figsize=(8, 5))
plt.plot(r[:len(vqe)], vqe, "b-o", ms=5, label="VQE")
if len(hf):
plt.plot(r[:len(hf)], hf, "orange", ls="--", label="HF")
if len(fci):
plt.plot(r[:len(fci)], fci, "g-", label="FCI (exact)")
plt.xlabel("O-H Bond Length (Å)")
plt.ylabel("Energy (Hartree)")
plt.title(
f"Potential Energy Surface — {self.mol_cfg.name.capitalize()}\n"
f"Basis: {self.mol_cfg.basis} | Ansatz: {self.vqe_cfg.ansatz_type}"
)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
if save:
fname = os.path.join(
self.vqe_cfg.results_dir,
f"pes_{self.mol_cfg.name}.png"
)
plt.savefig(fname, dpi=150, bbox_inches="tight")
log.info(f"PES plot saved: {fname}")
plt.show()
# ════════════════════════════════════════════════════════════
# 8. MAIN PIPELINE ORCHESTRATOR
# ════════════════════════════════════════════════════════════
def run_vqe_pipeline(
molecule_name: str = "water",
ansatz_type: str = "uccsd",
mapper: str = "jordan_wigner",
use_simulator: bool = True,
backend_name: str = "ibm_boston",
run_pes_scan: bool = False,
) -> dict:
"""
Full Milestone 2 VQE pipeline.
Steps:
1. Configure molecule and VQE settings
2. Connect to IBM backend or local simulator
3. Build qubit Hamiltonian
4. Construct ansatz circuit
5. Run VQE optimization
6. Benchmark against HF and FCI
7. Plot convergence and save results
8. (Optional) Run potential energy surface scan
"""
log.info("=" * 60)
log.info("QUICK Pipeline — Milestone 2: VQE Baseline")
log.info("=" * 60)
# ── 1. Configuration ─────────────────────────────────────
mol_cfg = MOLECULES.get(molecule_name)
if mol_cfg is None:
raise ValueError(
f"Unknown molecule: '{molecule_name}'. "
f"Choose from: {list(MOLECULES.keys())}"
)
vqe_cfg = VQEConfig(
ansatz_type = ansatz_type,
mapper = mapper,
optimizer = "COBYLA",
max_iterations = 1000,
shots = 8192,
use_simulator = use_simulator,
backend_name = backend_name,
results_dir = "results/milestone2/",
)
# ── 2. Backend connection ────────────────────────────────
backend_mgr = IBMBackendManager(vqe_cfg)
backend = backend_mgr.get_backend()
backend_info = backend_mgr.get_backend_info()
log.info(f"Backend info: {backend_info}")
# ── 3. Hamiltonian ───────────────────────────────────────
hb = HamiltonianBuilder(mol_cfg, vqe_cfg)
H = hb.build()
hf_energy = hb.get_hf_energy()
fci_energy = hb.get_fci_energy()
# ── 4. Ansatz ────────────────────────────────────────────
ab = AnsatzBuilder(hb, vqe_cfg)
ansatz = ab.build()
params = ab.get_initial_parameters()
f"Circuit | Qubits: {ansatz.num_qubits} | "
f"Depth: {ansatz.decompose().depth()} | "
f"Parameters: {ansatz.num_parameters}"
)
# ── 5. VQE optimization ──────────────────────────────────
runner = VQERunner(
hamiltonian = H,
ansatz = ansatz,
backend = backend,
cfg = vqe_cfg,
nuclear_repulsion = hb.nuclear_repulsion,
)
result = runner.run(params)
# ── 6. Analysis and benchmarking ─────────────────────────
analyzer = ResultsAnalyzer(
vqe_result = result,
mol_cfg = mol_cfg,
vqe_cfg = vqe_cfg,
hf_energy = hf_energy,
fci_energy = fci_energy,
)
analyzer.print_summary()
analyzer.validate_convergence()
analyzer.plot_convergence(save=True)
analyzer.save_results()
# ── 7. Optional PES scan ─────────────────────────────────
if run_pes_scan and molecule_name == "water":
log.info("Running potential energy surface scan ...")
bond_lengths = np.linspace(0.7, 2.5, 12)
scanner = BondLengthScanner(mol_cfg, vqe_cfg, backend)
scan_result = scanner.scan_oh_bond(bond_lengths)
scanner.plot_pes(scan_result)
log.info("Milestone 2 complete.")
return result
# ════════════════════════════════════════════════════════════
# 9. COMMAND LINE INTERFACE
# ════════════════════════════════════════════════════════════
def parse_args():
parser = argparse.ArgumentParser(
description="QUICK Pipeline — Milestone 2: VQE Baseline"
)
parser.add_argument(
"--molecule",
default = "water",
choices = list(MOLECULES.keys()),
help = "Molecule to simulate (default: water)"
)
parser.add_argument(
"--ansatz",
default = "uccsd",
choices = ["uccsd", "two_local", "efficient_su2"],
help = "Ansatz type (default: uccsd)"
)
parser.add_argument(
"--mapper",
default = "jordan_wigner",
choices = ["jordan_wigner", "parity", "bravyi_kitaev"],
help = "Qubit mapper (default: jordan_wigner)"
)
parser.add_argument(
"--backend",
default = "ibm_boston",
help = "IBM backend name (default: ibm_boston)"
)
parser.add_argument(
"--hardware",
action = "store_true",
help = "Use real IBM hardware (requires IBM_TOKEN env var)"
)
parser.add_argument(
"--pes-scan",
action = "store_true",
help = "Run potential energy surface scan"
)
return parser.parse_args()
if __name__ == "__main__":
args = parse_args()
run_vqe_pipeline(
molecule_name = args.molecule,
ansatz_type = args.ansatz,
mapper = args.mapper,
use_simulator = not args.hardware,
backend_name = args.backend,
run_pes_scan = args.pes_scan,
)
Report abuse