Ising model simulation using metropolis algorithm
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.