Magnetization vs Time (FIG 8.7)

The exact value of the magnetization of a material at a given point of time is uncertain due to the random nature of the spins. This can be demonstrated by plotting the Magnetization as a function of time for a varying range of temperatures. In this Ising Model simulation a 10 x 10 lattice configuration is used with 1000 Monte Carlo sweeps for each temperature point. The time is not real time but solely a product of the Monte Carlo method, this means that one unit of time is one sweep through the entire lattice or one Monte Carlo sweep. This gives each spin the opportunity to flip during each time step. [1]

The beginning temperature is set to 1.5 and the result of which can be seen in Figure 1. The magnetization of the lattice is in the ferromagnetic state which is to be expected at this low temperature. There are a couple of instances in time where the magnetization drops  but these are fairly rare occurrences and just a small fluctuation due to the random nature in the Monte Carlo calculation. For the most part however all the spins are aligned and the lattice is in a ferromagnetic state.

Figure 1. Magnetization vs Time for T = 1.5

When the temperature is slightly increased to 2 as can be seen in figure 2. The fluctuations have increased significantly but the material is still in a ferromagnetic state. The average value of the magnetism has dropped to around to a magnetization value of roughly 0.9. The importance of the increase in fluctuations as they signal that the model is approaching a second order phase transition known as the critical point [1]. This point is sensitive to minor changes in the magnetic field or temperature. This will be shown even clearer when the temperature is raised closer to the Curie temperature.

Figure 2. Magnetization vs Time for T = 2

When this temperature is raised further to 2.25 the fluctuations are much larger and the system fluctuates around ± 0.8. The Curie temperature is 2.27 so the value the simulation is at is very close to that resulting in changes in phase from ferromagnetic to paramagnetic at certain instances of time. This is shown below in figure 3.

Figure 3. Magnetization vs Time for T = 2.25

In figure 4 the temperature is raised to 4 and the fluctuations have decreased and are situated approximately at M = ±0.2. The temperature is now above 2.27 and the system is now in the paragmagnetic state.

Figure 4. Magnetization vs Time for T = 4

The results achieved validated up by the the results found in Giordanos book In figure 5 & 6. The fluctuations due to the random nature of the Ising Model show it is difficult to predict the exact magnetization at a given temperature. However it is possible to get a good idea of what level of magnetization at a given temperature.

Figure 5. Magnetization vs Time from Computational physics by Giordano for T = 1.5 & T = 2

Figure 6. Magnetization vs Time from Computational physics by Giordano for T = 2.25 & T = 4

Beowulf Cluster - The beginning of something big

Beowulf Cluster - The beginning of something big 

A Beowulf cluster is essentially a system consisting of at least 2 PCs working together to perform a given task. Each PC is known as a node with one PC set up to be the server node and the rest as client nodes. The server node controls the whole cluster and serves files to the client nodes. The client nodes are effectively “dumb nodes” and are given a IP address from the server node and can be operated remotely via the server node. The server node is the gateway node to the outside world and is the node that the user would be interfacing with [7]. The operating system that will be used for the cluster is Ubuntu which is a Linux distribution. There are two main types of software packages used to achieve the Beowulf cluster, PVM (Parallel Virtual Machine) and MPI (Message Passing Interface) libraries. The main differences between both libraries are:

• Process Control: PVM can start and stop tasks, to find out which tasks are running, and possibly where they are running. MPI is much more static, It contains functions to start a group of tasks and to send a kill signal to a group of tasks
• Resource Control: PVM is dynamic and can utilize nodes when required. This allows load balancing, task migration and fault tolerance. MPI is static and each node is utilized simultaneously. All potential communication paths are known at start-up. The advantage of this is higher performance but lacks in load balancing, task migration and fault tolerance.
• Message Passing operations: PVM provides very simple message passing capabilities. MPI has a much wider range of message passing capabilities MPI-2 has 248 functions for message-passing operations
• Fault Tolerance: In PVM If one node fails, the entire system adapts, and the application continues without the failed node. In MPI If one node fails, the entire application fails due to MPIs static nature. [8]

For the application in this project, MPI is the better of the two options due to its higher processing power due to its superior message passing capabilities. The application will not take advantage of PVMs many benefits. Due to the repeated nature of the application there is no need for task migration. Fault tolerance would be beneficial, but this is targeted more at high end parallel processing systems with hundreds of nodes running week/month long simulations. The probability of one of the 8 nodes failing is low but in the case that it does the simulation times for the application won’t be nearly long as the high-end systems causing little disruption.

