代写Ising model、代写phase transitions、代做Python、Python编程语言调试

日期：2019-11-09 11:27

The Ising model and phase transitions

Remarks on completing the module

This assignment is summatively assessed. It is imperative that you submit the notebook on time.

Ising model

The task for this assignment is to implement the Ising Model introduced in the lecture. The structure in terms of

a code skeleton provided below needs to be followed. Otherwise the automatic tests, which allow you to test

different parts of you implementation, will not work.

We consider an Ising model, in which the interaction energy of spin is calculated from

where the sum is over the 4 nearest neighbours of .

We will restrict the calculation to a 2 dimensional grid throughout, hence 4 nearest neighbours, ,

. Notice that the expression for is different from the form considered in the lecture.

To simplify the calculations, we will use the dimensionless interaction energy

where

in the following. Here, is Boltzmann's constant, and is the temperature.

Given all spin states, we calculate the ensenble-averaged macroscopic magnetization , as

The brackets denote the ensemble average. The parameter has the dimensions of energy, and

consequently is dimensionless.

Follow the numbered steps in the following cells.

The cells below describe how to proceed, step-by-step. Begin by reading through all steps, comparing the

instructions to the Ising model discussed in the lecture. Complete the assignment using the cells below.

Several cells allow you to test your implementation step by step.

1. Set up the regular grid

Set up a 2D grid in the form of a python class, called Grid . The class should contain

the spin array

the value of J

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 2/10

The spin array should be a 2D array of dimension , with . We will address a particular spin with its 2D

Cartesian coordinates , where and are the indices in the 2D array. So, for

example, spin refers to the spin located at vertex of the grid.

Initialize the spins on the grid randomly, meaning each spin can be either up, , or down, , with

equal probability.

When performing calculations on the grid, we will assume periodic boundary conditions

no marks

2. Calculate the energy

Write a method, energy , as part of the class Grid , which calculates the interaction energy of a given spin,

, by summing over its four nearest neighbours. The function should take the grid array, , and the cell

indices and as parameters. It should return a python tuple containing two dimensionless energies

corresponding to the energy of the current spin state of cell , and the energy of the flipped

spin state .

This means that for a cell with spin state , the method should return and

vice versa.

Remember to account for periodic boudnary conditions on the grid.

You can test the implementation of this method using the test cells provded. What are the interaction energies

of cells (6,6), (15,0) and (31, 17) of the assignment grid given below (please include the answer to this question

in the PDF you hand in).

3. Calculate the probability of flipping a spin

The probability that the spin at vertex is flipped depends on the spin states of its neighbours and the

value of as explained in the lecture.

Write a method prob_flip which calculates the probability that spin is flipped, given the (dimensionless)

interaction energies for the current state and the flipped state .

The probability for a flip is given by

You can test the implementation of this method using the test cells provided. What are the probabilities for cells

(12, 12), (18,0) and (31, 12) of the assignment grid given below (please include the answer to this question in

the PDF you hand in)?

2 marks

4. Calculate the macroscopic magnetisation,

Write a method which calculates the current macroscopic magnetisation of a given grid, and add it to the

Grid class. The function should take the grid-array as a parameter and return the mean, macroscopic

magnetisation,

M

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 3/10

You can test the implementation of this method using the test cells provded. Calculate the magnetisation of the

assignment grid. State the answer to 3 significant digits on the PDF file you hand in.

5. Red-black sweep

Write a method to sweep over all spins in turn, first in (say) and then in in a red-black pattern. Red-black

means, first loop over all the red cells on a red-black chessboard pattern (looping over them first, then ).

Once all the red cells are done, then update all the black cells. Add this method to the Grid class.

For each spin in turn, flip the spin or not following the criterion discussed in the lecture. This means that the

spin in each cell in turn should be flipped with a probability (as discussed in step 3).

You can use the methods implemented in step 2 and 3.

6. Thermalisation and magnetisation measurement

Starting from a random configuration, the system needs to be evolved over a number of full red-black sweeps

in order to reach thermal equilibrium. This thermalization is part of the method you develop in this step.

Write a method that starts by sweeping the grid times to allow for the system to thermalize before you

carry out any measurements.

Next, the method should perform a further sweeps, while in addition computing and recording the

