ising/
ising.rs

1use lattice_core::Graph;
2
3#[derive(Default)]
4struct Observable {
5    count: usize,
6    sum: f64,
7    sum_sq: f64,
8}
9
10impl Observable {
11    fn add(&mut self, value: f64) {
12        self.count += 1;
13        self.sum += value;
14        self.sum_sq += value * value;
15    }
16
17    fn mean(&self) -> f64 {
18        if self.count == 0 {
19            0.0
20        } else {
21            self.sum / self.count as f64
22        }
23    }
24
25    fn error(&self) -> f64 {
26        if self.count < 2 {
27            0.0
28        } else {
29            let n = self.count as f64;
30            let mean = self.mean();
31            let var = (self.sum_sq / n) - (mean * mean);
32            (var.max(0.0) / (n - 1.0)).sqrt()
33        }
34    }
35}
36
37struct Lcg {
38    state: u64,
39}
40
41impl Lcg {
42    fn new(seed: u64) -> Self {
43        Self { state: seed }
44    }
45
46    fn next_f64(&mut self) -> f64 {
47        self.state = self
48            .state
49            .wrapping_mul(6364136223846793005)
50            .wrapping_add(1);
51        let value = self.state >> 11;
52        (value as f64) / ((1u64 << 53) as f64)
53    }
54}
55
56fn main() {
57    println!("Metropolis Algorithm for Classical Ferromagnetic Ising Model");
58
59    let beta = 1.0 / 2.269;
60
61    let dim = 2usize;
62    let length = 32usize;
63    let graph = Graph::simple(dim, length);
64
65    let mut rng = Lcg::new(12345);
66    let mut spins = vec![1.0f64; graph.num_sites()];
67
68    let mut energy = Observable::default();
69    let mut mag2 = Observable::default();
70
71    let sweeps = 65536usize;
72    let therm = sweeps / 8;
73
74    for mcs in 0..(therm + sweeps) {
75        for site in 0..graph.num_sites() {
76            let mut diff = 0.0;
77            for k in 0..graph.num_neighbors(site) {
78                diff += 2.0 * spins[site] * spins[graph.neighbor(site, k)];
79            }
80            if rng.next_f64() < (-beta * diff).exp() {
81                spins[site] = -spins[site];
82            }
83        }
84
85        if mcs > therm {
86            let mut ene = 0.0;
87            for bond in 0..graph.num_bonds() {
88                ene -= spins[graph.source(bond)] * spins[graph.target(bond)];
89            }
90            energy.add(ene / graph.num_sites() as f64);
91
92            let mag = spins.iter().sum::<f64>() / graph.num_sites() as f64;
93            mag2.add(mag * mag);
94        }
95    }
96
97    println!("dimension = {dim}");
98    println!("length = {length}");
99    println!("beta = {beta}");
100    println!(
101        "energy density = {} +/- {}",
102        energy.mean(),
103        energy.error()
104    );
105    println!(
106        "magnetization density^2 = {} +/- {}",
107        mag2.mean(),
108        mag2.error()
109    );
110}
ising/ising.rs

ising/
ising.rs
