HMC harmonic.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 18:23, 19 February 2025 Werner (Talk | contribs) ← Previous diff |
Revision as of 18:24, 19 February 2025 Werner (Talk | contribs) Next diff → |
||
| Line 1: | Line 1: | ||
| ==Context== | ==Context== | ||
| - | This page is part of my public lectures on [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures Algorithms and Computations in theoretical physics] (see the announcement for a syllabus). For more information, see [[Oxford_Lectures_2025|the main page of the public lectures]]. | + | This page is part of my public lectures on [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures Algorithms and Computations in theoretical physics] at the University of Oxford (see the announcement for a syllabus). For more information, see [[Oxford_Lectures_2025|the main page of the public lectures]]. |
| ==Python program== | ==Python program== | ||
Revision as of 18:24, 19 February 2025
Contents |
Context
This page is part of my public lectures on Algorithms and Computations in theoretical physics at the University of Oxford (see the announcement for a syllabus). For more information, see the main page of the public lectures.
Python program
import math, random
def U(x):
U = 0.0
for k in range(N):
k_minus = (k - 1) % N
x_minus = x[k_minus]
if k == 0: x_minus -= L
U += (x[k] - x_minus) ** 2 / 2.0
return U
def grad_U(x, k):
k_plus = (k + 1) % N
k_minus = (k - 1) % N
x_plus = x[k_plus]
if k == N - 1: x_plus += L
x_minus = x[k_minus]
if k == 0: x_minus -= L
return 2.0 * x[k] - x_minus - x_plus
def LeapFrog(x, p, DeltaTau, I):
p = [p[k] - DeltaTau * grad_U(x, k) / 2.0 for k in range(N)]
for iter in range(I):
x = [x[k] + DeltaTau * p[k] for k in range(N)]
if iter != I - 1: p = [p[k] - DeltaTau * grad_U(x, k) for k in range(N)]
p = [p[k] - DeltaTau * grad_U(x, k) / 2.0 for k in range(N)]
return x, p
N = 8
L = 2 * N
DeltaTau = 0.1
I = 16
Iter = 100000
x = [L / N * k for k in range(N)]
U_old = U(x)
for step in range(Iter):
#
# Choose momenta and compute kinetic energy
#
p = [random.gauss(0.0, 1.0) for k in range(N)]
K_old = sum([y ** 2 / 2.0 for y in p])
#
# Do I sweeps of Leapfrog. CPU time is proportional to (2 * I + 1) * N
#
x_new, p_new = LeapFrog(x, p, DeltaTau, I)
#
# Do the Metropolis step
#
U_new = U(x_new)
K_new = sum([y ** 2 / 2.0 for y in p_new])
if random.uniform(0.0, 1.0) < math.exp(- U_new + U_old - K_new + K_old):
x = x_new[:]
U_old = U_new
Further information
References
Krauth, W. XXX
