Pennylance simulation for CU-graphene.py

# -*- coding: utf-8 -*-
"""Untitled16.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1yxYuwPDvT3d-fyjghv2uzIeUkarhbtEP

## 📦 Environment Setup

We begin by installing all necessary packages for the simulation:

- **PennyLane**: A hybrid quantum-classical platform that supports quantum chemistry workflows.
- **PennyLane-QChem**: Extends PennyLane with quantum chemistry tools, including Hamiltonian generation.
- **Basis Set Exchange (BSE)**: Used to dynamically fetch basis set data (e.g., STO-3G) for atoms such as Cu.

This setup enables the construction of real-space molecular Hamiltonians consistent with the SCF–QSVT methodology described in the paper by Lamane et al. (2023).
"""

!pip install basis-set-exchange --quiet

import pennylane as qml
from pennylane import numpy as np
from pennylane.qchem import Molecule, molecular_hamiltonian
import matplotlib.pyplot as plt
from scipy.linalg import expm
import pandas as pd

"""## Molecular Structure: Cu₂–C₆ Hybrid System

We define a minimal testbed system composed of a copper dimer (Cu₂) coupled to a hexagonal carbon ring (C₆), simulating Cu–graphene interaction.

This configuration serves as a tractable system to explore the self-consistent application of quantum signal processing techniques (e.g., QSVT) to approximate Fermi–Dirac filtering \( f(H) \), as proposed in the reference work.

We instantiate a `Molecule` object with atomic symbols and 3D coordinates. The `load_data=True` option ensures proper basis set retrieval from BSE.

"""

# 1. Define Cu2 + C6 structure
symbols = ['Cu', 'Cu', 'C', 'C', 'C', 'C', 'C', 'C']
coordinates = np.array([
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 2.2],
    [1.4, 0.0, -1.2],
    [0.7, 1.2, -1.2],
    [-0.7, 1.2, -1.2],
    [-1.4, 0.0, -1.2],
    [-0.7, -1.2, -1.2],
    [0.7, -1.2, -1.2],
])

# 2. Create Molecule object
mol = Molecule(
    symbols=symbols,
    coordinates=coordinates,
    charge=0,
    load_data = True
)

"""## Second-Quantized Molecular Hamiltonian

Using the `molecular_hamiltonian()` function, we construct the electronic Hamiltonian of the Cu₂–C₆ system in second quantization.

This step corresponds to the generation of the real-space Hamiltonian \( H[n] \), which is subsequently used in both classical and quantum approximations of the electron density.

We work in an active space of 4 electrons and 4 orbitals for tractability, consistent with the reduced models discussed in Section 3.1 of the paper.

"""

# 3. Build Hamiltonian with active space
H_base, n_qubits = molecular_hamiltonian(
    mol,
    active_electrons=4,
    active_orbitals=4
)

H_base = qml.utils.sparse_hamiltonian(H_base).toarray()
dim = H_base.shape[0]

"""## Chebyshev-Based Approximation of Fermi–Dirac Distribution

To simulate the quantum signal processing steps of QSVT, we construct an approximation to the Fermi–Dirac function \( f(H) = [1 + e^{\beta(H - \mu)}]^{-1} \)
using Chebyshev polynomial expansion:
\[
f(H) \approx \sum_k c_k T_k(H)
\]
This expansion, described in Section 4.2 of the paper, avoids matrix diagonalization and is efficiently implementable on quantum hardware via block-encoding.

In our classical setting, it allows us to mimic the quantum filtering used to extract the electron density \( n(r_j) \) directly from the Hamiltonian matrix.

"""

# 4. Chebyshev tools
def chebyshev_coefficients(f, N, a=-1, b=1):
    k = np.arange(N)
    x = np.cos(np.pi * (k + 0.5) / N)
    x = 0.5 * (x * (b - a) + (b + a))
    fx = f(x)
    c = (2 / N) * np.dot(fx, np.cos(np.outer(k, np.pi * (k + 0.5) / N)))
    c[0] /= 2
    return c

def apply_chebyshev(H_scaled, c):
    T0 = np.eye(H_scaled.shape[0])
    T1 = H_scaled
    result = c[0] * T0 + c[1] * T1
    for k in range(2, len(c)):
        T2 = 2 * H_scaled @ T1 - T0
        result += c[k] * T2
        T0, T1 = T1, T2
    return result

"""## Self-Consistent Field Loop via QSVT-Inspired Filtering

We implement a simplified SCF loop, iteratively updating the Hamiltonian \( H[n] \) based on the output density from the Chebyshev-approximated Fermi operator.

At each iteration:
1. \( f(H) \) is approximated via Chebyshev polynomials (QSVT-style),
2. The diagonal of \( f(H) \) is used as the density \( n^{(k)} \),
3. A simple density-dependent potential is added to form \( H[n^{(k)}] \),
4. Convergence is checked by \( \|n^{(k+1)} - n^{(k)}\| < \varepsilon \)

This mirrors the iterative real-space DFT refinement loop shown in **Figure 2** and **Algorithm 1** of the paper.

"""

