Examples ======== This page contains practical examples demonstrating various use cases of the Free Fermion Library. Basic Examples -------------- Simple Pfaffian Calculation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff # Create a 4x4 skew-symmetric matrix A = np.array([[0, 1, 2, 3], [-1, 0, 4, 5], [-2, -4, 0, 6], [-3, -5, -6, 0]]) # Compute pfaffian pf_value = ff.pf(A) print(f"Pfaffian: {pf_value}") # Verify: pf(A)^2 should equal det(A) det_value = np.linalg.det(A) print(f"det(A): {det_value}") print(f"pf(A)^2: {pf_value**2}") Two-Site System Analysis ~~~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff # Two-site system n_sites = 2 alphas = ff.jordan_wigner_alphas(n_sites) # Hopping Hamiltonian: H = -t(a†₀a₁ + a†₁a₀) t = 1.0 A = np.array([[0, -t], [-t, 0]]) H = ff.build_H(n_sites, A) print("Hamiltonian coefficient matrix:") ff._print(H) # Generate ground state rho = ff.generate_gaussian_state(n_sites, H, alphas) # Compute correlation matrix gamma = ff.compute_2corr_matrix(rho, n_sites, alphas) print("\nCorrelation matrix:") ff._print(gamma) Intermediate Examples --------------------- Kitaev Chain Model ~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff def kitaev_chain(n_sites, mu, t, delta): """ Create Kitaev chain Hamiltonian. H = -μΣᵢnᵢ - tΣᵢ(a†ᵢaᵢ₊₁ + h.c.) + ΔΣᵢ(aᵢaᵢ₊₁ + h.c.) """ # Chemical potential term A = -mu * np.eye(n_sites) # Hopping term for i in range(n_sites - 1): A[i, i+1] = -t A[i+1, i] = -t # Pairing term B = np.zeros((n_sites, n_sites)) for i in range(n_sites - 1): B[i, i+1] = delta B[i+1, i] = -delta return ff.build_H(n_sites, A, B) # Parameters n_sites = 6 mu = 0.5 # Chemical potential t = 1.0 # Hopping strength delta = 0.8 # Pairing strength # Build Hamiltonian H = kitaev_chain(n_sites, mu, t, delta) # Diagonalize eigenvals, eigenvecs = ff.eigh_sp(H) print("Energy eigenvalues:") energies = np.diag(eigenvals)[n_sites:] # Positive eigenvalues ff._print(energies) # Check for zero modes (topological phase) zero_modes = np.abs(energies) < 1e-10 print(f"\nNumber of zero modes: {np.sum(zero_modes)}") Random Matrix Ensemble ~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff import matplotlib.pyplot as plt def random_gaussian_ensemble(n_sites, num_samples=100): """Generate statistics for random Gaussian ensembles.""" eigenvalues = [] for _ in range(num_samples): # Random Hermitian matrix A = np.random.randn(n_sites, n_sites) A = A + A.T # Build Hamiltonian H = ff.build_H(n_sites, A) # Diagonalize evals, _ = ff.eigh_sp(H) eigenvalues.extend(np.diag(evals)[n_sites:]) return np.array(eigenvalues) # Generate ensemble n_sites = 4 eigenvals = random_gaussian_ensemble(n_sites, num_samples=200) print(f"Generated {len(eigenvals)} eigenvalues") print(f"Mean: {np.mean(eigenvals):.3f}") print(f"Std: {np.std(eigenvals):.3f}") Advanced Examples ----------------- Correlation Function Analysis ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff def analyze_correlations(n_sites, A_matrix): """Analyze correlation functions for a given system.""" alphas = ff.jordan_wigner_alphas(n_sites) H = ff.build_H(n_sites, A_matrix) rho = ff.generate_gaussian_state(n_sites, H, alphas) # Compute different correlation matrices gamma = ff.compute_2corr_matrix(rho, n_sites, alphas) cov = ff.compute_cov_matrix(rho, n_sites, alphas) # Extract physical quantities occupations = np.diag(gamma)[n_sites:] # ⟨a†ᵢaᵢ⟩ # Correlation lengths (simplified) correlations = [] for i in range(n_sites): for j in range(i+1, n_sites): corr = gamma[i, j + n_sites] # ⟨a†ᵢaⱼ⟩ correlations.append((j-i, abs(corr))) return { 'occupations': occupations, 'correlations': correlations, 'gamma': gamma, 'covariance': cov } # Example: Uniform hopping chain n_sites = 6 A = np.zeros((n_sites, n_sites)) for i in range(n_sites - 1): A[i, i+1] = A[i+1, i] = -1.0 results = analyze_correlations(n_sites, A) print("Site occupations:") ff._print(results['occupations']) print("\nDistance vs correlation strength:") for dist, corr in results['correlations'][:5]: print(f"Distance {dist}: {corr:.4f}") Perfect Matching in Graphs ~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff import networkx as nx def analyze_perfect_matchings(n_vertices): """Analyze perfect matchings using pfaffian method.""" if n_vertices % 2 != 0: print("Need even number of vertices for perfect matchings") return # Generate random planar graph G = ff.generate_random_planar_graph(n_vertices, seed=42) if G is None: print("Could not generate planar graph") return # Method 1: Brute force enumeration matchings_brute = ff.find_perfect_matchings(G) # Method 2: Pfaffian calculation pfo_matrix = ff.pfo_algorithm(G, verbose=False) pf_value = ff.pf(pfo_matrix) print(f"Graph with {n_vertices} vertices:") print(f"Edges: {G.number_of_edges()}") print(f"Perfect matchings (brute force): {len(matchings_brute)}") print(f"Perfect matchings (pfaffian): {int(abs(pf_value))}") # Verify they match if abs(len(matchings_brute) - abs(pf_value)) < 1e-10: print("✓ Methods agree!") else: print("✗ Methods disagree!") return G, matchings_brute, pfo_matrix # Test with different sizes for n in [4, 6, 8]: print(f"\n{'-'*40}") analyze_perfect_matchings(n) Symplectic Geometry Example ~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff def symplectic_transformation_example(): """Demonstrate symplectic transformations.""" n_sites = 3 # Create random Hamiltonian A = np.random.randn(n_sites, n_sites) A = A + A.T H = ff.build_H(n_sites, A) print("Original Hamiltonian:") ff._print(H) # Symplectic diagonalization L, U = ff.eigh_sp(H) print(f"\nSymplectic matrix U is symplectic: {ff.is_symp(U)}") print(f"Eigenvalue matrix in canonical form: {ff.check_canonical_form(L)}") # Verify diagonalization: U† H U = L H_diag = U.conj().T @ H @ U print("\nDiagonalized Hamiltonian:") ff._print(H_diag) print("\nExpected eigenvalue structure:") ff._print(L) # Check if they match if np.allclose(H_diag, L): print("✓ Diagonalization successful!") else: print("✗ Diagonalization failed!") return H, L, U symplectic_transformation_example() Utility Examples ---------------- Matrix Cleaning ~~~~~~~~~~~~~~~ :: import numpy as np import ff # Create matrix with numerical noise clean_matrix = np.array([[1.0, 0.0, 0.5], [0.0, 2.0, 0.0], [0.5, 0.0, 1.5]]) noise = 1e-12 * np.random.randn(3, 3) noisy_matrix = clean_matrix + noise print("Original clean matrix:") ff._print(clean_matrix) print("\nWith numerical noise:") ff._print(noisy_matrix, k=15) print("\nAfter cleaning:") cleaned = ff.clean(noisy_matrix, threshold=1e-10) ff._print(cleaned) Custom Printing ~~~~~~~~~~~~~~~ :: import numpy as np import ff # Complex matrix with varying magnitudes matrix = np.array([[1.23456789, 1e-8 + 2.3456j], [0.000123456, 9.87654321]]) print("Default NumPy printing:") print(matrix) print("\nCustom precision (3 decimal places):") ff._print(matrix, k=3) print("\nCustom precision (6 decimal places):") ff._print(matrix, k=6) Performance Examples -------------------- Large System Benchmark ~~~~~~~~~~~~~~~~~~~~~~ :: import numpy as np import ff import time def benchmark_large_system(n_sites): """Benchmark performance for large systems.""" print(f"Benchmarking {n_sites}-site system...") # Generate random Hamiltonian A = np.random.randn(n_sites, n_sites) A = A + A.T # Time Hamiltonian construction start = time.time() H = ff.build_H(n_sites, A) h_time = time.time() - start # Time diagonalization start = time.time() eigenvals, eigenvecs = ff.eigh_sp(H) diag_time = time.time() - start # Time state generation alphas = ff.jordan_wigner_alphas(n_sites) start = time.time() rho = ff.generate_gaussian_state(n_sites, H, alphas) state_time = time.time() - start print(f" Hamiltonian construction: {h_time:.4f}s") print(f" Diagonalization: {diag_time:.4f}s") print(f" State generation: {state_time:.4f}s") print(f" Total: {h_time + diag_time + state_time:.4f}s") # Benchmark different sizes for n in [5, 10, 15]: benchmark_large_system(n) print()