1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
use rand::{
distributions::{Bernoulli, Distribution},
rngs::SmallRng,
};
use crate::{algorithm::Algorithm, lattice::Lattice};
/// The Metropolis single spin Ising algorithm.
pub struct Metropolis {
/// The probabilities of flipping an up spin depending on the number
/// of neighbouring up spins.
up_n_up_flip_dists: [Bernoulli; 5],
/// Same as `up_n_up_flip_dists` but for a down spin.
down_n_up_flip_dists: [Bernoulli; 5],
}
impl Metropolis {
#[must_use]
pub fn new(temperature: f64) -> Self {
let mut slf = Self {
// Initialize all probabilities with 1.0
// unwrap: 1.0 is obviously between 0 and 1.
up_n_up_flip_dists: [Bernoulli::new(1.0).unwrap(); 5],
down_n_up_flip_dists: [Bernoulli::new(1.0).unwrap(); 5],
};
// Set probabilities that are not 1.0
slf.set_temperature(temperature);
slf
}
}
impl Algorithm for Metropolis {
// Updates the temperature dependent probabilities of flipping a spin.
fn set_temperature(&mut self, val: f64) {
let minus_beta = -1.0 / val;
// unwrap: Probability is between 0 and 1 for positive temperatures.
let flip_dist_with_energy_diff_4 = Bernoulli::new((minus_beta * 4.0).exp()).unwrap();
let flip_dist_with_energy_diff_8 = Bernoulli::new((minus_beta * 8.0).exp()).unwrap();
// Central spin up:
// Neighbouring up: 0; energy diff: 2 * (+1) * (-4) = -8 <= 0
// Neighbouring up: 1; energy diff: 2 * (+1) * (-2) = -4 <= 0
// Neighbouring up: 2; energy diff: 2 * (+1) * (0) = 0 <= 0
// Neighbouring up: 3; energy diff: 2 * (+1) * (+2) = 4 > 0
self.up_n_up_flip_dists[3] = flip_dist_with_energy_diff_4;
// Neighbouring up: 4; energy diff: 2 * (+1) * (+4) = 8 > 0
self.up_n_up_flip_dists[4] = flip_dist_with_energy_diff_8;
// Central spin down:
// Neighbouring up: 0; energy diff: 2 * (-1) * (-4) = 8 > 0
self.down_n_up_flip_dists[0] = flip_dist_with_energy_diff_8;
// Neighbouring up: 1; energy diff: 2 * (-1) * (-2) = 4 > 0
self.down_n_up_flip_dists[1] = flip_dist_with_energy_diff_4;
// Neighbouring up: 2; energy diff: 2 * (-1) * (0) = 0 <= 0
// Neighbouring up: 3; energy diff: 2 * (-1) * (+2) = -4 <= 0
// Neighbouring up: 4; energy diff: 2 * (-1) * (+4) = -8 <= 0
}
// Picks a random spin and flips it with some probability.
fn step(&mut self, lattice: &mut Lattice, rng: &mut SmallRng) {
// Pick a random spin.
let (picked_spin, picked_spin_value) = lattice.rand_spin_with_value(rng);
// Get neighbouring spins.
let neighbouring_spins = lattice.neighbouring_spins(&picked_spin);
// Count neighbouring up spins.
// Safety: Neighbouring spins of a valid spin are also valid.
let n_neighbouring_spins_up = neighbouring_spins
.iter()
.filter(|&spin| unsafe { lattice.value_unchecked(spin) })
.count();
// Get the flipping distribution.
// Safety: Since there are only 4 neighbours, the number of neighbouring
// up spins can not be higher that 4.
let flip_dist = if picked_spin_value {
unsafe {
self.up_n_up_flip_dists
.get_unchecked(n_neighbouring_spins_up)
}
} else {
unsafe {
self.down_n_up_flip_dists
.get_unchecked(n_neighbouring_spins_up)
}
};
// Flip depending on the probability.
if flip_dist.sample(rng) {
// Safety: Random spin is always in bounds for the same lattice.
unsafe { lattice.set_value_unchecked(&picked_spin, !picked_spin_value) };
}
}
}