Jing Qiang Goh

We can study the phase transition of a spin-glass based on the well-known Ising model and Metropolis algorithm. The flipping of a spin depends on its neighboring spin configurations, and periodic boundary conditions are imposed such that we can study an 'infinite' system using a finite sized spin system. Periodic boundary conditions are commonly found in the other physical problems as well.

Assuming that we have a L times L square lattice system with spin[i, j] indicating a spin positioned at i-th and j-th site of the system. Let us mark the position of spin along the i-dimension as 0, 1, 2, ..., L - 1, and the same applied to the j-dimension. A straight forward implementation of periodic boundary conditions is using an if-else loop, such as (Python example):

# di1 and di2, indices of spin to be manipulated
# L: The size of square lattice edge, spin: A 2D numpy array
if i == 0:
    di1 = L - 1
    di2 = 1
elif i == L - 1:
    di1 = L - 2
    di2 = 0
else:
    di1 = i - 1
    di2 = i + 1
 
# spin[di1, j] will be the previous neighbor of spin[i, j]
# spin[di2, j] will be the next neighbor of spin[i, j]

Similarly, we can implement another if-else loop to handle the index j.

If you refer to articles describing the Monte Carlo simulation for the phase transition of spins, you are likely to learn an elegant solution for the implementation of periodic boundary conditions - using the mathematical operator Modulo.

An example that illustrates the implementation of periodic boundary conditions using the Modulo operator:

# L: The size of a square lattice edge, spin: A 2D numpy array
# The previous neighbor of spin[i, j], along the i-dimension
spin[(i - 1) % L, j]
 
# The next neighbor of spin[i, j], along the i-dimension
spin[(i + 1) % L, j]

Why do we consider the Modulo operator? You may have noticed that a Modulo operator can allow us to consider periodic boundary condition for the nearby spins further away from the nearest neighbor of spin[i, j]. Considering a system with L=10, spin[(i + 7) % L, j] can evaluate the 7th-spin away from spin[i, j] such that the 7-th nearest neighbor to the right of i=7 will be i=4. We cannot easily implement the periodic boundary conditions for neighbors further away by relying on simple conditional, if-else statements.

It could be a good idea to avoid the Modulo operation which needs more computational time when L is a huge number. Perhaps we can initialize another array to mark the indices of of lattice site such that we can refer this array for implementing the periodic boundary conditions.

import numpy as np # Use numpy array
 
L = 4
idx = list(range(L))
idx.extend(list(range(L)))
idx.extend(list(range(L)))
# Convert the list to numpy array
idx = np.array(idx, dtype=int) 
 
# The initialized array idx has the sequence of
# 0, 1, 2 .. L - 1, 0, 1, 2, ..., L - 1, 0, 1, 2, ..., L - 1
# idx with L = 4: [0 1 2 3 0 1 2 3 0 1 2 3]
 
# The previous kth-neighbor of spin[i, j]
spin[idx[i + L - k], j]
# The next kth-neighbor of spin[i, j]
spin[idx[i + L + k], j]

The above example assumes that we only need to evaluate the farthest neighbor of spin[i, j] with the size of L. If we only consider the nearest neighbor of a single spin, an array of [(L - 1), 0, 1, ..., (L - 1), 0] is sufficient for this purpose.

I would like to have a quick idea on the elapsed times required for these different periodic boundary schemes. In the following examples, let us consider the flipping of every spin within a square lattice by applying the Metropolis algorithm. 100 Monte Carlo sweeps are performed with each sweep having all the spins flipped, one by one. The nearest neighbor interaction is assumed in these examples. Below has a snippet of code for the example of using the Modulo method:

start = time.time()
# Initialize the system
L = args.L
print(L)
spin = np.ones((L, L))  # 2D square lattice, spin up
T = 300  # 300 K, for temperature
 
# Method 2, using modulus method
random.seed(10)
for k in range(args.num):
    for i in range(L):
        for j in range(L):
            # eflip, the change in the energy of system if we flip the
            # spin[i, j]. eflip depends on the configuration of 4 neighboring
            # spins. For instance, with reference to spin[i, j], we should evaluate
            # eflip based on spin[i+1, j], spin[i-1, j], spin[i, j+1], spin[i, j-1]
            eflip = 2*spin[i, j]*(
                spin[((i - 1) % L), j] +  # -1 in i-dimension
                spin[((i + 1) % L), j] +  # +1 in i-dimension
                spin[i, ((j - 1) % L)] +  # -1 in j-dimension
                spin[i, ((j + 1) % L)]    # +1 in j-dimension
            )
            # Metropolis algorithm
            if eflip = 0.0:
                spin[i, j] = -1.0*spin[i, j]
            else:
                if (random.random() = math.exp(-1.0*eflip/T)):
                    spin[i, j] = -1.0*spin[i, j]
 
end = time.time()
print(spin)
print(end - start)

The following table summarizes the rough estimation of elapsed time from the test examples (since I am not performing a careful profiling analysis). The absolute values are not always consistent, nevertheless they can provide us the relative standings. You may download the Python codes, modify and test run the codes on your own computing platform.

It turns out that if we implement the periodic boundary schemes as shown by the snippet of codes (Modulo operator) mentioned, it can complete the calculations in the shortest time as we scale up the size of L. If I try to repeat the same test with Matlab codes instead, we arrive at the following results:

Based on these examples, we realize that Python codes require more time in accessing the numpy array element. Surprisingly, the Python examples with Modulo operation can perform calculations faster than the examples implementing if-else conditions. Referring to the Matlab examples, a simple if-else condition is sufficient to provide us the fastest calculation as we only consider the nearest-neighbor spin interaction. Modulo operation in Matlab requires significantly more time than the other schemes.

Note that these tests DO NOT aim to compare the performance of Python to Matlab. Edit: I believe it is the Just-in-time compilation routine can explain that the Matlab codes compute faster than Python at this level of programming interface. One could also argue that the Python codes are not optimized for the simulation purpose and this results in the example with Modulo operation as the "best".

Other than demonstrating the effective ways of implementing boundary condiitons, this post also demonstrates that due to the different backend features of the programming platforms, we can arrive at different conclusions even though the semantic of programming appears the same.