# 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 > 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 < 0:
s *= -1
Energy += delta_E
Mag += 2*s
elif random.random() < 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.
$$$$


## Naming variables

The usual rule of variable names in snake_case applies, i.e. energyFunctional would become energy_functional. Class names on the other hand should be written in CamelCase. I don't mind using single capital letters for matrices.

Why use a,b for discrete indices? I would use one of i,j,k,l,n,m,p,q,r.

Use descriptive names: e.g. energies instead of E.

## Merging conditions

if delta_E < 0:
s *= -1
Energy += delta_E
Mag += 2*s
elif random.random() < expBeta[delta_E]:
s *= -1
Energy += delta_E
Mag += 2*s


you could simply write

if delta_E < 0 or random.random() < expBeta[delta_E]:
s *= -1
Energy += delta_E
Mag += 2*s


## String formatting

Use the sweet f-strings.

sweep_time = int(time.time() - start_time)
print(f"{sweeps} monte carlo sweeps completed in {sweep_time} seconds")


## Logging warnings

Consider using the logging library. I'd write warnings to stderr, but it's up to you, see.

import sys

print("Warning: Number of sweeps exceeded 10000", file=sys.stderr)
print("         setting number of sweeps to 10000", file=sys.stderr)


Truncating it to a single line allows for easier parsing later.

print("Warning: Number of sweeps truncated to 10000.", file=sys.stderr)


## Refactorisation

If performance wasn't the primary goal, I'd introduce a few new functions. I'd separate the timing part anyway.

start_time = time.time()
energies, mags = monte_carlo(n, state, energy, mag, beta, sweeps)
elapsed_seconds = int(time.time() - start_time)
print(f"{sweeps} monte carlo sweeps completed in {elapsed_seconds} seconds")


## monte_carlo

Applying these ideas, the monteCarlo function becomes the following.

def monte_carlo(n, state, energy, mag, beta, sweeps):
if sweeps > 10000:
print("Warning: Number of sweeps truncated to 10000.", file=sys.stderr)
sweeps = 10000

exp_betas = np.exp(-beta*np.arange(0,9))

energies = np.zeros(sweeps)
mags = np.zeros(sweeps)
for t in range(sweeps):
for tt in range(n*n):
j = random.randint(0, n-1)
k = random.randint(0, n-1)
s = state[j,k]

neighbour_sum = (state[(j-1)%n, k] +
state[j, (k-1)%n] + state[j, (k+1)%n] +
state[(j+1)%n, k])
energy_diff = 2*s*neighbour_sum

if energy_diff < 0 or random.random() < exp_betas[energy_diff]:
s *= -1
energy += energy_diff
mag += 2*s

state[j, k] = s

energies[t], mags[t] = energy, mag

return energies, mags

• Albeit, numpy random is better for scientific and statistical analysis purposes. (stackoverflow.com/a/7030595/8081272) @Gucio I have gone through your answer and you're correct. The only problem I am facing right now is, if I run my earlier code (the one in question), the average magnetisation is zero above temp >2.269. This ain't case when I use numpy's random. What I will do is, I will try @Andrew 's code with MT as generator. – Kartik Chhajed May 15 at 9:04
• I have updated my question, pointing the problem faced with your suggestions wrt my original code. Though, your suggestion is very helpful for beginners. For example, as a newb coder, I earlier did not know when to use snake_case and when to use camelCase. – Kartik Chhajed May 15 at 11:31
• @KartikChhajed Can you post the exact code used in generating your plots? Pastebin or something would work. – Andrew May 15 at 11:48
• I have updated my question again with code. Sorry, I have not incorporated your suggestions in my code: like using snake_case. I was avoiding posting this complicated stuff earlier. Thank You for your advise given more first. – Kartik Chhajed May 15 at 12:07
• @KartikChhajed No problem. Generally it is best not to edit your original question as the answers arrive. It is also good practice to not accept answers too quickly. It seems that I have overlooked something regarding the batch-generation of these pseudo-random numbers. I've removed the relevant parts from my answer. I'll update it if I think of something useful. – Andrew May 15 at 13:30

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 :

1. numpy uses the PCG64 random generator which, according to numpy, has better statistical properties than legacy MT19937
2. numpy used to use MT19937
3. 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.

• I have tried importing from numpy.random import Generator. But, it shows import error: ImportError: cannot import name 'Generator' from 'numpy.random' (/usr/lib/python3/dist-packages/numpy/random/__init__.py)`. – Kartik Chhajed May 15 at 9:30
• Generator seems to have appeared in numpy v1.17.0, so see if your numpy version is up to date. I can run the snippet on numpy v1.18.1. – Gucio May 15 at 9:37
• An answer that elaborates on random number generation is always welcome. No need to delete yours. – Andrew May 15 at 10:27