From c892056ea68e52b4ec128d6e943071f3aff78ad3 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Thu, 17 Nov 2022 11:40:45 +0100 Subject: [PATCH] 06: Bart sent me the feedback by mail --- .../exercise_sheet_06.html | 1320 +++++++++++++++++ 1 file changed, 1320 insertions(+) create mode 100644 Exercise sheet 6/feedback/2022-10-18 20:21:28.640362 UTC/exercise_sheet_06.html diff --git a/Exercise sheet 6/feedback/2022-10-18 20:21:28.640362 UTC/exercise_sheet_06.html b/Exercise sheet 6/feedback/2022-10-18 20:21:28.640362 UTC/exercise_sheet_06.html new file mode 100644 index 0000000..a0912c2 --- /dev/null +++ b/Exercise sheet 6/feedback/2022-10-18 20:21:28.640362 UTC/exercise_sheet_06.html @@ -0,0 +1,1320 @@ + + + + + +exercise_sheet_06 + + + + + + + + + + + + + + + + + + +
+
+
+

exercise_sheet_06 (Score: 100.0 / 100.0)

+
+
    + + + + + + + + + + + + + + + + + + + + +
  1. Written response (Score: 25.0 / 25.0)
    2. + + + + + + + + + + +
  2. Test cell (Score: 15.0 / 15.0)
    3. + + + + +
  3. Test cell (Score: 5.0 / 5.0)
    4. + + + + + + + + + + +
  4. Test cell (Score: 10.0 / 10.0)
    5. + + + + +
  5. Coding free-response (Score: 10.0 / 10.0)
    6. + + +
  6. Comment
    7. + + + + + + + + + +
  7. Test cell (Score: 10.0 / 10.0)
    8. + + + + +
  8. Coding free-response (Score: 10.0 / 10.0)
    9. + + + + + + + + + + +
  9. Test cell (Score: 5.0 / 5.0)
    10. + + + + +
  10. Coding free-response (Score: 10.0 / 10.0)
    11. + + + + + + + +
+
+
+
+
+
+ + +
+
+ +
+
+
+
+

Exercise sheet

Some general remarks about the exercises:

+
    +
  • For your convenience functions from the lecture are included below. +Feel free to reuse them without copying to the exercise solution box.
    • +
  • For each part of the exercise a solution box has been added, but you + may insert additional boxes. Do not hesitate to add Markdown boxes for +textual or LaTeX answers (via Cell > Cell Type > Markdown). But make sure to replace any part that says YOUR CODE HERE or YOUR ANSWER HERE and remove the raise NotImplementedError().
    • +
  • Please make your code readable by humans (and not just by the Python + interpreter): choose informative function and variable names and use +consistent formatting. Feel free to check the PEP 8 Style Guide for Python for the widely adopted coding conventions or this guide for explanation.
    • +
  • Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting Kernel > Restart & Run All in the jupyter menu.
    • +
  • For some exercises test cases have been provided in a separate cell in the form of assert statements. When run, a successful test will give no output, whereas a failed test will display an error message.
    • +
  • Each sheet has 100 points worth of exercises. Note that only the +grades of sheets number 2, 4, 6, 8 count towards the course examination. + Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are +just for feedback.
    • +
+

Please fill in your name here:

+ +
+
+ +
+
+
In [1]:
+
+ 
NAME = "Kees van Kempen"
+NAMES_OF_COLLABORATORS = "Bart Steeman" # We did not finish it together.
+
+ +
+
+
+ +
+
+
+
+
+
+ +
+
+ +
+
+
+
+

Exercise sheet 6

+

Code from the lecture

