Jump to content

Rep:Mod:JPS1124

From ChemWiki

Third Year CMP Compulsory Experiment James Simpson (CID:00733493)

Introduction to the Ising Model

The Model

Table 1: The Relationship Between the Number of Dimensions and the Number of Neighbours
Number of Dimensions Number of Neighbors
1 2
2 4
3 6
D 2D

The Ising model is an physics model used in order to understand the behaviour of ferromagnets. Ferromagnets are materials in which the magnetic dipoles of the material align so that an overall magnetic dipole is exhibited by the material. This effect is due to the favourable energy minimisation due to the alignment of the dipoles or spins. However this will unfavourably decrease the entropy. In the model only spins between neighbouring lattice points interact. This interaction is defined as E=12JiNjneighbours(i)sisj, where J is a constant and si and sj are the spins of the lattice point and its neighbour. The third rule is that a cell at the edge of the lattice will interact with another cell at the other edge of the lattice; this is so that all possible neighbours will be interacted with. The number of neighbours in a particular number of dimensions is expressed in table 1.

The interaction energy of a number of particles in a particular number of dimensions can be expressed as E=DNJ. This can be shown using the initial equation in the lab script of E=12JiNjneighbours(i)sisj. From table 1 it is clear that the number of neighbours each lattice cell has is equal to twice the number of dimensions. The first step is to show what the lowest interaction energy must be; in the lowest energy all the magnetic spins are parallel and so sisj=1 this means that the jneighbours(i)sisj term in the equation will simply be equal to the number of neighbours which is 2-D. This is simply done N times as there as N number of particles and so iNjneighbours(i)sisj=2ND. This is then multiplied by the constant of 12J, giving the interaction energy to be DNJ. From this point it is possible to consider the multiplicity of the system. Ordinarily in chemistry the multiplicity of a system is given by the equation Multiplicity=2s+1, but this cannot be done in this case as the values of the magnetic spin are integers and so the lattice cells can be considered as quasi-boson particles. As the equation stated in the previous sentence is designed for electrons it must be ignored, as in the case of a one-lattice cell unit it would give a multiplicity of 3 where a value of 2 is logically expected. All the magnetic spins in a ferromagnetic material will be aligned so that the spins are all parallel, however the spins of each lattice cell can have a value of +1 or a value of -1. This means that the number of micro-states is 2 and so the multiplicity is 2. It then follows that the entropy of the system given by S=kbln(Ω), where Ω is the number of micro-states, in this case the multiplicity of the system. Therefore the entropy of the system is simply S=kbln(2) which is 9.570x10-24 JK-1.

Phase Transition

In the lowest energy configuration the interaction energy is -3000J. However, when one spin changes the interaction energy of the system will naturally increase. In order to consider how much the energy increases when a spin flips, the effect this will have on the system must be thought about. It is equivalent to removing the spin entirely and putting back into the system another spin, but with a direction opposite to the rest of the system. When a spin is removed the system will lose six interactions, one for each neighbour, and then when the opposite spin is added there will be six interactions created: however this will increase the interaction energy by 6J. So, the flipping of the spin destabilizes the energy by 12J in total as 6J of favourable interactions are lost and 6J of unfavourable interactions are gained. This leads to the interaction energy of this system being -2988J. The new entropy of the system will be S=kbln(2000) as the new spin can occur anywhere in the lattice and there are 1000 options for that, and in addition there will be 2 options for every point in the lattice where the opposite spin will be. This is because the majority of the system can a have a spin of +1, where the opposite spin will be -1. The other option is where most of the system will have a spin of -1 and in this case the opposite spin will have a value of +1. This means when compared to the lowest energy configuration the entropy will have increased by S=kbln(1000), which is 9.54x10-23 JK-1.

Figure 1: The 1-D and 2-D Lattices Given in the Lab Script used in this Exercise

