Rep:IL:COH15

Introduction

In this wikipage, Ising lattice model was built and simulations were performed to investigate impact of temperature on energy, magnetisation and heat capacity of the model lattice using a self-replicated Monte Carlo simulation.

Tasks

Task 1

Show that the lowest possible energy for the Ising model is $E = - D N J$ , where $D$ is the number of dimensions and $N$ is the total number of spins. What is the multiplicity of this state? Calculate its entropy.

In Ising lattice model, it is assumed that interactions between spins in adjacent lattice cells is the only energetic contribution when no external magnetic field is applied. The defined interaction energy of the lattice is $E = - \frac{1}{2} J \sum_{i}^{N} \sum_{j \in n e i g h b o u r s (i)} s_{i} s_{j}$ .

For the lattice to have lowest energy, product of $s_{i} s_{j}$ must be positive, so the spin of neighbouring cell must be identical throughout the whole lattice. As a result, there are 2 possible degenerate states (i.e. all +1 or -1 spin). Sum of neighbour spin is equivalent to 2 times the dimension of the lattice we are modelling, since $s_{i} s_{j}$ always returns as 1, the sum of interactions that all cell unit in the lattice has is equal to N. Therefore, giving the lowest possible energy $E = - D N J$ , after the $\frac{1}{2}$ cancels out with the factor of 2 in the dimension part.

Entropy is defined by $S = k_{B} l n W$ , where $W$ is the number of degenerate states in the lowest energy configuration. In this case, $W = 2$ and so the entropy is $k_{B} l n 2$ .

Task 2

Imagine that the system is in the lowest energy configuration. To move to a different state, one of the spins must spontaneously change direction ("flip"). What is the change in energy if this happens (D=3,N=1000)? How much entropy does the system gain by doing so?

In a 3D lattice with 1000 unit cells, flipping the spin of any one cell is essentially changing the interaction energies with its 6 neighbouring cells from negative to positive, giving an overall $E_{c h a n g e} = + 12 J$ .Each cell in the lattice has equal probability for spin change to occur, which means there are 1000 different equally probable configurations at this state (i.e. $W = 1000$ ). The entropy is now $k_{B} l n 1000$ , which has increased by $k_{B} l n 500$ compared to the lowest possible state.

Task 3

Calculate the magnetisation of the 1D and 2D lattices in figure 1. What magnetisation would you expect to observe for an Ising lattice with D = 3, N=1000 at absolute zero?\|-

Knowing that total magnetisation of the lattice system at low temperatures with large number of parallel spins can be assumed to be $M = \sum_{i} s_{i}$ .

The calculated magnetisation of the 1D lattice in figure 1 is $+ 1$ and same for the 2D lattice.

At absolute zero, a 3D lattice with 1000 unit cells would be at its lowest energy configuration with all parallel spins, which would give $M = + 1000 o r - 1000$ .

Task 4

