# Ising Model simulation in Numpy

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]

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


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.

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]

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