Commit e655603f authored by Valerio's avatar Valerio
Browse files

notebook ex 6

parent dcd57c15
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from timeit import default_timer as timer\n",
"import random\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#does the next level of binning\n",
"def binning_step(data):\n",
" newdata = []\n",
" for i in range(int(len(data)/2)):\n",
" newdata.append((data[2*i]+data[2*i+1])/2.)\n",
" stdev = np.std(newdata)/np.sqrt(len(newdata)-1)\n",
" return [stdev,newdata]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#perform a binning analysis of the data, up to the level where at least minsamples samples remain\n",
"def binning(data, minsamples=8):\n",
"\tlevels = []\n",
"\tstdevs = []\n",
"\tstdevs.append(np.std(data)/np.sqrt(len(data)-1))\n",
"\tlevels.append(0)\n",
"\tlevel = 1\n",
"\twhile(len(data) >= 2*minsamples):\n",
"\t\tstdev,data = binning_step(data)\n",
"\t\tstdevs.append(stdev)\n",
"\t\tlevels.append(level)\n",
"\t\tlevel += 1\n",
"\treturn [levels, stdevs]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#perform a Wolff cluster update step on sys of size H*L\n",
"#absolute value of return is the cluster size, twice the return is the change in magnetisation\n",
"def Wolff_step(sys,H,L,probability):\n",
"\tflipped = 1\n",
"\t#choose random spin\n",
"\tspin_v = random.randint(0,H-1)\n",
"\tspin_h = random.randint(0,L-1)\n",
"\t#flip the spin\n",
"\tsys[spin_v][spin_h] *= -1\n",
"\t#keep track of what we do to the magnetisation\n",
"\tflip_sign = sys[spin_v][spin_h]\n",
"\t#check neighbours recursively\n",
"\tspins_to_check = [[spin_v,spin_h]]\n",
"\twhile(len(spins_to_check) > 0):\n",
"\t\tindices = spins_to_check.pop()\n",
"\t\tspin_v = indices[0]\n",
"\t\tspin_h = indices[1]\n",
"\t\tcurrent_spin = sys[spin_v][spin_h]\n",
"\t\t#add neighbours to cluster and schedule them with correct probability\n",
"\t\tif(current_spin * sys[(spin_v+1)%H][spin_h] == -1 and random.random() < probability):\n",
"\t\t\tspins_to_check.append([(spin_v+1)%H,spin_h])\n",
"\t\t\tsys[(spin_v+1)%H][spin_h] *= -1\n",
"\t\t\tflipped += 1\n",
"\n",
"\t\tif(current_spin * sys[(spin_v-1)%H][spin_h] == -1 and random.random() < probability):\n",
"\t\t\tspins_to_check.append([(spin_v-1)%H,spin_h])\n",
"\t\t\tsys[(spin_v-1)%H][spin_h] *= -1\n",
"\t\t\tflipped += 1\n",
"\n",
"\t\tif(current_spin * sys[spin_v][(spin_h+1)%L] == -1 and random.random() < probability):\n",
"\t\t\tspins_to_check.append([spin_v,(spin_h+1)%L])\n",
"\t\t\tsys[spin_v][(spin_h+1)%L] *= -1\n",
"\t\t\tflipped += 1\n",
"\n",
"\t\tif(current_spin * sys[spin_v][(spin_h-1)%L] == -1 and random.random() < probability):\n",
"\t\t\tspins_to_check.append([spin_v,(spin_h-1)%L])\n",
"\t\t\tsys[spin_v][(spin_h-1)%L] *= -1\n",
"\t\t\tflipped += 1\n",
"\treturn flipped*flip_sign\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"#perform one single spin update on sys of extent H*L. \n",
"#returns change in magnetisation\n",
"def MCMC_step(sys,H,L,exptable):\n",
"\t#pick random spin\n",
"\tspin_v = random.randint(0,H-1)\n",
"\tspin_h = random.randint(0,L-1)\n",
"\t#calculate local field\n",
"\th = sys[spin_v][(spin_h-1)%L] + sys[spin_v][(spin_h+1)%L] + sys[(spin_v+1)%H][spin_h] + sys[(spin_v-1)%H][spin_h] \n",
"\t\n",
"\t#flip with metropolis probability\n",
"\tif(random.random() < exptable[sys[spin_v][spin_h]*h]):\n",
"\t\tif(sys[spin_v][spin_h] == 1):\n",
"\t\t\tsys[spin_v][spin_h] = -1\n",
"\t\t\treturn -2\n",
"\t\telse:\n",
"\t\t\tsys[spin_v][spin_h] = 1\n",
"\t\t\treturn 2\n",
"\telse:\n",
"\t\treturn 0\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"#returns the magnetisation of sys\n",
"def magnetisation(sys):\n",
"\tmagn = 0\n",
"\tL = len(sys[0])\n",
"\tH = len(sys)\n",
"\tfor a in sys:\n",
"\t\tfor s in a:\n",
"\t\t\tmagn += s\n",
"\treturn float(magn)/float(H*L)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"#run a simulation using cluster updates\n",
"def run_cluster_simulation(dT,Tmin=1.8,Tmax=3.2,meas=8192*2,equil=1000*2):\n",
"\n",
"\t#######################\n",
"\t###\tInitialisation\n",
"\t#######################\n",
"\n",
"\t#arrays to store info about correlations\n",
"\tstatistics_magnetisation = []\n",
"\tstatistics_magnetisationq = []\n",
"\n",
"\t#arrays to store the expectation values and their standard deviations for all temperatures\n",
"\tmagnetisations = []\n",
"\tmagnetisationsq = []\n",
"\tmagnetisations_std = []\n",
"\tmagnetisationsq_std = []\n",
"\ttemperatures = []\n",
"\n",
"\t#define system parameters\n",
"\tL = 16\n",
"\tH = 16\n",
"\tJ = 1.\n",
"\t#initialise the system\n",
"\tsys = []\n",
"\tfor i in range(H):\n",
"\t\ta = []\n",
"\t\tfor j in range(L):\n",
"\t\t\ta.append(1)\n",
"\t\tsys.append(a)\n",
"\t#initialise the magnetisation\n",
"\ttotal_magnetisation = L*H\n",
"\n",
"\tT = Tmin\n",
"\n",
"\t\n",
"\t#######################\n",
"\t###\tSimulation\n",
"\t#######################\n",
"\t#simulate for the specified temperatures\n",
"\twhile(T <= Tmax):\n",
"\t\t#arrays to store all samples\n",
"\t\tmagdata = [] #magnetisation\n",
"\t\tmagqdata = [] #magnetisation squared\n",
"\t\tT += dT\n",
"\t\tbeta = 1./T\n",
"\t\tprobability = 1. - np.exp(-2.*beta) #probability of accepting the addition of a site to the cluster\n",
"\n",
"\t\t#equilibrate\n",
"\t\tfor i in range(int(equil)):\n",
"\t\t\ttotal_magnetisation += 2*Wolff_step(sys,H,L,probability)\n",
"\t\t#measure\n",
"\t\tfor i in range(int(meas)):\n",
"\t\t\twolff_return = Wolff_step(sys,H,L,probability)\n",
"\t\t\tm = np.abs(wolff_return)/float(L*H)\n",
"\t\t\ttotal_magnetisation += 2*wolff_return\n",
"\t\t\tmagdata.append(m)\n",
"\t\t\tmagqdata.append((total_magnetisation/float(L*H))**2)\n",
"\n",
"\t\t#extract means and their standard deviations\n",
"\t\tmagnetisations.append(np.mean(magdata))\n",
"\t\tmagnetisationsq.append(np.mean(magqdata))\n",
"\t\tmagnetisations_std.append(np.std(magdata))\n",
"\t\tmagnetisationsq_std.append(np.std(magqdata))\n",
"\n",
"\t\t#save temperature\n",
"\t\ttemperatures.append(T)\n",
"\n",
"\t\t#do a binning analysis of the data and store it\n",
"\t\tstatistics_magnetisation.append(binning(magdata))\n",
"\t\tstatistics_magnetisationq.append(binning(magqdata))\n",
"\n",
"\t#calculate correlation times for all temperatures from the binning results. binning level 5 seems reasonable for most temperatures\n",
"\ttaus = []\n",
"\ttausq = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\ttaus.append(0.5*(statistics_magnetisation[i][1][5]**2/magnetisations_std[i]**2*(meas-1.)-1.))\n",
"\t\ttausq.append(0.5*(statistics_magnetisationq[i][1][5]**2/magnetisationsq_std[i]**2*(meas-1.)-1.))\n",
"\n",
"\n",
"\n",
"\t#######################\n",
"\t###\tPlotting\n",
"\t#######################\n",
"\n",
"\t#correlation time magnetisation\n",
"\tplt.plot(temperatures,taus)\n",
"\tplt.title(\"Correlation time Wolff, magnetization\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Correlation time\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_corr_times_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#correlation time magnetisation squared\n",
"\tplt.plot(temperatures,taus)\n",
"\tplt.title(\"Correlation time Wolff, magnetization$^2$\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Correlation time\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_corr_times_mq.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#binning analysis results\n",
"\tfor i in range(len(statistics_magnetisation)):\n",
"\t\tplt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1])\n",
"\tleg = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\tleg.append(\"T=\" + str(round(temperatures[i],3)))\n",
"\tplt.legend(leg)\n",
"\tplt.title(\"Binning analysis Wolff, magnetisation\")\n",
"\tplt.xlabel(\"Binning Level\")\n",
"\tplt.ylabel(\"Error\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_binning_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\tfor i in range(len(statistics_magnetisationq)):\n",
"\t\tplt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1])\n",
"\tleg = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\tleg.append(\"T=\" + str(round(temperatures[i],3)))\n",
"\tplt.legend(leg)\n",
"\tplt.title(\"Binning analysis Wolff, magnetisation$^2$\")\n",
"\tplt.xlabel(\"Binning Level\")\n",
"\tplt.ylabel(\"Error\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_binning_mq.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#calculate errorbars\n",
"\t#binning level 5 seems reasonable for most temperatures\n",
"\terrorbars_m = []\n",
"\terrorbars_mq = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\terrorbars_m.append(statistics_magnetisation[i][1][5])\n",
"\t\terrorbars_mq.append(statistics_magnetisationq[i][1][5])\n",
"\n",
"\t#plot magnetisation with errorbars\n",
"\tplt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1)\n",
"\tplt.title(\"Magnetisation, Wolff\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Magnetisation\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#plot magnetisation squared with errorbars\n",
"\tplt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1)\n",
"\tplt.title(\"Magnetisation$^2$, Wolff\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Magnetisation^2\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"wolff_mq.pdf\")\n",
"\tplt.clf()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"#run a simulation using single spin updates\n",
"def run_metropolis_simulation(dT,Tmin=1.8,Tmax=3.2,meas=131072*4,equil=8192*2):\n",
"\n",
"\t#######################\n",
"\t###\tInitialisation\n",
"\t#######################\n",
"\n",
"\t#arrays to store info about correlations\n",
"\tstatistics_magnetisation = []\n",
"\tstatistics_magnetisationq = []\n",
"\n",
"\t#arrays to store the expectation values and their standard deviations for all temperatures\n",
"\tmagnetisations = []\n",
"\tmagnetisationsq = []\n",
"\tmagnetisations_std = []\n",
"\tmagnetisationsq_std = []\n",
"\ttemperatures = []\n",
"\n",
"\t#define system parameters\n",
"\tL = 16\n",
"\tH = 16\n",
"\tJ = 1.\n",
"\t#initialise the system\n",
"\tsys = []\n",
"\tfor i in range(H):\n",
"\t\ta = []\n",
"\t\tfor j in range(L):\n",
"\t\t\ta.append(1)\n",
"\t\tsys.append(a)\n",
"\t#initialise the magnetisation\n",
"\ttotal_magnetisation = L*H\n",
"\n",
"\tT = Tmin\n",
"\n",
"\t\n",
"\t#######################\n",
"\t###\tSimulation\n",
"\t#######################\n",
"\t#simulate for the specified temperatures\n",
"\twhile(T <= Tmax):\n",
"\t\t#arrays to store all samples\n",
"\t\tmagdata = [] #magnetisation\n",
"\t\tmagqdata = [] #magnetisation squared\n",
"\t\tT += dT\n",
"\t\tbeta = 1./T\n",
"\t\texptable = [1.,-1.,np.exp(-4.*J*beta),-1.,np.exp(-8.*J*beta),1.,-1.,1.,-1.] #lookuptable for the exponents we may need\n",
"\n",
"\t\t#equilibrate\n",
"\t\tfor i in range(int(equil)):\n",
"\t\t\ttotal_magnetisation += MCMC_step(sys,H,L,exptable)\n",
"\t\t#measure\n",
"\t\tfor i in range(int(meas)):\n",
"\t\t\ttotal_magnetisation += MCMC_step(sys,H,L,exptable)\n",
"\t\t\tm = np.abs(total_magnetisation)/float(L*H)\n",
"\t\t\tmagdata.append(m)\n",
"\t\t\tmagqdata.append((m)**2)\n",
"\n",
"\t\t#extract means and their standard deviations\n",
"\t\tmagnetisations.append(np.mean(magdata))\n",
"\t\tmagnetisationsq.append(np.mean(magqdata))\n",
"\t\tmagnetisations_std.append(np.std(magdata))\n",
"\t\tmagnetisationsq_std.append(np.std(magqdata))\n",
"\n",
"\t\t#save temperature\n",
"\t\ttemperatures.append(T)\n",
"\n",
"\t\t#do a binning analysis of the data and store it\n",
"\t\tstatistics_magnetisation.append(binning(magdata))\n",
"\t\tstatistics_magnetisationq.append(binning(magqdata))\n",
"\n",
"\t#calculate correlation times for all temperatures from the binning results. binning level 14 seems okay, but still not good at certain temperatures (CSD)\n",
"\ttaus = []\n",
"\ttausq = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\ttaus.append(0.5*(statistics_magnetisation[i][1][14]**2/magnetisations_std[i]**2*(meas-1.)-1.))\n",
"\t\ttausq.append(0.5*(statistics_magnetisationq[i][1][14]**2/magnetisationsq_std[i]**2*(meas-1.)-1.))\n",
"\n",
"\n",
"\n",
"\t#######################\n",
"\t###\tPlotting\n",
"\t#######################\n",
"\n",
"\t#correlation time magnetisation\n",
"\tplt.plot(temperatures,taus)\n",
"\tplt.title(\"Correlation time Metropolis, magnetization\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Correlation time\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_corr_times_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#correlation time magnetisation squared\n",
"\tplt.plot(temperatures,taus)\n",
"\tplt.title(\"Correlation time Metropolis, magnetization$^2$\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Correlation time\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_corr_times_mq.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#binning analysis results\n",
"\tfor i in range(len(statistics_magnetisation)):\n",
"\t\tplt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1])\n",
"\tleg = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\tleg.append(\"T=\" + str(round(temperatures[i],3)))\n",
"\tplt.legend(leg)\n",
"\tplt.title(\"Binning analysis Metropolis, magnetisation\")\n",
"\tplt.xlabel(\"Binning Level\")\n",
"\tplt.ylabel(\"Error\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_binning_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\tfor i in range(len(statistics_magnetisationq)):\n",
"\t\tplt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1])\n",
"\tleg = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\tleg.append(\"T=\" + str(round(temperatures[i],3)))\n",
"\tplt.legend(leg)\n",
"\tplt.title(\"Binning analysis Metropolis, magnetisation$^2$\")\n",
"\tplt.xlabel(\"Binning Level\")\n",
"\tplt.ylabel(\"Error\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_binning_mq.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#calculate errorbars\n",
"\t#binning level 14 seems reasonable for most temperatures\n",
"\terrorbars_m = []\n",
"\terrorbars_mq = []\n",
"\tfor i in range(len(temperatures)):\n",
"\t\terrorbars_m.append(statistics_magnetisation[i][1][14])\n",
"\t\terrorbars_mq.append(statistics_magnetisationq[i][1][14])\n",
"\n",
"\t#plot magnetisation with errorbars\n",
"\tplt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1)\n",
"\tplt.title(\"Magnetisation, Metropolis\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Magnetisation\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_m.pdf\")\n",
"\tplt.clf()\n",
"\n",
"\t#plot magnetisation squared with errorbars\n",
"\tplt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1)\n",
"\tplt.title(\"Magnetisation$^2$, Metropolis\")\n",
"\tplt.xlabel(\"Temperature\")\n",
"\tplt.ylabel(\"Magnetisation^2\")\n",
"\t#plt.show()\n",
"\tplt.savefig(\"metropolis_mq.pdf\")\n",
"\tplt.clf()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Figure size 432x288 with 0 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#run\n",
"run_metropolis_simulation(0.1)\n",
"run_cluster_simulation(0.1)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}
%% Cell type:code id: tags:
``` python
import numpy as np
from timeit import default_timer as timer
import random
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
#does the next level of binning
def binning_step(data):
newdata = []
for i in range(int(len(data)/2)):
newdata.append((data[2*i]+data[2*i+1])/2.)
stdev = np.std(newdata)/np.sqrt(len(newdata)-1)
return [stdev,newdata]
```
%% Cell type:code id: tags:
``` python
#perform a binning analysis of the data, up to the level where at least minsamples samples remain
def binning(data, minsamples=8):
levels = []
stdevs = []
stdevs.append(np.std(data)/np.sqrt(len(data)-1))
levels.append(0)
level = 1
while(len(data) >= 2*minsamples):
stdev,data = binning_step(data)
stdevs.append(stdev)
levels.append(level)
level += 1
return [levels, stdevs]
```
%% Cell type:code id: tags:
``` python
#perform a Wolff cluster update step on sys of size H*L
#absolute value of return is the cluster size, twice the return is the change in magnetisation
def Wolff_step(sys,H,L,probability):
flipped = 1
#choose random spin
spin_v = random.randint(0,H-1)
spin_h = random.randint(0,L-1)
#flip the spin
sys[spin_v][spin_h] *= -1
#keep track of what we do to the magnetisation
flip_sign = sys[spin_v][spin_h]
#check neighbours recursively
spins_to_check = [[spin_v,spin_h]]
while(len(spins_to_check) > 0):
indices = spins_to_check.pop()
spin_v = indices[0]
spin_h = indices[1]
current_spin = sys[spin_v][spin_h]
#add neighbours to cluster and schedule them with correct probability
if(current_spin * sys[(spin_v+1)%H][spin_h] == -1 and random.random() < probability):
spins_to_check.append([(spin_v+1)%H,spin_h])
sys[(spin_v+1)%H][spin_h] *= -1
flipped += 1
if(current_spin * sys[(spin_v-1)%H][spin_h] == -1 and random.random() < probability):
spins_to_check.append([(spin_v-1)%H,spin_h])
sys[(spin_v-1)%H][spin_h] *= -1
flipped += 1
if(current_spin * sys[spin_v][(spin_h+1)%L] == -1 and random.random() < probability):
spins_to_check.append([spin_v,(spin_h+1)%L])
sys[spin_v][(spin_h+1)%L] *= -1
flipped += 1
if(current_spin * sys[spin_v][(spin_h-1)%L] == -1 and random.random() < probability):
spins_to_check.append([spin_v,(spin_h-1)%L])
sys[spin_v][(spin_h-1)%L] *= -1
flipped += 1
return flipped*flip_sign
```
%% Cell type:code id: tags:
``` python
#perform one single spin update on sys of extent H*L.
#returns change in magnetisation
def MCMC_step(sys,H,L,exptable):
#pick random spin
spin_v = random.randint(0,H-1)
spin_h = random.randint(0,L-1)
#calculate local field
h = sys[spin_v][(spin_h-1)%L] + sys[spin_v][(spin_h+1)%L] + sys[(spin_v+1)%H][spin_h] + sys[(spin_v-1)%H][spin_h]
#flip with metropolis probability
if(random.random() < exptable[sys[spin_v][spin_h]*h]):
if(sys[spin_v][spin_h] == 1):
sys[spin_v][spin_h] = -1
return -2
else:
sys[spin_v][spin_h] = 1
return 2
else:
return 0
```
%% Cell type:code id: tags:
``` python
#returns the magnetisation of sys
def magnetisation(sys):
magn = 0
L = len(sys[0])
H = len(sys)
for a in sys:
for s in a:
magn += s
return float(magn)/float(H*L)
```