Markov ising.py
From Werner KRAUTH
(Difference between revisions)
Revision as of 21:38, 22 September 2015 Werner (Talk | contribs) ← Previous diff |
Current revision Werner (Talk | contribs) (→Version) |
||
Line 1: | Line 1: | ||
- | This page presents the program markov_disks_box.py, a Markov-chain algorithm for four disks in a square box of sides 1. | + | This page presents the Python3 program markov_ising.py, a Markov-chain algorithm for the Ising model on an LXL square lattice in two dimensions, with periodic boundary conditions. The program uses the Metropolis algorithm. |
__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= | ||
- | |||
- | import random | ||
- | L = [[0.25, 0.25], [0.75, 0.25], [0.25, 0.75], [0.75, 0.75]] | ||
- | sigma = 0.15 | ||
- | sigma_sq = sigma ** 2 | ||
- | delta = 0.1 | ||
- | n_steps = 1000 | ||
- | for steps in range(n_steps): | ||
- | a = random.choice(L) | ||
- | b = [a[0] + random.uniform(-delta, delta), a[1] + random.uniform(-delta, delta)] | ||
- | min_dist = min((b[0] - c[0]) ** 2 + (b[1] - c[1]) ** 2 for c in L if c != a) | ||
- | box_cond = min(b[0], b[1]) < sigma or max(b[0], b[1]) > 1.0 - sigma | ||
- | if not (box_cond or min_dist < 4.0 * sigma ** 2): | ||
- | a[:] = b | ||
- | print L | ||
- | |||
- | =Version= | ||
- | See history for version information. | ||
- | |||
- | [[Category:Python]] | ||
- | |||
- | |||
- | |||
import random, math | import random, math | ||
- | L = 16 | + | L = 6 |
N = L * L | N = L * L | ||
nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N, | nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N, | ||
(i // L) * L + (i - 1) % L, (i - L) % N) \ | (i // L) * L + (i - 1) % L, (i - L) % N) \ | ||
for i in range(N)} | for i in range(N)} | ||
- | nsteps = 1000000 | + | nsteps = 10000000 |
- | T = 2.0 | + | T = 1.0 |
beta = 1.0 / T | beta = 1.0 / T | ||
S = [random.choice([1, -1]) for k in range(N)] | 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): | for step in range(nsteps): | ||
k = random.randint(0, N - 1) | k = random.randint(0, N - 1) | ||
- | delta_E = 2.0 * S[k] * sum(S[nn] for nn in nbr[k]) | + | 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): | if random.uniform(0.0, 1.0) < math.exp(-beta * delta_E): | ||
S[k] *= -1 | S[k] *= -1 | ||
- | print S, sum(S) | + | 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. | ||
+ | |||
+ | [[Category:Python]] [[Category:Oxford_2024]] [[Category:MOOC_SMAC]] |
Current revision
This page presents the Python3 program markov_ising.py, a Markov-chain algorithm for the Ising model on an LXL square lattice in two dimensions, with periodic boundary conditions. The program uses the Metropolis algorithm.
Contents |
[edit]
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.
[edit]
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)
[edit]
Version
See history for version information.