py_default_default.py

# -*- coding: UTF-8 -*-



import os
import sys
import codecs

if '--interactive' not in sys.argv[1:]:
    if sys.version_info[0] == 2:
        sys.stdout = codecs.getwriter('UTF-8')(sys.stdout, 'strict')
        sys.stderr = codecs.getwriter('UTF-8')(sys.stderr, 'strict')
    else:
        sys.stdout = codecs.getwriter('UTF-8')(sys.stdout.buffer, 'strict')
        sys.stderr = codecs.getwriter('UTF-8')(sys.stderr.buffer, 'strict')

if '/usr/share/texlive/texmf-dist/scripts/pythontex' and '/usr/share/texlive/texmf-dist/scripts/pythontex' not in sys.path:
    sys.path.append('/usr/share/texlive/texmf-dist/scripts/pythontex')
from pythontex_utils import PythonTeXUtils
pytex = PythonTeXUtils()

pytex.docdir = os.getcwd()
if os.path.isdir('.'):
    os.chdir('.')
    if os.getcwd() not in sys.path:
        sys.path.append(os.getcwd())
else:
    if len(sys.argv) < 2 or sys.argv[1] != '--manual':
        sys.exit('Cannot find directory .')
if pytex.docdir not in sys.path:
    sys.path.append(pytex.docdir)



pytex.id = 'py_default_default'
pytex.family = 'py'
pytex.session = 'default'
pytex.restart = 'default'

pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '0'
pytex.line = '96'

print('=>PYTHONTEX:STDOUT#0#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#0#code#\n')
pytex.before()

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.linalg import expm
import matplotlib
matplotlib.use('Agg')  # Use non-interactive backend

# Set random seed for reproducibility
np.random.seed(42)

# Parameters for quantum simulation
hbar = 1.0  # Reduced Planck constant (natural units)
m = 1.0     # Particle mass (natural units)
L = 10.0    # Box length
N = 1000    # Number of grid points
dx = L / N  # Grid spacing
dt = 0.001  # Time step

# Create spatial grid
x = np.linspace(0, L, N)

# Define potential: harmonic oscillator + barrier
omega = 1.0  # Oscillator frequency
V = 0.5 * m * omega**2 * (x - L/2)**2  # Harmonic potential
V += 2.0 * np.exp(-((x - L/3)/(L/20))**2)  # Gaussian barrier for tunneling

# Create kinetic energy operator (finite difference)
kinetic_diag = -2 * np.ones(N)
kinetic_off = np.ones(N-1)
T_matrix = -(hbar**2 / (2*m*dx**2)) * (
    diags([kinetic_off, kinetic_diag, kinetic_off], [-1, 0, 1], shape=(N, N))
)

# Hamiltonian matrix
H = T_matrix + diags(V, 0)

# Initial wavefunction: Gaussian wave packet
x0 = L/4  # Initial position
sigma = L/20  # Wave packet width
k0 = 5.0  # Initial momentum
psi_initial = np.exp(-(x - x0)**2 / (2*sigma**2)) * np.exp(1j * k0 * x)
psi_initial = psi_initial / np.sqrt(np.trapezoid(np.abs(psi_initial)**2, x))

# Time evolution operator
U = expm(-1j * H.toarray() * dt / hbar)

# Evolve the wavefunction
times = np.arange(0, 5.0, dt)
psi = psi_initial.copy()

# Store some snapshots
snapshot_times = [0, 1.0, 2.0, 3.0, 4.0]
snapshots = []
snapshot_indices = []

for i, t in enumerate(times):
    if any(abs(t - st) < dt/2 for st in snapshot_times):
        snapshots.append(psi.copy())
        snapshot_indices.append(i)
    psi = U @ psi

