# The problem

The partition function for the quantum harmonic oscillator can be written in the path integral formulation as $$Z\propto\int Dx(\tau)\exp\left(-\frac{S_E}{\hbar}\right)=\int Dx(\tau)\exp\left(-\frac{\int_0^{\beta \hbar}d\tau\left(\frac{1}{2}m\dot x^2(\tau)+\frac{1}{2}m\omega^2 x^2(\tau)\right)}{\hbar}\right)$$

To attack this problem numerically, we can make this quantity discrete by choosing a sequence of N points x1,...,xN on the path x(τ) with distance a between each other; rescaling these variables we can write $$\frac{S_E^{discrete}}{\hbar}=\sum_{j=1}^N\left(\frac{1}{2\eta}(y_{j+1}-y_j)^2+\frac{\eta}{2}y_j^2\right)$$ where η=aω. The continuum is then recovered in the limit of small η.

In this discrete form, the internal energy of the system transforms as follows: $$\frac{U^{cont}}{\hbar\omega}=\partial_\beta\log Z=\hbar\omega\left(\frac{1}{2}+\frac{1}{e^{N\eta}-1}\right) \ \Rightarrow \ \frac{U^{discrete}}{\hbar\omega}=\frac{1}{2\eta}-\frac{1}{2\eta^2}\langle \Delta y^2\rangle+\frac{1}{2}\langle y^2\rangle$$ where the averages are first taken over a single path, and then over all paths generated in the Monte Carlo simulation according to the distribution function $$p(y_1,...,y_N)=e^{-S_E^{discrete}/\hbar}.$$

My goal is to compute the internal energy and compare it to the theoretical value for various values of N, in function of the inverse temperature Nη.

# The algorithm

In order to sample paths that follow this distribution, we can employ a MCMC method based on the Metropolis algorithm:

1. Initialize the discrete path in a cold configuration (all zeros) or hot configuration (random numbers)

2. For each Monte Carlo sweep, loop on all j=1,...,N and propose to modify y_j with a uniformly sampled y_p in the interval [y_j-δ,y_j+δ] for some parameter δ.

3. The update is accepted with probability $$r=e^{-\Delta S_{E}^{discrete}},$$ where the difference in the action (between the original and the modified path) can be easily computed as
\begin{align*} \Delta S_{E}^{discrete}&=(y^p)^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)-\frac{y^p}{\eta}(y_{j_0+1}+y_{j_0-1})-y_{j_0}^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)\\&-\frac{1}{2}y_{j_0}(y_{j_0+1}+y_{j_0-1}). \end{align*}

# The code

As kindly reviewed by a user of this forum, I report here the code that I wrote of a MCMC simulation of the quantum harmonic oscillator in the path integral formulation:

import numpy as np
import matplotlib.pyplot as plt
from numpy.random import default_rng, Generator
import time as tm

def metropolis_ho(
path: np.ndarray,
rand: Generator,
eta: float,
delta: float,
ntimes: int,
equilibrium_start: int = 10_000,
) -> tuple[
np.ndarray,  # observables 1
np.ndarray,  # observables 2
]:
"""Monte Carlo Simulation"""
n = len(path)

# Initialize arrays of observables
obs1 = np.empty(ntimes)
obs2 = np.empty(ntimes)

# Useful constants
c1 = 1/eta
c2 = c1 + eta/2

# Iterate loop on all sites
for i in range(ntimes):
for j in range(n):
for _ in range(3):
# Set y as j-th point on path
y = path[j]

# Propose modification
y_p = rand.uniform(y - delta, y + delta)

# Calculate accept probability
force = path[(j + 1) % n] + path[j - 1]
p_ratio = c1*force*(y_p - y) + c2*(y*y - y_p*y_p)

# Accept-reject
if rand.random() < np.exp(p_ratio):
path[j] = y_p

# Average of y^2 on the path
obs1[i] = np.dot(path, path)/n

# Average of Delta y^2 on the path
diff = path - np.roll(path, -1)
obs2[i] = np.dot(diff, diff)/n

# Get rid of non-equilibrium states and decorrelate
n_term = equilibrium_start
obs1 = obs1[n_term:]
obs2 = obs2[n_term:]

return obs1, obs2

def U(obs1,obs2,eta):
"""Computes internal energy"""

c1=1/(2*eta)
c2=1/(2*eta**2)

av1=np.average(obs1)
av2=np.average(obs2)

return c1-c2*av2+av1/2

