Commit 6e2338b0 authored by Valerio's avatar Valerio
Browse files

solutions ex7

parent e655603f
{
"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": "markdown",
"metadata": {},
"source": [
"Running this notebook locally on your laptop may take long time. Order of 1h."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#perform a Wolff cluster update step on sys of size H*L, with anisotropic couplings\n",
"#absolute value of return is the cluster size, twice the return is the change in magnetisation\n",
"def Wolff_step_aniso(sys,H,L,probability_horizontal, probability_vertical):\n",
" flipped = 1\n",
" #choose random spin\n",
" spin_v = random.randint(0,H-1)\n",
" spin_h = random.randint(0,L-1)\n",
" #flip the spin\n",
" sys[spin_v][spin_h] *= -1\n",
" #keep track of what we do to the magnetisation\n",
" flip_sign = sys[spin_v][spin_h]\n",
" #check neighbours recursively\n",
" spins_to_check = [[spin_v,spin_h]]\n",
" while(len(spins_to_check) > 0):\n",
" indices = spins_to_check.pop()\n",
" spin_v = indices[0]\n",
" spin_h = indices[1]\n",
" current_spin = sys[spin_v][spin_h]\n",
" #add neighbours to cluster and schedule them with correct probability\n",
" if(current_spin * sys[(spin_v+1)%H][spin_h] == -1 and random.random() < probability_vertical):\n",
" spins_to_check.append([(spin_v+1)%H,spin_h])\n",
" sys[(spin_v+1)%H][spin_h] *= -1\n",
" flipped += 1\n",
"\n",
" if(current_spin * sys[(spin_v-1)%H][spin_h] == -1 and random.random() < probability_vertical):\n",
" spins_to_check.append([(spin_v-1)%H,spin_h])\n",
" sys[(spin_v-1)%H][spin_h] *= -1\n",
" flipped += 1\n",
"\n",
" if(current_spin * sys[spin_v][(spin_h+1)%L] == -1 and random.random() < probability_horizontal):\n",
" spins_to_check.append([spin_v,(spin_h+1)%L])\n",
" sys[spin_v][(spin_h+1)%L] *= -1\n",
" flipped += 1\n",
"\n",
" if(current_spin * sys[spin_v][(spin_h-1)%L] == -1 and random.random() < probability_horizontal):\n",
" spins_to_check.append([spin_v,(spin_h-1)%L])\n",
" sys[spin_v][(spin_h-1)%L] *= -1\n",
" flipped += 1\n",
" return flipped*flip_sign"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#returns the magnetisation of sys\n",
"def magnetisation(sys):\n",
" magn = 0\n",
" L = len(sys[0])\n",
" H = len(sys)\n",
" for a in sys:\n",
" for s in a:\n",
" magn += s\n",
" return float(magn)/float(H*L)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class data_cl:\n",
" def __init__(self,L):\n",
" self.magnetisations = []\n",
" self.magnetisationsq = []\n",
" self.binders2 = []\n",
" self.binders4 = []\n",
" self.chis = []\n",
" self.L = L"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#run a simulation using cluster updates with anisotropy\n",
"def run_aniso_cluster_simulation(Jmin=0.5,Jmax=2.5,meas=8192*2,equil=1000*2):\n",
"\n",
" #######################\n",
" ###\tInitialisation\n",
" #######################\n",
"\n",
" Js = np.linspace(Jmin,Jmax,30)\n",
" delta = 0.05\n",
" gamma = 1.\n",
"\n",
" datas = {8: data_cl(8), 12: data_cl(12), 16: data_cl(16), 20: data_cl(20), 24: data_cl(24), 28: data_cl(28)}\n",
"\n",
" \n",
" #define system parameters\n",
" for i in reversed([8,12,16,20,24,28]):\n",
" d = datas[i]\n",
" L = d.L\n",
" H = d.L\n",
" sys = []\n",
" for i in range(H):\n",
" a = []\n",
" for j in range(L):\n",
" a.append(1)\n",
" sys.append(a)\n",
" #initialise the magnetisation\n",
" total_magnetisation = L*H\n",
"\n",
" #######################\n",
" ###\tSimulation\n",
" #######################\n",
" #simulate for the specified temperatures\n",
" for J in Js:\n",
" print(\"{} {}\".format(L,J))\n",
" #arrays to store all samples\n",
" #magdata = [] #magnetisation\n",
" #magqdata = [] #magnetisation squared\n",
" #magqqdata = [] #magnetisation **4\n",
" mag = 0.\n",
" magq = 0.\n",
" magqq = 0.\n",
" prob_h = 1. - np.exp(-2.*delta*J)\n",
" prob_v = 1. - delta*gamma\n",
"\n",
" e_s = equil\n",
" m_s = meas\n",
" if(J < 1.5):\n",
" e_s = int(8*equil)\n",
" m_s = int(8*meas)\n",
" else:\n",
" e_s = int(equil*2.)\n",
" m_s = int(meas*2.)\n",
"\n",
" #equilibrate\n",
" for i in range(e_s):\n",
" total_magnetisation += 2*Wolff_step_aniso(sys,H,L,prob_h,prob_v)\n",
" #measure\n",
" for i in range(m_s):\n",
" wolff_return = Wolff_step_aniso(sys,H,L,prob_h,prob_v)\n",
" m = np.abs(wolff_return)/float(L*H)\n",
" total_magnetisation += 2*wolff_return\n",
" #magdata.append(m)\n",
" #magdata.append(np.abs(total_magnetisation)/float(L*H))\n",
" #magqdata.append((total_magnetisation/float(L*H))**2)\n",
" #magqqdata.append((total_magnetisation/float(L*H))**4)\n",
" mag += m\n",
" magq += (total_magnetisation/float(L*H))**2\n",
" magqq += (total_magnetisation/float(L*H))**4\n",
"\n",
" #extract means and their standard deviations\n",
" #d.magnetisations.append(np.mean(magdata))\n",
" #d.magnetisationsq.append(np.mean(magqdata))\n",
" #d.binders2.append(np.mean(magqdata)/np.mean(magdata)**2)\n",
" #d.binders4.append(np.mean(magqqdata)/np.mean(magqdata)**2)\n",
" #d.chis.append(np.mean(magqdata)-np.mean(magdata)**2)\n",
" magex = float(mag)/float(m_s)\n",
" magqex = float(magq)/float(m_s)\n",
" magqqex = float(magqq)/float(m_s)\n",
" d.magnetisations.append(magex)\n",
" d.magnetisationsq.append(magqqex)\n",
" d.binders2.append(magqex/magex**2)\n",
" d.binders4.append(magqqex/magqex**2)\n",
" d.chis.append(magqex-magex**2)\n",
" \n",
"\n",
"\n",
"\n",
"\n",
" #######################\n",
" ###\tPlotting\n",
" #######################\n",
"\n",
" #plot magnetisation\n",
" leg = []\n",
" for i in [8,12,16,20,24,28]:\n",
" plt.plot(Js,datas[i].magnetisations,'o-',ms=1)\n",
" leg.append(\"L=\" + str(i))\n",
" plt.legend(leg)\n",
" plt.title(\"Magnetisation\")\n",
" plt.xlabel(\"J/$\\Gamma$\")\n",
" plt.ylabel(\"Magnetisation\")\n",
" #plt.show()\n",
" plt.savefig(\"M.pdf\")\n",
" plt.clf()\n",
"\n",
" #plot second order Binder\n",
" leg = []\n",
" for i in [8,12,16,20,24,28]:\n",
" plt.plot(Js,datas[i].binders2,'o-',ms=1)\n",
" leg.append(\"L=\" + str(i))\n",
" plt.legend(leg)\n",
" plt.title(\"Second order Binder cumulant\")\n",
" plt.xlabel(\"J/$\\Gamma$\")\n",
" plt.ylabel(\"$U_2$\")\n",
" plt.savefig(\"U2.pdf\")\n",
" plt.clf()\n",
"\n",
" #plot fourth order Binder\n",
" leg = []\n",
" for i in [8,12,16,20,24,28]:\n",
" plt.plot(Js,datas[i].binders4,'o-',ms=1)\n",
" leg.append(\"L=\" + str(i))\n",
" plt.legend(leg)\n",
" plt.title(\"Fourth order Binder cumulant\")\n",
" plt.xlabel(\"J/$\\Gamma$\")\n",
" plt.ylabel(\"$U_4$\")\n",
" plt.savefig(\"U4.pdf\")\n",
" plt.clf()\n",
"\n",
" #plot susceptibility\n",
" leg = []\n",
" for i in [8,12,16,20,24,28]:\n",
" plt.plot(Js,datas[i].chis,'o-',ms=1)\n",
" leg.append(\"L=\" + str(i))\n",
" plt.legend(leg)\n",
" plt.title(\"Susceptibility\")\n",
" plt.xlabel(\"J/$\\Gamma$\")\n",
" plt.ylabel(\"$\\chi$\")\n",
" plt.savefig(\"Chi.pdf\")\n",
" plt.clf()"
]
}
],
"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:markdown id: tags:
Running this notebook locally on your laptop may take long time. Order of 1h.
%% Cell type:code id: tags:
``` python
#perform a Wolff cluster update step on sys of size H*L, with anisotropic couplings
#absolute value of return is the cluster size, twice the return is the change in magnetisation
def Wolff_step_aniso(sys,H,L,probability_horizontal, probability_vertical):
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_vertical):
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_vertical):
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_horizontal):
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_horizontal):
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
#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)
```
%% Cell type:code id: tags:
``` python
class data_cl:
def __init__(self,L):
self.magnetisations = []
self.magnetisationsq = []
self.binders2 = []
self.binders4 = []
self.chis = []
self.L = L
```
%% Cell type:code id: tags:
``` python
#run a simulation using cluster updates with anisotropy
def run_aniso_cluster_simulation(Jmin=0.5,Jmax=2.5,meas=8192*2,equil=1000*2):
#######################
### Initialisation
#######################
Js = np.linspace(Jmin,Jmax,30)
delta = 0.05
gamma = 1.
datas = {8: data_cl(8), 12: data_cl(12), 16: data_cl(16), 20: data_cl(20), 24: data_cl(24), 28: data_cl(28)}
#define system parameters
for i in reversed([8,12,16,20,24,28]):
d = datas[i]
L = d.L
H = d.L
sys = []
for i in range(H):
a = []
for j in range(L):
a.append(1)
sys.append(a)
#initialise the magnetisation
total_magnetisation = L*H
#######################
### Simulation
#######################
#simulate for the specified temperatures
for J in Js:
print("{} {}".format(L,J))
#arrays to store all samples
#magdata = [] #magnetisation
#magqdata = [] #magnetisation squared
#magqqdata = [] #magnetisation **4
mag = 0.
magq = 0.
magqq = 0.
prob_h = 1. - np.exp(-2.*delta*J)
prob_v = 1. - delta*gamma
e_s = equil
m_s = meas
if(J < 1.5):
e_s = int(8*equil)
m_s = int(8*meas)
else:
e_s = int(equil*2.)
m_s = int(meas*2.)
#equilibrate
for i in range(e_s):
total_magnetisation += 2*Wolff_step_aniso(sys,H,L,prob_h,prob_v)
#measure
for i in range(m_s):
wolff_return = Wolff_step_aniso(sys,H,L,prob_h,prob_v)
m = np.abs(wolff_return)/float(L*H)
total_magnetisation += 2*wolff_return
#magdata.append(m)
#magdata.append(np.abs(total_magnetisation)/float(L*H))
#magqdata.append((total_magnetisation/float(L*H))**2)
#magqqdata.append((total_magnetisation/float(L*H))**4)
mag += m
magq += (total_magnetisation/float(L*H))**2
magqq += (total_magnetisation/float(L*H))**4
#extract means and their standard deviations
#d.magnetisations.append(np.mean(magdata))
#d.magnetisationsq.append(np.mean(magqdata))
#d.binders2.append(np.mean(magqdata)/np.mean(magdata)**2)
#d.binders4.append(np.mean(magqqdata)/np.mean(magqdata)**2)
#d.chis.append(np.mean(magqdata)-np.mean(magdata)**2)
magex = float(mag)/float(m_s)
magqex = float(magq)/float(m_s)
magqqex = float(magqq)/float(m_s)
d.magnetisations.append(magex)
d.magnetisationsq.append(magqqex)
d.binders2.append(magqex/magex**2)
d.binders4.append(magqqex/magqex**2)
d.chis.append(magqex-magex**2)
#######################
### Plotting
#######################
#plot magnetisation
leg = []
for i in [8,12,16,20,24,28]:
plt.plot(Js,datas[i].magnetisations,'o-',ms=1)
leg.append("L=" + str(i))
plt.legend(leg)
plt.title("Magnetisation")
plt.xlabel("J/$\Gamma$")
plt.ylabel("Magnetisation")
#plt.show()
plt.savefig("M.pdf")
plt.clf()
#plot second order Binder
leg = []
for i in [8,12,16,20,24,28]:
plt.plot(Js,datas[i].binders2,'o-',ms=1)
leg.append("L=" + str(i))
plt.legend(leg)
plt.title("Second order Binder cumulant")
plt.xlabel("J/$\Gamma$")
plt.ylabel("$U_2$")
plt.savefig("U2.pdf")
plt.clf()
#plot fourth order Binder
leg = []
for i in [8,12,16,20,24,28]:
plt.plot(Js,datas[i].binders4,'o-',ms=1)
leg.append("L=" + str(i))
plt.legend(leg)
plt.title("Fourth order Binder cumulant")
plt.xlabel("J/$\Gamma$")
plt.ylabel("$U_4$")
plt.savefig("U4.pdf")
plt.clf()
#plot susceptibility
leg = []
for i in [8,12,16,20,24,28]:
plt.plot(Js,datas[i].chis,'o-',ms=1)
leg.append("L=" + str(i))
plt.legend(leg)
plt.title("Susceptibility")
plt.xlabel("J/$\Gamma$")
plt.ylabel("$\chi$")
plt.savefig("Chi.pdf")
plt.clf()
```
This diff is collapsed.
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment