Calculating the energy of the harmonic oscillator using a Monte Carlo method

Question

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?

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.