+ +
+
+ +
+
+
In [2]:
+
+ 
import numpy as np
+import matplotlib.pylab as plt
+
+rng = np.random.default_rng()
+%matplotlib inline
+
+def aligned_init_config(width):
+    '''Produce an all +1 configuration.'''
+    return np.ones((width,width),dtype=int)
+
+def plot_ising(config,ax,title):
+    '''Plot the configuration.'''
+    ax.matshow(config, vmin=-1, vmax=1, cmap=plt.cm.binary)
+    ax.title.set_text(title)
+    ax.set_yticklabels([])
+    ax.set_xticklabels([])
+    ax.set_yticks([])
+    ax.set_xticks([])
+
+from collections import deque
+
+def neighboring_sites(s,w):
+    '''Return the coordinates of the 4 sites adjacent to s on an w*w lattice.'''
+    return [((s[0]+1)%w,s[1]),((s[0]-1)%w,s[1]),(s[0],(s[1]+1)%w),(s[0],(s[1]-1)%w)]
+
+def cluster_flip(state,seed,p_add):
+    '''Perform a single Wolff cluster move with specified seed on the state with parameter p_add.'''
+    w = len(state)
+    spin = state[seed]
+    state[seed] = -spin  
+    cluster_size = 1
+    unvisited = deque([seed])   # use a deque to efficiently track the unvisited cluster sites
+    while unvisited:   # while unvisited sites remain
+        site = unvisited.pop()  # take one and remove from the unvisited list
+        for nbr in neighboring_sites(site,w):
+            if state[nbr] == spin and rng.uniform() < p_add:
+                state[nbr] = -spin
+                unvisited.appendleft(nbr)
+                cluster_size += 1
+    return cluster_size
+
+def wolff_cluster_move(state,p_add):
+    '''Perform a single Wolff cluster move on the state with addition probability p_add.'''
+    seed = tuple(rng.integers(0,len(state),2))
+    return cluster_flip(state,seed,p_add)
+
+def compute_magnetization(config):
+    '''Compute the magnetization M(s) of the state config.'''
+    return np.sum(config)
+
+def run_ising_wolff_mcmc(state,p_add,n):
+    '''Run n Wolff moves on state and return total number of spins flipped.'''
+    total = 0
+    for _ in range(n):
+        total += wolff_cluster_move(state,p_add)
+    return total
+
+def sample_autocovariance(x,tmax):
+    '''Compute the autocorrelation of the time series x for t = 0,1,...,tmax-1.'''
+    x_shifted = x - np.mean(x)
+    return np.array([np.dot(x_shifted[:len(x)-t],x_shifted[t:])/len(x) 
+                     for t in range(tmax)])
+
+def find_correlation_time(autocov):
+    '''Return the index of the first entry that is smaller than 
+    autocov[0]/e or the length of autocov if none are smaller.'''
+    smaller = np.where(autocov < np.exp(-1)*autocov[0])[0]
+    return smaller[0] if len(smaller) > 0 else len(autocov)
+
+ +
+
+
+ +
+
+
+
+
+

MCMC simulation of the XY model

(100 Points)

+

Goal of this exercise: Practice implementing MCMC simulation + of the XY spin model using an appropriate cluster algorithm and +analyzing the numerical effectiveness.

+

The XY model is a relative of the Ising model in which the discrete ±1 spin at each lattice site is replaced by a continuous 2-dimensional spin on the unit circle

+S1={(x,y)∈R2:x2+y2=1}.

To be precise, we consider a w×w lattice with periodic boundary conditions and a XY configuration s=(s1,…,sN)∈Γ=S1N, N=w2, with Hamiltonian that is very similar to the Ising model,

+HXY(s)=−J∑i∼jsi⋅sj.

Here, as in the Ising model, the sum runs over nearest neighbor pairs i and j and si⋅sj is the usual Euclidean inner product of the vectors si,sj∈S1. We will only consider the ferromagnetic XY model and set J=1 in the remainder of the exercise. Note that nowhere in the definition the x- and y-components + of the spins are related to the two directions of the lattice (one +could also have studied the XY model on a one-dimensional or +three-dimensional lattice and the spins would still have two +components). As usual we are interested in sampling configurations s∈Γ with distribution π(s) given by the Boltzmann distribution

+π(s)=1ZXYe−βHXY(s),β=1/T.

The XY model admits a (local) cluster algorithm that is very similar to the Wolff algorithm of the Ising model. It amounts to the following recipe:

+
    +
  1. Sample a uniform seed site iseed in 1,…,N and an independent uniform unit vector n^∈S1.
    2. +
  2. Grow a cluster C starting from the seed iseed consisting only of sites j whose spin sj is "aligned" with the seed, in the sense that sj⋅n^ has the same sign as siseed⋅n^, or (sj⋅n^)(siseed⋅n^)>0. + Like in the Ising model this is done iteratively by examining the +neighbors of sites that are already in the cluster, and adding those +that are aligned with appropriate probability. The difference with the +Ising model is that this probability depends on the spins si and sj that are linked (meaning that sj is an aligned neighbor of si) via the formula +padd(si,sj)=1−exp⁡(−2β(si⋅n^)(sj⋅n^)).
    3. +
  3. Once the cluster C is constructed, all of its spins are "flipped" in the sense that they are reflected in the plane perpendicular to n^, i.e. sj→sj−2(sj⋅n^)n^.
    4. +
