VQE for Quantum Chemistry

Variational Quantum Eigensolver for molecular energy calculations.

Overview

VQE (Variational Quantum Eigensolver) is a hybrid quantum-classical algorithm for finding the ground state energy of molecules. It's one of the most promising near-term applications of quantum computers.

This example demonstrates:

• Molecular Hamiltonian encoding
• Variational circuit design
• Classical optimization loop
• Energy convergence analysis

The Algorithm

┌─────────────────────────────────────────────────────────────┐
│                    VQE Algorithm                             │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│   ┌──────────────────┐                                      │
│   │  Ansatz Circuit  │                                      │
│   │  U(θ)            │──────┐                               │
│   └──────────────────┘      │                               │
│           ▲                 ▼                               │
│           │         ┌───────────────┐                       │
│   ┌───────┴───────┐ │   Measure     │                       │
│   │   Classical   │ │   ⟨H⟩         │                       │
│   │   Optimizer   │ └───────────────┘                       │
│   │   (COBYLA)    │         │                               │
│   └───────────────┘         │                               │
│           ▲                 │                               │
│           └─────────────────┘                               │
│                 E(θ)                                        │
│                                                             │
└─────────────────────────────────────────────────────────────┘

Complete Example: H₂ Molecule

"""
VQE for H2 Molecule Ground State Energy
========================================

Calculates the ground state energy of molecular hydrogen (H2)
using the Variational Quantum Eigensolver algorithm.

Expected result: ~-1.137 Hartree (near equilibrium bond length)
"""

import numpy as np
from scipy.optimize import minimize
from softqcos_sdk import QCOSClient, AsyncQCOSClient
import time

# ============================================================================
# MOLECULAR DATA
# ============================================================================

# H2 Hamiltonian in STO-3G basis (simplified, 2 qubits)
# H = c0*I + c1*Z0 + c2*Z1 + c3*Z0Z1 + c4*X0X1 + c5*Y0Y1
H2_COEFFICIENTS = {
    "I": -0.4804,
    "Z0": 0.3435,
    "Z1": -0.4347,
    "Z0Z1": 0.5716,
    "X0X1": 0.0910,
    "Y0Y1": 0.0910,
}

# Known exact ground state energy for validation
H2_EXACT_ENERGY = -1.1373  # Hartree


def create_ansatz_circuit(theta: list[float]) -> str:
    """
    Create parameterized ansatz circuit for H2 VQE.

    Uses RY-CNOT-RY ansatz (hardware efficient):

         ┌───────┐     ┌───────┐
    q0: ─┤ Ry(θ0)├──■──┤ Ry(θ2)├─
         └───────┘┌─┴─┐└───────┘
    q1: ─┤ Ry(θ1)├┤ X ├┤ Ry(θ3)├─
         └───────┘└───┘└───────┘
    """
    t0, t1, t2, t3 = theta

    circuit = f"""
OPENQASM 2.0;
include "qelib1.inc";

qreg q[2];
creg c[2];

// Initial rotation layer
ry({t0}) q[0];
ry({t1}) q[1];

// Entangling layer
cx q[0], q[1];

// Final rotation layer
ry({t2}) q[0];
ry({t3}) q[1];
"""
    return circuit


def measure_pauli_expectation(client: QCOSClient, base_circuit: str,
                               pauli_string: str, shots: int = 2048) -> float:
    """
    Measure expectation value of a Pauli string.

    Args:
        client: QCOSClient instance
        base_circuit: Base ansatz circuit (without measurements)
        pauli_string: e.g., "Z0", "Z0Z1", "X0X1"
        shots: Number of measurement shots

    Returns:
        Expectation value in [-1, 1]
    """
    # Add measurement basis rotation if needed
    circuit = base_circuit.strip()

    # Remove any existing measurements
    circuit_lines = [l for l in circuit.split('\n') if 'measure' not in l.lower()]

    # Add basis rotation for X and Y measurements
    if "X" in pauli_string:
        # X basis: apply H before measurement
        for i, char in enumerate(pauli_string):
            if char == "X":
                qubit = int(pauli_string[i+1])
                circuit_lines.append(f"h q[{qubit}];")

    if "Y" in pauli_string:
        # Y basis: apply S†H before measurement
        for i, char in enumerate(pauli_string):
            if char == "Y":
                qubit = int(pauli_string[i+1])
                circuit_lines.append(f"sdg q[{qubit}];")
                circuit_lines.append(f"h q[{qubit}];")

    # Add measurement
    circuit_lines.append("measure q -> c;")

    measurement_circuit = '\n'.join(circuit_lines)

    # Execute
    result = client.execute(qasm=measurement_circuit, shots=shots)

    # Calculate expectation value
    expectation = 0.0
    total = sum(result.counts.values())

    for bitstring, count in result.counts.items():
        # Calculate parity based on Pauli string
        parity = 1

        # Extract qubit indices from Pauli string
        for i, char in enumerate(pauli_string):
            if char in "XYZ":
                qubit_idx = int(pauli_string[i+1])
                # Bitstring is reversed (q0 is rightmost)
                bit = int(bitstring[-(qubit_idx+1)])
                if bit == 1:
                    parity *= -1

        expectation += parity * count / total

    return expectation