The Curie temperature is the temperature below which ferromagnetism will be exhibited. At temperatures below the Curie temperature the stabilisation energy will be large enough to compensate for the loss in entropy. However above the Curie temperature this is not the case and the entropy effect will dominate, leading to the material showing diamagnetism. Magnetisation is given as M=isi and so in all cases the magnetisation is simply the difference between the number of +1 spins and the number of -1 spins. In the case of the 1-D lattice there are three +1 spins and two -1 spins and so the magnetisation will be simply +1. For the 2-D case there are thirteen +1 spins and twelve -1 spins, meaning that the magnetisation again in that case is +1. Although both of the systems in figure 1 are not very magnetised, the Ising lattice in 3-D containing 1000 lattice cells at absolute zero will be highly magnetised. As the system is at absolute zero there will be no thermal energy available in order to overcome the spin flipping energy barriers, therefore it is expected that the system will adopt the lowest energy configuration, which is where all the spins align. This means that the value of the magnetisation will be equal to the number of lattice cells, in this case 1000. Therefore the magnetisation is either -1000 or +1000, but it cannot be known which of these two options is correct without further investigation.

Calculating the Energy and Magnetisation

Whenever the ipython programme was loaded the following two lines were run in order to start the session %load_ext autoreload and %autoreload 2. Firstly the files IsingLattice.py and ILcheck.py were extracted and stored in the H:Drive.

Modifying the Files

This section involved two parts: one of these was to find the magnetisation while the other was to find the energy of a random arrangement of a specific lattice size, determined by the user. The code used to find the magnetisation is shown below:

   def magnetisation(self):
       magnetisation = 0
       for i in range(self.n_rows):
           for j in range (self.n_cols):
               magnetisation += self.lattice [i,j]
       #Return the total magnetisation of the current lattice configuration.
       return magnetisation

The code used to determine the energy of the lattice is shown below:

   def energy(self):
       ener = 0.0
       for i in range(self.n_rows):
           for j in range (self.n_cols):
           # 1st row, 1st column corner
               if j == 0:
           # Last row, 1st column corner
                   if i == self.n_rows-1:
                       ener += self.lattice [i,j]*self.lattice [0,j]
                       ener += self.lattice [i,j]*self.lattice [i,j+1]
           # Rest of 1st column
                   else:
                       ener += self.lattice [i,j]*self.lattice [i+1,j]
                       ener += self.lattice [i,j]*self.lattice [i,j+1]
           # 1st row, last column corner
               elif j == self.n_cols-1:
           # Last row, last column corner
                   if i == self.n_rows-1:
                       ener += self.lattice [i,j]*self.lattice [0,j]
                       ener += self.lattice [i,j]*self.lattice [i,0]
           # Rest of last column
                   else:
                       ener += self.lattice [i,j]*self.lattice [i+1,j]
                       ener += self.lattice [i,j]*self.lattice [i,0]
               else:
           # Rest of last column
                   if i == self.n_rows-1:
                       ener += self.lattice [i,j]*self.lattice [0,j]
                       ener += self.lattice [i,j]*self.lattice [i,j+1]
           # Rest of lattice
                   else:
                       ener += self.lattice [i,j]*self.lattice [i+1,j]
                       ener += self.lattice [i,j]*self.lattice [i,j+1]
       energy = ener*-1
       "Return the total energy of the current lattice configuration."
       return energy

Testing the Code

This code was then tested using the file Ilcheck.py which created three lattices and checked the energy and magnetisation of these lattices. One configuration corresponded to the energy minimum, one to the energy maximum, and one to an random intermediate state: this was found to work as expected. The result of this and a checkpoint, while that did not work as expected, are shown below.

Table 2: The Correct and an Incorrect Checkpoint Files
Figure 2: A Incorrect Checkpoint File Figure 3: The Correct Checkpoint File

Introduction to the Monte Carlo simulation

Average Energy and Magnetisation