All the components are now in including the switch, 8 PCs and all the peripherals. And I  have just started wiping them and installing Ubuntu.

Figure 1. Pile of Junk

Magnetization and energy versus temperature (FIG 8.8 & 8.9)

The main ingredient of the Ising model is the energy of a pair of spins. The spin can only contain 2 possible values, up or down (+1 or -1). The Ising model tells us on each bond between spins the energy of the bonds is minimized when the 2 spins are aligned. It gives ability to display the phase transition of the material from ferromagnetic to paramagnetic. [1]

Figure 1 – Spin arrangement of ferromagnetic material (left) paramagnetic material (right)

Theory

For example, if we take Iron, at low temperatures (low depending on the characteristics of the material) the spins will be pointing in the same direction and the sample is said to be ferromagnetic. As the temperature is increased some spins begin to flip until the critical temperature is reached which is known as the curie temperature (which differs for different materials). At this point the spin orientation is random and the sample is now paramagnetic. A spin si can only be one of two values ±1. Each of these spins interacts with its nearest neighbour. A simple Ising model assumes an interaction only between nearest neighbours so that the energy of the system can be described as: [1, 4, 5]

Equation 1. energy of the system

<ij> is the sum over all pairs of nearest neighbours spins and j is known as the exchange constant and is assumed to be positive, The probability of finding the system in any particular state is:

Equation 2. he probability of finding the system in any particular state

Where  kB = Boltzmann constant, Eα is the energy of state α, and T is the temperature. The measure of Magnetization of the system is M_α which is the sum of values for sj multiplied by Pα

Equation 3. Magnetization of the system

Using these equations, it is possible to simulate a basic 2D Ising model. This model works by cycling through a number of defined temperature points within a set range. At each temperature point the Monte Carlo method is applied. This means that at every temperature point the simulation is run a set amount of times. The data from the simulation is then added together and divided by the number of simulations ran at that temperature point. This is done as the phase spins are random in nature as it is a quantum phenomenon. The simulation takes an initialised 10x10 lattice (2D Array) that is initially ferromagnetic (+1 in each cell) and cycles through each cell. At each cell Equation 1 is applied. If the value of the equation 1 is less than 0 then the spin is flipped (+1 to -1 or vice versa). If the value is greater than 0 than Equation 2 is applied. If the value of Equation 2 is greater than a randomly generated +1 or -1 then the spin is flipped. [1,4]

Equation 4. Energy per spin using Monte Carlo 

The spontaneous magnetization is calculated by applying Equation 3 for each MC sweep. And again, dividing by the number of MC sweeps × number of cells in array.

Results

Figure 2 below shows the change of spins with an increase in temperature for a 100 x 100 lattice. A snapshot is taken at a given at each temperature from 1 to 5 in steps of 0.1. The temperature at each stage is below on the x - axis individual squares represent a spin with Red = +1, Blue = -1

Figure 2. 100 x 100 lattice 

It can be seen at T = 1 all spins are +1 meaning that all spins are aligned and the material is ferromagnetic. As the temperature rises between 1-2 a small amount of random spins begin to flip around the lattice but revert back to +1. once T > 2 the amount of spins that begin to flip begin increasing until it reaches the curie temperature of roughly T = 2.3 where there is a dramatic increase. From T > 2.3 the spins are randomly flipped in groups.  The number of +1 and -1 spins throughout the material are now roughly equal and the material is now paramagnetic.

The simulations for magnetism and energy can be seen below. The simulations where run twice to show the random variations. The first simulation is shown in blue and the second simulation shown in red.

The Magnetization of the system can be seen in Figure 3. At T = 1 the magnetization is 1 and slowly drops off until it reaches the curie temperature. There is a clear drop off in magnetization at the curie temperature region which shows the phase change. After this point the magnetization hovers around zero. This can be seen in Figure 3 and is validated in Figure 4.

Figure 3. Simulated spontaneous magnetization as a function of temperature

Figure 4. Spontaneous magnetization as a function of temperature from Giordano (Fig 8.8 ) [1]

The energy of the system can be seen to rise as the phase changes, with the transition at the curie temperature. The Simulated thermal average of the energy per spin versus temperature can be seen in figure 5. and validated in figure 6.