Complete the functions energy( ) and magnetisation( ), which should return the energy of the lattice and the total magnetisation, respectively. In the energy( ) function you may assume that $J = 1.0$ at all times (in fact, we are working in reduced units in which $J = k_{B}$ , but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment.

The energy( ) function:

   def energy(self):
       "Return the total energy of the current lattice configuration."
       row_energy=0
       for row in self.lattice:
           en_r=0
           for i in range(self.n_rows):
               if row[i]==row[i-1]:
                   en_r = en_r-1
               else:
                   en_r = en_r +1
           row_energy=row_energy+en_r
       col_energy=0
       for j in range(self.n_cols):
           en_c=0
           for i in range(self.n_rows):
               if self.lattice[i,j]==self.lattice[i-1,j]:
                   en_c=en_c-1
               else:
                   en_c=en_c+1
           col_energy=col_energy+en_c
       energy =col_energy+row_energy
       return energy

The magnetisation( ) function:

   def magnetisation(self):
       "Return the total magnetisation of the current lattice configuration."
       magnetisation =np.sum(self.lattice)    
       return magnetisation

Task 5

Run the ILcheck.py script from the IPython Qt console using the command

%run ILcheck.py

The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.

Fig. 1 Image of ILcheck.py* output*

Task 6

How many configurations are available to a system with 100 spins? To evaluate these expressions, we have to calculate the energy and magnetisation for each of these configurations, then perform the sum. Let's be very, very, generous, and say that we can analyse $1 \times 1 0^{9}$ configurations per second with our computer. How long will it take to evaluate a single value of ${⟨ M ⟩}_{T}$ ?

Number of configurations available $= 2^{100}$

Time needed to complete analysis $= 4.01 \times 1 0^{13} y e a r s$

Task 7

Implement a single cycle of the above algorithm in the montecarlocycle(T) function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. You may assume that the energy returned by your energy( ) function is in units of $k_{B}$ ! Complete the statistics( ) function. This should return the following quantities whenever it is called: $⟨ E ⟩$ , $⟨ E^{2} ⟩$ , $⟨ M ⟩$ , $⟨ M^{2} ⟩$ , and the number of Monte Carlo steps that have elapsed.

def montecarlostep(self, T):
       # complete this function so that it performs a single Monte Carlo step
       self.n_cycles=self.n_cycles+1
       energy_0 = self.energy()
       lattice_0=self.lattice.copy()
       #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))
       self.lattice[random_i,random_j]=self.lattice[random_i,random_j]*(-1)
       energy_1=self.energy()
       energy_diff=energy_1-energy_0
       #the following line will choose a random number in the range[0,1) for you
       random_number = np.random.random()
       prob=np.exp(energy_diff/(-1*T))
       if energy_diff>0:
           if prob<random_number:
               self.lattice=lattice_0
       
       M_1=self.magnetisation()
       E_1=self.energy()
       self.E=self.E+E_1
       self.M=self.M+M_1
       self.E2=self.E2+(E_1**2)
       self.M2=self.M2+(M_1**2)
       return E_1,M_1

   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
       return np.array([self.E,self.E2,self.M, self.M2])/self.n_cycles

Task 8

If $T < T_{C}$ , do you expect a spontaneous magnetisation (i.e. do you expect $⟨ M ⟩ \neq 0$ )? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to PNG and attach this to your report. You should also include the output from your statistics( ) function.

When $T < T_{C}$ , spontaneous magnetisation occurs as the system is dominated by the enthalpy. This is due to the fact that entropy gained from spin flipping is not great enough to overcome the increase in enthalpy. However, as temperature increases, effect of entropic drive becomes more dominant and spontaneous magnetisation occur readily.

Fig. 2 Output of ILanim.py

Task 9

Use the script ILtimetrial.py to record how long your current version of IsingLattice.py takes to perform 2000 Monte Carlo steps. This will vary, depending on what else the computer happens to be doing, so perform repeats and report the error in your average!

Result Output: Took 2.907136615182344s

Task 10

Look at the documentation for the NumPy sum function. You should be able to modify your magnetisation( ) function so that it uses this to evaluate M. The energy is a little trickier. Familiarise yourself with the NumPy roll and multiply functions, and use these to replace your energy double loop (you will need to call roll and multiply twice!).

Optimised code for energy( ) and magnetisation( ) functions

   def energy(self):
       "Return the total energy of the current lattice configuration."
       y=self.lattice
       energy=-1*(np.sum(np.multiply(y,np.roll(y,1,axis=0)))+np.sum(np.multiply(y,np.roll(y,1,axis=1))))
       return energy

   def magnetisation(self):
       "Return the total magnetisation of the current lattice configuration."
       magnetisation =np.sum(self.lattice)    
       return magnetisation

Task 11

Use the script ILtimetrial.py to record how long your new version of IsingLattice.py takes to perform 2000 Monte Carlo steps. This will vary, depending on what else the computer happens to be doing, so perform repeats and report the error in your average!

Time recorded after updating code (10 repeats): $0.32691241028169621 \pm 0.0078827643826809655 s$

Task 12, 13 & 14

The script ILfinalframe.py runs for a given number of cycles at a given temperature, then plots a depiction of the final lattice state as well as graphs of the energy and magnetisation as a function of cycle number. This is much quicker than animating every frame! Experiment with different temperature and lattice sizes. How many cycles are typically needed for the system to go from its random starting position to the equilibrium state? Modify your statistics( ) and montecarlostep( ) functions so that the first N cycles of the simulation are ignored when calculating the averages. You should state in your report what period you chose to ignore, and include graphs from ILfinalframe.py to illustrate your motivation in choosing this figure.
Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, with error bars, for an 8 X 8 lattice. Use your intuition and results from the script ILfinalframe.py to estimate how many cycles each simulation should be. The temperature range 0.25 to 5.0 is sufficient. Use as many temperature points as you feel necessary to illustrate the trend, but do not use a temperature spacing larger than 0.5. The NumPy function savetxt() stores your array of output data on disk — you will need it later. Save the file as 8x8.dat so that you know which lattice size it came from.
Repeat the final task of the previous section for the following lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32. Make sure that you name each datafile that your produce after the corresponding lattice size! Write a Python script to make a plot showing the energy per spin versus temperature for each of your lattice sizes. Hint: the NumPy loadtxt function is the reverse of the savetxt function, and can be used to read your previously saved files into the script. Repeat this for the magnetisation. As before, use the plot controls to save your a PNG image of your plot and attach this to the report. How big a lattice do you think is big enough to capture the long range fluctuations?

2 X 2 Lattice (**cycles ignored**= 80)
Fig. 5 T=0.5	Fig. 6 T=1.0	Fig. 7 T=3.0

4 X 4 Lattice (**cycles ignored**= 320)
Fig. 5 T=0.5	Fig. 6 T=1.0	Fig. 7 T=3.0

8 X 8 Lattice (**cycles ignored**= 1280)
Fig. 5 T=0.5	Fig. 6 T=1.0	Fig. 7 T=3.0

16 X 16 Lattice (**cycles ignored**= 5120)
Fig. 5 T=0.5	Fig. 6 T=1.0	Fig. 7 T=3.0

32 X 32 Lattice (**cycles ignored**= 20480)
Fig. 5 T=0.5	Fig. 6 T=1.0	Fig. 7 T=3.0

Task 15

By definition, $C = \frac{\partial ⟨ E ⟩}{\partial T}$ . From this, show that $C = \frac{V a r [E]}{k_{B} T^{2}}$ (where $V a r [E]$ is the variance in $E$ )

$C = \frac{\partial ⟨ E ⟩}{\partial T} = \frac{(⟨ E^{2} ⟩ - {⟨ E ⟩}^{2})}{k_{B} T^{2}} = \frac{(Δ E)^{2}}{k_{B} T^{2}} = \frac{V a r [E]}{k_{B} T^{2}}$

Task 16 & 17

Write a Python script to make a plot showing the heat capacity versus temperature for each of your lattice sizes from the previous section. You may need to do some research to recall the connection between the variance of a variable, \mathrm{Var}[X], the mean of its square \left\langle X^2\right\rangle, and its squared mean \left\langle X\right\rangle^2. You may find that the data around the peak is very noisy — this is normal, and is a result of being in the critical region. As before, use the plot controls to save your a PNG image of your plot and attach this to the report.
A C++ program has been used to run some much longer simulations than would be possible on the college computers in Python. You can view its source code here if you are interested. Each file contains six columns: T, E, E^2, M, M^2, C (the final five quantities are per spin), and you can read them with the NumPy loadtxt function as before. For each lattice size, plot the C++ data against your data. For one lattice size, save a PNG of this comparison and add it to your report — add a legend to the graph to label which is which. To do this, you will need to pass the label="..." keyword to the plot function, then call the legend() function of the axis object (documentation here).

Dimension	2x2	4x4	8x8	16x16	32x32
Runtime	5000	50000	100000	350000	600000
Number of cycles ignored	80	320	1280	5120	20480

The conditions set for simulation were listed in the table above. These conditions were set based on results obtained in Task 12. In Task 15, the relationship between $C$ and $V a r [E]$ is shown. Using this equation, heat capacity values were calculated across the temperature range from 0.25 to 0.50.

2 X 2 Lattice (**cycles ignored**= 80)
Energy Plot
Magnetisation Plot
Heat Capacity Plot

The curvature at low temperature indicated that the variance value is greater than it should be, hinting that the number of cycles was insufficient to allow equilibrium state to be reached. Increasing runtime could help improve experimental data to match with C++ data better. In general, experimental data agreed with C++ data.

4 X 4 Lattice (**cycles ignored**= 320)
Energy Plot
Magnetisation Plot
Heat Capacity Plot

The curvature at low temperature was due to similar reason as mentioned for 2x2 lattice above.

8 X 8 Lattice (**cycles ignored**= 1280)
Energy Plot
Magnetisation Plot
Heat Capacity Plot

Large error bar at low temperature indicates that runtime was insufficient.

16 X 16 Lattice (**cycles ignored**= 5120)
Energy Plot
Magnetisation Plot
Heat Capacity Plot

32 X 32 Lattice (**cycles ignored**= 20480)
Energy Plot
Magnetisation Plot
Heat Capacity Plot

Observations:

In general, all curve above matches with C++ data. The heat capacity per spin increases as temperature increases until a peak value is reached. Then, heat capacity per spin decreases as temperature further increases. The maximum heat capacity per spin in each lattice lie in similar temperature range (approx. from 2.0 to 2.8). From figure , maximum heat capacity of lattice increase as lattice size (L) increases. Increase in lattice size (L) also narrows peak width of the plot, as the increment increases more rapidly in bigger lattice.

Fig. 2 Heat Capacity per spin VS Temperature of Each Lattice Sizes

Task 18 & 19

write a script to read the data from a particular file, and plot C vs T, as well as a fitted polynomial. Try changing the degree of the polynomial to improve the fit — in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.
Modify your script from the previous section. You should still plot the whole temperature range, but fit the polynomial only to the peak of the heat capacity! You should find it easier to get a good fit when restricted to this region.

Task 20

Find the temperature at which the maximum in C occurs for each datafile that you were given. Make a text file containing two colums: the lattice side length (2,4,8, etc.), and the temperature at which C is a maximum. This is your estimate of $T_{C}$ for that side length. Make a plot that uses the scaling relation given above to determine $T_{C, \infty}$ . By doing a little research online, you should be able to find the theoretical exact Curie temperature for the infinite 2D Ising lattice. How does your value compare to this? Are you surprised by how good/bad the agreement is? Attach a PNG of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.

Theoretical value of $T_{C, \infty}$ is 2.269. By plotting $T_{C}$ values collected against $\frac{1}{L}$ , a linear best fit line was generated using Polyfit and return a value of 2.277 for $T_{C, \infty}$ with a percentage difference of 0.353%. This experimental value of $T_{C, \infty}$ matches with the theoretical value from reference journal. Although the points around the best fit line might seem to be quite far from the best fit line, the standard error value is only 0.00251. This shows that the influence of experiment error is quite small. Further simulations could be run on lattice with bigger dimensions to extend and give more data points to the plot, which could increase the reliability of the $T_{C, \infty}$ , as well as performing repeats.

Crashed Image:

https://wiki.ch.ic.ac.uk/wiki/images/2/2e/32x32_fit4_1.8_2.8_600000.png
Order:4; T=1.8-2.9, Curie T:2.301, C++ T=2.30

https://wiki.ch.ic.ac.uk/wiki/images/a/a9/16x16_fit3_2.0_2.8.png
Order:3; T=2.0-2.8; Curie T:2.287, C++ T=2.25

https://wiki.ch.ic.ac.uk/wiki/images/7/7f/8x8_fit6_15_28.png
Order:6; T=1.5-2.8; Curie T:2.308, C++ T=2.31

https://wiki.ch.ic.ac.uk/wiki/index.php?title=File:4x4_fit4_20_28.png
Order:4; T=2.0-2.8; Curie T:2.457, C++ T=2.40

https://wiki.ch.ic.ac.uk/wiki/index.php?title=File:32x32_fit3_2.0_2.8_600000.png
Order:3; T=2.0-2.8; Curie T:2.494, C++ T=2.48