A system that contains 100 lattice cells with each lattice cell being allowed to be one of two states, spin up or spin down, has a certain number of states available to it. Using the equation for the number of micro-states Ω=nN, where n is the number of energy levels available and N is the number of particles, it is found that the number of micro-states will be Ω=2100 or 1.27x1030 states. If the computer can analyse 109 configurations in a second then it would take 1.27x1021 seconds to analyse all the configurations. As this number is too large to be meaningful it has been converted to years, and in years it would take 4.02x1013 years to analyse all the configurations, longer than the age of the universe!

Importance Sampling.

The code used for the Monte Carlo simulation is shown below:

   def montecarlostep(self, T):
       # complete this function so that it performs a single Monte Carlo step
       energy = self.energy()
       magnetisation = self.magnetisation()
       #the following two lines will select the coordinates of the random spin for you
       random_i = np.random.choice(range(0, self.n_rows))
       random_j = np.random.choice(range(0, self.n_cols))
       #the following line will choose a random number in the rang e[0,1) for you
       random_number = np.random.random()
       # Randomly change a spin
       if self.lattice [random_i, random_j] == 1:
       # From 1 go to -1
           self.lattice [random_i, random_j] = -1
       else:
       # From -1 go to 1
           self.lattice [random_i, random_j] = 1
       energy1 = self.energy()
       magnetisation1 = self.magnetisation()
       deltaenergy = energy1 - energy
       if deltaenergy < 0:
       # Energy goes down
           energy = energy1
           magnetisation = magnetisation1
       elif random_number <= exp(- deltaenergy/ T):
       # Energy goes up but smaller than random number
           energy = energy1
           magnetisation = magnetisation1
       else:
       # Energy goes up and larger than random number
           self.lattice [random_i, random_j] = - self.lattice [random_i, random_j]
       self.E += energy
       self.E2 += energy**2
       self.M += magnetisation
       self.M2 += magnetisation**2
       self.n_cycles += 1
       return energy, magnetisation
   def statistics(self):
       # complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
       aveE = self.E/self.n_cycles
       aveE2 = self.E2/self.n_cycles
       aveM = self.M/self.n_cycles
       aveM2 = self.M2/self.n_cycles
       return aveE, aveE2, aveM, aveM2
Table 3: Monte Carlo Simulation Results
Figure 4: Example 1 of the Minimum energy
Figure 5: Example 2 of the Minimum energy
Average Energy -1.47164536741 -1.4658836689
Average Magnetisation -0.616646698616 0.606508668904
Average Energy2 2.34781017039 2.34409081376
Average Magnetisation2 0.475055536142 0.471877403174

The code was then tested using the file ILanim.py; this ran a Monte Carlo simulation of an eight by eight lattice and displayed the output at a temperature of 1 temperature unit. The results of two runs of this testing, including the display output, can be seen opposite. In one case all the spins in the system became +1, while in the other case the spins in the system became -1. As has been stated earlier, below the Curie temperature the lattice will be ferromagnetic and so the spins will align. This means that it is expected that there will be spontaneous magnetisation.

Accelerating the Code

In order to check if the code has been accelerated it was necessary to check how quick the original code was, as an accelerated code must be quicker. The speed of the original code was tested using the file ILtimetrial.py, which ran 2000 steps of the Monte Carlo simulation. This was done 5 times and then averaged, as each time it ran there were slightly different timings. The times of each run and the average can be seen from table 4 below.

Table 4: The Time Taken, in Seconds, for the Original Code to Perform 2000 Monte Carlo Steps
Attempt 1 Attempt 2 Attempt 3 Attempt 4 Attempt 5 Average Time Taken
6.5321323358111965 6.5576305262353145 6.548430656233478 6.566359750713659 6.561720323517164 6.553254719

The standard error of the sample was found to be 0.006049722. Using the numpy.sum function, a new code for the magnetisation was made. This code is shown below:

def magnetisation(self):

       magnetisation = np.sum(self.lattice)
       return magnetisation