print(f"Quantum simulation completed: {len(times)} time steps")
print(f"Stored {len(snapshots)} wavefunction snapshots")
print(f"Conservation check - final norm: {np.trapezoid(np.abs(psi)**2, x):.6f}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '1'
pytex.line = '165'

print('=>PYTHONTEX:STDOUT#1#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#1#code#\n')
pytex.before()

# Create quantum mechanics visualization
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

# Plot potential
ax1.plot(x, V, 'k-', linewidth=2, label='Potential V(x)')
ax1.set_ylabel('Energy (natural units)', fontsize=12)
ax1.set_title('Quantum Tunneling Through Potential Barrier', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot wavefunction snapshots
colors = plt.cm.viridis(np.linspace(0, 1, len(snapshots)))
for i, (psi_snap, color) in enumerate(zip(snapshots, colors)):
    probability = np.abs(psi_snap)**2
    ax2.fill_between(x, 0, probability, alpha=0.6, color=color,
                     label=f't = {snapshot_times[i]:.1f}')

ax2.set_xlabel('Position x (natural units)', fontsize=12)
ax2.set_ylabel('Probability Density |ψ(x,t)|²', fontsize=12)
ax2.set_title('Wavefunction Evolution: Quantum Tunneling Dynamics', fontsize=14, fontweight='bold')
ax2.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/quantum_tunneling.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Quantum tunneling visualization saved to assets/quantum_tunneling.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '2'
pytex.line = '215'

print('=>PYTHONTEX:STDOUT#2#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#2#code#\n')
pytex.before()

from scipy.special import hermite
from scipy.special import factorial

# Quantum harmonic oscillator parameters
omega = 1.0  # Angular frequency
m = 1.0      # Mass
hbar = 1.0   # Reduced Planck constant

# Length scale
x0 = np.sqrt(hbar / (m * omega))

# Create position grid
x_osc = np.linspace(-4*x0, 4*x0, 1000)

# Calculate first 6 energy eigenstates
n_states = 6
wavefunctions = []
energies = []

for n in range(n_states):
    # Energy eigenvalue
    E_n = hbar * omega * (n + 0.5)
    energies.append(E_n)

    # Normalization factor
    norm = 1.0 / np.sqrt(2**n * factorial(n)) * (m*omega/(np.pi*hbar))**(1/4)

    # Hermite polynomial
    H_n = hermite(n)

    # Wavefunction
    xi = x_osc / x0  # Dimensionless coordinate
    psi_n = norm * np.exp(-xi**2 / 2) * H_n(xi)
    wavefunctions.append(psi_n)

print(f"Calculated {n_states} harmonic oscillator eigenstates")
print(f"Energy levels: {[f'{E:.2f}' for E in energies[:4]]}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '3'
pytex.line = '255'

print('=>PYTHONTEX:STDOUT#3#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#3#code#\n')
pytex.before()

# Visualize harmonic oscillator eigenstates
fig, ax = plt.subplots(figsize=(12, 10))

# Plot potential
V_osc = 0.5 * m * omega**2 * x_osc**2
ax.plot(x_osc/x0, V_osc/(hbar*omega), 'k-', linewidth=3, label='Potential V(x)')

# Plot energy levels and wavefunctions
colors = plt.cm.Set1(np.linspace(0, 1, n_states))
for n, (psi, E, color) in enumerate(zip(wavefunctions, energies, colors)):
    # Energy level
    ax.axhline(E/(hbar*omega), color=color, linestyle='--', alpha=0.7)

    # Wavefunction (scaled and shifted)
    psi_scaled = psi * 2 + E/(hbar*omega)  # Scale and shift
    ax.plot(x_osc/x0, psi_scaled, color=color, linewidth=2,
            label=f'n={n}, E={E/(hbar*omega):.1f}ℏω')

ax.set_xlabel('Position x/x₀', fontsize=14)
ax.set_ylabel('Energy / ℏω', fontsize=14)
ax.set_title('Quantum Harmonic Oscillator: Energy Eigenstates and Wavefunctions',
             fontsize=16, fontweight='bold')
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
ax.set_ylim(0, 6)

plt.tight_layout()
plt.savefig('assets/harmonic_oscillator.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Harmonic oscillator visualization saved to assets/harmonic_oscillator.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '4'
pytex.line = '307'

print('=>PYTHONTEX:STDOUT#4#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#4#code#\n')
pytex.before()

# Particle in a box simulation
L_box = 1.0  # Box length
m_box = 1.0  # Particle mass
hbar = 1.0   # Reduced Planck constant

# Position grid
x_box = np.linspace(0, L_box, 1000)

# Calculate first 5 energy levels and wavefunctions
n_levels = 5
box_energies = []
box_wavefunctions = []

for n in range(1, n_levels + 1):
    # Energy eigenvalue
    E_n = (n**2 * np.pi**2 * hbar**2) / (2 * m_box * L_box**2)
    box_energies.append(E_n)

    # Wavefunction
    psi_n = np.sqrt(2/L_box) * np.sin(n * np.pi * x_box / L_box)
    box_wavefunctions.append(psi_n)

# Calculate probability current for superposition state
# Create superposition of n=1 and n=2 states
c1, c2 = 1/np.sqrt(2), 1/np.sqrt(2)  # Equal superposition
psi_superposition = c1 * box_wavefunctions[0] + c2 * box_wavefunctions[1]

print(f"Particle in box: calculated {n_levels} energy levels")
print(f"Energy ratios E_n/E_1: {[E/box_energies[0] for E in box_energies]}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '5'
pytex.line = '339'

print('=>PYTHONTEX:STDOUT#5#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#5#code#\n')
pytex.before()

# Visualize particle in a box
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))

# Left panel: Energy levels and wavefunctions
colors = plt.cm.viridis(np.linspace(0, 1, n_levels))
for n, (psi, E, color) in enumerate(zip(box_wavefunctions, box_energies, colors)):
    # Energy level
    ax1.axhline(E, color=color, linestyle='--', alpha=0.7, linewidth=1)

    # Wavefunction (scaled and shifted)
    psi_scaled = psi * 0.3 + E  # Scale and shift
    ax1.plot(x_box, psi_scaled, color=color, linewidth=2,
             label=f'n={n+1}, E={E:.2f}')

# Draw box walls
ax1.axvline(0, color='black', linewidth=4)
ax1.axvline(L_box, color='black', linewidth=4)
ax1.axhline(0, color='black', linewidth=2)

ax1.set_xlabel('Position x', fontsize=12)
ax1.set_ylabel('Energy', fontsize=12)
ax1.set_title('Particle in a Box: Quantum Energy Levels', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Right panel: Superposition state and probability density
ax2.plot(x_box, psi_superposition, 'b-', linewidth=2, label='ψ(x) = (ψ₁ + ψ₂)/√2')
ax2.fill_between(x_box, 0, np.abs(psi_superposition)**2, alpha=0.5, color='red',
                 label='|ψ(x)|²')

ax2.axvline(0, color='black', linewidth=4)
ax2.axvline(L_box, color='black', linewidth=4)
ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Amplitude / Probability Density', fontsize=12)
ax2.set_title('Quantum Superposition: Interference Pattern', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/particle_in_box.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Particle in box visualization saved to assets/particle_in_box.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '6'
pytex.line = '407'

print('=>PYTHONTEX:STDOUT#6#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#6#code#\n')
pytex.before()

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

# Set random seed for reproducibility
np.random.seed(42)

@jit(nopython=True)
def ising_energy(spins, J, h):
    """Calculate total energy of Ising model configuration"""
    N = spins.shape[0]
    energy = 0.0

    # Nearest neighbor interactions
    for i in range(N):
        for j in range(N):
            # Periodic boundary conditions
            right = spins[i, (j+1) % N]
            down = spins[(i+1) % N, j]
            energy -= J * spins[i, j] * (right + down)

    # External field
    energy -= h * np.sum(spins)
    return energy

@jit(nopython=True)
def metropolis_step(spins, beta, J, h):
    """Single Metropolis Monte Carlo step"""
    N = spins.shape[0]

    for _ in range(N * N):
        # Random site
        i, j = np.random.randint(0, N, 2)

        # Calculate energy change for spin flip
        neighbors = (spins[i, (j+1) % N] + spins[i, (j-1) % N] +
                    spins[(i+1) % N, j] + spins[(i-1) % N, j])

        delta_E = 2 * J * spins[i, j] * neighbors + 2 * h * spins[i, j]

        # Metropolis acceptance
        if delta_E < 0 or np.random.random() < np.exp(-beta * delta_E):
            spins[i, j] *= -1

@jit(nopython=True)
def magnetization(spins):
    """Calculate magnetization per site"""
    return np.mean(spins)

# Ising model parameters
N = 64  # Lattice size
J = 1.0  # Ferromagnetic coupling
h = 0.0  # No external field

# Temperature range for phase transition study
T_min, T_max = 1.0, 4.0
n_temps = 20
temperatures = np.linspace(T_min, T_max, n_temps)

# Critical temperature (analytical result for 2D Ising)
T_c = 2 * J / np.log(1 + np.sqrt(2))  # ≈ 2.269

print(f"2D Ising model simulation: {N}×{N} lattice")
print(f"Theoretical critical temperature: T_c = {T_c:.3f}")

# Monte Carlo simulation
n_steps = 5000
n_equilibration = 1000

magnetizations = []
susceptibilities = []
energies = []
heat_capacities = []

for T in temperatures:
    beta = 1.0 / T

    # Initialize random configuration
    spins = np.random.choice([-1, 1], size=(N, N))

    # Equilibration
    for step in range(n_equilibration):
        metropolis_step(spins, beta, J, h)

    # Measurement phase
    mag_samples = []
    energy_samples = []

    for step in range(n_steps):
        metropolis_step(spins, beta, J, h)

        if step % 10 == 0:  # Sample every 10 steps
            mag_samples.append(abs(magnetization(spins)))
            energy_samples.append(ising_energy(spins, J, h) / (N*N))

    # Calculate thermodynamic quantities
    mag_mean = np.mean(mag_samples)
    mag_var = np.var(mag_samples)
    energy_mean = np.mean(energy_samples)
    energy_var = np.var(energy_samples)

    magnetizations.append(mag_mean)
    susceptibilities.append(beta * mag_var * N * N)
    energies.append(energy_mean)
    heat_capacities.append(beta**2 * energy_var * N * N)

# Store final configuration for visualization
final_spins = spins.copy()

print(f"Monte Carlo simulation completed: {len(temperatures)} temperatures")
print(f"Final magnetization at T={T:.2f}: {magnetizations[-1]:.4f}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '7'
pytex.line = '521'

print('=>PYTHONTEX:STDOUT#7#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#7#code#\n')
pytex.before()

# Visualize Ising model results
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12))

# Magnetization vs temperature
ax1.plot(temperatures, magnetizations, 'bo-', linewidth=2, markersize=6)
ax1.axvline(T_c, color='red', linestyle='--', linewidth=2, label=f'T_c = {T_c:.3f}')
ax1.set_xlabel('Temperature T', fontsize=12)
ax1.set_ylabel('Magnetization |⟨m⟩|', fontsize=12)
ax1.set_title('Ising Model: Magnetic Phase Transition', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Magnetic susceptibility
ax2.plot(temperatures, susceptibilities, 'ro-', linewidth=2, markersize=6)
ax2.axvline(T_c, color='red', linestyle='--', linewidth=2, label=f'T_c = {T_c:.3f}')
ax2.set_xlabel('Temperature T', fontsize=12)
ax2.set_ylabel('Susceptibility χ', fontsize=12)
ax2.set_title('Magnetic Susceptibility: Critical Behavior', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Energy vs temperature
ax3.plot(temperatures, energies, 'go-', linewidth=2, markersize=6)
ax3.axvline(T_c, color='red', linestyle='--', linewidth=2, label=f'T_c = {T_c:.3f}')
ax3.set_xlabel('Temperature T', fontsize=12)
ax3.set_ylabel('Energy per site ⟨E⟩/N²', fontsize=12)
ax3.set_title('Internal Energy: Thermodynamic Behavior', fontsize=14, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Spin configuration
im = ax4.imshow(final_spins, cmap='RdBu', vmin=-1, vmax=1, interpolation='nearest')
ax4.set_title(f'Spin Configuration at T = {temperatures[-1]:.2f}', fontsize=14, fontweight='bold')
ax4.set_xlabel('x', fontsize=12)
ax4.set_ylabel('y', fontsize=12)

# Add colorbar
cbar = plt.colorbar(im, ax=ax4, shrink=0.8)
cbar.set_label('Spin Value', fontsize=12)
cbar.set_ticks([-1, 0, 1])
cbar.set_ticklabels(['↓ (-1)', '0', '↑ (+1)'])

plt.tight_layout()
plt.savefig('assets/ising_model.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Ising model visualization saved to assets/ising_model.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '8'
pytex.line = '583'

print('=>PYTHONTEX:STDOUT#8#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#8#code#\n')
pytex.before()

# 1D Ising model exact solution
def ising_1d_exact(T, J, h, N):
    """Exact solution for 1D Ising model using transfer matrix"""
    beta = 1.0 / T

    # Transfer matrix elements
    T11 = np.exp(beta * (J + h))    # ↑↑
    T12 = np.exp(beta * (-J))       # ↑↓
    T21 = np.exp(beta * (-J))       # ↓↑
    T22 = np.exp(beta * (J - h))    # ↓↓

    # Transfer matrix
    T_matrix = np.array([[T11, T12], [T21, T22]])

    # Eigenvalues
    eigenvals = np.linalg.eigvals(T_matrix)
    lambda_max = np.max(eigenvals)

    # Thermodynamic quantities (thermodynamic limit)
    f_per_site = -T * np.log(lambda_max)  # Free energy per site

    # For finite system
    Z = np.trace(np.linalg.matrix_power(T_matrix, N))  # Partition function
    F = -T * np.log(Z)  # Free energy

    return F, f_per_site, Z

# Calculate exact 1D results
T_range_1d = np.linspace(0.5, 5.0, 100)
J_1d = 1.0
h_1d = 0.1  # Small external field
N_1d = 50

free_energies = []
partition_functions = []

for T in T_range_1d:
    F, f_per_site, Z = ising_1d_exact(T, J_1d, h_1d, N_1d)
    free_energies.append(f_per_site)
    partition_functions.append(Z)

print(f"1D Ising model exact calculation: {len(T_range_1d)} temperature points")
print(f"System size: N = {N_1d}, coupling J = {J_1d}, field h = {h_1d}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '9'
pytex.line = '629'

print('=>PYTHONTEX:STDOUT#9#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#9#code#\n')
pytex.before()

# Visualize thermodynamic properties
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# Free energy per site
ax1.plot(T_range_1d, free_energies, 'b-', linewidth=2)
ax1.set_xlabel('Temperature T', fontsize=12)
ax1.set_ylabel('Free Energy per Site f', fontsize=12)
ax1.set_title('1D Ising Model: Free Energy (Exact Solution)', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Partition function (log scale)
ax2.semilogy(T_range_1d, partition_functions, 'r-', linewidth=2)
ax2.set_xlabel('Temperature T', fontsize=12)
ax2.set_ylabel('Partition Function Z', fontsize=12)
ax2.set_title('Partition Function: Statistical Weight', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/thermodynamics.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Thermodynamics visualization saved to assets/thermodynamics.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '10'
pytex.line = '674'

print('=>PYTHONTEX:STDOUT#10#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#10#code#\n')
pytex.before()

# 2D square lattice tight-binding model
def tight_binding_2d(kx, ky, t, epsilon_0):
    """Energy dispersion for 2D square lattice"""
    return epsilon_0 - 2*t*(np.cos(kx) + np.cos(ky))

def graphene_dispersion(kx, ky, t):
    """Energy dispersion for graphene (honeycomb lattice)"""
    # Nearest neighbor phase factors
    f = 1 + np.exp(1j * kx) + np.exp(1j * (kx/2 + np.sqrt(3)*ky/2))
    return t * np.abs(f), -t * np.abs(f)  # Two bands

# Parameters
t = 1.0  # Hopping parameter
epsilon_0 = 0.0  # Onsite energy

# Create 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 band structure
E_square = tight_binding_2d(KX, KY, t, epsilon_0)
E_graphene_plus, E_graphene_minus = graphene_dispersion(KX, KY, t)

# Density of states calculation
def calculate_dos(energies, n_bins=200):
    """Calculate density of states"""
    E_flat = energies.flatten()
    E_min, E_max = E_flat.min(), E_flat.max()
    E_bins = np.linspace(E_min, E_max, n_bins)
    hist, bin_edges = np.histogram(E_flat, bins=E_bins)
    E_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
    dos = hist / (E_bins[1] - E_bins[0]) / len(E_flat)
    return E_centers, dos

E_dos_square, dos_square = calculate_dos(E_square)
E_dos_graphene, dos_graphene = calculate_dos(np.concatenate([E_graphene_plus.flatten(),
                                                           E_graphene_minus.flatten()]))

print(f"Band structure calculation completed: {k_points}×{k_points} k-points")
print(f"Square lattice bandwidth: {E_square.max() - E_square.min():.3f}")
print(f"Graphene bandwidth: {E_graphene_plus.max() - E_graphene_minus.min():.3f}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '11'
pytex.line = '720'

print('=>PYTHONTEX:STDOUT#11#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#11#code#\n')
pytex.before()

# Visualize band structure
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12))

# Square lattice band structure
im1 = ax1.contourf(KX, KY, E_square, levels=50, cmap='viridis')
ax1.set_xlabel('kₓ', fontsize=12)
ax1.set_ylabel('kᵧ', fontsize=12)
ax1.set_title('Square Lattice: Energy Dispersion E(k)', fontsize=14, fontweight='bold')
plt.colorbar(im1, ax=ax1, label='Energy E/t')

# Graphene band structure (valence band)
im2 = ax2.contourf(KX, KY, E_graphene_minus, levels=50, cmap='plasma')
ax2.set_xlabel('kₓ', fontsize=12)
ax2.set_ylabel('kᵧ', fontsize=12)
ax2.set_title('Graphene: Valence Band Dispersion', fontsize=14, fontweight='bold')
plt.colorbar(im2, ax=ax2, label='Energy E/t')

# Density of states - square lattice
ax3.plot(E_dos_square, dos_square, 'b-', linewidth=2, label='Square lattice')
ax3.fill_between(E_dos_square, 0, dos_square, alpha=0.3)
ax3.set_xlabel('Energy E/t', fontsize=12)
ax3.set_ylabel('Density of States', fontsize=12)
ax3.set_title('Electronic Density of States: Square Lattice', fontsize=14, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Density of states - graphene
ax4.plot(E_dos_graphene, dos_graphene, 'r-', linewidth=2, label='Graphene')
ax4.fill_between(E_dos_graphene, 0, dos_graphene, alpha=0.3, color='red')
ax4.axvline(0, color='black', linestyle='--', alpha=0.7, label='Dirac point')
ax4.set_xlabel('Energy E/t', fontsize=12)
ax4.set_ylabel('Density of States', fontsize=12)
ax4.set_title('Electronic Density of States: Graphene', fontsize=14, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/band_structure.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Band structure visualization saved to assets/band_structure.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '12'
pytex.line = '776'

print('=>PYTHONTEX:STDOUT#12#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#12#code#\n')
pytex.before()

# Fermi surface calculation
def fermi_surface_2d(kx, ky, t, mu, T=0.01):
    """Calculate Fermi surface for 2D system"""
    E = tight_binding_2d(kx, ky, t, 0)

    # Fermi-Dirac distribution
    if T > 0:
        f = 1.0 / (1.0 + np.exp((E - mu) / T))
    else:
        f = (E < mu).astype(float)

    return E, f

# Parameters for different filling levels
t_fermi = 1.0
mu_values = [-3.0, -2.0, 0.0, 2.0]  # Different chemical potentials
filling_names = ['Low filling', 'Quarter filling', 'Half filling', 'High filling']

# Higher resolution for Fermi surface
k_fermi = 200
kx_fermi = np.linspace(-np.pi, np.pi, k_fermi)
ky_fermi = np.linspace(-np.pi, np.pi, k_fermi)
KX_fermi, KY_fermi = np.meshgrid(kx_fermi, ky_fermi)

# Calculate Fermi surfaces
fermi_data = []
for mu in mu_values:
    E, f = fermi_surface_2d(KX_fermi, KY_fermi, t_fermi, mu, T=0.1)
    fermi_data.append((E, f))

print(f"Fermi surface calculation: {k_fermi}×{k_fermi} k-points")
print(f"Chemical potentials: {mu_values}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '13'
pytex.line = '811'

print('=>PYTHONTEX:STDOUT#13#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#13#code#\n')
pytex.before()

# Visualize Fermi surfaces
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
axes = axes.flatten()

for i, ((E, f), mu, name) in enumerate(zip(fermi_data, mu_values, filling_names)):
    ax = axes[i]

    # Plot Fermi surface as contour
    contour = ax.contour(KX_fermi, KY_fermi, E, levels=[mu], colors='red', linewidths=3)

    # Fill occupied states
    filled = ax.contourf(KX_fermi, KY_fermi, f, levels=50, cmap='Blues', alpha=0.7)

    ax.set_xlabel('kₓ/π', fontsize=12)
    ax.set_ylabel('kᵧ/π', fontsize=12)
    ax.set_title(f'{name}: μ = {mu:.1f}t', fontsize=14, fontweight='bold')
    ax.set_xlim(-np.pi, np.pi)
    ax.set_ylim(-np.pi, np.pi)

    # Add Brillouin zone boundary
    ax.plot([-np.pi, np.pi, np.pi, -np.pi, -np.pi],
            [-np.pi, -np.pi, np.pi, np.pi, -np.pi],
            'k--', alpha=0.5, linewidth=1)

    # Colorbar for occupation
    if i == 1:  # Add colorbar to one plot
        cbar = plt.colorbar(filled, ax=ax, shrink=0.8)
        cbar.set_label('Occupation f(k)', fontsize=10)

plt.tight_layout()
plt.savefig('assets/fermi_surface.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Fermi surface visualization saved to assets/fermi_surface.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '14'
pytex.line = '860'

print('=>PYTHONTEX:STDOUT#14#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#14#code#\n')
pytex.before()

# 1D monoatomic chain phonon dispersion
def phonon_monoatomic(k, K, M):
    """Phonon frequency for 1D monoatomic chain"""
    return np.sqrt(4*K/M * np.sin(k/2)**2)

# 1D diatomic chain phonon dispersion
def phonon_diatomic(k, K, M1, M2):
    """Phonon frequencies for 1D diatomic chain"""
    M_sum = M1 + M2
    M_prod = M1 * M2

    # Discriminant
    disc = (M_sum)**2 - 4*M_prod*np.cos(k)**2

    # Optical and acoustic branches
    omega_plus = np.sqrt(K * (M_sum + np.sqrt(disc)) / (2*M_prod))
    omega_minus = np.sqrt(K * (M_sum - np.sqrt(disc)) / (2*M_prod))

    return omega_plus, omega_minus

# Parameters
K = 1.0  # Spring constant
M = 1.0  # Mass (monoatomic)
M1, M2 = 1.0, 2.0  # Masses (diatomic)

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

# Calculate dispersions
omega_mono = phonon_monoatomic(k_points_1d, K, M)
omega_optical, omega_acoustic = phonon_diatomic(k_points_1d, K, M1, M2)

# Density of phonon states
def phonon_dos(k_vals, omega_vals, n_bins=100):
    """Calculate phonon density of states"""
    # Use only positive k values to avoid double counting
    k_pos = k_vals[k_vals >= 0]
    omega_pos = omega_vals[k_vals >= 0]

    omega_min, omega_max = 0, omega_pos.max()
    omega_bins = np.linspace(omega_min, omega_max, n_bins)

    hist, bin_edges = np.histogram(omega_pos, bins=omega_bins)
    omega_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
    dos = hist / (omega_bins[1] - omega_bins[0]) / len(omega_pos)

    return omega_centers, dos

# Calculate phonon DOS
omega_dos_mono, dos_mono = phonon_dos(k_points_1d, omega_mono)
omega_dos_opt, dos_opt = phonon_dos(k_points_1d, omega_optical)
omega_dos_ac, dos_ac = phonon_dos(k_points_1d, omega_acoustic)

print(f"Phonon dispersion calculation: {len(k_points_1d)} k-points")
print(f"Monoatomic max frequency: {omega_mono.max():.3f}")
print(f"Diatomic optical max: {omega_optical.max():.3f}")
print(f"Diatomic acoustic max: {omega_acoustic.max():.3f}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '15'
pytex.line = '920'

print('=>PYTHONTEX:STDOUT#15#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#15#code#\n')
pytex.before()

# Visualize phonon dispersion
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))

# Dispersion relations
ax1.plot(k_points_1d/np.pi, omega_mono, 'b-', linewidth=2, label='Monoatomic chain')
ax1.plot(k_points_1d/np.pi, omega_optical, 'r-', linewidth=2, label='Optical branch')
ax1.plot(k_points_1d/np.pi, omega_acoustic, 'g-', linewidth=2, label='Acoustic branch')

ax1.set_xlabel('Wave vector k/π', fontsize=12)
ax1.set_ylabel('Frequency ω√(M/K)', fontsize=12)
ax1.set_title('Phonon Dispersion Relations', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-1, 1)

# Density of states
ax2.plot(omega_dos_mono, dos_mono, 'b-', linewidth=2, label='Monoatomic')
ax2.fill_between(omega_dos_mono, 0, dos_mono, alpha=0.3, color='blue')

ax2.plot(omega_dos_opt, dos_opt, 'r-', linewidth=2, label='Optical')
ax2.fill_between(omega_dos_opt, 0, dos_opt, alpha=0.3, color='red')

ax2.plot(omega_dos_ac, dos_ac, 'g-', linewidth=2, label='Acoustic')
ax2.fill_between(omega_dos_ac, 0, dos_ac, alpha=0.3, color='green')

ax2.set_xlabel('Frequency ω√(M/K)', fontsize=12)
ax2.set_ylabel('Density of States', fontsize=12)
ax2.set_title('Phonon Density of States', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/phonon_dispersion.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Phonon dispersion visualization saved to assets/phonon_dispersion.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '16'
pytex.line = '979'

print('=>PYTHONTEX:STDOUT#16#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#16#code#\n')
pytex.before()

# Simple mean-field solution of Hubbard model
def hubbard_mean_field(t, U, n, T=0.1, max_iter=1000, tol=1e-6):
    """
    Mean-field solution of Hubbard model on 2D square lattice
    Returns self-consistent chemical potential and magnetization
    """
    # Initial guess
    mu = U * n / 2  # Chemical potential
    m = 0.1  # Magnetization (small initial value)

    for iteration in range(max_iter):
        # Mean-field Hamiltonian eigenvalues
        # For 2D square lattice: epsilon_k = -2t(cos(kx) + cos(ky))
        # With mean-field correction: E_k^± = epsilon_k ± U*m/2

        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)

        epsilon_k = -2*t*(np.cos(KX) + np.cos(KY))
        E_plus = epsilon_k - mu + U*m/2   # Spin-up band
        E_minus = epsilon_k - mu - U*m/2  # Spin-down band

        # Fermi-Dirac distribution
        f_plus = 1.0 / (1.0 + np.exp(E_plus / T))
        f_minus = 1.0 / (1.0 + np.exp(E_minus / T))

        # New density and magnetization
        n_new = np.mean(f_plus + f_minus)
        m_new = np.mean(f_plus - f_minus)

        # Check convergence
        if abs(m_new - m) < tol:
            break

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

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

    return mu, m, f_plus, f_minus, epsilon_k

# Parameters for different interaction strengths
t_hub = 1.0
density = 0.8  # Filling factor
U_values = [0.0, 2.0, 4.0, 6.0]  # Interaction strengths

results = []
for U in U_values:
    mu, m, f_up, f_dn, eps_k = hubbard_mean_field(t_hub, U, density)
    results.append((U, mu, m, f_up, f_dn, eps_k))
    print(f"U/t = {U/t_hub:.1f}: μ = {mu:.3f}, m = {m:.4f}")

print(f"Hubbard model mean-field calculation completed")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '17'
pytex.line = '1038'

print('=>PYTHONTEX:STDOUT#17#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#17#code#\n')
pytex.before()

# Plot Hubbard model results
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12))

# Magnetization vs interaction strength
U_plot = [result[0] for result in results]
m_plot = [abs(result[2]) for result in results]  # Take absolute value

ax1.plot(U_plot, m_plot, 'ro-', linewidth=2, markersize=8)
ax1.set_xlabel('Interaction U/t', fontsize=12)
ax1.set_ylabel('Magnetization |m|', fontsize=12)
ax1.set_title('Hubbard Model: Magnetic Instability', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Chemical potential vs interaction
mu_plot = [result[1] for result in results]
ax2.plot(U_plot, mu_plot, 'bo-', linewidth=2, markersize=8)
ax2.set_xlabel('Interaction U/t', fontsize=12)
ax2.set_ylabel('Chemical Potential μ/t', fontsize=12)
ax2.set_title('Chemical Potential: Correlation Effects', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Spin-resolved occupations for strong coupling case
U_strong = results[-1]  # Strongest coupling
f_up_strong = U_strong[3]
f_dn_strong = U_strong[4]

im3 = ax3.contourf(f_up_strong, levels=50, cmap='Reds')
ax3.set_title(f'Spin-Up Occupation (U = {U_strong[0]:.1f}t)', fontsize=14, fontweight='bold')
ax3.set_xlabel('kₓ index', fontsize=12)
ax3.set_ylabel('kᵧ index', fontsize=12)
plt.colorbar(im3, ax=ax3, shrink=0.8)

im4 = ax4.contourf(f_dn_strong, levels=50, cmap='Blues')
ax4.set_title(f'Spin-Down Occupation (U = {U_strong[0]:.1f}t)', fontsize=14, fontweight='bold')
ax4.set_xlabel('kₓ index', fontsize=12)
ax4.set_ylabel('kᵧ index', fontsize=12)
plt.colorbar(im4, ax=ax4, shrink=0.8)

plt.tight_layout()
plt.savefig('assets/hubbard_model.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Hubbard model visualization saved to assets/hubbard_model.pdf")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '18'
pytex.line = '1097'

print('=>PYTHONTEX:STDOUT#18#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#18#code#\n')
pytex.before()

# Demonstrate numerical precision considerations
def convergence_study():
    """Study convergence of numerical methods"""

    # Test case: numerical integration of oscillatory function
    def integrand(x):
        return np.sin(10*x) * np.exp(-x**2)

    # Analytical result (approximate)
    analytical = 0.4995  # Pre-computed high-precision result

    # Different grid sizes
    N_values = [10, 20, 50, 100, 200, 500, 1000]
    errors_trapz = []
    errors_simpson = []

    for N in N_values:
        x = np.linspace(0, 3, N)
        y = integrand(x)

        # Trapezoidal rule
        integral_trapz = np.trapz(y, x)
        error_trapz = abs(integral_trapz - analytical)
        errors_trapz.append(error_trapz)

        # Simpson's rule (if N is odd)
        if N % 2 == 1:
            from scipy.integrate import simpson
            integral_simpson = simpson(y, x)
            error_simpson = abs(integral_simpson - analytical)
        else:
            error_simpson = np.nan
        errors_simpson.append(error_simpson)

    return N_values, errors_trapz, errors_simpson

N_vals, err_trapz, err_simp = convergence_study()

print("Numerical integration convergence study completed")
print(f"Grid sizes tested: {N_vals}")


pytex.after()
pytex.command = 'code'
pytex.set_context('')
pytex.args = ''
pytex.instance = '19'
pytex.line = '1140'

print('=>PYTHONTEX:STDOUT#19#code#')
sys.stderr.write('=>PYTHONTEX:STDERR#19#code#\n')
pytex.before()

# Plot convergence study
fig, ax = plt.subplots(figsize=(12, 8))

ax.loglog(N_vals, err_trapz, 'bo-', linewidth=2, markersize=6, label='Trapezoidal Rule')

# Filter out NaN values for Simpson's rule
N_simp = [N for N, err in zip(N_vals, err_simp) if not np.isnan(err)]
err_simp_clean = [err for err in err_simp if not np.isnan(err)]
ax.loglog(N_simp, err_simp_clean, 'rs-', linewidth=2, markersize=6, label="Simpson's Rule")

# Theoretical scaling lines
N_theory = np.array(N_vals)
ax.loglog(N_theory, 1e-1 * (N_theory[0]/N_theory)**2, 'k--', alpha=0.7, label='N⁻² scaling')
ax.loglog(N_theory, 1e-3 * (N_theory[0]/N_theory)**4, 'k:', alpha=0.7, label='N⁻⁴ scaling')

ax.set_xlabel('Number of Grid Points N', fontsize=12)
ax.set_ylabel('Absolute Error', fontsize=12)
ax.set_title('Numerical Integration: Convergence Analysis', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/convergence_study.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Convergence study visualization saved to assets/convergence_study.pdf")


pytex.after()


pytex.cleanup()