claude-scientific-skills/scientific-skills/pennylane/references/quantum_chemistry.md

Quantum Chemistry with PennyLane

Table of Contents

1. Molecular Hamiltonians
2. Variational Quantum Eigensolver (VQE)
3. Molecular Structure
4. Basis Sets and Mapping
5. Excited States
6. Quantum Chemistry Workflows

Molecular Hamiltonians

Building Molecular Hamiltonians

import pennylane as qml
from pennylane import qchem
import numpy as np

# Define molecule
symbols = ['H', 'H']
coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.74])  # Angstroms

# Generate Hamiltonian
hamiltonian, n_qubits = qchem.molecular_hamiltonian(
    symbols,
    coordinates,
    charge=0,
    mult=1,  # Spin multiplicity
    basis='sto-3g',
    method='dhf'  # Dirac-Hartree-Fock
)

print(f"Hamiltonian: {hamiltonian}")
print(f"Number of qubits needed: {n_qubits}")

Jordan-Wigner Transformation

# Hamiltonian is automatically in qubit form via Jordan-Wigner
# Manual transformation:
from pennylane import fermi

# Fermionic operators
a_0 = fermi.FermiC(0)  # Creation operator
a_1 = fermi.FermiA(1)  # Annihilation operator

# Convert to qubits
qubit_op = qml.qchem.jordan_wigner(a_0 * a_1)

Bravyi-Kitaev Transformation

# Alternative mapping (more efficient for some systems)
from pennylane.qchem import bravyi_kitaev

# Build Hamiltonian with Bravyi-Kitaev
hamiltonian, n_qubits = qchem.molecular_hamiltonian(
    symbols,
    coordinates,
    mapping='bravyi_kitaev'
)

Custom Hamiltonians

# Build Hamiltonian from coefficients and operators
coeffs = [0.2, -0.8, 0.5]
obs = [
    qml.PauliZ(0),
    qml.PauliZ(0) @ qml.PauliZ(1),
    qml.PauliX(0) @ qml.PauliX(1)
]

H = qml.Hamiltonian(coeffs, obs)

# Or use simplified syntax
H = 0.2 * qml.PauliZ(0) - 0.8 * qml.PauliZ(0) @ qml.PauliZ(1) + 0.5 * qml.PauliX(0) @ qml.PauliX(1)

Variational Quantum Eigensolver (VQE)

Basic VQE Implementation

from pennylane import numpy as np

# Define device
dev = qml.device('default.qubit', wires=n_qubits)

# Hartree-Fock state preparation
hf_state = qchem.hf_state(electrons=2, orbitals=n_qubits)

def ansatz(params, wires):
    """Variational ansatz."""
    qml.BasisState(hf_state, wires=wires)

    for i in range(len(wires)):
        qml.RY(params[i], wires=i)

    for i in range(len(wires)-1):
        qml.CNOT(wires=[i, i+1])

@qml.qnode(dev)
def vqe_circuit(params):
    ansatz(params, wires=range(n_qubits))
    return qml.expval(hamiltonian)

# Optimize
opt = qml.GradientDescentOptimizer(stepsize=0.4)
params = np.random.normal(0, np.pi, n_qubits, requires_grad=True)

for n in range(100):
    params, energy = opt.step_and_cost(vqe_circuit, params)

    if n % 20 == 0:
        print(f"Step {n}: Energy = {energy:.8f} Ha")

print(f"Final ground state energy: {energy:.8f} Ha")

UCCSD Ansatz

from pennylane.qchem import UCCSD

# Singles and doubles excitations
singles, doubles = qchem.excitations(electrons=2, orbitals=n_qubits)

@qml.qnode(dev)
def uccsd_circuit(params):
    # Hartree-Fock reference
    qml.BasisState(hf_state, wires=range(n_qubits))

    # UCCSD ansatz
    UCCSD(params, wires=range(n_qubits), s_wires=singles, d_wires=doubles)

    return qml.expval(hamiltonian)

