Heat-bath+algorithm+for+the+Ising+model

toc =Introduction=
 * Class Session 08: Heat-bath algorithm for the Ising model **

We discuss a simple algorithm for the Ising model at equilibrium with a bath at inverse temperature //β// (see section 5.2.2 in SMAC). The system is described by a set of //N// spins {//σ// 1 ,...,//σ N //} each of them taking the value +1 or –1. The energy of the system is proportional to the Boltzman weight of energy math E[\{\sigma\}]=-\displaystyle \sum_{\langle ij\rangle} \sigma_i \sigma_j math where the sum is over neighboring sites.

=The heat bath algorithm= During the lecture we have discussed the Metropolis algorithm for this system. Rather than flipping a spin at a random site, we now thermalize this spin with its local environment: the state of a spin after a flip does not depend on its state before the flip, but only on its molecular field //h// (defined by the sum of the spins of its neighbors). Show that, in the presence of a molecular field //h// at site //k//, the spins point up and down with probabilities: math \displaystyle \pi_h^+ = \frac 1{1+e^{-2\beta h}} \qquad \qquad \pi_h^- = \frac 1{1+e^{+2\beta h}} math Those define the transition probabilities at each step of the heat bath algorithm (for a uniformly chosen spin to become respectively up or down, irrespective of its initial state). Show that detailed balance is verified with respect to the Boltzmann-Gibbs weight at inverse temperature //β.// The following program implements this algorithm in two dimension on a square lattice of lenght //L// with periodic boundary conditions: code format="python" import random,math,pylab,numpy L = 100 Nstep = 100000 beta = 0.3 neighbors = [(0,-1),(0,1),(-1,0),(1,0)] s = 1 for i in range(L)] for j in range(L)] energy = -2. * L**2 toten = 0. for step in range(Nstep):    kx = random.randrange(L)    ky = random.randrange(L)    h = sum( s[(kx+dx)%L][(ky+dy)%L] for (dx,dy) in neighbors)    olds = s[kx][ky]    Upsilon = random.random    s[kx][ky] = -1    if Upsilon < 1./(1+math.exp(-2*beta*h)): s[kx][ky] = 1    if olds != s[kx][ky]: energy-=2*h*s[kx][ky]    toten += energy print toten / Nstep / L**2 print sum(sum(s[i][j] for i in range(L)) for j in range(L) )/L**2. pylab.matshow(numpy.array(s)) pylab.show [[code
 * 1) initialisation
 * 1) evolution

=Half order and coupling algorithms=



Let us first introduce the concept of half order for the configurations in the Ising model (see figure 1). We define that, for a site, an up spin is larger than a down spin, and a whole configuration {//σ// 1 ,...,//σ N //} of spins is larger than an other configuration {//σ//' 1 ,...,//σ//'// N //} if the spins on each sites //k// satisfy //σ k //≥//σ//'// k //. This defines an //half-order// among the configurations (//half//, because not all configurations can be compared)..
 * Show that the half order is preserved at each iteration of the heat bath algorithm.
 * In the forward coupling, show that it is sufficient to focus on the coupling of two particular configurations. Which ones?

The same remark applies for the backward coupling, and allows to draw a perfect sample. Here is an implementation: code format="python" import random,math,pylab,numpy L=10 beta=0.01 neighbors=[(0,-1),(0,1),(-1,0),(1,0)] splus0= 1 for i in range(L)] for j in range(L)] sminus0=[[ -1 for i in range(L)] for j in range(L)] splus,sminus=splus0[:],sminus0[:] changes=[] while splus!=sminus :   splus,sminus=splus0[:],sminus0[:]    kx=random.randrange(L)    ky=random.randrange(L)    Upsilon=random.random    changes.insert(0,[kx,ky,Upsilon])    for [kx,ky,Upsilon] in changes:        for s in (splus,sminus):            h=sum( s[(kx+dx)%L][(ky+dy)%L] for (dx,dy) in neighbors)            s[kx][ky]=-1            if Upsilon<1./(1+math.exp(-2*beta*h)): s[kx][ky]=1 print "here is a perfect sample", splus print "steps to find it", len(changes) pylab.matshow(numpy.array(splus)) pylab.show [[code What is the advantage of this method compared to the naive implementation of the backward coupling?
 * 1) initialisation
 * 1) evolution

url}?f=print| [Print this page]