statistical_mechanics.py

#!/usr/bin/env python3
"""
Statistical Mechanics Monte Carlo Simulations
==============================================

This module implements Monte Carlo simulations for statistical mechanics systems:
- 2D Ising model with Metropolis algorithm
- Exact 1D Ising model solutions
- Thermodynamic property calculations
- Phase transition analysis

Keywords: ising model monte carlo, statistical mechanics python, metropolis algorithm,
phase transition simulation, partition function calculation, magnetic susceptibility

Author: CoCalc Scientific Templates
License: MIT
"""

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

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

class IsingModel2D:
    """
    2D Ising model Monte Carlo simulation using Metropolis algorithm.

    Hamiltonian: H = -J Σ_{<i,j>} S_i S_j - h Σ_i S_i
    where S_i = ±1 are spin variables.
    """

    def __init__(self, N, J=1.0, h=0.0):
        """
        Parameters:
        -----------
        N : int
            Linear lattice size (N×N total spins)
        J : float
            Nearest-neighbor coupling constant
        h : float
            External magnetic field
        """
        self.N = N
        self.J = J
        self.h = h
        self.spins = np.random.choice([-1, 1], size=(N, N))

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

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

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

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

        for _ in range(N * N):
            # Choose 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 criterion
            if delta_E < 0 or np.random.random() < np.exp(-beta * delta_E):
                spins[i, j] *= -1

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

    def simulate(self, temperature, n_steps=5000, n_equilibration=1000, sample_interval=10):
        """
        Run Monte Carlo simulation at given temperature.

        Parameters:
        -----------
        temperature : float
            Temperature T
        n_steps : int
            Number of Monte Carlo steps after equilibration
        n_equilibration : int
            Number of equilibration steps
        sample_interval : int
            Sampling interval for measurements

        Returns:
        --------
        observables : dict
            Dictionary containing measured observables
        """
        beta = 1.0 / temperature

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

        # Measurement phase
        magnetizations = []
        energies = []

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

            if step % sample_interval == 0:
                mag = abs(self.magnetization(self.spins))
                energy = self.total_energy(self.spins, self.J, self.h) / (self.N**2)

                magnetizations.append(mag)
                energies.append(energy)

        # Calculate thermodynamic quantities
        mag_mean = np.mean(magnetizations)
        mag_var = np.var(magnetizations)
        energy_mean = np.mean(energies)
        energy_var = np.var(energies)

        observables = {
            'magnetization': mag_mean,
            'susceptibility': beta * mag_var * self.N**2,
            'energy': energy_mean,
            'heat_capacity': beta**2 * energy_var * self.N**2,
            'final_configuration': self.spins.copy()
        }

        return observables

class IsingModel1D:
    """
    1D Ising model exact solution using transfer matrix method.

    The transfer matrix approach provides exact thermodynamic properties
    for comparison with Monte Carlo results and validation.
    """

    def __init__(self, J=1.0, h=0.0):
        """
        Parameters:
        -----------
        J : float
            Nearest-neighbor coupling constant
        h : float
            External magnetic field
        """
        self.J = J
        self.h = h

    def transfer_matrix(self, temperature):
        """
        Construct transfer matrix for given temperature.

        Returns:
        --------
        T_matrix : ndarray
            2×2 transfer matrix
        """
        beta = 1.0 / temperature

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

        return np.array([[T11, T12], [T21, T22]])

    def exact_properties(self, temperature, N=None):
        """
        Calculate exact thermodynamic properties.

        Parameters:
        -----------
        temperature : float
            Temperature T
        N : int, optional
            System size for finite-size calculations

        Returns:
        --------
        properties : dict
            Exact thermodynamic properties
        """
        T_matrix = self.transfer_matrix(temperature)
        eigenvals = np.linalg.eigvals(T_matrix)
        lambda_max = np.max(eigenvals)

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

        properties = {
            'free_energy_per_site': f_per_site,
            'eigenvalues': eigenvals
        }

        # Finite system properties
        if N is not None:
            Z = np.trace(np.linalg.matrix_power(T_matrix, N))
            F_total = -temperature * np.log(Z)
            properties['partition_function'] = Z
            properties['free_energy_total'] = F_total

        return properties

def critical_temperature_2d():
    """
    Return the exact critical temperature for 2D Ising model.

    T_c = 2J / ln(1 + √2) ≈ 2.269 J
    """
    return 2.0 / np.log(1 + np.sqrt(2))