# Initialize parameters
n_params = len(singles) + len(doubles)
params = np.zeros(n_params, requires_grad=True)

# Optimize
opt = qml.AdamOptimizer(stepsize=0.1)
for n in range(100):
    params, energy = opt.step_and_cost(uccsd_circuit, params)

Adaptive VQE

def adaptive_vqe(hamiltonian, n_qubits, max_gates=10):
    """Adaptive VQE: Grow ansatz iteratively."""
    dev = qml.device('default.qubit', wires=n_qubits)

    # Start with HF state
    operations = []
    params = []

    hf_state = qchem.hf_state(electrons=2, orbitals=n_qubits)

    @qml.qnode(dev)
    def circuit(p):
        qml.BasisState(hf_state, wires=range(n_qubits))

        for op, param in zip(operations, p):
            op(param)

        return qml.expval(hamiltonian)

    # Iteratively add gates
    for _ in range(max_gates):
        # Find best gate to add
        best_op = None
        best_improvement = 0

        for candidate_op in generate_candidates():
            # Test adding this operation
            test_ops = operations + [candidate_op]
            test_params = params + [0.0]

            improvement = evaluate_improvement(test_ops, test_params)

            if improvement > best_improvement:
                best_improvement = improvement
                best_op = candidate_op

        if best_improvement < threshold:
            break

        operations.append(best_op)
        params.append(0.0)

        # Optimize current ansatz
        opt = qml.AdamOptimizer(stepsize=0.1)
        for _ in range(50):
            params = opt.step(circuit, params)

    return circuit, params

Molecular Structure

Defining Molecules

# Simple diatomic
h2_symbols = ['H', 'H']
h2_coords = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.74])

# Water molecule
h2o_symbols = ['O', 'H', 'H']
h2o_coords = np.array([
    0.0, 0.0, 0.0,      # O
    0.757, 0.586, 0.0,  # H
   -0.757, 0.586, 0.0   # H
])

# From XYZ format
molecule = qchem.read_structure('molecule.xyz')
symbols, coords = molecule

Geometry Optimization

def optimize_geometry(symbols, initial_coords, basis='sto-3g'):
    """Optimize molecular geometry."""

    def energy_surface(coords):
        H, n_qubits = qchem.molecular_hamiltonian(
            symbols, coords, basis=basis
        )

        # Run VQE to get energy
        energy = run_vqe(H, n_qubits)
        return energy

    # Classical optimization of nuclear coordinates
    from scipy.optimize import minimize

    result = minimize(
        energy_surface,
        initial_coords,
        method='BFGS',
        options={'gtol': 1e-5}
    )

    return result.x, result.fun

optimized_coords, min_energy = optimize_geometry(h2_symbols, h2_coords)
print(f"Optimized geometry: {optimized_coords}")
print(f"Energy: {min_energy} Ha")

Bond Dissociation Curves

def dissociation_curve(symbols, axis=2, distances=None):
    """Calculate potential energy surface."""

    if distances is None:
        distances = np.linspace(0.5, 3.0, 20)

    energies = []

    for d in distances:
        coords = np.zeros(6)
        coords[axis] = d  # Set bond length

        H, n_qubits = qchem.molecular_hamiltonian(
            symbols, coords, basis='sto-3g'
        )

        energy = run_vqe(H, n_qubits)
        energies.append(energy)

        print(f"Distance: {d:.2f} Å, Energy: {energy:.6f} Ha")

    return distances, energies

# H2 dissociation
distances, energies = dissociation_curve(['H', 'H'])

import matplotlib.pyplot as plt
plt.plot(distances, energies)
plt.xlabel('Bond length (Å)')
plt.ylabel('Energy (Ha)')
plt.title('H2 Dissociation Curve')
plt.show()

Basis Sets and Mapping

Basis Set Selection