The new code for determining the energy is show below. This was done using the numpy multiply and numpy roll functions.

   def energy(self):
       np.multiply(self.lattice, np.roll(self.lattice,1,axis=0))
       np.multiply(self.lattice, np.roll(self.lattice,1,axis=1))
       return energy

The code still gave the expected result when tested using the file ILcheck.py, this can be seen to the below:

Figure 6: The Checkpoint File of the Re-optimised Code
Figure 6: The Checkpoint File of the Re-optimised Code

The file ILtimetrial was used again in order to find how long the new code would take to perform 2000 Monte Carlo steps. As can be seen below, the new code was much faster than the original, 0.387094820036961 seconds compared to 6.553254719 seconds. The standard error was also decreased, 0.000177801 compared to 0.006049722. This means that the new code worked more quickly than the original and was more consistent in the time taken to perform 2000 Monte Carlo steps.


Table 5: The Time Taken, in Seconds, for the Accelerated Code to Perform 2000 Monte Carlo Steps
Attempt 1 Attempt 2 Attempt 3 Attempt 4 Attempt 5 Average Time Taken
0.384887314998096 0.38467463684389713 0.3845507255513354 0.3855771603227396 0.3848634022924955 0.387094820036961

The Effect of Temperature

Correcting the Averaging Code

The behaviour of the lattice using the Ising model can now be tested in order to probe further into the Curie temperature, the area where the change of domination between the enthalpic and entropic terms takes place. As it takes time for the system to reach the equilibrium state this will affect the outcome of any experiment taking place, so it will be necessary to ignore the first few Monte Carlo steps until the equilibrium state is reached; as has been seen earlier, the energy sharply decreased before the minimum energy was reached. The energy and magnetisation should only be averaged after equilibrium has been reached. Lattices of different sizes and different temperature had 150000 Monte Carlo steps performed on them using the file Il.finalframe.py and the number of steps required to reach equilibrium are shown below.

Table 6: The Number of Monte Carlo Steps Needed for the Lattice to reach the Minimum Energy
Lattice Size Temperature Monte Carlo Steps Needed Picture
8x8 1 656
Figure 7: 8x8 Lattice, T=1, Run 1
Figure 7: 8x8 Lattice, T=1, Run 1
8x8 1 459
Figure 8: 8x8 Lattice, T=1, Run 2
Figure 8: 8x8 Lattice, T=1, Run 2
8x8 1 478
Figure 9: 8x8 Lattice, T=1, Run 3
Figure 9: 8x8 Lattice, T=1, Run 3
8x8 1 571
Figure 10: 8x8 Lattice, T=1, Run 4
Figure 10: 8x8 Lattice, T=1, Run 4
8x8 1.5 -
Figure 11: 8x8 Lattice, T=1.5
Figure 11: 8x8 Lattice, T=1.5
4x4 1 19
Figure 12: 4x4 Lattice, T=1, Run 1
Figure 12: 4x4 Lattice, T=1, Run 1
4x4 1 59
Figure 13: 4x4 Lattice, T=1, Run 2
Figure 13: 4x4 Lattice, T=1, Run 2
4x4 1 47
Figure 14: 4x4 Lattice, T=1, Run 3
Figure 14: 4x4 Lattice, T=1, Run 3
4x4 1 47
Figure 15: 4x4 Lattice, T=1, Run 4
Figure 15: 4x4 Lattice, T=1, Run 4
4x4 1.5 -
Figure 16: 4x4 Lattice, T=1.5
Figure 16: 4x4 Lattice, T=1.5

