Friday, 2 February 2018
Monday, 29 January 2018
2D Ising model v1.1.2
This updated version now has accurate figures which match those of the results obtained in Giordano's book. The simulation has been run twice and the first simulation is represented by the color red and the second simulation is represented by the color blue.
"""
2D Ising model v1.1.2
Modified
Version 1.0: Thu Nov 9 17:34:37 2017
Version 1.1: Sun Jan 21 18:54:34 2018
Version 1.1.1 Mon Jan 22 20:34:28 2018
Version 1.1.2 Mon Jan 29 11:08:00 2018
@author: jameslowe1995
"""
#libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import rand
import numpy
#Adjustable Varibles
H = 2 # Initializng lattice: 1 = all +1s, 2 = all -1s, 3 = random +1s and -1s
Tfinal = 5 #Final temperature
Tinterval = .1 # Intervals between temperature points
T = 1 #Starting temperature
L = 10 #Lattice array L*L
MonteCarloSweeps = 1000 #Monte Carlo Sweeps
#Fixed variables
i = 1
j = 1
r = 0
m = 0
Kb = 1.38064852e-23 #Boltzmanns constant
if H == 1:
#Initialize lattice all +1s
lattice = numpy.zeros((L, L))+1
if H == 2:
#Initialize lattice all -1s
lattice = numpy.zeros((L, L))-1
if H ==3:
#Initialize lattice random +1s and -1s
lattice = 2*np.random.randint(2, size=(L,L))-1
#plot intial lattice
plt.figure(4,figsize=(10, 10), dpi=80)
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel("Initial configuration")
while T <= Tfinal:
sweeps = 0
x = 0
y = 0
z = 0
energy = 0
magnetism =0
E = 0
M = 0
S = 0
latticetotal = numpy.zeros((L, L))
while sweeps < MonteCarloSweeps:
while z == 0:
if x < (L):
if y < (L):
s = lattice[x, y]
#boundary condition 1
if y == 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(L-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 2
if y == 0:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 3
if y == 0:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 4
if y > 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(L-1)%L,y]
#boundary condition 5
if y > 0:
if y != L-1:
if x == L-1:
Eflip = lattice[(x-1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(0)%L,y]
#boundary condition 6
if y == L-1:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(L-1)%L,y]+lattice[x,(y-1)%L]
#boundary condition 7
if y == L-1:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 8
if y == L-1:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(L-1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 9
if y < L-1:
if y > 0:
if x < L-1:
if x > 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
Eflip2 = 2*s*Eflip
if Eflip2 < 0:
s = s*-1 #flip the value
elif rand() < np.exp(-Eflip2*(1/T)):
s = s*-1 #flip the value
lattice[x,y] = s
if y <= L-1:
x = x + 1
if x == L:
x = 0
y = y + 1
if y == L:
x = 0
y = 0
sweeps = sweeps + 1
value = lattice[x,y]
energy += Eflip*(-value)
magnetism = np.sum(lattice)
E = E + energy
E = E/3
M = M + magnetism
# S = S + (magnetism*magnetism)
if sweeps == MonteCarloSweeps:
z = 1
E = (E/(MonteCarloSweeps*L*L))
M = (M/(L*L*MonteCarloSweeps*100))
# S = (((S/(MonteCarloSweeps*L*L)))-((E*E)/(MonteCarloSweeps*L*L)))*(1/T)*(1/T)
print("Energy = ", E, " - magnetism = ",abs(M), " at temperature = ", T)
plot the lattice
plt.figure(figsize=(10, 10), dpi=80);
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel(T)
plt.figure(2,figsize=(10, 10))
plt.ylim([-2,0])
plt.xlim([0,5])
plt.plot(T, E, 'ro', )
plt.xlabel("Temperature",)
plt.ylabel("Energy per Spin ")
plt.title("Thermal average of the energy versus temperature")
plt.figure(3,figsize=(10, 10))
plt.ylim([-0.1,1.1])
plt.xlim([0,5])
#plt.plot([0,0], [0,6], 'k-', lw=2)
plt.plot(T, abs(M), 'bo')
plt.xlabel("Temperature")
plt.ylabel("Magnetization ")
plt.title("Spontaneous magnetization as a function of temperature")
plt.figure(8,figsize=(10, 10))
plt.plot(T, S, 'ro', )
plt.xlabel("",)
plt.ylabel(" ")
plt.title("")
T = T + Tinterval
"""
2D Ising model v1.1.2
Modified
Version 1.0: Thu Nov 9 17:34:37 2017
Version 1.1: Sun Jan 21 18:54:34 2018
Version 1.1.1 Mon Jan 22 20:34:28 2018
Version 1.1.2 Mon Jan 29 11:08:00 2018
@author: jameslowe1995
"""
#libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import rand
import numpy
#Adjustable Varibles
H = 2 # Initializng lattice: 1 = all +1s, 2 = all -1s, 3 = random +1s and -1s
Tfinal = 5 #Final temperature
Tinterval = .1 # Intervals between temperature points
T = 1 #Starting temperature
L = 10 #Lattice array L*L
MonteCarloSweeps = 1000 #Monte Carlo Sweeps
#Fixed variables
i = 1
j = 1
r = 0
m = 0
Kb = 1.38064852e-23 #Boltzmanns constant
if H == 1:
#Initialize lattice all +1s
lattice = numpy.zeros((L, L))+1
if H == 2:
#Initialize lattice all -1s
lattice = numpy.zeros((L, L))-1
if H ==3:
#Initialize lattice random +1s and -1s
lattice = 2*np.random.randint(2, size=(L,L))-1
#plot intial lattice
plt.figure(4,figsize=(10, 10), dpi=80)
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel("Initial configuration")
while T <= Tfinal:
sweeps = 0
x = 0
y = 0
z = 0
energy = 0
magnetism =0
E = 0
M = 0
S = 0
latticetotal = numpy.zeros((L, L))
while sweeps < MonteCarloSweeps:
while z == 0:
if x < (L):
if y < (L):
s = lattice[x, y]
#boundary condition 1
if y == 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(L-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 2
if y == 0:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 3
if y == 0:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 4
if y > 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(L-1)%L,y]
#boundary condition 5
if y > 0:
if y != L-1:
if x == L-1:
Eflip = lattice[(x-1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(0)%L,y]
#boundary condition 6
if y == L-1:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(L-1)%L,y]+lattice[x,(y-1)%L]
#boundary condition 7
if y == L-1:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 8
if y == L-1:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(L-1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 9
if y < L-1:
if y > 0:
if x < L-1:
if x > 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
Eflip2 = 2*s*Eflip
if Eflip2 < 0:
s = s*-1 #flip the value
elif rand() < np.exp(-Eflip2*(1/T)):
s = s*-1 #flip the value
lattice[x,y] = s
if y <= L-1:
x = x + 1
if x == L:
x = 0
y = y + 1
if y == L:
x = 0
y = 0
sweeps = sweeps + 1
value = lattice[x,y]
energy += Eflip*(-value)
magnetism = np.sum(lattice)
E = E + energy
E = E/3
M = M + magnetism
# S = S + (magnetism*magnetism)
if sweeps == MonteCarloSweeps:
z = 1
E = (E/(MonteCarloSweeps*L*L))
M = (M/(L*L*MonteCarloSweeps*100))
# S = (((S/(MonteCarloSweeps*L*L)))-((E*E)/(MonteCarloSweeps*L*L)))*(1/T)*(1/T)
print("Energy = ", E, " - magnetism = ",abs(M), " at temperature = ", T)
plot the lattice
plt.figure(figsize=(10, 10), dpi=80);
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel(T)
plt.figure(2,figsize=(10, 10))
plt.ylim([-2,0])
plt.xlim([0,5])
plt.plot(T, E, 'ro', )
plt.xlabel("Temperature",)
plt.ylabel("Energy per Spin ")
plt.title("Thermal average of the energy versus temperature")
plt.figure(3,figsize=(10, 10))
plt.ylim([-0.1,1.1])
plt.xlim([0,5])
#plt.plot([0,0], [0,6], 'k-', lw=2)
plt.plot(T, abs(M), 'bo')
plt.xlabel("Temperature")
plt.ylabel("Magnetization ")
plt.title("Spontaneous magnetization as a function of temperature")
plt.figure(8,figsize=(10, 10))
plt.plot(T, S, 'ro', )
plt.xlabel("",)
plt.ylabel(" ")
plt.title("")
T = T + Tinterval
Monday, 22 January 2018
2D Ising model v1.1.1
This updated version now has an accurate Y-axis scale on magnetism and Energy graphs. A 10x10 lattice is simulated as to compare with the graphs from Giordano s book. The temperature interval is set to .1 as to achieve a quick simulation.
"""
2D Ising model v1.1.1
Modified
Version 1.0: Thu Nov 9 17:34:37 2017
Version 1.1: Sun Jan 21 18:54:34 2018
Version 1.1.1 Mon Jan 22 20:34:28 2018
@author: jameslowe1995
"""
#libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from random import randint
from numpy.random import rand
import numpy
#Adjustable Varibles
H = 2 # Initializng lattice: 1 = all +1s, 2 = all -1s, 3 = random +1s and -1s
Tfinal = 4 #Final temperature
Tinterval = .1 # Intervals between temperature points
T = 1 #Starting temperature
L = 10 #Lattice array L*L
MonteCarloSweeps = 1000 #Monte Carlo Sweeps
#Fixed variables
i = 1
j = 1
r = 0
m = 0
Kb = 1.38064852e-23 #Boltzmanns constant
if H == 1:
#Initialize lattice all +1s
lattice = numpy.zeros((L, L))+1
if H == 2:
#Initialize lattice all -1s
lattice = numpy.zeros((L, L))-1
if H ==3:
#Initialize lattice random +1s and -1s
lattice = 2*np.random.randint(2, size=(L,L))-1
#plot intial lattice
plt.figure(4,figsize=(10, 10), dpi=80)
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel("Initial configuration")
while T <= Tfinal:
sweeps = 0
x = 0
y = 0
z = 0
energy = 0
magnetism =0
E = 0
M = 0
latticetotal = numpy.zeros((L, L))
while sweeps < MonteCarloSweeps:
while z == 0:
if x < (L):
if y < (L):
s = lattice[x, y]
#boundary condition 1
if y == 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(L-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 2
if y == 0:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 3
if y == 0:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 4
if y > 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(L-1)%L,y]
#boundary condition 5
if y > 0:
if y != L-1:
if x == L-1:
Eflip = lattice[(x-1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(0)%L,y]
#boundary condition 6
if y == L-1:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(L-1)%L,y]+lattice[x,(y-1)%L]
#boundary condition 7
if y == L-1:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 8
if y == L-1:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(L-1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 9
if y < L-1:
if y > 0:
if x < L-1:
if x > 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
Eflip2 = 2*s*Eflip
if Eflip2 < 0:
s = s*-1 #flip the value
elif rand() < np.exp(-Eflip2*(1/T)):
s = s*-1 #flip the value
lattice[x,y] = s
if y <= L-1:
x = x + 1
if x == L:
x = 0
y = y + 1
if y == L:
x = 0
y = 0
sweeps = sweeps + 1
value = lattice[x,y]
energy += Eflip*(-value)
magnetism = np.sum(lattice)
E = E + energy
E = E/4
M = M + magnetism
if sweeps == MonteCarloSweeps:
z = 1
E = (E/(MonteCarloSweeps*L*L))
M = (M/(L*L*MonteCarloSweeps*100))
print("Energy = ", E, " - magnetism = ",abs(M), " at temperature = ", T)
#plot the lattice
plt.figure(figsize=(10, 10), dpi=80);
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel(T)
plt.figure(2,figsize=(15, 10))
plt.plot(T, E, 'ro', )
plt.xlabel("Temperature (T)",)
plt.ylabel("Energy ")
plt.figure(3,figsize=(15, 10))
plt.plot(T, abs(M), 'bo')
plt.xlabel("Temperature (T)")
plt.ylabel("Magnetization ")
T = T + Tinterval
"""
2D Ising model v1.1.1
Modified
Version 1.0: Thu Nov 9 17:34:37 2017
Version 1.1: Sun Jan 21 18:54:34 2018
Version 1.1.1 Mon Jan 22 20:34:28 2018
@author: jameslowe1995
"""
#libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from random import randint
from numpy.random import rand
import numpy
#Adjustable Varibles
H = 2 # Initializng lattice: 1 = all +1s, 2 = all -1s, 3 = random +1s and -1s
Tfinal = 4 #Final temperature
Tinterval = .1 # Intervals between temperature points
T = 1 #Starting temperature
L = 10 #Lattice array L*L
MonteCarloSweeps = 1000 #Monte Carlo Sweeps
#Fixed variables
i = 1
j = 1
r = 0
m = 0
Kb = 1.38064852e-23 #Boltzmanns constant
if H == 1:
#Initialize lattice all +1s
lattice = numpy.zeros((L, L))+1
if H == 2:
#Initialize lattice all -1s
lattice = numpy.zeros((L, L))-1
if H ==3:
#Initialize lattice random +1s and -1s
lattice = 2*np.random.randint(2, size=(L,L))-1
#plot intial lattice
plt.figure(4,figsize=(10, 10), dpi=80)
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel("Initial configuration")
while T <= Tfinal:
sweeps = 0
x = 0
y = 0
z = 0
energy = 0
magnetism =0
E = 0
M = 0
latticetotal = numpy.zeros((L, L))
while sweeps < MonteCarloSweeps:
while z == 0:
if x < (L):
if y < (L):
s = lattice[x, y]
#boundary condition 1
if y == 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(L-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 2
if y == 0:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 3
if y == 0:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y] + lattice[x,(L-1)%L]
#boundary condition 4
if y > 0:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(L-1)%L,y]
#boundary condition 5
if y > 0:
if y != L-1:
if x == L-1:
Eflip = lattice[(x-1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+lattice[(0)%L,y]
#boundary condition 6
if y == L-1:
if x == 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(L-1)%L,y]+lattice[x,(y-1)%L]
#boundary condition 7
if y == L-1:
if x < L-1:
if x != 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(0)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 8
if y == L-1:
if x == L-1:
Eflip = lattice[(0)%L,y] + lattice[x,(L-1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
#boundary condition 9
if y < L-1:
if y > 0:
if x < L-1:
if x > 0:
Eflip = lattice[(x+1)%L,y] + lattice[x,(y+1)%L] + lattice[(x-1)%L,y]+ lattice[x,(y-1)%L]
Eflip2 = 2*s*Eflip
if Eflip2 < 0:
s = s*-1 #flip the value
elif rand() < np.exp(-Eflip2*(1/T)):
s = s*-1 #flip the value
lattice[x,y] = s
if y <= L-1:
x = x + 1
if x == L:
x = 0
y = y + 1
if y == L:
x = 0
y = 0
sweeps = sweeps + 1
value = lattice[x,y]
energy += Eflip*(-value)
magnetism = np.sum(lattice)
E = E + energy
E = E/4
M = M + magnetism
if sweeps == MonteCarloSweeps:
z = 1
E = (E/(MonteCarloSweeps*L*L))
M = (M/(L*L*MonteCarloSweeps*100))
print("Energy = ", E, " - magnetism = ",abs(M), " at temperature = ", T)
#plot the lattice
plt.figure(figsize=(10, 10), dpi=80);
X, Y = np.meshgrid(range(L), range(L))
plt.pcolormesh(X, Y, lattice, cmap=plt.cm.RdBu);
plt.xlabel(T)
plt.figure(2,figsize=(15, 10))
plt.plot(T, E, 'ro', )
plt.xlabel("Temperature (T)",)
plt.ylabel("Energy ")
plt.figure(3,figsize=(15, 10))
plt.plot(T, abs(M), 'bo')
plt.xlabel("Temperature (T)")
plt.ylabel("Magnetization ")
T = T + Tinterval
On a side note by setting the Monte Carlo sweeps to 1 it is possible to see the relationship between each temperature point. Below is a 100x100 lattice with a temperature interval of 0.1
Subscribe to:
Posts (Atom)