Speeding Up Your Code (especially python)#

To speed up your python code, what you want to do is avoid having a for-loop over steps; instead you should work with many of the spins at once.

Heat Bath#

Before we do this, though, we need to make a minor change to the Monte Carlo move we are making (otherwise we accidentally get a pathological Markov Chain)

Our current move is

  • pick a random spin

  • consider flipping it

We are going to switch to another move - the heat bath move.

  • Pick a random spin (say spin \(i\))

  • Consider

    • making it spin-up with probability \(\frac{\exp[-\beta E(...\uparrow_i,...)]}{\exp[-\beta E(...\uparrow_i,...)]+\exp[-\beta E(...\downarrow_i,...)}\)

    • making it spin-down with probability \(\frac{\exp[-\beta E(...\downarrow_i,...)]}{\exp[-\beta E(...\uparrow_i,...)]+\exp[-\beta E(...\downarrow_i,...)}\)

Notice that this move doesn’t actually care what the current spin was! This means sometimes you are going to replace a configuration with the same configuration. (If you think for a minute you can convince yourself that if you’re considering doing this, you accept it with probability one).

When we do this we need to be a little careful with the acceptance probability. When we were flipping a spin, \(T(c\rightarrow c') = T(c' \rightarrow c) =1/N\) and the ratio of the \(\frac{T(c' \rightarrow c)}{T(c \rightarrow c')}=1\). Therefore it didn’t come into the acceptance probability.

In the new move, let’s assume we are moving from

  • our current configuration \(c = \{...\downarrow_i...\}\) to

  • our new configuration \(c' = \{...\uparrow_i...\}\)

Then we find that

  • \(T(c \rightarrow c')=(1/N) \frac{\exp[-\beta E(...\uparrow_i,...)]}{\exp[-\beta E(...\uparrow_i,...)]+\exp[-\beta E(...\downarrow_i,...)}\) (since this is the probability of choosing spin-up) and

  • \(T(c' \rightarrow c)=(1/N) \frac{\exp[-\beta E(...\downarrow_i,...)]}{\exp[-\beta E(...\uparrow_i,...)]+\exp[-\beta E(...\downarrow_i,...)}\).

Therefore, the ratio

\[\frac{T(c' \rightarrow c)}{T(c \rightarrow c')} = \frac{\exp[-\beta E(...\downarrow_i,...)]}{\exp[-\beta E(...\uparrow_i,...)]}\]

But our ratio of Boltzmann factors gives us $\(\frac{\pi(c')}{\pi(c)} = \frac{\pi(\{...\uparrow_i...\})}{\pi(\{...\downarrow_i...\})} = \frac{\exp[-\beta E(...\uparrow_i,...)]}{\exp[-\beta E(...\downarrow_i,...)]}\)$

Therefore the total acceptance ratio (the product of these two terms)

\[\frac{T(c' \rightarrow c)}{T(c \rightarrow c')} \frac{\pi(c')}{\pi(c)} =1 \]

Therefore, we accept 100% of the moves that we consider (although this isn’t as good as it sounds because some of that time we are accepting a move that doesn’t do anything).

If you do a little bit of algebra, you can rewrite this entire move so that the probability you flip your spin (i.e. the probability you compare to a random number) is \((1+1/\exp(-\beta \Delta E))^{-1}\) where \(\Delta E\) is the energy difference between your current configuration and the one with the spin-flipped. Replace your acceptance criterion with this new criterion. This should be the only change you need. Verify you get the same answers.

Many Spins At Once#

Let’s consider the \(7 \times 7\) lattice. Which spins can we flip simultaneously?

graph.png

We can consider flipping all of the red squares at once, because the red squares aren’t neighbors of any other red squares (even with the periodic boundary conditions) and so if we flip them all at one time, there isn’t going to be any conflicts. The same is true of the blue, yellow, and green colors (the yellow and green exist to deal with the boundary conditions and the odd-size lattice).

In python, we can access all the colors by doing

# This is for a L x L lattice
x,y=np.indices((L,L))
red=np.logical_and((x+y) % 2==0   , np.logical_and(x<L-1,y<L-1))
red[L-1,L-1]=True
blue=np.logical_and((x+y) % 2==1   , np.logical_and(x<L-1,y<L-1))
green=np.logical_and((x+y) % 2==0   , np.logical_or(x==L-1,y==L-1))
green[L-1,L-1]=False
yellow=np.logical_and((x+y) % 2==1   , np.logical_or(x==L-1,y==L-1))

For each color (say red as an example)

  • write a function that computes the \(\Delta E\) for all the red spins without using any for-loops. Two things you might find useful here are that spins[red] will access all the red spins. In addition, np.roll(spins,1,axis=0) will return an array where all the spins are rolled over to the right.

  • After calling this function, write one line of python that returns an array which says which red spins should be flipped.

  • Write one line of python which flips just those spins.

Just to give a sense, this is how long my version of the naive python code and the efficient python code take for 1000 sweeps:

Slow Code Wall time: 1min 20s
Efficient Code Wall time: 1.06 s