2D Ising model v1

# -*- coding: utf-8 -*-

"""
2D Ising model v1.0

Created on Wed Nov  9 17:34:37 2017

@author: jameslowe1995
"""

#libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import numpy


#Adjustable Varibles
H = 1                                # Initializing 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 = 10 # 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


while T <= Tfinal:

    
    sweeps = 1
    x = 0
    y = 0
    z = 0
        
    while sweeps < MonteCarloSweeps:

      while z == 0:
            if x < (L):
                if y < (L):
                    value = lattice[x,y]
                    
# =======================================================================
#                 
#                   boundary conditions are required to calculate the energy
#                   of each nearest neighbor, the following boundary conditions 
#                   show the neighbors of a 3x3 lattice                 
#                     
#                     boundary condition 1 = 2 + 3 + 4 + 7
#                     boundary condition 2 = 1 + 3 + 5 + 8
#                     boundary condition 3 = 1 + 2 + 6 + 9
#                     boundary condition 4 = 1 + 5 + 6 + 7
#                     boundary condition 5 = 2 + 4 + 6 + 8
#                     boundary condition 6 = 3 + 4 + 5 + 9
#                     boundary condition 7 = 1 + 4 + 8 + 9 
#                     boundary condition 8 = 2 + 5 + 7 + 9
#                     boundary condition 9 = 3 + 6 + 7 + 8
#                     
#                     xxxxx xxxxx xxxxx 
#                     x      x x      x x      x 
#                     x  1  x x  2  x x  3  x 
#                     x      x x      x x      x 
#                     xxxxx xxxxx xxxxx   
#                     
#                     xxxxx xxxxx xxxxx 
#                     x      x x      x x      x 
#                     x  4  x x  5  x x  6  x 
#                     x      x x      x x      x 
#                     xxxxx xxxxx xxxxx 
#                     
#                     xxxxx xxxxx xxxxx 
#                     x      x x      x x      x 
#                     x  7  x x  8  x x  9  x 
#                     x      x x      x x      x 
#                     xxxxx xxxxx xxxxx 
#                     
# ======================================================================
        
                    #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]
                        
                    if Eflip <= 0:
                        value = value*-1 #flip the value
                    if Eflip > 0:
                            r = randint(0,1) # random number 0-1
                            if r <= exp(-Eflip/Kb*(T)):
                                value = value*-1 #flip the value
                    lattice[x,y] = value
                    
                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
                
                if sweeps == MonteCarloSweeps:
                    z = 1
    
    
    #Calculate energy and magnetism    
    energy = 0
    magnetism =0
    x = 1
    y = 1
    value = lattice[x,y]
    if x <= (L):
        if y <= (L):
            value = lattice[x,y]
            NeigbourTotal = lattice[(x+1), y] + lattice[x,(y+1)] + lattice[(x-1), y] + lattice[x,(y-1)]
            energy += -NeigbourTotal*value
    magnetism = np.sum(lattice)
    print("Energy = ", energy, " - magnetism = ",abs(magnetism), " at temperature = ", T)
    
    #plot the lattice               
    plt.figure()
    plt.imshow(lattice)
    plt.xlabel(T)
    
    #plot energy
    plt.figure(2)
    plt.plot(T, energy, 'ro', )
    plt.xlabel("Temperature (T)",)
    plt.ylabel("Energy ")

    #plot magnetization         
    plt.figure(3)
    plt.plot(T, abs(magnetism), 'bo')
    plt.xlabel("Temperature (T)")
    plt.ylabel("Magnetization ")

    
    T = T + Tinterval

#END             



The code above outputs the following lattice array configurations at varying temperatures along with the initial configuration. At this stage in the development of the code the lattice array is not behaving as expected. The resulting arrays tend to change randomly at each temperature value. The graphs are also showing the random values. However, the base code is now in place.