I am trying to calculate the energy, magnetization and specific heat of a two dimensional lattice using the metropolis monte carlo algorithm.

import numpy as np
import random

#creating the initial array
def init_spin_array(rows, cols):
    return np.ones((rows, cols))

#calcuating the nearest neighbours
def find_neighbors(spin_array, lattice, x, y):
    left   = (x, y - 1)
    right  = (x, (y + 1) % lattice)
    top    = (x - 1, y)
    bottom = ((x + 1) % lattice, y)

    return [spin_array[left[0], left[1]],
            spin_array[right[0], right[1]],
            spin_array[top[0], top[1]],
            spin_array[bottom[0], bottom[1]]]

#calculating the energy of the configuration
def energy(spin_array, lattice, x ,y):
    return 2 * spin_array[x, y] * sum(find_neighbors(spin_array, lattice, x, y))

#main code
def main10():
    #defining the number of initial sweeps, the lattice size, and number of monte carlo sweeps
    lattice = 10
    sweeps = 1000
    e1= e0 = 0
    for temperature in np.arange(0.1, 4.0, 0.2):
        #setting up initial variables
        spin_array = init_spin_array(lattice, lattice)
        mag = np.zeros(sweeps + RELAX_SWEEPS)
        spec = np.zeros(sweeps + RELAX_SWEEPS)
        Energy = np.zeros(sweeps + RELAX_SWEEPS)
        # the Monte Carlo
        for sweep in range(sweeps + RELAX_SWEEPS):
            for i in range(lattice):
                for j in range(lattice):
                    e = energy(spin_array, lattice, i, j)
                    if e <= 0:
                        spin_array[i, j] *= -1
                    elif np.exp((-1.0 * e)/temperature) > random.random():
                        spin_array[i, j] *= -1

            #Thermodynamic Variables 

            mag[sweep] = abs(sum(sum(spin_array))) / (lattice ** 2)

            Energy[sweep] = energy(spin_array,lattice,i,j)/ (lattice ** 2)

            #Specific Heat
            e0 = e0 + energy(spin_array,lattice,i,j)               
            e1 = e1 + energy(spin_array,lattice,i,j) *energy(spin_array,lattice,i,j)
            spec[sweep]=((e1/(sweeps*lattice) - e0*e0/(sweeps*sweeps*lattice*lattice)) / (temperature * temperature))

        #Printing the thermodynamic variables    

        print(temperature,sum(Energy[RELAX_SWEEPS:]) / sweeps, sum(mag[RELAX_SWEEPS:]) / sweeps,sum(spec[RELAX_SWEEPS:]) / sweeps)


I appear to be calculating the Magnetization correctly, which makes me convinced that my monte carlo algorithm is correct. But I appear to have some errors when calculating the Total Energy and the Specific Heat. I am trying to think of possible errors for them but I am genuinely struggling!

The graphs of them are shown below.



Specific Heat

Any help would be much appreciated! :)


1 Answer 1


Your specific heat is indeed not correct. You should get a peak centered on the critical temperature $T_c\simeq 2.27$. The specific heat is $$C=\big[\langle E^2\rangle-\langle E\rangle^2\big]/k_BT^2$$ where $\langle\ldots\rangle$ denotes the average over thermal fluctuations. In a Monte Carlo simulation, this average becomes $$\langle E^n\rangle\simeq {1\over\rm sweeps}\sum_{{\rm sweep}=1}^{\rm sweeps}[E({\rm sweep})]^n$$ to keep your notation. The specific heat can only be computed after all the sweeps have been performed. In your Python code, you compute the specific heat for each iteration. The calculation should be outside the loop over 'sweep'. You should accumulate $E$ and $E^2$ at each iteration (what it is done correctly with $e_0$ and $e_1$) and then, at the end of the loop, normalize them by dividing them by sweeps ($e_1=e_1/{\rm sweeps}$) and finally compute $\langle E^2\rangle-\langle E\rangle^2$ as $e_1-e_0^2$.

  • $\begingroup$ Thank you so much for you help! With your corrections it's now working :) $\endgroup$
    – P.Blah
    Nov 12, 2017 at 14:18
  • 1
    $\begingroup$ Happy to have been helpful! Consider upvoting the answer... Thanks! $\endgroup$
    – Christophe
    Nov 12, 2017 at 18:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.