+

(a) Verify by a calculation, to be included using markdown and LaTeX below, that the probabilities padd(si,sj) are the appropriate ones to ensure detailed balance for the Boltzmann distribution π(s). Hint: Follow the same reasoning as for the Ising model. Compare the probabilities involved in producing the cluster C in state s and state s′. Why do the probabilities only differ at the boundary edges in the cluster C? (25 pts)

+ +
+
+ +
+
+
+
Student's answer + Score: 25.0 / 25.0 (Top) +
+
+
+

The condition for detailed balance is P(s→s′)P(s′→s)=π(s′)π(s). +Let's work to that expression by following the idea of the proof for the Ising model proof.

+

Starting from the left-hand side, consider the same cluster from the same seed and unit vector n^ but between different states s and s′. Flipping between s and s′ is done by the same cluster C, + but from opposite sign. The probability of this happening is given by +the boundaries of the clusters, namely that they are not added to the +cluster, with probability 1−padd per boundary site. This gives a ratio P(s→s′)P(s′→s)=∏<i,j>+(1−padd(si,sj))∏<i,j>−(1−padd(si,sj)) for pairs <i,j>± with si∈C and sj neighbouring si just outside the cluster with ±(si⋅n^)(sj⋅n^)>0. Substituting the formula for padd in, we find P(s→s′)P(s′→s)=∏<i,j>+exp⁡(−2β(si⋅n^)(sj⋅n^))∏<i,j>−exp⁡(2β(si⋅n^)(sj⋅n^))=∏<i,j>exp⁡(−2β(si⋅n^)(sj⋅n^))=exp⁡(−2β∑<i,j>(si⋅n^)(sj⋅n^)) with <i,j> describing both types of pairs.

+

For the right-hand side, π(s′)π(s)=eβ[HXY(s)−HXY(s′)] just from the definition of π. +This difference in energy we can write as +HXY(s)−HXY(s′)=∑<i,j>[si⋅sj−si′⋅sj]=2∑<i,j>(si⋅n^)(sj⋅n^) +by summing over all pairs <i,j> + as above, as pairs not across the boundary do not give rise to an +energy change. The last step is due to the flips from both signs to be +accounted for. +This shows P(s→s′)P(s′→s)=π(s′)π(s), thus that detailed balance holds.

+ +
+
+ +
+ +
+
+
+
+

(b) In order to implement the cluster update described above, we take the state to be described by a Numpy array of dimension (w,w,2), for which we have already provided a function xy_aligned_init_config to generate an all-aligned initial state. Write the function xy_cluster_flip, that grows and flips a cluster starting from the given seed site and n^ and returns the cluster size, and xy_cluster_move, that performs the previous function to a random seed and direction n^ and also returns the cluster size. (20 pts)