def phase_transition_study(lattice_size=64, T_min=1.0, T_max=4.0, n_temps=20):
    """
    Study the magnetic phase transition in 2D Ising model.

    Parameters:
    -----------
    lattice_size : int
        Linear size of the lattice
    T_min, T_max : float
        Temperature range
    n_temps : int
        Number of temperature points

    Returns:
    --------
    temperatures : array
        Temperature points
    results : list
        List of simulation results at each temperature
    """
    print(f"Phase transition study: {lattice_size}×{lattice_size} lattice")

    temperatures = np.linspace(T_min, T_max, n_temps)
    T_c = critical_temperature_2d()
    print(f"Theoretical critical temperature: T_c = {T_c:.3f}")

    ising = IsingModel2D(N=lattice_size, J=1.0, h=0.0)
    results = []

    for i, T in enumerate(temperatures):
        print(f"  Temperature {i+1}/{n_temps}: T = {T:.2f}")
        observables = ising.simulate(T, n_steps=5000, n_equilibration=1000)
        results.append(observables)

    return temperatures, results, T_c

def exact_vs_monte_carlo_comparison():
    """
    Compare 1D Ising model exact solution with Monte Carlo simulation.
    """
    print("=== 1D Ising Model: Exact vs Monte Carlo ===")

    # Parameters
    J = 1.0
    h = 0.1  # Small field to break symmetry
    temperatures = np.linspace(0.5, 5.0, 20)

    # Exact solution
    exact_model = IsingModel1D(J=J, h=h)
    exact_results = []

    for T in temperatures:
        props = exact_model.exact_properties(T, N=100)
        exact_results.append(props)

    print(f"Exact 1D calculation completed: {len(temperatures)} temperatures")
    return temperatures, exact_results

def thermodynamic_integration_demo():
    """
    Demonstrate calculation of thermodynamic properties.
    """
    print("=== Thermodynamic Integration Demo ===")

    # Calculate free energy as function of temperature
    temperatures = np.linspace(0.5, 5.0, 100)
    model_1d = IsingModel1D(J=1.0, h=0.1)

    free_energies = []
    partition_functions = []

    for T in temperatures:
        props = model_1d.exact_properties(T, N=50)
        free_energies.append(props['free_energy_per_site'])
        partition_functions.append(props['partition_function'])

    return temperatures, free_energies, partition_functions

def extract_critical_exponents(temperatures, magnetizations, T_c):
    """
    Extract critical exponents from Monte Carlo data.

    Near T_c: m ∝ |T - T_c|^β where β ≈ 1/8 for 2D Ising model
    """
    # Focus on temperatures near T_c
    mask = np.abs(temperatures - T_c) < 0.5
    T_fit = temperatures[mask]
    m_fit = np.array(magnetizations)[mask]

    # Only use temperatures below T_c where magnetization is non-zero
    below_Tc = T_fit < T_c
    if np.sum(below_Tc) > 3:
        reduced_temp = (T_c - T_fit[below_Tc]) / T_c
        log_m = np.log(m_fit[below_Tc] + 1e-10)  # Avoid log(0)
        log_t = np.log(reduced_temp + 1e-10)

        # Linear fit in log-log space
        coeffs = np.polyfit(log_t, log_m, 1)
        beta_critical = coeffs[0]

        print(f"Estimated critical exponent β = {beta_critical:.3f}")
        print(f"Theoretical value β = 0.125")

        return beta_critical
    else:
        print("Insufficient data points below T_c for critical exponent extraction")
        return None

if __name__ == "__main__":
    # Run demonstrations
    temperatures, results, T_c = phase_transition_study(lattice_size=32, n_temps=10)

    # Extract observables
    magnetizations = [r['magnetization'] for r in results]
    energies = [r['energy'] for r in results]
    susceptibilities = [r['susceptibility'] for r in results]

    print(f"\nFinal results:")
    print(f"Lowest T magnetization: {magnetizations[0]:.4f}")
    print(f"Highest T magnetization: {magnetizations[-1]:.4f}")

    # Critical exponent analysis
    extract_critical_exponents(temperatures, magnetizations, T_c)

    # Run other demonstrations
    exact_vs_monte_carlo_comparison()
    thermodynamic_integration_demo()