Ising model simulation using metropolis algorithm

Question

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.
```

Answer 1

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

Instead of

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

which is easier to read.

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

• Kartik Chhajed: 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: 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.
• Eman Yalpsid: @KartikChhajed Can you post the exact code used in generating your plots? Pastebin or something would work.
• Kartik Chhajed: 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.
• Eman Yalpsid: @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.

Answer 2

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 :

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.

• Kartik Chhajed: 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).
• Gucio: 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.
• Eman Yalpsid: An answer that elaborates on random number generation is always welcome. No need to delete yours.