"""
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
