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:

Pick a random element of the matrix
calculate the energy change of the system if that element is flipped (turn 1 to -1 vice versa)
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:

the main matrix (the one subjected to the function),
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

Reinderien · Accepted Answer · 2022-03-15 15:25:51Z

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

Stack Exchange Network

Ising Model simulation in Numpy

1 Answer 1

Suggested

You must log in to answer this question.

Hot Network Questions

Ising Model simulation in Numpy

1 Answer 1

Suggested

You must log in to answer this question.

Related

Hot Network Questions