def exact_U(n,eta):
""""Theoretical expectation for internal energy"""

return 0.5+1./(np.exp(n*eta)-1)

begin=tm.time()

#Initialize arrays to plot U(n)
n_array=np.asarray([5,10,15,20,25,30,40,50])
U_array=[]

for n in n_array:

#Initialize path
hot = True
rand: Generator = default_rng(seed=0)

if hot:
start = rand.uniform(low=-1, high=1, size=n)
else:
start = np.zeros(n)

#Run MCMC
obs1, obs2 = metropolis_ho(
path=start,
rand=rand,
eta=0.01,
delta=1,
ntimes=80000,
equilibrium_start=1000,
)

#Compute energy
U_array.append(U(obs1,obs2,eta=0.01))

plt.plot(np.dot(n_array,0.01),U_array,linestyle='None',marker='s',color='k')

#Set arrays for plotting exact U
n_array=np.arange(0.01,300,0.1)
exact_U_array=[]

for n in n_array:

#Compute theoretical value
exact_U_array.append(exact_U(n,eta))

plt.plot(np.dot(eta,n_array),exact_U_array,'b')
plt.ylim(0,12)
plt.xlim(0,3)
plt.show()

end=tm.time()
print(f"Total runtime: {end-begin}.")


# Results

Even without error bars, it is pretty evident that my results are just not that accurate. Can the code be improved or fixed in order to obtain a closer fit?

• Can you explain why markov-chain is relevant? It's not obvious from the code description. Commented Oct 17, 2022 at 16:32
• @TobySpeight Sorry, I hoped it would be clear from the 'theory' part. The generated paths constitute a Markov chain, as each path is obtained after a Monte Carlo sweep from the previous one. The goal is to reach an equilibrium where paths are generated according to the specified distribution. Commented Oct 17, 2022 at 17:19
• Thanks - it wasn't obvious to me, but I think some of the theory went over my head. It's been a quarter-century or more since I last had contact with Markov processes! Commented Oct 17, 2022 at 19:34
• Which part of the algorithmic description corresponds to the loop over range(3)? Commented Oct 17, 2022 at 22:03
• @Reinderien Good catch, that was not described there. It is actually a technique from a paper that is supposed to make everything better, but it's quite fuzzy and may probably be removed without consequence. Commented Oct 17, 2022 at 22:59

You will probably get feedback more specific to the physics problem you are solving at physics.stackexchange.com
[Thanks to Reinderien in the comments for clarification.]

As for the code itself, there are some issues that may make debugging harder:

• Bad naming - names should explain what they are actually naming
# Useful constants
c1 = 1/eta
c2 = c1 + eta/2


I can't tell what c1 means by looking at its name.

def U(obs1,obs2,eta):
"""Computes internal energy"""


Should be named get_internal_energy().

• Comments should explain why, not what
# Set y as j-th point on path
y = path[j]


I can see what it does by reading the code. The question is why do we need this line?

def U(obs1,obs2,eta):
"""Computes internal energy"""


You write in the docstring what should be in the name of the method.

• Inconsistent spacing

Your spacing is all over the place, though this issue is minor and can be fixed easily. Reformat the file (e.g. Ctrl+Alt+L in PyCharm) with whatever linter you're using, it inreases readability at zero cost.

• Hi, thanks for your answer. Are you sure this is not the right site? I am not looking to fix the algorithm, which I believe is correct (although I may be proven wrong). Rather, I'm looking to improve the accuracy of the results. Anyway, I will make the fixes you suggested to the code. Commented Oct 17, 2022 at 17:19
• @MyCodeisaFlyingCircus what does the accuracy of the results depends on if not the algorithm? Commented Oct 17, 2022 at 17:31
• Sorry, I thought you meant changing the algorithm altogether. You may be right, I'm new to SE and don't really know the best site for this kind of questions. I posted this here because it seemed one of the most active and responsive of these sites. Commented Oct 17, 2022 at 17:45
• It's fine here (asterisk). CRSE frequently welcomes issues of performance but not correctness. Performance may be time efficiency etc., or "numerical performance" as is the case here. I consider this to be on topic since the algorithm basically works; it just doesn't perform well. Commented Oct 17, 2022 at 21:49
• All of that said, there is a case to be made that this question (especially now that it has excellent theoretical descriptions), with some reworking, would get feedback more specific to the numerical physics problem at physics.stackexchange.com. I think both sites offer different perspectives and it would be reasonable to cross post there. Commented Oct 17, 2022 at 21:51