As can be seen from table 6, the amount of time required for the energy to be minimised varied depending on a number of factors, the temperature of the system and the size of the lattice. Decreasing the number of lattice points from an 8x8 lattice to a 4x4 lattice, a decrease of 75%, resulted, on average, in the number of Monte Carlo steps needed for the energy to be minimised being decreased by a factor of 10. The number of Monte Carlo steps needed when the temperature was increased from 1 to 1.5 is not shown. The reason for this is clear if the diagrams within table 6 are examined: at a higher temperature more high level energy levels will be populated and so the minimum energy will be less easily observed, indeed has not been observed in either the 8x8 or 4x4 lattice at 1.5. The number of Monte Carlo steps needed was also different in each run as can be seen from table 6. It is better to overestimate the number of steps that will be required than to underestimate. All the values of Monte Carlo steps required for a 4x4 lattice were less than 100 and for a 8x8 lattice they were all less than 1000: therefore the first 100 steps should be ignored for a 4x4 lattice and the first 1000 steps should be ignored for an 8x8 lattice. The file ILfinalframe.py was then modified with the following addition and alterations - a new variable was made n_ignore, which was the number of Monte Carlo steps that were to be ignored.

       if self.n_cycles >= self.n_ignore:
           self.E += energy
           self.E2 += energy**2
           self.M += magnetisation
           self.M2 += magnetisation**2
           self.n_cycles += 1
       else:
           pass
       return energy, magnetisation
  
   def statistics(self):
       # complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
       aveE = self.E/(self.n_cycles-self.n_ignore)
       aveE2 = self.E2/(self.n_cycles-self.n_ignore)
       aveM = self.M/(self.n_cycles-self.n_ignore)
       aveM2 = self.M2/(self.n_cycles-self.n_ignore)
       return aveE, aveE2, aveM, aveM2, self.n_cycles

This code was altered so that, if the count was below the number of Monte Carlo steps needed before the lattice that reached the minimum energy, then the energy and magnetisation were not added to. The second change was that the average energy, energy squared, magnetisation and magnetisation squared were altered so that they only averaged for the number of counts that they had had added. In table 7. it is possible to see the number of ignored Monte Carlo steps for different size lattice.

Table 7: The Number of Monte Carlo Steps Ignored for Different Lattice Sizes
Lattice Size Monte Carlo Steps Ignored Picture
2x2 0
Figure 17: 2x2 Lattice
Figure 17: 2x2 Lattice
4x4 100
Figure 18: 4x4 Lattice
Figure 18: 4x4 Lattice
8x8 1000
Figure 19: 8x8 Lattice
Figure 19: 8x8 Lattice
16x16 5000
Figure 20: 16x16 Lattice
Figure 20: 16x16 Lattice
32x32 100000
Figure 21: 32x32 Lattice
Figure 21: 32x32 Lattice

Running Over a Range of Temperatures

The following code was changed in the file IsingLattice.py and the empty lists El and Ml were made:

   def statistics(self):
       # complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
       aveE = self.E/(self.n_cycles-self.n_ignore)
       aveE2 = self.E2/(self.n_cycles-self.n_ignore)
       aveM = self.M/(self.n_cycles-self.n_ignore)
       aveM2 = self.M2/(self.n_cycles-self.n_ignore)
       sdeve = np.std(self.El)
       sdevm = np.std(self.Ml)
       serre = sdeve / sqrt(self.n_cycles-self.n_ignore)
       serrm = sdevm / sqrt(self.n_cycles-self.n_ignore)
       return aveE, aveE2, aveM, aveM2, self.n_cycles, serre, serrm