def calculate_energy(theta: list[float], client: QCOSClient,
                     shots: int = 2048, verbose: bool = False) -> float:
    """
    Calculate the expectation value of H2 Hamiltonian.

    E(θ) = Σ cᵢ ⟨ψ(θ)|Pᵢ|ψ(θ)⟩
    """
    base_circuit = create_ansatz_circuit(theta)

    energy = H2_COEFFICIENTS["I"]  # Identity term

    # Measure each Pauli term
    pauli_terms = ["Z0", "Z1", "Z0Z1", "X0X1", "Y0Y1"]

    for pauli in pauli_terms:
        expectation = measure_pauli_expectation(client, base_circuit, pauli, shots)
        energy += H2_COEFFICIENTS[pauli] * expectation

        if verbose:
            print(f"   ⟨{pauli}⟩ = {expectation:+.4f}")

    return energy


def run_vqe():
    """Run complete VQE optimization for H2."""

    print("⚛️  VQE for H₂ Molecule")
    print("=" * 60)

    client = QCOSClient()

    # Initial parameters (random)
    np.random.seed(42)
    theta_init = np.random.uniform(-np.pi, np.pi, 4)

    print(f"\n📝 Configuration:")
    print(f"   Qubits: 2")
    print(f"   Parameters: 4")
    print(f"   Ansatz: RY-CNOT-RY")
    print(f"   Optimizer: COBYLA")

    # Track optimization history
    history = []
    iteration = [0]

    def objective(theta):
        energy = calculate_energy(theta, client, shots=1024)
        history.append(energy)
        iteration[0] += 1

        if iteration[0] % 5 == 0:
            print(f"   Iteration {iteration[0]:3d}: E = {energy:.6f} Ha")

        return energy

    print(f"\n⏳ Starting optimization...")
    print(f"   Initial energy: {objective(theta_init):.6f} Ha")
    print(f"   Exact energy:   {H2_EXACT_ENERGY:.6f} Ha")
    print()

    start_time = time.time()

    # Optimize
    result = minimize(
        objective,
        theta_init,
        method='COBYLA',
        options={'maxiter': 50, 'rhobeg': 0.5}
    )

    elapsed = time.time() - start_time

    # Final high-precision measurement
    final_energy = calculate_energy(result.x, client, shots=4096)

    print(f"\n✅ Optimization Complete")
    print(f"   Total iterations: {iteration[0]}")
    print(f"   Total time: {elapsed:.1f}s")

    print(f"\n📊 Results:")
    print(f"   VQE Energy:   {final_energy:.6f} Ha")
    print(f"   Exact Energy: {H2_EXACT_ENERGY:.6f} Ha")
    print(f"   Error:        {abs(final_energy - H2_EXACT_ENERGY):.6f} Ha")
    print(f"   Accuracy:     {abs(final_energy - H2_EXACT_ENERGY) * 627.5:.2f} kcal/mol")

    # Chemical accuracy is ~1 kcal/mol
    if abs(final_energy - H2_EXACT_ENERGY) * 627.5 < 1.0:
        print("\n🎉 Chemical accuracy achieved!")
    elif abs(final_energy - H2_EXACT_ENERGY) * 627.5 < 5.0:
        print("\n✅ Good accuracy (within 5 kcal/mol)")
    else:
        print("\n⚠️  Consider more iterations or different ansatz")

    print(f"\n📈 Optimal Parameters:")
    for i, t in enumerate(result.x):
        print(f"   θ{i} = {t:.4f} rad ({np.degrees(t):.1f}°)")

    return {
        "energy": final_energy,
        "exact_energy": H2_EXACT_ENERGY,
        "parameters": result.x,
        "history": history,
        "iterations": iteration[0],
    }


def main():
    result = run_vqe()

    # Plot convergence (if matplotlib available)
    try:
        import matplotlib.pyplot as plt

        plt.figure(figsize=(10, 6))
        plt.plot(result['history'], 'b-', linewidth=2, label='VQE Energy')
        plt.axhline(y=result['exact_energy'], color='r', linestyle='--',
                   label=f'Exact ({result["exact_energy"]:.4f} Ha)')
        plt.xlabel('Iteration')
        plt.ylabel('Energy (Hartree)')
        plt.title('VQE Convergence for H₂')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.savefig('vqe_convergence.png', dpi=150)
        print("\n📈 Convergence plot saved to: vqe_convergence.png")
    except ImportError:
        pass

    return result


if __name__ == "__main__":
    main()

Expected Output

⚛️  VQE for H₂ Molecule
============================================================

📝 Configuration:
   Qubits: 2
   Parameters: 4
   Ansatz: RY-CNOT-RY
   Optimizer: COBYLA