value of after every sweep. Use the method you developed in step 4 and the sweep implementaton of step

5.

and are input parameters of the method. The method should return a numpy array of length

, containing the magentisations measured after each measurement sweep.

Add this method to the Grid class.

7. Thermalisation

Plot the magnetisation over time for 1000 full mesh sweeps for and (include the thermalisation

period in the plot). Include this plot in the PDF you hand in. Save the plot to a file 'Thermalization.pdf'

4 marks

β = 0.1, 0.8 1.6

8. Ensemble-averaged magnetisation, , as a function of

Once the system has thermalized (a good choice is thermalisation sweeps), measure the timeaveraged

magnetisation over 1000 sweeps. From this, estimate the ensemble-averaged magnetisation, . Plot

as a function of for . Include this plot in the PDF you

hand in. Save the plot to a file 'Magnetisation.pdf'.

M¯ β

Ntherm = 400

M¯

M¯ 1/β β = 1.6, 1.3, 1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.4, 0.1

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 4/10

p g p

4 marks

9. Critical temperature

The critical temperature, , in the mean field approximation, is given by

where h is the number of nearest-neighbours, as discussed in the lecture. Use this to calculate the critical value

which corresponds to the critical temperature and mark it in the plot produced in step 8 (If you could not

produce the plot in step 8 you may also mention the numerical value of on the solution you hand in).

2 marks

10. Mean field approximation

The ensemble-averaged value of the mean magnitization in the mean field approximation, is the solution of

as discussed in the lecture. This equation can not be solved analytically for , for given .

Rewrite the equation in terms of using the relation between and derived before, and solve the resulting

formula numerically for . You may want to use the root finding

methods implemented in previous exercises. Redo the plot of step 8, but on top of the numerical result, over

plot the mean field approximation. Add a legend to the plot which allows to distinguish the solution obtained in

step 8 from the mean field approximation.

Save the plot as 'MeanField.pdf'

4 marks

β = 1.6, 1.3, 1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.4, 0.1

Solution

Use the cells below to complete the assignment. Use the test cells provided to make sure you are on the right

track.

In [1]:

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

import sys

import os

sys.path.append(os.getcwd())

from scipy.interpolate import CubicSpline

import pickle

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 5/10

In [ ]:

You can use the cells below to test your implementation

class Grid:

def __init__(self, size, beta):

'''This function initialises the grid, i.e. it sets the

grid size, the value of beta and initialises the cells of the

grid with randomly chosen 'plus' (1) or 'minus' (-1) states.'''

# set self.size, self.beta, and self.cells

# self.cells is a 2D array, so that self.cells[i,j]=+1 or -1, the spin in gr

# YOUR CODE HERE

raise NotImplementedError()

def energy(self, i, j, beta, grid):

'''This function calculates the energies 'e_plus' and 'e_minus'

corresponding to the two possible states 'plus' and 'minus' for the spin at

of a given grid with a given value of beta.

returns: the two energy states 'e_current' and 'e_flip' as a tuple.

'e_current' is the energy of the spin (i,j) in its current spin state

'e_flip' is the energy of the spin (i,j) if you were to flip its spin

'''

# YOUR CODE HERE

raise NotImplementedError()

return e_current, e_flip

def prob_flip(self, e_current, e_flip):

'''This function calculates the probability of a spin flip

for a given spin, given the energies e_current and e_flip

of the current and the flipped state for the cell.

returns: the probability for the flip'''

# YOUR CODE HERE

raise NotImplementedError()

return probability_for_flip

def sweep(self):

'''This function carries out a single red-black sweep.

returns: nothing. For each spin in turn, it compute the probability for

flipping the spin, using the prob_flip function. Comparing the probablity,

it draws a random number to decide whether or not the flip the spin '''

# YOUR CODE HERE

raise NotImplementedError()

def magnetisation(self, grid):

'''This function calculates the mean magnetisation of all the spin in the gr

returns: the mean magnetisation M'''

# YOUR CODE HERE

raise NotImplementedError()

return M

def do_sweeps(self, n_therm, n_measure):

'''This function carries out n_therm thermalisation sweeps and n_measure mea

At the end of each measurement sweep the average magnetisation is computed a

returns: an array of length 'n_measure' containing the recorded magnetisatio

It uses the sweep function, and the magnitization function'''