From the file ILtemperaturerange.py the code was changed to:

   from IsingLattice import *
   from matplotlib import pylab as pl
   import numpy as np
   n_rows = 4
   n_cols = 4
   n_ignore = 100
   il = IsingLattice(n_rows, n_cols)
   il.lattice = np.ones((n_rows, n_cols))
   spins = n_rows*n_cols
   runtime = 10000
   times = range(runtime)
   temps = np.arange(0.25, 5.0, 0.25)
   energies = []
   magnetisations = []
   energysq = []
   magnetisationsq = []
   serrma = []
   serren = []
   for t in temps:
       for i in times:
           if i % 100 == 0:
               print(t, i)
           energy, magnetisation = il.montecarlostep(t)
       aveE, aveE2, aveM, aveM2, n_cycles, serre, serrm = il.statistics()
       energies.append(aveE)
       energysq.append(aveE2)
       magnetisations.append(aveM)
       magnetisationsq.append(aveM2)
       serrma.append(serrm)
       serren.append(serre)
       #reset the IL object for the next cycle
       il.E = 0.0
       il.E2 = 0.0
       il.M = 0.0
       il.M2 = 0.0
       il.n_cycles = 0
   fig = pl.figure()
   enerax = fig.add_subplot(2,1,1)
   enerax.set_ylabel("Energy per spin")
   enerax.set_xlabel("Temperature")
   enerax.errorbar(temps, np.array(energies)/spins, xerr=0, yerr= np.array(serren)/spins)
   enerax.set_ylim([-2.1, 2.1])
   magax = fig.add_subplot(2,1,2)
   magax.set_ylabel("Magnetisation per spin")
   magax.set_xlabel("Temperature")
   magax.errorbar(temps, np.array(magnetisations)/spins, xerr=0, yerr= np.array(serrma)/spins)
   magax.set_ylim([-1.1, 1.1])
   enerax.plot(temps, np.array(energies)/spins)
   magax.plot(temps, np.array(magnetisations)/spins)
   pl.show()
   final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))
   np.savetxt("8x8.dat", final_data)
Table 8: The Change in Energy and Magnetisation of Lattices with Temperature
Figure 22: 8x8 Lattice
Figure 22: 8x8 Lattice
Figure 23: 4x4 Lattice
Figure 23: 4x4 Lattice
8x8 4x4

The Effect of System Size