⏳ Starting optimization...
   Initial energy: -0.234521 Ha
   Exact energy:   -1.137300 Ha

   Iteration   5: E = -0.892341 Ha
   Iteration  10: E = -1.054123 Ha
   Iteration  15: E = -1.098456 Ha
   Iteration  20: E = -1.125789 Ha
   Iteration  25: E = -1.134567 Ha

✅ Optimization Complete
   Total iterations: 28
   Total time: 45.2s

📊 Results:
   VQE Energy:   -1.135842 Ha
   Exact Energy: -1.137300 Ha
   Error:        0.001458 Ha
   Accuracy:     0.91 kcal/mol

🎉 Chemical accuracy achieved!

📈 Optimal Parameters:
   θ0 = 0.0523 rad (3.0°)
   θ1 = 3.1416 rad (180.0°)
   θ2 = 0.0000 rad (0.0°)
   θ3 = 0.0000 rad (0.0°)

Larger Molecules: LiH

"""VQE for Lithium Hydride (LiH) - 4 qubits"""

# LiH Hamiltonian coefficients (simplified)
LIH_COEFFICIENTS = {
    "I": -7.4983,
    "Z0": 0.2252,
    "Z1": 0.2252,
    "Z2": -0.4815,
    "Z3": -0.4815,
    "Z0Z1": 0.1209,
    "Z0Z2": 0.1653,
    "Z0Z3": 0.1740,
    "Z1Z2": 0.1740,
    "Z1Z3": 0.1653,
    "Z2Z3": 0.1761,
    "X0X1": 0.0453,
    "X2X3": 0.0453,
}

LIH_EXACT_ENERGY = -7.8823  # Hartree


def create_lih_ansatz(theta: list[float]) -> str:
    """4-qubit hardware-efficient ansatz for LiH."""

    circuit = f"""
OPENQASM 2.0;
include "qelib1.inc";

qreg q[4];
creg c[4];

// Layer 1: Single qubit rotations
ry({theta[0]}) q[0];
ry({theta[1]}) q[1];
ry({theta[2]}) q[2];
ry({theta[3]}) q[3];

// Layer 1: Entangling
cx q[0], q[1];
cx q[2], q[3];

// Layer 2: Single qubit rotations
ry({theta[4]}) q[0];
ry({theta[5]}) q[1];
ry({theta[6]}) q[2];
ry({theta[7]}) q[3];

// Layer 2: Entangling
cx q[1], q[2];

// Layer 3: Final rotations
ry({theta[8]}) q[0];
ry({theta[9]}) q[1];
ry({theta[10]}) q[2];
ry({theta[11]}) q[3];
"""
    return circuit

Optimization with SoftQCOS

"""Use SoftQCOS circuit optimization for VQE."""

from softqcos_sdk import QCOSClient

client = QCOSClient()

# Create ansatz
theta = [0.1, 0.2, 0.3, 0.4]
circuit = create_ansatz_circuit(theta)

# Optimize circuit for hardware
optimized = client.optimize(
    circuit=circuit,
    target="ionq",
    objective="max_fidelity"
)

print(f"Original gates: {optimized.original_gate_count}")
print(f"Optimized gates: {optimized.optimized_gate_count}")
print(f"Gate reduction: {optimized.improvements['gate_reduction']}")

# Use optimized circuit in VQE loop
# This can significantly improve results on noisy hardware

Best Practices

1. Shot Budgeting

# More shots for final energy, fewer during optimization
def smart_energy_calculation(theta, iteration, max_iterations):
    if iteration < max_iterations * 0.7:
        shots = 512  # Fast exploration
    elif iteration < max_iterations * 0.9:
        shots = 1024  # Refinement
    else:
        shots = 4096  # Final precision

    return calculate_energy(theta, client, shots=shots)

2. Gradient-based Optimization

# Use parameter-shift rule for gradients
def compute_gradient(theta, client, shift=np.pi/2):
    gradient = np.zeros_like(theta)

    for i in range(len(theta)):
        theta_plus = theta.copy()
        theta_minus = theta.copy()
        theta_plus[i] += shift
        theta_minus[i] -= shift

        E_plus = calculate_energy(theta_plus, client)
        E_minus = calculate_energy(theta_minus, client)

        gradient[i] = (E_plus - E_minus) / (2 * np.sin(shift))

    return gradient

3. Error Mitigation

# Zero-noise extrapolation
def mitigated_energy(theta, client, noise_scales=[1.0, 1.5, 2.0]):
    energies = []

    for scale in noise_scales:
        # Execute with amplified noise (circuit folding)
        energy = calculate_energy(theta, client, noise_scale=scale)
        energies.append(energy)

    # Extrapolate to zero noise
    coeffs = np.polyfit(noise_scales, energies, deg=1)
    mitigated = coeffs[1]  # y-intercept

    return mitigated

Next Steps

• QAOA Optimization - Solve combinatorial problems
• Batch Processing - Run VQE efficiently at scale
• Azure Quantum - Execute VQE on real quantum hardware