Figure 5. Simulated thermal average of the energy per spin versus temperature

Figure 6. thermal average of the energy per spin versus temperature from Giordano (Fig 8.9) [1]

Specific heat per spin (FIG 8.10)

The specific heat per spin is calculated according to the fluctuation-dissipation theorem:

Equation 1. fluctuation-dissipation theorem

Where T is the temperature at a given point, Kb is the Boltzmann’s constant and (ΔE)^2 is the total energy variance and can be calculated as follows: 

Equation 2. total energy variance

Where < E >^2 is the average energy level and < E^2 > is the average squared energy level:

Figure 1. Specific heat simulated

Figure 2. Specific heat computational physics by Giordano (Fig 8.10) [1]

correlation functions (8.11 )

The correlation function shows the relationship between a specific a spin S0 and its neighbors Si the function can be described as :

 The program cycles through each spin in a 20 x 20 lattice. for each spin at a given point S0, the average values of its neighbors Si are multiplied together. To start Si is = 1 so the neighbors are the immediate spins nearest to s0. This goes through all the way to i = 10/ The results of which can be seen below for varying temperatures. 

Figure 1 - Correlation Function simulated

Figure 2 - Correlation Function from Computational Physics by Giordano (Fig 8.11) [1]

Interim Report

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

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 

Supervisor Meeting outcome 22/01/2018

Objectives from meeting:

• Adjust y axis resolution of 2D model so that the energy and magnetism values match those from Giordanos book. 
• Write up what has been achieved thus far for interim report 
• Purchase Ethernet switch online for cluster
• Time permitting reproduce Landaus 1D ising model   
•  Application for ICHEC - Proposal and Abstract now required can be taken from interim report.
• Look into summer of HPC (closing date).

Supervisor Meeting outcome 18/01/2018

After a brief hiatus due to exams a meeting was held with my project supervisor to discuss the next steps of the project to be completed for Monday 22/01/2017. These outcomes include:

• Simulation of a 1D Ising model from the book computational physics (Landau).
• Once this is working, integrate the working code with the existing 2D model I had produced before Christmas.
• Write up a grant application to the Irish Centre for High-End computing for the use of their supercomputer.

On a side note the PCs for the Beowulf cluster are now in my possession and i will begin to slowly piece the cluster together over time. The main focus of the project as it stands will be the operation of the model.

References

[1] N.J. Giordano. “Statistical Mechanics, Phase Transitions and the Ising Model” in Computational Physics, 1^st ed., Alison Reeves, Ed. New Jersey: Pretence Hall 1997, pp.204-231.

[2] IEEE-CS/ACM Joint Task Force - Software Engineering Ethics and Professional Practices, [online] Available: https://www.computer.org/web/education/code-of-ethics

[3] The Institution of Engineers of Ireland. - Code of Ethics, [online] Available: http://www.engineersireland.ie/about/code-of-ethics-and-bye-laws.aspx

[4] J. Block, P. Virnau and T. Preis, “Multi-GPU accelerated multi-spin Monte Carlo simulations of the 2D Ising model”, Computer Physics Communications, Volume 181, Issue 9, pp. 1549-1556, Sep 2010.

[5] A. Züleyha, M. Ziya, Y. Selçuk, Ö.M. Kemal and T. Mesut, “Simulation of glioblastoma multiforme (GBM) tumor cells using ising model on the Creutz Cellular Automaton” Physica A: Statistical Mechanics and its Applications, Volume 486, pp 901-907, Mar 2017.

[6] S.D. Mostovoy, O.V. Pavlovsky, “Critical Casimir effects in 2D Ising model with curved defect lines”, Physics Letters A, Volume 382, Issue 5, pp 276-282, 2018.

 [7] D. Ridge, D. Becker, P. Merkey and T. Sterling, "Beowulf: harnessing the power of parallelism in a pile-of-PCs," 1997 IEEE Aerospace Conference, Snowmass at Aspen, CO, 1997, pp. 79-91 vol.2.

[8] S. Sampath, B. B. Sagar and B. R. Nanjesh, "Performance evaluation and comparison of MPI and PVM using a cluster based parallel computing architecture," 2013 International Conference on Circuits, Power and Computing Technologies (ICCPCT), Nagercoil, 2013, pp. 1253-1258.

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.