# Minimal basis (fastest, least accurate)
H_sto3g, n_qubits = qchem.molecular_hamiltonian(
    symbols, coords, basis='sto-3g'
)

# Double-zeta basis
H_631g, n_qubits = qchem.molecular_hamiltonian(
    symbols, coords, basis='6-31g'
)

# Large basis (slower, more accurate)
H_ccpvdz, n_qubits = qchem.molecular_hamiltonian(
    symbols, coords, basis='cc-pvdz'
)

Active Space Selection

# Select active orbitals
active_electrons = 2
active_orbitals = 2

H_active, n_qubits = qchem.molecular_hamiltonian(
    symbols,
    coords,
    active_electrons=active_electrons,
    active_orbitals=active_orbitals
)

print(f"Full system: {len(symbols)} electrons")
print(f"Active space: {active_electrons} electrons in {active_orbitals} orbitals")
print(f"Qubits needed: {n_qubits}")

Fermion-to-Qubit Mappings

# Jordan-Wigner (default)
H_jw, n_q_jw = qchem.molecular_hamiltonian(
    symbols, coords, mapping='jordan_wigner'
)

# Bravyi-Kitaev
H_bk, n_q_bk = qchem.molecular_hamiltonian(
    symbols, coords, mapping='bravyi_kitaev'
)

# Parity
H_parity, n_q_parity = qchem.molecular_hamiltonian(
    symbols, coords, mapping='parity'
)

print(f"Jordan-Wigner terms: {len(H_jw.ops)}")
print(f"Bravyi-Kitaev terms: {len(H_bk.ops)}")

Excited States

Quantum Subspace Expansion

def quantum_subspace_expansion(hamiltonian, ground_state_params, excitations):
    """Calculate excited states via subspace expansion."""

    @qml.qnode(dev)
    def ground_state():
        ansatz(ground_state_params, wires=range(n_qubits))
        return qml.state()

    # Get ground state
    psi_0 = ground_state()

    # Generate excited state basis
    basis = [psi_0]

    for exc in excitations:
        @qml.qnode(dev)
        def excited_state():
            ansatz(ground_state_params, wires=range(n_qubits))
            # Apply excitation
            apply_excitation(exc)
            return qml.state()

        psi_exc = excited_state()
        basis.append(psi_exc)

    # Build Hamiltonian matrix in subspace
    n_basis = len(basis)
    H_matrix = np.zeros((n_basis, n_basis))

    for i in range(n_basis):
        for j in range(n_basis):
            H_matrix[i, j] = np.vdot(basis[i], hamiltonian @ basis[j])

    # Diagonalize
    eigenvalues, eigenvectors = np.linalg.eigh(H_matrix)

    return eigenvalues, eigenvectors

SSVQE (Subspace-Search VQE)

def ssvqe(hamiltonian, n_states=3):
    """Calculate multiple states simultaneously."""

    def cost_function(params):
        states = []

        for i in range(n_states):
            @qml.qnode(dev)
            def state_i():
                ansatz(params[i], wires=range(n_qubits))
                return qml.state()

            states.append(state_i())

        # Energy expectation
        energies = [np.vdot(s, hamiltonian @ s) for s in states]

        # Orthogonality penalty
        penalty = 0
        for i in range(n_states):
            for j in range(i+1, n_states):
                overlap = np.abs(np.vdot(states[i], states[j]))
                penalty += overlap ** 2

        return sum(energies) + 1000 * penalty

    # Initialize parameters for all states
    params = [np.random.random(n_params) for _ in range(n_states)]

    opt = qml.AdamOptimizer(stepsize=0.01)
    for _ in range(100):
        params = opt.step(cost_function, params)

    return params

Quantum Chemistry Workflows

Full VQE Workflow

