Exercise sheet

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

NAME = ""
NAMES_OF_COLLABORATORS = ""

Exercise sheet 6

Code from the lecture

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)

6.1. 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 $\pm 1$ spin at each lattice site is replaced by a continuous $2$-dimensional spin on the unit circle

\[ S_1 = \{(x,y)\in\mathbb{R}^2: x^2+y^2=1\}. \]

To be precise, we consider a $w\times w$ lattice with periodic boundary conditions and a XY configuration $s = (s_1,\ldots,s_N) \in \Gamma = S_1^N$, $N=w^2$, with Hamiltonian that is very similar to the Ising model,

\[ H_{XY}(s) = - J \sum_{i\sim j} s_i \cdot s_j. \]

Here, as in the Ising model, the sum runs over nearest neighbor pairs $i$ and $j$ and $s_i \cdot s_j$ is the usual Euclidean inner product of the vectors $s_i,s_j \in S_1$. 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\in \Gamma$ with distribution $\pi(s)$ given by the Boltzmann distribution

\[ \pi(s) = \frac{1}{Z_{XY}} e^{-\beta H_{XY}(s)}, \qquad \beta = 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:

Sample a uniform seed site $i_{\text{seed}}$ in $1,\ldots,N$ and an independent uniform unit vector $\hat{n} \in S_1$.
Grow a cluster $C$ starting from the seed $i_{\text{seed}}$ consisting only of sites $j$ whose spin $s_j$ is “aligned” with the seed, in the sense that $s_j\cdot\hat{n}$ has the same sign as $s_{i_{\text{seed}}}\cdot \hat{n}$, or $(s_j\cdot\hat{n})(s_{i_{\text{seed}}}\cdot \hat{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 $s_i$ and $s_j$ that are linked (meaning that $s_j$ is an aligned neighbor of $s_i$) via the formula $$ p_{\text{add}}(s_i,s_j) = 1 - \exp\big( -2\beta\,(s_i\cdot\hat{n})(s_j\cdot\hat{n})\big).$$
Once the cluster $C$ is constructed, all of its spins are “flipped” in the sense that they are reflected in the plane perpendicular to $\hat{n}$, i.e. $s_j \to s_j - 2(s_j\cdot\hat{n})\hat{n}$.

(a) Verify by a calculation, to be included using markdown and LaTeX below, that the probabilities $p_{\text{add}}(s_i,s_j)$ are the appropriate ones to ensure detailed balance for the Boltzmann distribution $\pi(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)

YOUR ANSWER HERE

(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 $\hat{n}$ and returns the cluster size, and xy_cluster_move, that performs the previous function to a random seed and direction $\hat{n}$ and also returns the cluster size. (20 pts)

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)
    # YOUR CODE HERE
    raise NotImplementedError()
    return cluster_size

def xy_cluster_move(state,beta):
    '''Perform a single Wolff cluster move on the state with addition probability p_add.'''
    # YOUR CODE HERE
    raise NotImplementedError()

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)

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,\ldots,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)

temperatures = np.linspace(0.5,1.5,11)
width = 25
equilibration_moves = 400
measurement_moves = 1000

# YOUR CODE HERE
raise NotImplementedError()

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)

# Plotting
# YOUR CODE HERE
raise NotImplementedError()

(d) Make an MCMC estimate (and plot!) of the mean square magnetization per spin $\mathbb{E}[ m^2(s) ]$ for the same set of temperatures, where

\[ m^2(s) = \left(\frac{1}{N}\sum_{i=1}^N s_i\right)\cdot\left(\frac{1}{N}\sum_{i=1}^N s_i\right). \]

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 = w^2$ updates to spins. Then use 100 equilibration sweeps and 200 measurements of $m^2(s)$, with 2 sweeps between each measurement. Store the measured values of $m^2(s)$ in the data set "square-magn" of dimension $(11,200)$ in xy_data.hdf5. Then read the data and plot estimates for $\mathbb{E}[ m^2(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)

measurements = 200
equil_sweeps = 100
measure_sweeps = 2

# YOUR CODE HERE
raise NotImplementedError()

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)

# Plotting
# YOUR CODE HERE
raise NotImplementedError()

(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)

width = 25
state = xy_aligned_init_config(width)
equil_sweeps = 200

# YOUR CODE HERE
raise NotImplementedError()

with h5py.File('xy_data.hdf5','r') as f:
    assert f["states"][()].shape == (11, 25, 25, 2)

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
# YOUR CODE HERE
raise NotImplementedError()

Exercise sheet

Contents

6. Exercise sheet#

6.1. MCMC simulation of the XY model#