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.   

Project Selection

Supervisor meeting 19/10/17

After a meeting with my supervisor over a cup of coffee, the following outcomes were reached:

• The project will be on the Ising model.
• The models will be run using python.
• If feasible, a Beowulf cluster will be put together using old PCs from DIT to run later simulations (will deal with later on in the project).

Computational Complexity

When looking at modeling big data, certain circumstances must be accounted for mainly concerning a computers processing power. My current device is Microsoft surface pro 4 which has 4 GB's of RAM and is powered by an Intel Core I5-6300U CPU @ 2.40GHz. The latest PC's in the Dublin Institute of Technology are the OptiPlex 5040 by Dell which has 8GB'S and powered by an Intel Core I3-6100 4GB 1600MHz.