Project report: Isotta Magistrali, Gabriele Mura¶
In [75]:
import numba
import random
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pickle
import os
import random
from joblib import Parallel, delayed
import multiprocessing
np.random.seed(42)
3D Ising Model: simulation of magnetization VS temperature for diffrent grid sizes¶
1. Introduction to the problem¶
The Ising model is a physical model which consists in a simplified version of ferromagnetic materials. In the model we have a square grid on which N atoms, described by the value of their spins $S_i = \{+1,-1\}$, are displayed.
The system's energy is described by the following formula: $E(S) = -J\sum_{(i,j)}S_i \cdot S_j$
The probability of finding the system in a certain configuration is given by the Boltzmann distribution:
$P_B(S) = \frac{1}{Z}e^{-\frac{1}{k \cdot T}E(S)}$
For computational purposes we can rescale the temperature $\frac{k \cdot T}{J}$ -> $T$ so the expression above becomes simply:
$P_B(S) = \frac{1}{Z}e^{\frac{1}{T} \sum_{(i,j)}S_i \cdot S_j}$
In order to describe the system at equilibrium we might want to compute some thermodynamic averages like the average magnetization:
$\langle M\rangle = \sum_S \frac{1}{N} \sum_i S_i P_B(S)$
and see how they behave as a function of temperature (i.e. visualize the phase transition).
The sum over S runs over all the possible configurations of the systems, which is $2^N$. Thus, for N very large, it becomes impossible to compute those quantities analytically and instead we use Monte Carlo simulations: those allow us to sample configurations with the desired probability (the Boltzmann distribution) in order to be able to compute an "empirical average magnetization":
For a value $t$ of the temperature we generate $P$ possible configurations and compute the average magnetization as:
$ M_{emp} = \frac{1}{P}\sum_{p=1}^P \frac{1}{N} \sum_i S_i $
While for dimension $ \leq 2$ there is an analytical solution and we know that we can see a phase transition only for $D \geq 2$, for the 3D Ising model no analytical solution has been found as of today.
In this project we investigate the phase transition of the magnetization as a function of temperature in the 3D case, using Monte Carlo simulations with Metropolis algorithm for different grid sizes. In all our simulations we will consider spins located in a cubic grid and we will use the side length $L$ of the cube as a parameter for the size, with corresponding number of spins $N = L^3$.
2. State representation and problem visualization¶
We will start by implementing functions to compute the energy and the magnetization. Moreover, we will also implement a function computing the energy difference between two configurations differing by the value of one spin (it will be useful to speed up computations during later simulations).
In [21]:
def energy(state, side):
E = 0
for i in range(side):
for j in range(side):
for k in range(side):
E -= state[i, j, k] * (state[(i + 1) % side, j, k] +
state[i, (j + 1) % side, k] +
state[i, j, (k + 1) % side])
return E
def delta_energy(state, side, i, j, k):
## Change in energy after flipping spin in position i, j, k :
## all terns in the sum remain the same except for the one that include the flipped spin
## dE = Ef - Ei = s * neighbours_sum - (-s) * neighbours_sum) = 2 * s * neighbours_sum
neighbours_sum = (state[(i - 1) % side, j, k] +
state[(i + 1) % side, j, k] +
state[i, (j - 1) % side, k] +
state[i, (j + 1) % side, k] +
state[i, j, (k - 1) % side] +
state[i, j, (k + 1) % side])
return 2 * state[i, j, k] * neighbours_sum
In [22]:
def magnetization(state, side):
return np.sum(state) / (side**3)
Here we generate a random state and show it
In [73]:
side = 4
state = np.random.choice([-1, 1], size=(side, side, side))
# show the initial state with seaborn heatmap for each layer
fig, axs = plt.subplots(1, side)
for i in range(side):
sns.heatmap(state[:, :, i], ax=axs[i], cbar=False, square=True, xticklabels=False, yticklabels=False)
print(f"Energy: {energy(state, side)}")
print(f"Magnetization: {magnetization(state, side)}")
print('\n\n\t\tRandom initial state layer by layer')
Energy: 28 Magnetization: 0.03125 Random initial state layer by layer
In [74]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = np.meshgrid(np.arange(side), np.arange(side), np.arange(side))
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
ax.scatter(X, Y, Z, c=state.flatten(), s=50, marker='o')
plt.title('Random initial state')
plt.show()
3. Simulations¶
As stated above, the interest of our project is to show the behaviour of the magnetization as a function of temperature. In order to do so, we will implement various simulation algorithms that allow us to sample from the Boltzmann distribution. This will enable us to compute the empirical magnetization value which we can expect to be close to the true one.
Metropolis¶
In [76]:
@numba.jit(nopython=True)
def metropolis(state, side, t):
## spin to flip
i = np.random.randint(0, side)
j = np.random.randint(0, side)
k = np.random.randint(0, side)
s = state[i, j, k]
## change in energy
dE = delta_energy(state, side, i, j, k)
if dE < 0: ## move accepted
s = -s
elif np.random.random() < np.exp(- dE / t): ## move accepted
s = -s
state[i, j, k] = s
In [77]:
def plainMC(state, side, t, discard_steps = 100, nstep = 1000, stride = 100):
## m will be the average magnetization
m = 0
effective_steps = 0
# First steps to get into the right regime with burn-in
for _ in range(discard_steps):
metropolis(state, side, t)
# actual simulation steps
for istep in range(nstep):
metropolis(state, side, t)
# adding distance between samples
if istep % stride == 0:
m += magnetization(state, side)
effective_steps += 1
return m / effective_steps
Here we run the simulation for $L = \{5, 10, 20, 30\}$ and temperatures ranging from 2.5 to 5.5 with 0.25 intervals. In order to have good results we tried different parameters for the burn-in steps, stride and total steps of the Monte Carlo simulation. We faced a tradeoff between quality of results, obtained with higher values for the three parameters, and computational complexity, which increases very fast with larger values (especially for large $L$).
In order to keep the same quality level among the different simulations and achieve good results both for low and large $L$ values, we decided to model the number of steps as a function of the number of spins.
In [78]:
# running the simulation
sides = [5, 10, 20, 30]
temperatures = np.linspace(2.5, 5.5, 121)
def simulation(side, temperature):
print(f'Side: {side}, temperature: {round(temperature, 2)}')
state = np.random.choice([-1,1], size=(side, side, side))
# check if the simulation has already been run
if os.path.exists(f'/stoch_proj_data/magnetization_side_{side}_temperature_{temperature}.pkl'):
m = pickle.load(open(f'/stoch_proj_data/magnetization_side_{side}_temperature_{temperature}.pkl', 'rb'))
return m
# otherwise run the simulation
m = plainMC(state, side, temperature, discard_steps=10*(side**2), nstep=100*(side**2) if side<=20 else 20*(side**2), stride=side**2)
pickle.dump(m, open(f'/stoch_proj_data/magnetization_side_{side}_temperature_{temperature}.pkl', 'wb'))
return m
num_cores = multiprocessing.cpu_count()
results = Parallel(n_jobs=num_cores)(delayed(simulation)(side, temperature) for side in sides for temperature in temperatures)
Side: 5, temperature: 2.6 Side: 5, temperature: 2.52 Side: 5, temperature: 2.55 Side: 5, temperature: 2.58 Side: 5, temperature: 2.5 Side: 5, temperature: 2.65 Side: 5, temperature: 2.68 Side: 5, temperature: 2.62 Side: 5, temperature: 2.7 Side: 5, temperature: 2.72 Side: 5, temperature: 2.75 Side: 5, temperature: 2.78 Side: 5, temperature: 2.8 Side: 5, temperature: 2.82 Side: 5, temperature: 2.85 Side: 5, temperature: 2.88 Side: 5, temperature: 2.9 Side: 5, temperature: 2.92 Side: 5, temperature: 2.95 Side: 5, temperature: 2.98 Side: 5, temperature: 3.0 Side: 5, temperature: 3.02 Side: 5, temperature: 3.05 Side: 5, temperature: 3.08 Side: 5, temperature: 3.1 Side: 5, temperature: 3.12 Side: 5, temperature: 3.15 Side: 5, temperature: 3.18 Side: 5, temperature: 3.2 Side: 5, temperature: 3.22 Side: 5, temperature: 3.25 Side: 5, temperature: 3.28 Side: 5, temperature: 3.3 Side: 5, temperature: 3.32 Side: 5, temperature: 3.35 Side: 5, temperature: 3.38 Side: 5, temperature: 3.4 Side: 5, temperature: 3.42 Side: 5, temperature: 3.45 Side: 5, temperature: 3.48 Side: 5, temperature: 3.5 Side: 5, temperature: 3.53 Side: 5, temperature: 3.55 Side: 5, temperature: 3.58 Side: 5, temperature: 3.6 Side: 5, temperature: 3.62 Side: 5, temperature: 3.65 Side: 5, temperature: 3.68 Side: 5, temperature: 3.7 Side: 5, temperature: 3.72 Side: 5, temperature: 3.75 Side: 5, temperature: 3.78 Side: 5, temperature: 3.8 Side: 5, temperature: 3.82 Side: 5, temperature: 3.85 Side: 5, temperature: 3.88 Side: 5, temperature: 3.9 Side: 5, temperature: 3.92 Side: 5, temperature: 3.95 Side: 5, temperature: 3.98 Side: 5, temperature: 4.0 Side: 5, temperature: 4.03 Side: 5, temperature: 4.05 Side: 5, temperature: 4.08 Side: 5, temperature: 4.1 Side: 5, temperature: 4.12 Side: 5, temperature: 4.15 Side: 5, temperature: 4.18 Side: 5, temperature: 4.2 Side: 5, temperature: 4.22 Side: 5, temperature: 4.25 Side: 5, temperature: 4.28 Side: 5, temperature: 4.3 Side: 5, temperature: 4.32 Side: 5, temperature: 4.35 Side: 5, temperature: 4.38 Side: 5, temperature: 4.4 Side: 5, temperature: 4.42 Side: 5, temperature: 4.45 Side: 5, temperature: 4.47 Side: 5, temperature: 4.5 Side: 5, temperature: 4.53 Side: 5, temperature: 4.55 Side: 5, temperature: 4.58 Side: 5, temperature: 4.6 Side: 5, temperature: 4.62 Side: 5, temperature: 4.65 Side: 5, temperature: 4.68 Side: 5, temperature: 4.7 Side: 5, temperature: 4.72 Side: 5, temperature: 4.75 Side: 5, temperature: 4.78 Side: 5, temperature: 4.8 Side: 5, temperature: 4.82 Side: 5, temperature: 4.85 Side: 5, temperature: 4.88 Side: 5, temperature: 4.9 Side: 5, temperature: 4.93 Side: 5, temperature: 4.95 Side: 5, temperature: 4.97 Side: 5, temperature: 5.0 Side: 5, temperature: 5.03 Side: 5, temperature: 5.05 Side: 5, temperature: 5.08 Side: 5, temperature: 5.1 Side: 5, temperature: 5.12 Side: 5, temperature: 5.15 Side: 5, temperature: 5.18 Side: 5, temperature: 5.2 Side: 5, temperature: 5.22 Side: 5, temperature: 5.25 Side: 5, temperature: 5.28 Side: 5, temperature: 5.3 Side: 5, temperature: 5.32 Side: 5, temperature: 5.35 Side: 5, temperature: 5.38 Side: 5, temperature: 5.4 Side: 5, temperature: 5.43 Side: 5, temperature: 5.45 Side: 5, temperature: 5.48 Side: 5, temperature: 5.5 Side: 10, temperature: 2.5 Side: 10, temperature: 2.52 Side: 10, temperature: 2.55 Side: 10, temperature: 2.58 Side: 10, temperature: 2.6 Side: 10, temperature: 2.62 Side: 10, temperature: 2.65 Side: 10, temperature: 2.68 Side: 10, temperature: 2.7 Side: 10, temperature: 2.72 Side: 10, temperature: 2.75 Side: 10, temperature: 2.78 Side: 10, temperature: 2.8 Side: 10, temperature: 2.82 Side: 10, temperature: 2.85 Side: 10, temperature: 2.88 Side: 10, temperature: 2.9 Side: 10, temperature: 2.98 Side: 10, temperature: 2.92 Side: 10, temperature: 3.0 Side: 10, temperature: 2.95 Side: 10, temperature: 3.02 Side: 10, temperature: 3.05 Side: 10, temperature: 3.08 Side: 10, temperature: 3.1 Side: 10, temperature: 3.12 Side: 10, temperature: 3.15 Side: 10, temperature: 3.18 Side: 10, temperature: 3.2 Side: 10, temperature: 3.22 Side: 10, temperature: 3.28 Side: 10, temperature: 3.3 Side: 10, temperature: 3.25 Side: 10, temperature: 3.32 Side: 10, temperature: 3.35 Side: 10, temperature: 3.38 Side: 10, temperature: 3.4Side: 10, temperature: 3.48 Side: 10, temperature: 3.5 Side: 10, temperature: 3.42 Side: 10, temperature: 3.53 Side: 10, temperature: 3.45 Side: 10, temperature: 3.55 Side: 10, temperature: 3.68 Side: 10, temperature: 3.58 Side: 10, temperature: 3.7 Side: 10, temperature: 3.6 Side: 10, temperature: 3.72 Side: 10, temperature: 3.62 Side: 10, temperature: 3.75 Side: 10, temperature: 3.65 Side: 10, temperature: 3.78 Side: 10, temperature: 3.8 Side: 10, temperature: 3.88 Side: 10, temperature: 3.9 Side: 10, temperature: 3.82 Side: 10, temperature: 3.92 Side: 10, temperature: 3.85 Side: 10, temperature: 3.95 Side: 10, temperature: 3.98 Side: 10, temperature: 4.0 Side: 10, temperature: 4.08 Side: 10, temperature: 4.1 Side: 10, temperature: 4.03 Side: 10, temperature: 4.05 Side: 10, temperature: 4.12 Side: 10, temperature: 4.15 Side: 10, temperature: 4.28 Side: 10, temperature: 4.18 Side: 10, temperature: 4.2 Side: 10, temperature: 4.3 Side: 10, temperature: 4.32 Side: 10, temperature: 4.22 Side: 10, temperature: 4.35 Side: 10, temperature: 4.38 Side: 10, temperature: 4.25 Side: 10, temperature: 4.4 Side: 10, temperature: 4.42 Side: 10, temperature: 4.47 Side: 10, temperature: 4.45 Side: 10, temperature: 4.5 Side: 10, temperature: 4.53 Side: 10, temperature: 4.55 Side: 10, temperature: 4.58 Side: 10, temperature: 4.6 Side: 10, temperature: 4.62 Side: 10, temperature: 4.68 Side: 10, temperature: 4.65 Side: 10, temperature: 4.7 Side: 10, temperature: 4.72 Side: 10, temperature: 4.75 Side: 10, temperature: 4.78 Side: 10, temperature: 4.8 Side: 10, temperature: 4.82 Side: 10, temperature: 4.85 Side: 10, temperature: 4.88 Side: 10, temperature: 4.9 Side: 10, temperature: 4.93 Side: 10, temperature: 4.95 Side: 10, temperature: 4.97 Side: 10, temperature: 5.0 Side: 10, temperature: 5.08 Side: 10, temperature: 5.03 Side: 10, temperature: 5.05 Side: 10, temperature: 5.1 Side: 10, temperature: 5.12 Side: 10, temperature: 5.15 Side: 10, temperature: 5.18 Side: 10, temperature: 5.28 Side: 10, temperature: 5.2 Side: 10, temperature: 5.3 Side: 10, temperature: 5.32Side: 10, temperature: 5.22 Side: 10, temperature: 5.35 Side: 10, temperature: 5.38 Side: 10, temperature: 5.25 Side: 10, temperature: 5.4 Side: 10, temperature: 5.43 Side: 10, temperature: 5.48 Side: 10, temperature: 5.45 Side: 10, temperature: 5.5 Side: 20, temperature: 2.5 Side: 20, temperature: 2.52 Side: 20, temperature: 2.55 Side: 20, temperature: 2.65 Side: 20, temperature: 2.58Side: 20, temperature: 2.68 Side: 20, temperature: 2.7 Side: 20, temperature: 2.6 Side: 20, temperature: 2.62 Side: 20, temperature: 2.72 Side: 20, temperature: 2.75 Side: 20, temperature: 2.85 Side: 20, temperature: 2.88 Side: 20, temperature: 2.78 Side: 20, temperature: 2.8 Side: 20, temperature: 2.9 Side: 20, temperature: 2.92 Side: 20, temperature: 2.82 Side: 20, temperature: 2.95 Side: 20, temperature: 2.98 Side: 20, temperature: 3.05 Side: 20, temperature: 3.08 Side: 20, temperature: 3.0 Side: 20, temperature: 3.02 Side: 20, temperature: 3.1 Side: 20, temperature: 3.12 Side: 20, temperature: 3.15 Side: 20, temperature: 3.18 Side: 20, temperature: 3.2 Side: 20, temperature: 3.25 Side: 20, temperature: 3.22 Side: 20, temperature: 3.28 Side: 20, temperature: 3.3 Side: 20, temperature: 3.32 Side: 20, temperature: 3.35 Side: 20, temperature: 3.38 Side: 20, temperature: 3.4 Side: 20, temperature: 3.65 Side: 20, temperature: 3.42 Side: 20, temperature: 3.45 Side: 20, temperature: 3.68 Side: 20, temperature: 3.48 Side: 20, temperature: 3.7 Side: 20, temperature: 3.5 Side: 20, temperature: 4.05 Side: 20, temperature: 3.53 Side: 20, temperature: 3.72 Side: 20, temperature: 3.55 Side: 20, temperature: 4.45 Side: 20, temperature: 3.58 Side: 20, temperature: 3.75 Side: 20, temperature: 3.6 Side: 20, temperature: 4.47 Side: 20, temperature: 3.62 Side: 20, temperature: 4.08 Side: 20, temperature: 4.5 Side: 20, temperature: 3.78 Side: 20, temperature: 4.53 Side: 20, temperature: 4.1 Side: 20, temperature: 3.8 Side: 20, temperature: 4.85 Side: 20, temperature: 4.12 Side: 20, temperature: 4.55 Side: 20, temperature: 3.82 Side: 20, temperature: 4.15 Side: 20, temperature: 3.85 Side: 20, temperature: 4.88 Side: 20, temperature: 4.58 Side: 20, temperature: 3.88 Side: 20, temperature: 4.9 Side: 20, temperature: 4.6 Side: 20, temperature: 4.18 Side: 20, temperature: 3.9 Side: 20, temperature: 4.2 Side: 20, temperature: 4.93 Side: 20, temperature: 4.62 Side: 20, temperature: 4.95 Side: 20, temperature: 5.25 Side: 20, temperature: 4.65Side: 20, temperature: 3.92 Side: 20, temperature: 4.22 Side: 20, temperature: 5.28 Side: 20, temperature: 4.97Side: 20, temperature: 3.95 Side: 20, temperature: 4.25 Side: 20, temperature: 5.3 Side: 20, temperature: 4.28Side: 20, temperature: 5.0 Side: 20, temperature: 3.98Side: 20, temperature: 4.68 Side: 20, temperature: 5.32 Side: 20, temperature: 4.3 Side: 20, temperature: 4.7Side: 20, temperature: 5.03 Side: 20, temperature: 4.0 Side: 20, temperature: 4.72 Side: 20, temperature: 4.32 Side: 20, temperature: 5.05 Side: 20, temperature: 5.35 Side: 20, temperature: 4.03 Side: 20, temperature: 4.75 Side: 20, temperature: 4.35 Side: 20, temperature: 5.38 Side: 20, temperature: 5.08 Side: 20, temperature: 4.78 Side: 20, temperature: 4.38 Side: 30, temperature: 2.62 Side: 20, temperature: 4.8 Side: 20, temperature: 5.1 Side: 20, temperature: 5.4 Side: 20, temperature: 4.4 Side: 20, temperature: 4.82 Side: 20, temperature: 5.43 Side: 20, temperature: 5.12 Side: 20, temperature: 5.45Side: 30, temperature: 2.65 Side: 20, temperature: 4.42 Side: 20, temperature: 5.48 Side: 20, temperature: 5.15 Side: 30, temperature: 2.68 Side: 20, temperature: 5.18 Side: 20, temperature: 5.5 Side: 30, temperature: 3.02 Side: 30, temperature: 2.5 Side: 20, temperature: 5.2 Side: 30, temperature: 2.7 Side: 20, temperature: 5.22 Side: 30, temperature: 2.52 Side: 30, temperature: 2.72 Side: 30, temperature: 3.05 Side: 30, temperature: 2.55 Side: 30, temperature: 2.75 Side: 30, temperature: 3.08 Side: 30, temperature: 2.58 Side: 30, temperature: 2.78 Side: 30, temperature: 2.6 Side: 30, temperature: 2.8 Side: 30, temperature: 3.42 Side: 30, temperature: 3.1 Side: 30, temperature: 2.82 Side: 30, temperature: 2.85 Side: 30, temperature: 3.12 Side: 30, temperature: 2.88 Side: 30, temperature: 3.15 Side: 30, temperature: 3.45 Side: 30, temperature: 2.9 Side: 30, temperature: 3.18 Side: 30, temperature: 2.92 Side: 30, temperature: 3.2 Side: 30, temperature: 2.95 Side: 30, temperature: 2.98 Side: 30, temperature: 3.22 Side: 30, temperature: 3.48 Side: 30, temperature: 3.25 Side: 30, temperature: 3.0 Side: 30, temperature: 3.28 Side: 30, temperature: 3.3 Side: 30, temperature: 3.5 Side: 30, temperature: 3.32 Side: 30, temperature: 3.35 Side: 30, temperature: 3.38 Side: 30, temperature: 3.4 Side: 30, temperature: 3.53 Side: 30, temperature: 3.55 Side: 30, temperature: 3.58 Side: 30, temperature: 3.6 Side: 30, temperature: 3.62 Side: 30, temperature: 3.65 Side: 30, temperature: 3.68 Side: 30, temperature: 3.7 Side: 30, temperature: 3.72 Side: 30, temperature: 3.75 Side: 30, temperature: 3.78 Side: 30, temperature: 3.8 Side: 30, temperature: 3.82 Side: 30, temperature: 3.85 Side: 30, temperature: 3.88 Side: 30, temperature: 3.9 Side: 30, temperature: 3.92 Side: 30, temperature: 3.95 Side: 30, temperature: 3.98 Side: 30, temperature: 4.0 Side: 30, temperature: 4.03 Side: 30, temperature: 4.05 Side: 30, temperature: 4.08 Side: 30, temperature: 4.1 Side: 30, temperature: 4.12 Side: 30, temperature: 4.15 Side: 30, temperature: 4.18 Side: 30, temperature: 4.2 Side: 30, temperature: 4.22 Side: 30, temperature: 4.25 Side: 30, temperature: 4.28 Side: 30, temperature: 4.3 Side: 30, temperature: 4.32 Side: 30, temperature: 4.35 Side: 30, temperature: 4.38 Side: 30, temperature: 4.4 Side: 30, temperature: 4.42 Side: 30, temperature: 4.45 Side: 30, temperature: 4.47 Side: 30, temperature: 4.5 Side: 30, temperature: 4.53 Side: 30, temperature: 4.55 Side: 30, temperature: 4.58 Side: 30, temperature: 4.6 Side: 30, temperature: 4.62 Side: 30, temperature: 4.65 Side: 30, temperature: 4.68 Side: 30, temperature: 4.7 Side: 30, temperature: 4.72 Side: 30, temperature: 4.75 Side: 30, temperature: 4.78 Side: 30, temperature: 4.8 Side: 30, temperature: 4.82 Side: 30, temperature: 4.85 Side: 30, temperature: 4.88 Side: 30, temperature: 4.9 Side: 30, temperature: 4.93 Side: 30, temperature: 4.95 Side: 30, temperature: 4.97 Side: 30, temperature: 5.0 Side: 30, temperature: 5.03 Side: 30, temperature: 5.05 Side: 30, temperature: 5.08 Side: 30, temperature: 5.1 Side: 30, temperature: 5.12 Side: 30, temperature: 5.15 Side: 30, temperature: 5.18 Side: 30, temperature: 5.2 Side: 30, temperature: 5.22 Side: 30, temperature: 5.25 Side: 30, temperature: 5.28 Side: 30, temperature: 5.3 Side: 30, temperature: 5.32 Side: 30, temperature: 5.35 Side: 30, temperature: 5.38 Side: 30, temperature: 5.4 Side: 30, temperature: 5.43 Side: 30, temperature: 5.45 Side: 30, temperature: 5.48 Side: 30, temperature: 5.5
In [87]:
# plotting the results
results = []
for side in sides:
magnetization_per_t = []
for temperature in temperatures:
m = pickle.load(open(f'/stoch_proj_data/magnetization_side_{side}_temperature_{temperature}.pkl', 'rb'))
magnetization_per_t.append(m)
results.append(magnetization_per_t)
fig, ax = plt.subplots(2, 2, figsize=(12, 8))
for i, side in enumerate(sides):
ax[i//2, i%2].scatter(temperatures, results[i])
ax[i//2, i%2].set_title(f'Side = {side}')
ax[i//2, i%2].set_xlabel('Temperature')
ax[i//2, i%2].set_ylabel('Magnetization')
plt.tight_layout()
plt.show()
From the simulations above we can extrapolate two final takeaways:
- as expected the critical temperature at which there is a phase transition is $T_c$ ~ $4.5$
- for low $L$ values, the deviations from the theoretical predictions are likely due to finite-size effects. Simulations get more and more precise as we increase the number of spins N i.e. the side of the grid. Indeed, analytically, we always consider $\lim_{N\to\infty}$
We can see that for $L = 30$ we have some outliers, but this is due to the tradeoff between accuracy of results and computational complexity that we described above, indicating that a larger number of simulation steps was needed to achieve a perfect result.
Further implementations¶
We also decided to implement different algorithms that were proven to improve convergence, but we did not achieve results as good as the ones shown above. Here, we just state the names and our implementations.
In [ ]:
# Wolff algorithm
def wolff(state, side, t):
i = np.random.randint(0, side)
j = np.random.randint(0, side)
k = np.random.randint(0, side)
cluster = [(i, j, k)]
s = state[i, j, k]
stack = [(i, j, k)]
while len(stack) > 0:
i, j, k = stack.pop()
i_plus = (i + 1) if (i + 1) < side else 0
j_plus = (j + 1) if (j + 1) < side else 0
k_plus = (k + 1) if (k + 1) < side else 0
i_minus = (i - 1) if (i - 1) >= 0 else (side - 1)
j_minus = (j - 1) if (j - 1) >= 0 else (side - 1)
k_minus = (k - 1) if (k - 1) >= 0 else (side - 1)
neighbours = [(i_plus, j, k), (i, j_plus, k), (i, j, k_plus), (i_minus, j, k), (i, j_minus, k), (i, j, k_minus)]
for n in neighbours:
if state[n] == s and n not in cluster and np.random.random() < 1 - np.exp(-2/t):
stack.append(n)
cluster.append(n)
for i, j, k in cluster:
state[i, j, k] = -state[i, j, k]
def wolff_simulation(state, side, t, discard_steps = 100, nstep = 1000, stride = 100):
m = 0
effective_steps = 0
for _ in range(discard_steps):
wolff(state, side, t)
for istep in range(nstep):
wolff(state, side, t)
if istep % stride == 0:
m += magnetization(state, side)
effective_steps += 1
return m / effective_steps
In [ ]:
# Simulated tempering
def metropolis_t(e, temperatures, g, it_old, it_new):
delta_beta = 1./temperatures[it_new] - 1./temperatures[it_new]
if e * delta_beta - g[it_new] + g[it_old] < 0: return True
else:
if np.random.random() < np.exp( -delta_beta * e + g[it_new] - g[it_old]): return True
return False
class MonteCarlo3D_temp:
def __init__(self, L, initial_configuration, temperatures, discard_steps = 500, nstep = 5000, stride=1):
self.L = L
self.discard_steps = discard_steps
self.nstep = nstep
self.volume = self.L ** 3
self.stride = stride
self.temperatures = temperatures
self.initial_configuration = initial_configuration
# simulated tempering function
def tempering(self, g=None, verbose=True):
self.g = np.zeros(len(self.temperatures)) if g is None else g
state = self.initial_configuration.copy()
e_old = energy(state)
e_history = []
state_history = []
temperature_history = []
accepted = 0
accepted_t = 0
n_temp = len(self.temperatures)
it_old = n_temp - 1 # present temperature index
# loop over steps
for istep in range(self.nstep):
# choose movetype (0=spin, 1=temperature)
move_type = np.random.randint(2)
# attempt move
if move_type == 0: # move spin
outcome = metropolis(state, self.L, self.temperatures[it_old])
e_new = energy(state)
it_new = it_old
else: # move temperature
di = np.random.randint( 3 ) - 1
it_new = it_old + di
if it_new < 0: it_new = 0
if it_new > n_temp - 1: it_new = n_temp - 1
e_new = e_old
outcome = metropolis_t(e_old, self.temperatures, self.g, it_old, it_new)
if outcome:
e_old = e_new
it_old = it_new
if move_type == 0: accepted = accepted + 1
else: accepted_t = accepted_t + 1
# output
e_history.append( e_old )
if verbose:
if istep % self.stride == 0: print( istep, e_old, self.temperatures[it_old] )
if istep % self.stride == 0:
state_history.append( state.copy() ) # .copy otherwise it stores the pointer to the same array
temperature_history.append( self.temperatures[it_old] )
print("Accepted spin = "+str(accepted))
print("Accepted temp = "+str(accepted_t))
m_per_t = []
for my_t in self.temperatures:
my_states = [ state for t, state in zip(temperature_history, state_history) if t == my_t ]
print(f"n of states selected per temp {round(my_t,2)} = {len(my_states)}")
m_states = [magnetization(state) for state in my_states]
m_per_t.append(np.mean(m_states))
return m_per_t
# running the simulation
side = 5
temperatures = np.linspace(2.5,5.5,61)
temperature_sets = dict()
m_final = dict()
for i,j in zip(range(0,12,3),range(4)):
temperature_sets[i] = temperatures[j:j+17:4]
temperature_sets[i+1] = temperatures[j+20:j+20+17:4]
temperature_sets[i+2] = temperatures[j+40:j+40+17:4]
initial_configuration = np.random.choice([-1,1], size=(side, side, side))
for i, t_set in temperature_sets.items():
print(f'Temperature set {i+1}/{len(temperature_sets)}')
ising = MonteCarlo3D_temp(side, initial_configuration, t_set, discard_steps=500*(side**3), nstep=5000*(side**3), stride=side**3)
m_per_t = ising.tempering(verbose=False)
m_final[i] = m_per_t
In [ ]:
# parallel tempering
class IsingModel3D:
def __init__(self, L, discard_steps = 1000, nstep = 5000, stride=1):
self.L = L
self.discard_steps = discard_steps
self.nstep = nstep
self.volume = self.L ** 3
self.stride = stride
self.t_energy = dict()
self.t_state = dict()
def _initial_configuration(self):
state = np.random.choice([-1,1], size=(self.L, self.L, self.L))
return state
def parallel_tempering(self, energies, temperatures, spin_configs):
ret_energies = energies.copy()
ret_spins = spin_configs.copy()
for i in range(len(temperatures)-1):
delta = (1/temperatures[i] - 1/temperatures[i+1])*(energies[i] - energies[i+1])
if np.random.random() < np.exp(delta):
temp_e = ret_energies[i]
temp_s = ret_spins[i]
ret_energies[i] = ret_energies[i+1]
ret_energies[i+1] = temp_e
ret_spins[i] = ret_spins[i+1]
ret_spins[i+1] = temp_s
return ret_energies, ret_spins
# Monte Carlo CPU version
def plainMC(self, t, state):
# Initialize the lattice randomly
e_history = []
state_history = []
# Monte Carlo steps
for istep in range(self.nstep):
for _ in range(side ** 2):
metropolis(state, self.L, t)
E = energy(state)
e_history.append(E)
state_history.append(state.copy())
return e_history, state_history
def simulation(self, temperatures, ptsteps=100):
states = [self._initial_configuration() for _ in range(len(temperatures))]
for t in temperatures:
self.t_energy[t] = []
self.t_state[t] = []
for i, t in enumerate(temperatures):
for _ in range(self.discard_steps):
for _ in range(side ** 2):
metropolis(states[i], self.L, t)
for mc in range(ptsteps):
print(f'Parallel tempering {mc+1}/{ptsteps}')
for i, t in enumerate(temperatures):
print(f' Temperature: {t}')
e_history, state_history = self.plainMC(t, states[i])
self.t_energy[t] += e_history
self.t_state[t] += state_history
energies = [e[-1] for e in list(self.t_energy.values())]
spin_configs = [s[-1] for s in list(self.t_state.values())]
changed_order_e, changed_order_s = self.parallel_tempering(energies, temperatures, spin_configs)
states = changed_order_s
for i,t in enumerate(temperatures):
self.t_energy[t][-1] = changed_order_e[i]
self.t_state[t][-1] = changed_order_s[i].copy()
m_per_t = []
for t in temperatures:
states_at_t = self.t_state[t]
m_per_state = [magnetization(state) for state in states_at_t]
m_per_t.append(np.mean(m_per_state))
return m_per_t