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"""
    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)

#Initialize arrays to plot U(n)

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)
        start = np.zeros(n)
    #Run MCMC
    obs1, obs2 = metropolis_ho(
    #Compute energy


#Set arrays for plotting exact U

for n in n_array:
    #Compute theoretical value


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


enter image description here

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?

  • 2
    \$\begingroup\$ Can you explain why markov-chain is relevant? It's not obvious from the code description. \$\endgroup\$ Commented Oct 17, 2022 at 16:32
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ Commented Oct 17, 2022 at 17:19
  • 1
    \$\begingroup\$ 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! \$\endgroup\$ Commented Oct 17, 2022 at 19:34
  • 1
    \$\begingroup\$ Which part of the algorithmic description corresponds to the loop over range(3)? \$\endgroup\$
    – Reinderien
    Commented Oct 17, 2022 at 22:03
  • \$\begingroup\$ @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. \$\endgroup\$ Commented Oct 17, 2022 at 22:59

1 Answer 1


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.

  • \$\begingroup\$ 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. \$\endgroup\$ Commented Oct 17, 2022 at 17:19
  • \$\begingroup\$ @MyCodeisaFlyingCircus what does the accuracy of the results depends on if not the algorithm? \$\endgroup\$ Commented Oct 17, 2022 at 17:31
  • \$\begingroup\$ 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. \$\endgroup\$ Commented Oct 17, 2022 at 17:45
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – Reinderien
    Commented Oct 17, 2022 at 21:49
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – Reinderien
    Commented Oct 17, 2022 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.