From the previous section the energies and spins of different lattice sizes at different temperatures were found so that the onset of phase transition could be seen. It is possible to show all the data sets on the same graph and by doing this the minimum lattice size needed to accurately model the long range fluctuations that occur within the system. The code used in order to do this is shown below:

   def temp(file):
       data = np.loadtxt(file)
       temp = data[:,0]
       return temp
   
   def energy(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energ = data[:,1]
       energy = energ/num
       return energy
   
   def magn(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       mag = data[:,3]
       magn = mag/num
       return magn
   
   def plot(a,b,c,d,e):
       fig = pl.figure()
       enerax = fig.add_subplot(2,1,1)
       enerax.set_ylabel("Energy per spin")
       enerax.set_xlabel("Temperature")
       enerax.set_ylim([-2.1, 0])
       magax = fig.add_subplot(2,1,2)
       magax.set_ylabel("Magnetisation per spin")
       magax.set_xlabel("Temperature")
       magax.set_ylim([-1.1, 1.1])
       enerax.plot(temp(a), energy(a),temp(b), energy(b),temp(c), energy(c),temp(d), energy(d),temp(e), energy(e))
       enerax.legend(["2x2", "4x4", "8x8", "16x16", "32x32"],bbox_to_anchor=(0., 1.02, 1., .102), loc=3,ncol=5, mode="expand", borderaxespad=0.)
       magax.plot(temp(a), magn(a),temp(b), magn(b),temp(c), magn(c),temp(d), magn(d),temp(e), magn(e))
       pl.legend()
   pl.show()
Figure 24: The Graph Showing the Energies and Magnetisations with Different Lattice Sizes
Figure 24: The Graph Showing the Energies and Magnetisations with Different Lattice Sizes

The graph made by this code is shown on figure 24, to the right and it is clear that an 8x8 lattice is the minimum size necessary to observe the long term fluctuations. From the energy part of the graph the energies are almost the same from an 8x8 lattice and larger and so it is pointless to compute for a larger lattice as it will not improve the results.

Determining the Heat Capacity

Increasing the temperature above the Curie temperature induces a phase transition. This means that the magnetisation of the system will rapidly drop and from this the heat capacity of the system can be found using the relationship C=ET=Var[E]kBT2. It is known that the heat capacity should become very strongly peaked at the phase transition temperature and the code used to plot a graph showing the heat capacity versus temperature for each of lattice size is:

   def temp(file):
       data = np.loadtxt(file)
       temp = data[:,0]
       return temp
   
   def temp2(file):
       data = np.loadtxt(file)
       temp = data[:,0]
       temp2 = temp*temp
       return temp2
   
   def energy(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energy = data[:,1]
       return energy
   
   def energy2(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energy2 = data[:,2]
       return energy2
   
   def var(file):
       var = energy2(file) - (energy(file)*energy(file))
       return var
   
   def heat(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       heat = (var(file)/temp2(file))/num
       return heat
   
   def plot(a,b,c,d,e):
       fig = pl.figure()
       pl.ylabel("Heat Capacity")
       pl.xlabel("Temperature")
       pl.ylim([0, 1.2])
       pl.plot(temp(a), heat(a),temp(b), heat(b),temp(c), heat(c),temp(d), heat(d),temp(e), heat(e))
       pl.legend(["2x2", "4x4", "8x8", "16x16", "32x32"],bbox_to_anchor=(0., 1.02, 1., .102), loc=3,ncol=5, mode="expand", borderaxespad=0.)
       pl.legend()
       pl.show()

The result of this code can be seen below:

Figure 25: The Heat Capacity versus Temperature

Locating the Curie Temperature

Comparison of Python Data with C++ Data

Figure 26: Graph Showing the Difference Between the Python and C++ Data
Figure 26: Graph Showing the Difference Between the Python and C++ Data

It is clear from the previous section that the heat capacity becomes strongly peaked in the vicinity of the critical temperature around 2 to 2.5. The peak of the heat capacity became more sharp as the lattice size was increased. If there was an infinite size lattice then the critical temperature would diverge at the Curie temperature. Obviously this is not possible and in fact, not only does the heat capacity not diverge with different lattice sizes but the Curie temperature also changes. However the temperature at which the maximum heat capacity is found is modeled using the equation TC,L=AL+TC,inf, TC,L is the Curie temperature of an LxL lattice and TC, inf is the Curie temperature of an infinity large lattice; A is a constant which is not important. It is possible to compare a C++ program, that has run much longer simulations, to the data computed using Python. The code used for this is shown below and the comparison for an 8x8 lattice is shown to the right:

   def temp(file):
       data = np.loadtxt(file)
       temp = data[:,0]
       return temp
   
   def temp2(file):
       data = np.loadtxt(file)
       temp = data[:,0]
       temp2 = temp*temp
       return temp2
   
   def energy(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energy = data[:,1]
       return energy
   
   def energy2(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energy2 = data[:,2]
       return energy2
   
   def energy2(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       data = np.loadtxt(file)
       energy2 = data[:,2]
       return energy2
   
   def var(file):
       var = energy2(file) - (energy(file)*energy(file))
       return var
   
   def heat(file):
       filenum = int(file.split('x')[0])
       num = filenum*filenum
       heat = (var(file)/temp2(file))/num
       return heat
   
   def cap(file):
       data = np.loadtxt(file)
       cap = data[:,5]
       return cap
   
   def plot(a,b):
       fig = pl.figure()
       pl.ylabel("Heat Capacity")
       pl.xlabel("Temperature")
       pl.ylim([0, 1.3])
       pl.plot(temp(a), heat(a),temp(b), cap(b))
       pl.legend(["My data", "C++ Data"])
       pl.legend()
       pl.show()

Polynomial fitting

In order to find where the heat capacity is at a maximum, the data will be fitted to a polynomial. This is done using the polyfit and polyval functions from NumPy. The result of this can be seen below with a 11th order polynomial fit along with the code used.

   def plot(file):
       data = np.loadtxt(file) #assume data is now a 2D array containing two columns, T and C
       T = data[:,0] #get the first column
       C = data[:,5] # get the last column
       fit = np.polyfit(T, C, 101) # fit a third order polynomial
       T_min = np.min(T)
       T_max = np.max(T)
       T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
       fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
       pl.ylabel("Heat Capacity")
       pl.xlabel("Temperature")
       pl.ylim([0, 1.3])
       pl.plot(T, C, T_range, fitted_C_values)
       pl.legend(["C++ data", "Polynomial"])
       pl.show()
Figure 27: Graph Showing the C++ Data and an 11th Order Polynomial
Figure 27: Graph Showing the C++ Data and an 11th Order Polynomial

Fitting in a particular temperature range

This code could then be modified so that it only fitted the data in the region required. The polynomial used to fit the data was a 3rd order one. The change is shown below:

   def plot(file):
       data = np.loadtxt(file) #assume data is now a 2D array containing two columns, T and C
       T = data[:,0] #get the first column
       C = data[:,5] # get the last column
       T_min = 2.0
       T_max = 2.6
       T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
       selection = np.logical_and(T > T_min, T < T_max) #choose only those rows where both conditions are true
       peak_T_values = T[selection]
       peak_C_values = C[selection]
       fit = np.polyfit(peak_T_values, peak_C_values, 3)
       fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
       pl.ylabel("Heat Capacity")
       pl.xlabel("Temperature")
       pl.ylim([0, 1.3])
       pl.plot(T, C, T_range, fitted_C_values)
       pl.legend(["C++ data", "Polynomial"], bbox_to_anchor=(0., 1.02, 1., .102), loc=3,ncol=5, mode="expand", borderaxespad=0.)
       pl.show()
Figure 28: Graph Showing the C++ Data and an 3rd Order Polynomial
Figure 28: Graph Showing the C++ Data and an 3rd Order Polynomial

Finding the peak in C

Table 8: The Values of Tmax and Cmax for Different Lattice Sizes
Lattice Size Tmax Cmax
2x2 2.52252252 0.414795316278
4x4 2.45345345 0.816615090356
8x8 2.36276276 1.19327965652
16x16 2.31931932 1.56571895874
32x32 2.29479479 1.83746564575

The code was then altered in order to find the maximum values of the heat capacity, and the temperature this occurred at, for different lattice sizes: from this the Curie temperature could be found. The values of Tmax and Cmax for the different size lattices can be seen on the right in table 8 (for 2x2, 4x4, and 8x8 a 3rd order polynomial was just, for a 16x16 a 5th order polynomial was used and for the 32x32 lattice a 7th order polynomial was used). The result from changing the code can be seen result below:

   def plot(file):
       data = np.loadtxt(file) #assume data is now a 2D array containing two columns, T and C
       T = data[:,0] #get the first column
       C = data[:,5] # get the last column
       T_min = 2.0
       T_max = 2.5
       T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
       selection = np.logical_and(T > T_min, T < T_max) #choose only those rows where both conditions are true
       peak_T_values = T[selection]
       peak_C_values = C[selection]
       fit = np.polyfit(peak_T_values, peak_C_values, 7)
       fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
       Cmax = np.max(fitted_C_values)
       Tmax = T_range[fitted_C_values == Cmax]
       pl.ylabel("Heat Capacity")
       pl.xlabel("Temperature")
       pl.ylim([0, 2])
       pl.plot(T, C, T_range, fitted_C_values)
       pl.legend(["C++ data", "Polynomial"], bbox_to_anchor=(0., 1.02, 1., .102), loc=3,ncol=5, mode="expand", borderaxespad=0.)
       pl.show()
       print (Cmax)
       print (Tmax)
Figure 29: Finding TC,inf
Figure 29: Finding TC,inf

These results could be ploted using the equation TC,L=AL+TC,inf, this was done on Microsoft Excel for simplicity and this can be seen from figure 29. The equation of the trend-line of this data is also shown and gives the value of TC,inf to be 2.2962. This compares quite well to a literature value of TC,inf of 2.2691. This could be improved by using a higher order polynomial in order to more accurately find Cmax and therefore more accurately find Cinf.

Reference

1. J. Kotze, Introduction to Monte Carlo methods for an Ising Model of a Ferromagnet, http://arxiv.org/pdf/0803.0217.pdf