+ +
+
+ +
+
+
In [3]:
+
Student's answer(Top)
+
+
+ 
def xy_aligned_init_config(width):
+    '''Return an array of dimension (width,width,2) representing aligned spins in x-direction.'''
+    return np.dstack((np.ones((width,width)),np.zeros((width,width))))
+
+def xy_cluster_flip(state,seed,nhat,beta):
+    '''Perform a cluster move with specified seed and vector nhat on the state at temperature beta.'''
+    w = len(state)
+    
+    # Let's flip along the way, starting with the seed.
+    state[seed] -= 2*state[seed]@nhat*nhat
+    # Use a double-ended queue to keep track of what elements to try to expand from.
+    unvisited = deque([seed])
+    cluster_size = 1
+    while unvisited:
+        site = unvisited.pop()
+        # As the site is already added to the queue, it has already been flipped
+        # so we need to take a minus into account.
+        s_i_n = -state[site]@nhat
+        for nbr in neighboring_sites(site,w):
+            s_j_n = state[nbr]@nhat
+            one_minus_p_add = np.exp(-2*beta*s_i_n*s_j_n)
+            if np.sign(s_i_n) == np.sign(s_j_n) and rng.uniform() > one_minus_p_add:
+                state[nbr] -= 2*s_j_n*nhat
+                unvisited.appendleft(nbr)
+                cluster_size += 1
+    return cluster_size
+
+def xy_cluster_move(state,beta):
+    '''Perform a single Wolff cluster move on the state with addition probability p_add.'''
+    seed = tuple(rng.integers(0,len(state),2))
+    phi = 2*np.pi*rng.uniform()
+    nhat = np.array([np.cos(phi), np.sin(phi)])
+    return xy_cluster_flip(state,seed,nhat,beta)
+
+ +
+
+ +
+
+ +
+
+
+
In [4]:
+
Grade cell: cell-2e1ebc198d10ea5b + + Score: 15.0 / 15.0 (Top) +
+
+
+ 
from nose.tools import assert_almost_equal
+assert 1 <= xy_cluster_flip(xy_aligned_init_config(4),(0,0),np.array([np.cos(0.5),np.sin(0.5)]),0.5) <= 16
+assert_almost_equal(np.mean([xy_cluster_flip(xy_aligned_init_config(3),(0,0),
+                                             np.array([np.cos(0.5),np.sin(0.5)]),0.3) 
+                             for _ in range(200)]),5.3,delta=0.7)
+assert_almost_equal(np.mean([xy_cluster_flip(xy_aligned_init_config(3),(1,2),
+                                             np.array([np.cos(0.2),np.sin(0.2)]),0.2) 
+                             for _ in range(200)]),4.3,delta=0.6)
+
+ +
+
+ +
+
+ +
+
+
+
In [5]:
+
Grade cell: cell-9445d8ff6dc2a16f + + Score: 5.0 / 5.0 (Top) +
+
+
+ 
assert 1 <= xy_cluster_move(xy_aligned_init_config(4),0.5) <= 16
+assert_almost_equal(np.mean([xy_cluster_move(xy_aligned_init_config(3),0.3) 
+                             for _ in range(200)]),3.6,delta=0.75)
+assert_almost_equal(np.mean([xy_cluster_move(xy_aligned_init_config(3),0.9) 
+                             for _ in range(200)]),6.3,delta=0.75)
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(c) Estimate and plot the average cluster size in equilibrium for a 25x25 lattice (w=25) for the range of temperatures T=0.5,0.6,…,1.5. + It is not necessary first to estimate the equilibration time: you may +start in a fully aligned state and use 400 moves for equilibration and +1000 for estimating the average cluster size. It is not necessary to +estimate errors for this average. Store your averages in the data set "cluster-size" (an array of size 11) in the HDF5-file xy_data.hdf5, just like you did in Exercise sheet 5. Then read the data from file and produce a plot. (20 pts)

+ +
+
+ +
+
+
In [6]:
+
Student's answer(Top)
+
+
+ 
temperatures = np.linspace(0.5,1.5,11)
+width = 25
+equilibration_moves = 400
+measurement_moves = 1000
+
+import h5py
+import time
+
+t_0 = time.time()
+def toc():
+    """Return how many seconds the function was last called."""
+    
+    if not 't_0' in globals():
+        return -np.inf
+    
+    global t_0
+    dt = time.time() - t_0
+    t_0 = time.time()
+    return dt
+
+with h5py.File("xy_data.hdf5", "a") as f:
+    if not "cluster-size" in f:
+        state = xy_aligned_init_config(width)
+        cluster_sizes = np.zeros(len(temperatures))
+        
+        for idx, T in enumerate(temperatures):
+            toc()
+            print("T = {:.1f}".format(T))
+            beta = 1/T
+            
+            # Equilibrate
+            for _ in range(equilibration_moves):
+                xy_cluster_move(state,beta)
+            print("\tequilibrating took {:2.3f} seconds".format(toc()))
+            
+            # Measure
+            for _ in range(measurement_moves):
+                cluster_sizes[idx] += xy_cluster_move(state,beta)
+            cluster_sizes[idx] /= measurement_moves
+            print("\tmeasuring     took {:2.3f} seconds".format(toc()))
+            print()
+            
+        f.create_dataset("cluster-size",data=cluster_sizes)
+
+ +
+
+ +
+
+ +
+
+
+
In [7]:
+
Grade cell: cell-69f351f86b801f63 + + Score: 10.0 / 10.0 (Top) +
+
+
+ 
with h5py.File('xy_data.hdf5','r') as f:
+    assert f["cluster-size"][()].shape == (11,)
+    assert_almost_equal(f["cluster-size"][4],225,delta=40)
+    assert_almost_equal(f["cluster-size"][10],8,delta=8)
+
+ +
+
+ +
+
+ +
+
+
+
In [8]:
+
Student's answer + Score: 10.0 / 10.0 (Top) +
+
+
+ 
# Plotting
+with h5py.File("xy_data.hdf5", "r") as f:
+    plt.figure()
+    # My girlfriend told me to use this marker.
+    plt.plot(temperatures, np.array(f["cluster-size"]), marker="4", markersize=20, linestyle="dashed", linewidth=.5)
+    plt.xlabel("temperature $T$")
+    plt.ylabel("average cluster size")
+    plt.title("Average cluster size for the $w \\times w$ XY model \n"
+              "(w = {}, #equilibrations = {}, #measurements = {})"
+              .format(width, equilibration_moves, measurement_moves))
+    plt.show()
+
+ +
+
+
+
+ Comments:

