From Werner KRAUTH

Revision as of 22:34, 4 March 2024 Werner (Talk | contribs) ← Previous diff		Current revision Werner (Talk | contribs)
Line 1:		Line 1:
-	import random, math	+	This page presents the Python3 program heatbath_ising.py, a Markov-chain algorithm for the Ising model on an L x L square lattice in two dimensions with periodic boundary conditions. This program uses the heatbath algorithm.
		+	__FORCETOC__
		+	=Description=
		+	The program is described in my 2024 Oxford Lecture No 8, and also in my book. This version of the program estimates the energy per particle, and the specific heat.
		+	=Program=
		+
		+
		+	import random, math
		+
	L = 6		L = 6
	N = L * L		N = L * L
Line 28:		Line 36:
	c_V = beta ** 2 * (E2_av - E_av ** 2) / float(N)		c_V = beta ** 2 * (E2_av - E_av ** 2) / float(N)
	print(E_av / N,c_V)		print(E_av / N,c_V)
		+
		+	=Version=
		+	See history for version information.
		+
		+	[[Category:Python]] [[Category:Oxford_2024]] [[Category:MOOC_SMAC]]

---

Current revision

This page presents the Python3 program heatbath_ising.py, a Markov-chain algorithm for the Ising model on an L x L square lattice in two dimensions with periodic boundary conditions. This program uses the heatbath algorithm.

Contents 1 Description 2 Program 3 Version

Description

The program is described in my 2024 Oxford Lecture No 8, and also in my book. This version of the program estimates the energy per particle, and the specific heat.

Program

import random, math

L = 6
N = L * L
nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N,
            (i // L) * L + (i - 1) % L, (i - L) % N) \
                                    for i in range(N)}
nsteps = 10000000
beta = 1.0
S = [random.choice([-1, 1]) for site in range(N)]
E = -0.5 * sum(S[k] * sum(S[nn] for nn in nbr[k]) \
                                for k in range(N))
E_tot, E2_tot = 0.0, 0.0
for step in range(nsteps):
    k = random.randint(0, N - 1)
    Upsilon = random.uniform(0.0, 1.0)
    h = sum(S[nn] for nn in nbr[k])
    Sk_old = S[k]
    S[k] = -1
    if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * h)):
        S[k] = 1
    if S[k] != Sk_old:
        E -= 2.0 * h * S[k]
    E_tot += E
    E2_tot += E ** 2
E_av  = E_tot / float(nsteps)
E2_av = E2_tot / float(nsteps)
c_V = beta ** 2 * (E2_av - E_av ** 2) / float(N)
print(E_av / N,c_V)

Version

See history for version information.