# YOUR CODE HERE

raise NotImplementedError()

return magnetisation

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 6/10

Test and assignment grids

The cell below loads the test grids, test_grid_1 and test_grid_2 , and their corresponding values for

, test_beta_1 and test_beta_2 .

In addition, it loads assignement_grid and assignment_beta .

Use the first two grids to test your implementation. Use the third grid to answer the assignment questions.

β

In [3]:

2. Interaction energy calculation

The cell below allows you to test your interaction energy calculation method. If it does not return an error your

implementation might be correct.

In [ ]:

grid 1 loaded, size= 32 , beta= 1.6

grid 2 loaded, size= 32 , beta= 0.8

grid 3 loaded, size= 32 , beta= 0.1

filename = 'test_grid_1.pickle'

f = open(filename, 'rb')

(test_grid_1, test_beta_1) = pickle.load(f)

f.close()

filename = 'test_grid_2.pickle'

f = open(filename, 'rb')

(test_grid_2, test_beta_2) = pickle.load(f)

f.close()

filename = 'assignment_grid.pickle'

f = open(filename, 'rb')

(assignment_grid, assignment_beta) = pickle.load(f)

f.close()

print(" grid 1 loaded, size=", len(test_grid_1)," , beta= ", test_beta_1)

print(" grid 2 loaded, size=", len(test_grid_2)," , beta= ", test_beta_2)

print(" grid 3 loaded, size=", len(assignment_grid)," , beta= ", assignment_beta)

g1 = Grid(len(test_grid_1), test_beta_1)

g2 = Grid(len(test_grid_2), test_beta_2)

cells = [(6,6), (15,0), (31,17)]

energies_1 = [(0.0, 6.4), (0.0, 6.4), (0.0, 6.4)]

energies_2 = [(0.8, 2.4000000000000004), (1.6, 1.6), (2.4, 0.8)]

for c, cell in enumerate(cells):

i = cell[0]

j = cell[1]

e_1 = g1.energy(i, j, test_beta_1, test_grid_1)

e_2 = g2.energy(i, j, test_beta_2, test_grid_2)

assert(np.isclose(e_1, energies_1[c]).all())

assert(np.isclose(e_2, energies_2[c]).all())

print("Your implementation might be correct!")

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 7/10

3. Probability calculation

The cell below allows you to test your probability calculation method. If it does not return an error your

implementation might be correct.

In [ ]:

4. Magnetisation calculation

The cell below allows you to test your magnetisation method. If it does not return an error your implementation

might be correct.

In [ ]:

The following hidden cell uses the assignment grid to test

the calculation of energies 2 marks

the calculation of the probabilities 2 marks

the calculation of the magnetization 2 marks

5 Implement the red-black sweep

6 Implement the thermalization step

7. Thermalisation

Plot the magnetisation over time for 1000 full mesh sweeps for and (include the thermalisation

period in the plot). Include this plot in the PDF you hand in. Save the plot as file 'Thermalisation.pdf'

4 marks

β = 0.1, 0.8 1.6

g1 = Grid(len(test_grid_1), test_beta_1)

energies = [(0.1, 0.3), (0.2, 0.2), (0.3, 0.1), (1.5, 1.6), (0.1, 1.6), (0.8, 2.4)]

