Errors in Quantum Monte-Carlo for the harmonic oscillator
Post your questions on this hwk in the blog at the end of the page.
If you hand your Hwk by e-mail, PLEASE write the file in the format HW6_YOUR-NAME and send it to M. Civelli

Introduction

In Quantum Monte Carlo (QMC), errors have two origins:

the sampling and the statistical error associated to it;
the discretization of the path.

Let us first note that the discretization error is the same in QMC and in the matrix squaring method. In this homework we will treat the case of the harmonic oscillator, for which the matrix squaring can be solved analytically (see Lecture), and we will discuss in details the problem of discretization. We also present methods to measure the free energy and the internal energy of the system.

-A- The central density

In the last lecture, we discussed the exact recursion for the matrix squaring procedure of the harmonic oscillator (see SMAC equations (3.33)-(3.38)). We use this calculation to compute the density π(x = 0) at inverse temperature β = 8, for different discretizations (number of slices).

Familiarize yourself with SMAC section 3.2.2, which applies matrix squaring to the harmonic oscillator, especially with the "initial condition" (β → 0) and the recursion relations for the functions f(β) and g(β) on page 146.
Write a (trivial) Python program for evaluating numerically the density matrix with N = 2n slices. That is: start from the initial condition (3.34) with inverse temperature β equal to β / N and iterate n times the recursion relations for f and g so as to arrive at inverse temperature β. We denote the discretization step by Δτ = β / N.
Plot the corresponding central density π(x = 0) at inverse temperature β = 8 as a function of Δτ = β / N. Does it converge to the exact value (3.38) as N increases? Don't forget that the probability density π(x) is normalized.
Determine numerically the scaling of the discretization error for the central density: with which power of Δτ does it converge to the exact value?
Explain this scaling of the discretization error from what we discussed about the Trotter approximation of SMAC eqs (3.30) and (3.31).

-B- The free energy

We now study the convergence of the free energy of the harmonic oscillator.

Modify the above calculation by including the recursion for the normalization constant c(β). This allows you to determine numerically the partition function Z(β), which depends on the discretization step Δτ.
Determine numerically, for the free energy F(β) (defined from Z(β) = exp( − β F(β)) ), the scaling of the error between the exact result (derived from (3.40)) and its numerical evaluation, as a function Δτ.

-C- The internal energy

-1- Expression of the internal energy E

Explain in details why the "internal energy" (i.e. the mean energy) of the harmonic oscillator is given by
$\displaystyle\langle E \rangle = \frac{\int dx H_x \rho(x,x'\to x,\beta)}{\int dx \rho(x,x,\beta)} \qquad (\star)$
where Hx and x→x' mean that the Hamiltonian H acts only on the first argument (x) of the density matrix, and that we then take the limit x→x'. You can have a look to equation (3.42).

-2- An estimator for the internal energy

A QMC simulation of a particle in a harmonic trap with N slices, using the Trotter approximation, generates configurations (x0,x1,...,xN−1)

1. Show in details (this is not straightforward) that

$\displaystyle \langle E \rangle = \Big\langle \frac{ H_{x_0} \rho(x_0,x_1,\beta/N)}{\rho(x_0,x_1,\beta/N)} \Big\rangle$

where the average on the right hand side is over paths (x0,x1,...,xN−1).
Hint: start from expression (٭) above and reexpress the unnormalized density matrix ρ(x,x',β) using its "path-integral" formulation (with N slices).

2. Show furthermore that (writing t instead of β / N)

$\displaystyle \langle E \rangle =\Big\langle x_0^2/2 + \frac{4 t + 2 t^3 - t^4 x_0^2 - 4 t^2 x_0(x_0-x_1) - 4(x_0-x_1)^2}{8 t^2} \Big\rangle \quad (\star\star)$

You might use (and justify why you use) the high-temperature asymptotic expression of the density matrix.

3. Implement this estimator in the below naive simulation program for a particle in a harmonic potential, the cousin of last week's program for the Morse potential:

from math import pi, exp, sqrt
from random import uniform as ran, randint as nran
def rho_free(x,xp,beta):
   return exp( - (x-xp)**2 / (2.*beta) )
N, beta, delta = 4, .5, 0.5
del_tau = beta/N
x = [0. for k in range(N)]
for iter in range(1000000):
   k = nran(0,N-1)
   x_old = x[k]
   x_new = x_old + ran(-delta,delta)
   xp = x[(k+1)%N]
   xm = x[k-1]
   pi_old = rho_free(xm,x_old,del_tau)*rho_free(x_old,xp,del_tau) * \
      exp(-del_tau*x_old**2/2.)
   pi_new = rho_free(xm,x_new,del_tau)*rho_free(x_new,xp,del_tau) * \
      exp(-del_tau*x_new**2/2.)
   Upsilon = pi_new/pi_old
   if (ran(0,1) < Upsilon): x[k] = x_new
   if iter%10000 == 0:
      print 'Path: ', x

4. Evaluate numerically the internal energy for different temperatures and discretizations. Check that the average corresponds to the analytical result. Study the error (using the bunching algorithm).
5. Generalize the above estimator (٭٭) for an arbitrary smooth (i.e. at least twice differentiable) potential.

[Print this page]

Homework6