Canonical Ensemble in Statistical Mechanics¶

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

Oscillatore Armonico Classico monodimensionale¶

Si considera un oscillatore armonico classico monodimensionale descritto dall'Hamiltoniana $$\large H(q,p) = \frac{p^2}{2} + \frac{q^2}{2}, $$

nell'insieme canonico a temperatura fissata.
Si adottano unità ridotte tali che $$\large k_B = 1, \qquad T = 1, \qquad \beta = \frac{1}{T} = 1. $$

Lo spazio delle fasi è continuo ma troncato a valori finiti di $q$ e $p$, al fine di rendere numericamente trattabile l'integrale di partizione.

Energia interna¶

Il valore numerico ottenuto è $$\large E \simeq 1.001. $$

Per l'oscillatore armonico classico, il teorema di equipartizione dell'energia implica che ogni termine quadratico nell'Hamiltoniana contribuisce con $\frac{1}{2}k_B T$. Poiché l'Hamiltoniana contiene due termini quadratici, $$\large \langle E \rangle = \frac{1}{2}T + \frac{1}{2}T = T = 1. $$

Il risultato numerico è quindi in perfetto accordo con la previsione teorica.

Funzione di partizione¶

La funzione di partizione classica è definita come $$\large Z(\beta) = \int dq\,dp\; e^{-\beta H(q,p)}. $$

Per l'oscillatore armonico classico non troncato, l'integrale può essere calcolato analiticamente e fornisce $$\large Z_\infty = \frac{2\pi}{\beta}. $$

Per $\beta = 1$ si ottiene $$\large Z_\infty = 2\pi \approx 6.283. $$

Il valore numerico $$\large Z \simeq 6.249 $$ risulta leggermente inferiore a $2\pi$ a causa del troncamento dello spazio delle fasi, che esclude le code gaussiane della distribuzione.

Energia libera¶

L'energia libera di Helmholtz è definita come $$\large F = -\frac{1}{\beta} \ln Z. $$

Inserendo il valore numerico di $Z$ si ottiene $$\large F \simeq -1.83, $$

in accordo con il valore teorico $$\large F_\infty = -\ln(2\pi). $$

L'energia libera non rappresenta un'energia meccanica, ma quantifica l'equilibrio tra energia interna ed entropia, misurando il numero effettivo di microstati accessibili al sistema.

Entropia¶

Nell'insieme canonico l'entropia è data da $$\large S = \frac{E - F}{T}. $$

Poiché $T=1$, si ottiene $$\large S \simeq 2.83. $$

Nel contesto classico continuo, l'entropia dipende dal volume dello spazio delle fasi accessibile e non è definita in senso assoluto: essa presenta ambiguità additive che dipendono dalla scelta delle unità e dal troncamento dello spazio delle fasi.

Interpretazione fisica complessiva¶

I risultati ottenuti possono essere così interpretati:

  • l'energia interna verifica il teorema di equipartizione;
  • la funzione di partizione misura il volume di fase accessibile;
  • l'energia libera quantifica la ``facilità'' termodinamica dello stato;
  • l'entropia rappresenta il logaritmo del numero di microstati compatibili.

Il calcolo numerico realizza quindi una misura diretta del volume dello spazio delle fasi pesato dalla distribuzione di Boltzmann.

In [6]:
import numpy as np

def boltzmann_weight(H, beta):
    '''
    Calculate the Boltzmann probability distribution for a given Hamiltonian.
        :param H: Hamiltonian of the system.
        :param beta: Inverse temperature (1/kT).
    :return: 
        Probability according to Boltzmann distribution.
    '''
    return np.exp(-beta * H)

def example_hamiltonian(q, p):
    '''
    Example Hamiltonian function for demonstration purposes.
        :param q: Generalized coordinate.
        :param p: Generalized momentum.
    :return: 
        Hamiltonian value for 1D harmonic oscillatore.
    '''
    return 0.5 * (q**2 + p**2)

def canonical_monte_carlo(hamiltonian, beta, num_samples, qmax=5.0, pmax=5.0):
    '''
    Approximate the partition function using Monte Carlo sampling.
        :param hamiltonian_func: Hamiltonian function of the system.
        :param beta: Inverse temperature (1/kT).
        :param num_samples: Number of samples for approximation.
    :return: 
        Approximation of partition function Z, free Energy.
    '''
    weights = []
    energies = []

    for _ in range(num_samples):
        q = np.random.uniform(-qmax, qmax)
        p = np.random.uniform(-pmax, pmax)
        H = hamiltonian(q, p)
        w = boltzmann_weight(H, beta)

        weights.append(w)
        energies.append(H)

    weights = np.array(weights)
    energies = np.array(energies)

    phase_space_volume = (2*qmax) * (2*pmax)

    Z = phase_space_volume * np.mean(weights)
    E = np.sum(energies * weights) / np.sum(weights)

    return Z, E

def free_energy(Z, beta):
    '''
    Compute the free energy from the partition function.
        :param Z: Partition function.
        :param beta: Inverse temperature (1/kT).
    '''
    return - (1/beta) * np.log(Z)

def entropy(E, F, beta):
    '''
    Compute the entropy from the free energy.
        :param E: Internal energy.
        :param F: Free energy.
        :param beta: Inverse temperature (1/kT).
    :return: 
        Entropy S.
    '''
    T = 1 / beta
    return (E - F) / T


# ===== Example =====
beta = 1.0
num_samples = 200_000

Z, E = canonical_monte_carlo(example_hamiltonian, beta, num_samples)
F = free_energy(Z, beta)
S = entropy(E, F, beta)

print("Z =", Z)
print("E =", E)
print("F =", F)
print("S =", S)

Z = 6.248663432626204
E = 1.0010993376920274
F = -1.8323675900991643
S = 2.833466927791192
In [ ]: