condensed_matter.py

#!/usr/bin/env python3
"""
Condensed Matter Physics Simulations
=====================================

This module implements fundamental condensed matter physics calculations:
- Electronic band structure using tight-binding models
- Fermi surface calculations and visualization
- Phonon dispersion relations
- Many-body systems (Hubbard model)

Keywords: band structure calculation python, tight binding model, fermi surface calculation,
phonon dispersion python, hubbard model simulation, electronic properties condensed matter

Author: CoCalc Scientific Templates
License: MIT
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# Set reproducible random seed
np.random.seed(42)

class TightBindingModel:
    """
    Tight-binding model for electronic band structure calculations.

    Supports various lattice geometries and can calculate:
    - Energy dispersion relations E(k)
    - Density of states
    - Fermi surfaces
    """

    def __init__(self, lattice_type='square'):
        """
        Parameters:
        -----------
        lattice_type : str
            Type of lattice ('square', 'triangular', 'honeycomb')
        """
        self.lattice_type = lattice_type

    def dispersion_square_2d(self, kx, ky, t=1.0, epsilon_0=0.0):
        """
        Energy dispersion for 2D square lattice.

        E(k) = ε₀ - 2t(cos(kₓ) + cos(kᵧ))

        Parameters:
        -----------
        kx, ky : array
            k-space coordinates
        t : float
            Hopping parameter
        epsilon_0 : float
            On-site energy
        """
        return epsilon_0 - 2*t*(np.cos(kx) + np.cos(ky))

    def dispersion_graphene(self, kx, ky, t=1.0):
        """
        Energy dispersion for graphene (honeycomb lattice).

        Returns both conduction and valence bands.

        Parameters:
        -----------
        kx, ky : array
            k-space coordinates
        t : float
            Hopping parameter
        """
        # Phase factor for nearest neighbors
        f = 1 + np.exp(1j * kx) + np.exp(1j * (kx/2 + np.sqrt(3)*ky/2))
        energy_magnitude = t * np.abs(f)

        return energy_magnitude, -energy_magnitude  # Conduction, valence bands

    def density_of_states(self, energies, n_bins=200, broadening=0.01):
        """
        Calculate density of states with optional Gaussian broadening.

        Parameters:
        -----------
        energies : array
            Energy values from band calculation
        n_bins : int
            Number of energy bins
        broadening : float
            Gaussian broadening width
        """
        E_flat = energies.flatten()
        E_min, E_max = E_flat.min() - 3*broadening, E_flat.max() + 3*broadening
        E_grid = np.linspace(E_min, E_max, n_bins)

        if broadening > 0:
            # Gaussian broadened DOS
            dos = np.zeros_like(E_grid)
            for E in E_flat:
                dos += np.exp(-(E_grid - E)**2 / (2*broadening**2))
            dos /= (np.sqrt(2*np.pi) * broadening * len(E_flat))
        else:
            # Histogram DOS
            hist, bin_edges = np.histogram(E_flat, bins=E_grid)
            dos = hist / (E_grid[1] - E_grid[0]) / len(E_flat)
            E_grid = (bin_edges[:-1] + bin_edges[1:]) / 2

        return E_grid, dos

    def fermi_surface(self, kx, ky, mu, temperature=0.01):
        """
        Calculate Fermi surface and electron occupation.

        Parameters:
        -----------
        kx, ky : array
            k-space grid
        mu : float
            Chemical potential
        temperature : float
            Temperature for Fermi-Dirac distribution

        Returns:
        --------
        energy : array
            Energy dispersion
        occupation : array
            Fermi-Dirac occupation function
        """
        energy = self.dispersion_square_2d(kx, ky)

        # Fermi-Dirac distribution
        if temperature > 0:
            occupation = 1.0 / (1.0 + np.exp((energy - mu) / temperature))
        else:
            occupation = (energy < mu).astype(float)

        return energy, occupation

class PhononModel:
    """
    Phonon dispersion calculations for various lattice models.

    Implements:
    - 1D monoatomic chain
    - 1D diatomic chain
    - Density of phonon states
    """

    @staticmethod
    def monoatomic_chain_1d(k, K=1.0, M=1.0):
        """
        Phonon frequency for 1D monoatomic chain.

        ω(k) = √(4K/M) |sin(k/2)|

        Parameters:
        -----------
        k : array
            Wave vector
        K : float
            Spring constant
        M : float
            Atomic mass
        """
        return np.sqrt(4*K/M) * np.abs(np.sin(k/2))

    @staticmethod
    def diatomic_chain_1d(k, K=1.0, M1=1.0, M2=2.0):
        """
        Phonon frequencies for 1D diatomic chain.

        Returns both optical and acoustic branches.

        Parameters:
        -----------
        k : array
            Wave vector
        K : float
            Spring constant
        M1, M2 : float
            Masses of the two atoms
        """
        M_sum = M1 + M2
        M_prod = M1 * M2

        # Discriminant for eigenvalue equation
        discriminant = M_sum**2 - 4*M_prod*np.cos(k)**2

        # Optical and acoustic branches
        omega_optical = np.sqrt(K * (M_sum + np.sqrt(discriminant)) / (2*M_prod))
        omega_acoustic = np.sqrt(K * (M_sum - np.sqrt(discriminant)) / (2*M_prod))

        return omega_optical, omega_acoustic

class HubbardModel:
    """
    Mean-field solution of the Hubbard model.

    H = -t Σ⟨i,j⟩σ (c†ᵢσ cⱼσ + h.c.) + U Σᵢ nᵢ↑ nᵢ↓

    Implements self-consistent mean-field approximation to study
    magnetic instabilities and correlation effects.
    """

    def __init__(self, lattice_type='square_2d'):
        """
        Parameters:
        -----------
        lattice_type : str
            Lattice geometry
        """
        self.lattice_type = lattice_type

    def mean_field_bands(self, kx, ky, t, mu, U, m):
        """
        Mean-field energy bands for 2D square lattice.

        Parameters:
        -----------
        kx, ky : array
            k-space coordinates
        t : float
            Hopping parameter
        mu : float
            Chemical potential
        U : float
            On-site interaction
        m : float
            Magnetization (order parameter)
        """
        # Non-interacting dispersion
        epsilon_k = -2*t*(np.cos(kx) + np.cos(ky))

        # Mean-field bands
        E_up = epsilon_k - mu + U*m/2     # Spin-up band
        E_down = epsilon_k - mu - U*m/2   # Spin-down band

        return E_up, E_down

    def self_consistent_solution(self, t, U, n_target, T=0.1, max_iter=1000, tol=1e-6):
        """
        Solve mean-field equations self-consistently.

        Parameters:
        -----------
        t : float
            Hopping parameter
        U : float
            Interaction strength
        n_target : float
            Target electron density
        T : float
            Temperature
        max_iter : int
            Maximum iterations
        tol : float
            Convergence tolerance

        Returns:
        --------
        mu : float
            Converged chemical potential
        m : float
            Converged magnetization
        converged : bool
            Whether solution converged
        """
        # Initial guess
        mu = U * n_target / 2
        m = 0.1  # Small initial magnetization

        # k-space grid
        k_mesh = 64
        kx = np.linspace(-np.pi, np.pi, k_mesh)
        ky = np.linspace(-np.pi, np.pi, k_mesh)
        KX, KY = np.meshgrid(kx, ky)

        for iteration in range(max_iter):
            # Calculate mean-field bands
            E_up, E_down = self.mean_field_bands(KX, KY, t, mu, U, m)

            # Fermi-Dirac distributions
            f_up = 1.0 / (1.0 + np.exp(E_up / T))
            f_down = 1.0 / (1.0 + np.exp(E_down / T))

            # New density and magnetization
            n_new = np.mean(f_up + f_down)
            m_new = np.mean(f_up - f_down)

            # Check convergence
            if abs(m_new - m) < tol:
                return mu, m_new, True

            # Update with mixing for stability
            m = 0.7 * m + 0.3 * m_new

            # Adjust chemical potential to maintain target density
            mu = mu + 0.1 * (n_target - n_new)

        return mu, m, False

def band_structure_comparison():
    """
    Compare band structures of different lattice types.
    """
    print("=== Band Structure Comparison ===")

    tb = TightBindingModel()

    # k-space grid
    k_points = 100
    kx = np.linspace(-np.pi, np.pi, k_points)
    ky = np.linspace(-np.pi, np.pi, k_points)
    KX, KY = np.meshgrid(kx, ky)

    # Calculate dispersions
    E_square = tb.dispersion_square_2d(KX, KY, t=1.0)
    E_graphene_c, E_graphene_v = tb.dispersion_graphene(KX, KY, t=1.0)

    # Calculate density of states
    E_dos_square, dos_square = tb.density_of_states(E_square, broadening=0.05)
    E_dos_graphene, dos_graphene = tb.density_of_states(
        np.concatenate([E_graphene_c.flatten(), E_graphene_v.flatten()]),
        broadening=0.05
    )

    print(f"Square lattice bandwidth: {E_square.max() - E_square.min():.3f}")
    print(f"Graphene bandwidth: {E_graphene_c.max() - E_graphene_v.min():.3f}")

    return {
        'square': {'energy': E_square, 'dos_E': E_dos_square, 'dos': dos_square},
        'graphene': {'energy_c': E_graphene_c, 'energy_v': E_graphene_v,
                    'dos_E': E_dos_graphene, 'dos': dos_graphene},
        'kx': KX, 'ky': KY
    }

def fermi_surface_evolution():
    """
    Study evolution of Fermi surface with chemical potential.
    """
    print("=== Fermi Surface Evolution ===")

    tb = TightBindingModel()

    # High-resolution k-space grid
    k_res = 200
    kx = np.linspace(-np.pi, np.pi, k_res)
    ky = np.linspace(-np.pi, np.pi, k_res)
    KX, KY = np.meshgrid(kx, ky)

    # Different chemical potentials (filling levels)
    mu_values = [-3.0, -2.0, 0.0, 2.0]
    fermi_data = []

    for mu in mu_values:
        energy, occupation = tb.fermi_surface(KX, KY, mu, temperature=0.1)
        fermi_data.append({'mu': mu, 'energy': energy, 'occupation': occupation})

    print(f"Calculated Fermi surfaces for μ = {mu_values}")
    return fermi_data, KX, KY

def phonon_comparison():
    """
    Compare phonon dispersions for different chain types.
    """
    print("=== Phonon Dispersion Comparison ===")

    phonon = PhononModel()

    # k-point path
    k = np.linspace(-np.pi, np.pi, 500)

    # Calculate dispersions
    omega_mono = phonon.monoatomic_chain_1d(k, K=1.0, M=1.0)
    omega_opt, omega_ac = phonon.diatomic_chain_1d(k, K=1.0, M1=1.0, M2=2.0)

    print(f"Monoatomic max frequency: {omega_mono.max():.3f}")
    print(f"Diatomic optical max: {omega_opt.max():.3f}")
    print(f"Diatomic acoustic max: {omega_ac.max():.3f}")

    return {
        'k': k,
        'monoatomic': omega_mono,
        'optical': omega_opt,
        'acoustic': omega_ac
    }

def hubbard_magnetic_transition():
    """
    Study magnetic transition in Hubbard model.
    """
    print("=== Hubbard Model Magnetic Transition ===")

    hubbard = HubbardModel()

    # Parameters
    t = 1.0
    density = 0.8
    U_values = [0.0, 2.0, 4.0, 6.0]

    results = []
    for U in U_values:
        mu, m, converged = hubbard.self_consistent_solution(t, U, density, T=0.1)

        result = {
            'U': U,
            'mu': mu,
            'magnetization': abs(m),  # Take absolute value
            'converged': converged
        }
        results.append(result)

        status = "converged" if converged else "not converged"
        print(f"U/t = {U/t:.1f}: μ = {mu:.3f}, |m| = {abs(m):.4f} ({status})")

    return results

def calculate_transport_properties(band_data):
    """
    Calculate basic transport properties from band structure.

    This is a simplified demonstration - real transport calculations
    require more sophisticated methods.
    """
    print("=== Transport Properties Calculation ===")

    # Extract band data
    energy = band_data['square']['energy']
    kx = band_data['kx']
    ky = band_data['ky']

    # Calculate group velocity (simplified)
    dk = kx[0, 1] - kx[0, 0]  # Grid spacing

    # Numerical derivatives for group velocity
    vx = -np.gradient(energy, dk, axis=1)  # vₓ = -∂E/∂kₓ
    vy = -np.gradient(energy, dk, axis=0)  # vᵧ = -∂E/∂kᵧ

    # Speed
    speed = np.sqrt(vx**2 + vy**2)

    print(f"Maximum group velocity: {speed.max():.3f}")
    print(f"Average group velocity: {np.mean(speed):.3f}")

    return {'vx': vx, 'vy': vy, 'speed': speed}

if __name__ == "__main__":
    # Run all demonstrations
    band_data = band_structure_comparison()
    fermi_data, KX, KY = fermi_surface_evolution()
    phonon_data = phonon_comparison()
    hubbard_results = hubbard_magnetic_transition()
    transport = calculate_transport_properties(band_data)

    print("\n=== Summary ===")
    print(f"Completed all condensed matter physics calculations")
    print(f"Band structure: {band_data['square']['energy'].shape} k-points")
    print(f"Fermi surfaces: {len(fermi_data)} chemical potentials")
    print(f"Phonon modes: {len(phonon_data['k'])} k-points")
    print(f"Hubbard model: {len(hubbard_results)} interaction strengths")