probabilities = [0.8187307530779818, 1, 1, 0.9048374180359595, 0.22313016014842982,

this_prob = []

for i, e in enumerate(energies):

this_prob.append(g1.prob_flip(e[0], e[1]))

assert(np.isclose(this_prob, probabilities).all())

print("Your implementation might be correct!")

g1 = Grid(len(test_grid_1), test_beta_1)

assert(np.isclose(g1.magnetisation(test_grid_1),0.193359375))

g2 = Grid(len(test_grid_2), test_beta_2)

assert(np.isclose(g2.magnetisation(test_grid_2),-0.3203125))

print("Your implementation might be correct!")

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 8/10

In [ ]:

8. Ensemble-averaged magnetisation, , as a function of

Once the system has thermalized (a good choice is thermalisation sweeps), measure the timeaveraged

magnetisation over 1000 sweeps. From this, estimate the ensemble-averaged magnetisation, . Plot

as a function of for . Include this plot in the PDF you

hand in. Save the plot to a file 'Magnetisation.pdf'.

Perform the following steps:

for each of the listed values of beta:

create a random grid

sweep the for N=400 initial 'thermalization' steps

then sweep for another 1000 steps, calculating after every sweep

use this to compute the ensemble-averaged magnitization, for that value of

plot as a function of

4 marks

# for each value of beta=[0.1, 0.8, 1.6]

# generate a grid of size=32, with the given value of beta

# perform N=1000 red-black sweeps (each sweep runs over the full 32x32 grid)

# calculate the mean magnetization, M for each sweep

# Plot M as a function of sweep number.

# You may want to use some of the plotting commands below.

# set-up the figure

fig, ax = plt.subplots(1,1, figsize = (8, 5))

file = "Thermalisation.pdf"

# caculate mag, the average magnetization, for N=1000 sweeps

# pass the value of beta into the plot command to generate the label, as in

# ax.plot(np.arange(len(mag[b])), mag[b], label='beta=%.2f'%beta)

# YOUR CODE HERE

raise NotImplementedError()

print("Calculation finished")

# plot the result, annotate the file, and save the file

ax.set_xlabel(r'$N_{steps}$')

ax.set_ylabel(r'$M$')

ax.set_ylim([-1.05, 1.05])

ax.legend()

plt.savefig(file)

fig.show()

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 9/10

In [ ]:

9. Critical temperature

The critical temperature, , in the mean field approximation, is given by

where h is the number of nearest-neighbours, as discussed in the lecture. Use this to calculate the critical value

which corresponds to the critical temperature and mark it in the plot produced in step 8 (If you could not

produce the plot in step 8 you may also mention the numerical value of on the solution you hand in).

Enter your answer in the box below. If you don't see a box, execute the hidden cell below.

2 marks

10. Mean field approximation

The ensemble-averaged value of the mean magnitization in the mean field approximation, is the solution of

as discussed in the lecture. This equation can not be solved analytically for , for given .

Enter your analytically calculated value of beta_c to 3 sig figs

beta_c= 0

# Step 8: Magnetisation

# set-up the figure

file = "Magnetisation.pdf"

fig, ax = plt.subplots(1,1, figsize = (8, 5))

# the range of values of beta

betas = [1.6,1.3,1.1,1.0,0.9,0.8,0.7,0.6,0.4,0.1]

# Loop over values of beta, computing the ensemble-averaged M for each beta

# name the resulting ensemble-averaged M mean_mags

# It is an array with the same dimension as betas

# YOUR CODE HERE

raise NotImplementedError()

# make the plot

ax.set_xlabel(r'$\beta^{-1}$')

ax.set_ylabel(r'$\bar M$')

ax.set_ylim([-1.05, 1.05])

ax.set_xlim(0,4)

plt.plot(1./np.array(betas), mean_mags)

plt.savefig(file)

fig.show()

2019/11/7 phase_transition

https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Phase_Transition_19/phase_transition.ipynb 10/10

Rewrite the equation in terms of using the relation between and derived before, and solve the resulting

formula numerically for . You may want to use the root finding

methods implemented in previous exercises. Redo the plot of step 8, but on top of the numerical result, over

plot the mean field approximation. Add a legend to the plot which allows to distinguish the solution obtained in

step 8 from the mean field approximation.

The numerical value of versus shows that the transition from at hight to at low is not

infinitely sharp. To quantify where the transition occurs, it might be useful to compute , the value of where

. Calculate this for both the numerical and mean field approximation, and indicate the point

on the plot.

Save the plot as 'MeanField.pdf'

4 marks

β T β

β = 1.6, 1.3, 1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.4, 0.1

M¯ β M¯ = 0 T M¯ > 0 T

β1/2 β

M¯ = 1/2

(β, M) = ( , 1/2)

¯ β1/2

In [ ]:

In [ ]:

# Implement the mean field calculation here: calculate mean_mag_MF, the mean field a

# for a given value of beta. Also implement the calculation of the critical value, b

# MFA

# YOUR CODE HERE

raise NotImplementedError()

# YOUR CODE HERE

raise NotImplementedError()