Fancy marker! +Good labeling

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

(d) Make an MCMC estimate (and plot!) of the mean square magnetization per spin E[m2(s)] for the same set of temperatures, where

+m2(s)=(1N∑i=1Nsi)⋅(1N∑i=1Nsi).

To choose the equilibration time and time between +measurement, use the average cluster size from (c) to estimate how many +moves correspond to 1 sweep, i.e. roughly N=w2 updates to spins. Then use 100 equilibration sweeps and 200 measurements of m2(s), with 2 sweeps between each measurement. Store the measured values of m2(s) in the data set "square-magn" of dimension (11,200) in xy_data.hdf5. +Then read the data and plot estimates for E[m2(s)] + including errors (based on batching or jackknife). If the errors are +too small to see, you may multiply them by some number and indicate this + in the title of the plot. (20 pts)

+ +
+
+ +
+
+
In [9]:
+
Student's answer(Top)
+
+
+ 
measurements = 200
+equil_sweeps = 100
+measure_sweeps = 2
+
+with h5py.File("xy_data.hdf5", "a") as f:
+    if not "square-magn" in f:
+        state = xy_aligned_init_config(width)
+        square_magn = np.zeros((len(temperatures), measurements))
+        
+        for idx, T in enumerate(temperatures):
+            toc()
+            print("T = {:.1f}".format(T))
+            beta = 1/T
+            # We let the sweep size depend on the width and the average cluster size
+            # from before:
+            sweep_size = np.ceil(width**2/f["cluster-size"][idx]).astype(int)
+            
+            # Equilibrate
+            for j in range(equil_sweeps*sweep_size):
+                xy_cluster_move(state,beta)
+                print("\tequilibrating... [{}/{}]".format(j, equil_sweeps*sweep_size), end='\r')
+            print("\tequilibrating took {:2.3f} seconds".format(toc()))
+            
+            # Measure
+            for j in range(measurements):
+                xy_cluster_move(state,beta)
+                square_magn[idx][j] = np.linalg.norm(np.mean(state, axis=(0,1)))**2
+                print("\tmeasuring... [{}/{}]".format(j, measurements), end='\r')
+                
+                # Sweeps between measurements
+                for _ in range(measure_sweeps*sweep_size):
+                    xy_cluster_move(state,beta)
+            print("\tmeasuring     took {:2.3f} seconds".format(toc()))
+            print()
+            
+        f.create_dataset("square-magn",data=square_magn)
+
+ +
+
+ +
+
+ +
+
+
+
In [10]:
+
Grade cell: cell-5d3a2f49615613b2 + + Score: 10.0 / 10.0 (Top) +
+
+
+ 
with h5py.File('xy_data.hdf5','r') as f:
+    assert f["square-magn"][()].shape == (11, 200)
+    assert_almost_equal(np.mean(f["square-magn"][4]),0.456,delta=0.02)
+    assert_almost_equal(np.mean(f["square-magn"][9]),0.023,delta=0.01)
+
+ +
+
+ +
+
+ +
+
+
+
In [11]:
+
Student's answer + Score: 10.0 / 10.0 (Top) +
+
+
+ 
# Plotting
+
+# As copied from last week:
+def jackknife_batch_estimate(data,observable,k):
+    '''Divide data into k batches and apply the function observable to each 
+    collection of all but one batches. Returns the mean and corrected 
+    standard error.'''
+    batches = np.reshape(data,(k,-1))
+    values = [observable(np.delete(batches,i,0).flatten()) for i in range(k)]
+    return np.mean(values), np.sqrt(k-1)*np.std(values)
+
+with h5py.File('xy_data.hdf5','r') as f:
+    k = 100
+    mean_sq, std_sq = np.array([jackknife_batch_estimate(data, np.mean, k) for data in f["square-magn"]]).T
+    plt.figure()
+    # I do plot a line to visualize the trend.
+    plt.errorbar(temperatures, mean_sq, yerr=10*std_sq,
+                 fmt='.--', markersize=10, capsize=4)
+    plt.xlabel("temperature $T$")
+    plt.ylabel("mean square magnetization $\\overline{m^2}$")
+    plt.title("Mean square magnetization for the $w \\times w$ XY model \n"
+              "(w = {}, #measurements = {}, $10\\sigma$ errors)"
+              .format(width, measurements))
+    plt.show()
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

