Ising model simulation using metropolis algorithm I am new to this community; I have tried my best to respect the policy of the community. I have written the Monte Carlo metropolis algorithm for the ising model. I want to optimize the code. I have tried my best. I want to optimize it further. The following is the code:

(I have used tricks like finding exponential only once, careful generation of random number, etc.)

import numpy as np
import time
import random

def monteCarlo(N,state,Energy,Mag,Beta,sweeps):
    if sweeps &gt; 10000:
        print("Warning:Number of sweeps exceeded 10000\n\tsetting number of sweeps to 10000")
        sweeps = 10000
    start_time = time.time()
    expBeta = np.exp(-Beta*np.arange(0,9))
    E = np.zeros(sweeps)
    M = np.zeros(sweeps)
    for t in range(sweeps):
        for tt in range(N*N):
            a = random.randint(0, N-1)
            b = random.randint(0, N-1)
            s = state[a,b]
            delta_E = 2*s*(state[(a+1)%N,b] + state[a,(b+1)%N] + state[(a-1)%N,b] + state[a,(b-1)%N])
            if delta_E &lt; 0:
                s *= -1
                Energy += delta_E
                Mag += 2*s
            elif random.random() &lt; expBeta[delta_E]:
                s *= -1
                Energy += delta_E
                Mag += 2*s
            state[a, b] = s
        E[t] = Energy
        M[t] = Mag
    print("%d monte carlo sweeps completed in %d seconds" %(sweeps,time.time()-start_time))
    return E,M #returning list of Energy and Magnetization set

#####lattice config#####
"""N is lattice size
nt is number of Temperature points
sweeps are number of mc steps per spin"""
print("Starting Ising Model Simulation")
N = int(input("Enter lattice size : "))
startTime = time.time()
nt = 10
N2 = N*N
sweeps = 10000 #mc steps per spin
"""we will plot the following wrt temperature, T"""
T = np.linspace(2, 3, nt) #you can set temperature range
"""preparing lattice with all spins up"""
state = np.ones((N,N),dtype="int")
Energy = -N2
Mag = N2
#temperature loop
#for k in tqdm_gui(range(nt)):
for k in range(nt):
    temp = T[k]
    Beta=1/temp
    print("____________________________________\nTemperature is %0.2f, time is %d" %(temp,time.time()-startTime))
    E, M = monteCarlo(N,state,Energy,Mag,Beta,sweeps) #list of Energy and Magnetization
    Energy = E[-1]
    Mag = M[-1]

#further code is finding magnetization, autocorrelation, specific heat, autocorrelation, etc.
```
I cannot yet comment due to my reputation, so I'll write this as an answer to your comment on Andrew's answer, and delete if someone comments this information or Andrew updates his answer.

Saying that numpy's random is not a good random number generator does not seem right to me.
From numpy's reference : 


  By default, Generator uses bits provided by PCG64 which has better statistical properties than the legacy MT19937 used in RandomState.


So to me it seems that :


numpy uses the PCG64 random generator which, according to numpy, has better statistical properties than legacy MT19937
numpy used to use MT19937
you can still chose to use the MT19937 random number generator 


from numpy.random import Generator, MT19937
rg = Generator(MT19937(12345))
rg.random()


Maybe there is something I'm missing, and maybe it is in this part of your comment 


  It produces 53-bit precision floats and has a period of 2**19937-1


If so, I'd be interested to know how numpy's random would still be flawed for a Monte-Carlo analysis.