def full_chemistry_workflow(symbols, coords, basis='sto-3g'):
    """Complete quantum chemistry calculation."""

    print("1. Building molecular Hamiltonian...")
    H, n_qubits = qchem.molecular_hamiltonian(
        symbols, coords, basis=basis
    )

    print(f"   Molecule: {' '.join(symbols)}")
    print(f"   Qubits: {n_qubits}")
    print(f"   Hamiltonian terms: {len(H.ops)}")

    print("\n2. Preparing Hartree-Fock state...")
    n_electrons = sum(qchem.atomic_numbers[s] for s in symbols)
    hf_state = qchem.hf_state(n_electrons, n_qubits)

    print("\n3. Running VQE...")
    energy, params = run_vqe(H, n_qubits, hf_state)

    print(f"\n4. Results:")
    print(f"   Ground state energy: {energy:.8f} Ha")

    print("\n5. Computing properties...")
    dipole = compute_dipole_moment(symbols, coords, params)
    print(f"   Dipole moment: {dipole:.4f} D")

    return {
        'energy': energy,
        'params': params,
        'dipole': dipole
    }

results = full_chemistry_workflow(['H', 'H'], h2_coords)

Molecular Property Calculation

def compute_molecular_properties(symbols, coords, vqe_params):
    """Calculate molecular properties from VQE solution."""

    # Energy
    H, n_qubits = qchem.molecular_hamiltonian(symbols, coords)
    energy = vqe_circuit(vqe_params)

    # Dipole moment
    dipole_obs = qchem.dipole_moment(symbols, coords)

    @qml.qnode(dev)
    def dipole_circuit(axis):
        ansatz(vqe_params, wires=range(n_qubits))
        return qml.expval(dipole_obs[axis])

    dipole = [dipole_circuit(i) for i in range(3)]
    dipole_magnitude = np.linalg.norm(dipole)

    # Particle number (sanity check)
    @qml.qnode(dev)
    def particle_number():
        ansatz(vqe_params, wires=range(n_qubits))
        N_op = qchem.particle_number(n_qubits)
        return qml.expval(N_op)

    n_particles = particle_number()

    return {
        'energy': energy,
        'dipole_moment': dipole_magnitude,
        'dipole_vector': dipole,
        'particle_number': n_particles
    }

Reaction Energy Calculation

def reaction_energy(reactants, products):
    """Calculate energy of chemical reaction."""

    # Calculate energies of reactants
    E_reactants = 0
    for molecule in reactants:
        symbols, coords = molecule
        H, n_qubits = qchem.molecular_hamiltonian(symbols, coords)
        E_reactants += run_vqe(H, n_qubits)

    # Calculate energies of products
    E_products = 0
    for molecule in products:
        symbols, coords = molecule
        H, n_qubits = qchem.molecular_hamiltonian(symbols, coords)
        E_products += run_vqe(H, n_qubits)

    # Reaction energy
    delta_E = E_products - E_reactants

    print(f"Reactant energy: {E_reactants:.6f} Ha")
    print(f"Product energy: {E_products:.6f} Ha")
    print(f"Reaction energy: {delta_E:.6f} Ha ({delta_E * 627.5:.2f} kcal/mol)")

    return delta_E

# Example: H2 dissociation
reactants = [((['H', 'H'], h2_coords_bonded))]
products = [((['H'], [0, 0, 0]), (['H'], [10, 0, 0]))]  # Separated atoms

delta_E = reaction_energy(reactants, products)

Best Practices

1. Start with small basis sets - Use STO-3G for testing, upgrade for production
2. Use active space - Reduce qubits by selecting relevant orbitals
3. Choose appropriate mapping - Bravyi-Kitaev often reduces circuit depth
4. Initialize with HF - Start VQE from Hartree-Fock state
5. Validate results - Compare with classical methods (FCI, CCSD)
6. Consider symmetries - Exploit molecular symmetries to reduce complexity
7. Use UCCSD for accuracy - UCCSD ansatz is chemically motivated
8. Monitor convergence - Check gradient norms and energy variance
9. Account for correlation - Ensure ansatz captures electron correlation
10. Benchmark thoroughly - Test on known systems before novel molecules