(e) Produce a single equilibrated state for each temperature and store them in the data set "states" of dimension (11,25,25,2) in xy_data.hdf5. Then read them and produce a table of plots using the provided function plot_xy, which shows colors based on the angle of the spin, each with a title to indicate the temperature. Can you observe the Kosterlitz–Thouless transition of the XY model? (15 pts)

+ +
+
+ +
+
+
In [12]:
+
Student's answer(Top)
+
+
+ 
width = 25
+state = xy_aligned_init_config(width)
+equil_sweeps = 200
+
+with h5py.File("xy_data.hdf5", "a") as f:
+    if not "states" in f:
+        states = np.zeros(temperatures.shape + state.shape)
+        for idx, T in enumerate(temperatures):
+            toc()
+            print("T = {:.1f}".format(T))
+            beta = 1/T
+            sweep_size = np.ceil(width**2/f["cluster-size"][idx]).astype(int)
+            
+            # Equilibrate
+            for j in range(equil_sweeps*sweep_size):
+                xy_cluster_move(state,beta)
+                print("\tequilibrating... [{}/{}]".format(j, equil_sweeps*sweep_size), end='\r')
+            print("\tequilibrating took {:2.3f} seconds".format(toc()))
+            states[idx] = state.copy()
+            print()
+        
+        f.create_dataset("states", data=states)
+
+ +
+
+ +
+
+ +
+
+
+
In [13]:
+
Grade cell: cell-8968d7e317e522f7 + + Score: 5.0 / 5.0 (Top) +
+
+
+ 
with h5py.File('xy_data.hdf5','r') as f:
+    assert f["states"][()].shape == (11, 25, 25, 2)
+
+ +
+
+ +
+
+ +
+
+
+
In [14]:
+
Student's answer + Score: 10.0 / 10.0 (Top) +
+
+
+ 
def plot_xy(state,ax,title=""):
+    '''Plot the XY configuration given by state. Takes an Axes object ax from 
+    matplotlib to draw to, and adds the specified title.'''
+    angles = np.arctan2(*np.transpose(state,axes=(2,0,1)))
+    ax.matshow(angles, vmin=-np.pi, vmax=np.pi, cmap=plt.cm.hsv)
+    ax.title.set_text(title)
+    ax.set_yticklabels([])
+    ax.set_xticklabels([])
+    ax.set_yticks([])
+    ax.set_xticks([])
+
+# Make a table of plots
+with h5py.File("xy_data.hdf5", "a") as f:
+    fig, axs = plt.subplots(3, 4, figsize=(10, 10))
+    axs = axs.ravel()
+    for idx, T in enumerate(temperatures):
+        plot_xy(f["states"][idx], axs[idx],"T = {:.1f}".format(T))
+    
+    # Hide the axes of unused plots.
+    [ax.axis('off') for ax in axs[len(temperatures):]]
+    plt.tight_layout()
+    plt.show()
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

The Kosterlitz-Thouless transition is vaguely visible. For T=1.0 and lower, vortices and anti-vortices can be recognized. Above this temperature, this structure vanishes.

+ +
+
+ +
+
+ + +
+
+
+
+
+ + + \ No newline at end of file