Markov ising.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 22:24, 4 March 2024 Werner (Talk | contribs) (→Program) ← Previous diff |
Revision as of 22:27, 4 March 2024 Werner (Talk | contribs) Next diff → |
||
| Line 1: | Line 1: | ||
| - | This page presents the program markov_ising.py, a Markov-chain algorithm for the Ising model on an LXL square lattice in two dimensions. | + | This page presents the program markov_ising.py, a Markov-chain algorithm for the Ising model on an LXL square lattice in two dimensions, with periodic boundary conditions. |
| + | |||
| __FORCETOC__ | __FORCETOC__ | ||
| =Description= | =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= | =Program= | ||
| Line 32: | Line 33: | ||
| E2_av = E2_tot / float(nsteps) | E2_av = E2_tot / float(nsteps) | ||
| 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= | =Version= | ||
Revision as of 22:27, 4 March 2024
This page presents the program markov_ising.py, a Markov-chain algorithm for the Ising model on an LXL square lattice in two dimensions, with periodic boundary conditions.
Contents |
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
T = 1.0
beta = 1.0 / T
S = [random.choice([1, -1]) for k 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)
h = sum(S[nn] for nn in nbr[k])
Sk_old = S[k]
delta_E = 2.0 * S[k] * h
if random.uniform(0.0, 1.0) < math.exp(-beta * delta_E):
S[k] *= -1
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.