# 5. SCF Loop
def scf_loop(H0, beta=5.0, mu=0.0, max_iter=20, tol=1e-4):
    H = H0.copy()
    densities = []
    fermi = lambda x: 1 / (1 + np.exp(beta * (x - mu)))
    λ_max, λ_min = np.max(np.linalg.eigvalsh(H0)), np.min(np.linalg.eigvalsh(H0))
    scale = 2 / (λ_max - λ_min)
    shift = -(λ_max + λ_min) / (λ_max - λ_min)

    for i in range(max_iter):
        H_scaled = scale * H + shift * np.eye(H.shape[0])
        c = chebyshev_coefficients(fermi, N=30)
        fH = apply_chebyshev(H_scaled, c)
        n_diag = np.real(np.diag(fH))
        densities.append(n_diag)

        if i > 0 and np.linalg.norm(densities[-1] - densities[-2]) < tol:
            print(f"SCF converged at iteration {i}")
            break

        V_n = np.diag(n_diag * 0.1)  # Potential from density
        H = H0 + V_n

    return densities[-1], H

n_final, H_final = scf_loop(H_base)

"""## QCBM: Generative Quantum Circuit for Density Sampling

We define a variational quantum circuit known as a Quantum Circuit Born Machine (QCBM). This acts as a quantum generative model to propose candidate electron density vectors.

While Lamane et al. focus on applying QSVT to propagate densities through filtering, future extensions may involve integrating generative models (e.g., VQEs, QCBMs) into geometry or charge density optimization workflows.

This cell produces a candidate density vector \( \tilde{n}_j \sim \text{QCBM}(\theta) \) for comparison against the SCF result.

"""

# 6. QCBM Generator (Quantum Circuit Born Machine)
dev = qml.device("default.qubit", wires=4)

@qml.qnode(dev)
def qcbm_circuit(params):
    for i in range(4):
        qml.RY(params[i], wires=i)
    for i in range(3):
        qml.CNOT(wires=[i, i+1])
    return [qml.expval(qml.PauliZ(i)) for i in range(4)]

def generate_qcbm_density(params):
    raw = qcbm_circuit(params)
    return 0.5 * (1 - np.array(raw))  # Z→[0,1]

params = np.array([0.1, 1.2, 2.0, 0.7], requires_grad=True)
qcbm_density = generate_qcbm_density(params)

"""## Density of States and Conductivity Estimation

We estimate the electrical conductivity using the density of states at the Fermi level, \( \sigma \propto g(E_F) \), consistent with the Kubo–Greenwood formalism.

This approach is justified in the paper (Section 4.4) where the authors note that post-converged densities allow access to spectral features of \( H \), enabling conductivity inference.

We diagonalize the final SCF-updated Hamiltonian and estimate:
\[
\sigma \approx g(E_F) = \left.\frac{dN}{dE}\right|_{E=E_F}
\]

"""

# 7. DOS and conductivity
E_F = 0.0
eigvals = np.linalg.eigvalsh(H_final)
hist, bins = np.histogram(eigvals, bins=100, density=True)
dos_E = 0.5 * (bins[:-1] + bins[1:])
g_EF = np.interp(E_F, dos_E, hist)

"""## Visualization and Quantitative Comparison

We conclude with two visual diagnostics:

- The **density of states (DOS)** near the Fermi level
- A comparison of the **QSVT-SCF-derived electron density** versus the **QCBM-generated density**

These outputs allow us to qualitatively assess the convergence of the SCF loop and the potential of generative quantum models to propose density profiles that are physically meaningful.

Such comparison is essential for extending the framework in Lamane et al. to hybrid SCF–QCBM feedback systems.

"""

# 8. Visualizations
plt.figure(figsize=(8, 4))
plt.plot(dos_E, hist, label='DOS(E)')
plt.axvline(E_F, color='red', linestyle='--', label='Fermi Level')
plt.title("Final DOS after SCF")
plt.xlabel("Energy")
plt.ylabel("DOS(E)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

plt.figure(figsize=(8, 4))
plt.plot(n_final, 'o-', label='SCF Final Density')
plt.plot(np.pad(qcbm_density, (0, len(n_final)-4), constant_values=0), 'x--', label='QCBM Candidate')
plt.title("SCF Electron Density vs QCBM")
plt.xlabel("Orbital index")
plt.ylabel("n(r)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

df = pd.DataFrame({
    "Orbital": np.arange(len(n_final)),
    "SCF Density": np.round(n_final, 4),
    "QCBM Candidate": np.round(np.pad(qcbm_density, (0, len(n_final)-4), constant_values=0), 4)
})
print(df)
print(f"\nEstimated conductivity from g(E_F): σ ≈ {g_EF:.4f} [normalized units]")