Ising Model simulation in Numpy

Question

I've written this code to simulate the Ising Model. The idea is a system represented as a matrix with each element having value 1 or -1. The input is an initial matrix of 1 and -1; it is then being subjected to this algorithm:

1. Pick a random element of the matrix
2. calculate the energy change of the system if that element is flipped (turn 1 to -1 vice versa)
3. if the flipping reduces total energy (after/before < 0): do it, else: flip it with probability

$$ P = \exp(-\text{constant}*\text{energy difference}) $$ $0ドル \le P \le 1$$

(I added an extra term in the code to modify its dynamics.)

Energy is calculated by observing the spin of an element with its neighbor, so I use two matrices for this calculation:

1. the main matrix (the one subjected to the function),
2. an adjacency matrix (representing connection between each elements from the main matrix: 1 if connected, else 0).

The idea is to multiply the row from the adjacency matrix to all elements of the main matrix.

The code works but since I have to apply the function for thousands of times and it uses two (rather large) matrices, the calculation takes a lot of time. Perhaps ideas for speed-ups (or just for writing a better code in general) and possibility to use numba which currently gives me complicated error messages.

import numpy as np
n = 64
n2 = n**2
main_matrix = 2*np.random.randint(2, size=(n, n))-1
adj_matrix = np.random.randint(2, size=(n2, n2)) # Assume random network is used
def mcmove(beta):
 a = np.random.randint(0, n)
 b = np.random.randint(0, n)
 s = main_matrix[a, b]
 energy_difference = np.sum(np.multiply(adj_matrix[a*n + b], main_matrix.reshape(1,-1)[0]))/np.sum(adj_matrix[a*n + b])
 cost = 2*s*(energy_difference - s*abs(np.sum(main_matrix))/n2)
 if cost < 0:
 s *= -1
 elif np.random.rand() < np.exp(-cost*beta):
 s *= -1
 main_matrix[a, b] = s
 return main_matrix
# How this function is used
total_spin = []
for i in range(1000):
 for i in range(n2):
 mcmove(0.001*i)
 total_spin.append(np.sum(main_matrix)) # sum all spins
total_spin

Answer 1

Stop using np.random.randint and family; they're deprecated in favour of the new generators.

You can assign a and b in one go by generating an integer array of length 2.

Add PEP484 type hints.

Don't np.multiply; in this case use np.dot since effectively you are calculating a dot product.

Don't reshape; your reshape is equivalent to a flatten().

Use np.abs instead of abs.

Don't elif; this is equivalent to an or sub-expression.

Don't return main_matrix; that's ignored and you're mutating in-place.

Use a different letter for your inner for i iteration, probably j.

Rather than 0.001 * i, find this value directly by iterating over a call to np.linspace.

Rather than total_spin.append, pre-allocate a properly typed, empty numpy array.

Other than the above there isn't much you can improve in terms of performance since the iterations are not independent: an evolved main_matrix depends on the value of main_matrix from the previous iteration. At that point you would drop down to a better language like C and calls to cblas or equivalent.

Suggested

import numpy as np
n = 64
n2 = n**2
rand = np.random.default_rng(seed=0)
main_matrix = rand.choice((-1, 1), size=(n, n))
adj_matrix = rand.integers(2, size=(n2, n2)) # Assume random network is used
def mcmove(beta: float) -> None:
 a, b = rand.integers(low=0, high=n, size=2)
 s = main_matrix[a, b]
 adj_row = adj_matrix[a*n + b, :]
 energy_difference = np.dot(adj_row, main_matrix.flatten()) / adj_row.sum()
 cost = 2 * s * (energy_difference - s * np.abs(main_matrix.sum()) / n2)
 if cost < 0 or rand.random() < np.exp(-cost * beta):
 main_matrix[a, b] = -s
def main() -> None:
 # How this function is used
 N = 1000
 total_spin = np.empty(shape=N, dtype=int)
 for i in range(N):
 for j in np.linspace(0, (n2 - 1)/N, n2):
 mcmove(j)
 total_spin[i] = main_matrix.sum() # sum all spins
 print(total_spin)
if __name__ == '__main__':
 main()