Tutorial: Writing a Path Integral Code in Python
From PIMC++
An Odyssey in Ergodicity
In this tutorial, we will write a Path Integral Code from scratch in python. At the end of this tutorial we will have a simple path integral code that will be able to calculate a properties of two (unpolarized) electrons in a quantum dot (which we will model as a double quantum harmonic well).
Throughout this tutorial we will work in atomic units (i.e. hartree and bohr radii). The relevant parameters we will work with
for electrons and inverse temperature 
(k=1.0 in our units)
Contents |
[edit] Preliminaries
To setup your environment, please run
source /afs/ictp/public/shared/smr1929/login.cshrc
Again verify that this is working by typing
TestPlot
which should produce a plot of 
.
Now call
cd /tmp/ setup_pimc
[edit] Formalism
The many-particle density matrix is defined as
where 
is a 
-dimensional vector representing the positions of all 
particles in the system.
The density matrix has the property that the thermal average of any observable quantity can be expressed as
Many properties we are interested in (the energy, pair correlation functions, etc.) are diagonal in position space. This allows us to compute properties from the diagonal of the density matrix, 
.
In PIMC, we expand the diagonal density matrix by dividing the exponential into 
pieces and inserting the resolution of unity in position space between each piece, yielding
where 
. Note that each density matrix, 
, in some sense, "connects" particle positions at time slice 
to the adjacent slice at slice 
. Thus, we can think of the progression of particle positions from time slice to time slice as a discrete path.
Haven't we just made more work for ourselves, increasing the dimensions over which we have to integrate? In some sense, we have, but in Monte Carlo integration, the computational demand does not increase rapidly with dimensionality as in standard numerical quadrature. For small 
(high temperature), we can write down very good approximations to 
. Physically, this is because particle motion is uncorrelated at very high temperature. By using Metropolis Monte Carlo to evaluate the above integral, we can then sample the density matrix at any lower temperature we want by choosing the number of time-slices, .
We can then compute thermal averages by writing
Above, we have divided in the integrand into a sampling probability and a value to average.
Let us label the sampling probability 
, where 
represents a particular set of paths for the particles. Let 
represent the value to be averaged for that path. Note that the first position in the path, 
is also the last position in the path. This indicates that the path closes upon itself, resulting in objects that resemble ring polymers.
The entire PIMC method can then be summarized by the following algorithm:
1. Move the path variables,, compute
and add to the running average. 3. Repeat until averages have sufficiently small error bar.
It turns out that it is much easier to work with the logarithm of the short-time density matrix instead of the density matrix itself. In particular, define the action, 
, as
We divide this action into two pieces: one resulting from the potential energy and one from the kinetic energy.
[edit] The potential action
We will work in a simple approximation for the potential action known as the primitive approximation.
In this approximation, the potential action is given simply by 
times the potential energy at each time slice.
In our two-electron problem, we will have an external potential representing the confining field of the dot and a modified Coulomb interaction between the particles.
[edit] The kinetic action
The kinetic energy operator for a single particle is simply 
, where 
. We can analytically work out that, in three dimensions
Note that we have used small 
and 
here to denote a single particle at two adjacent time slices. Written this way, the kinetic action looks like a spring potential acting between the time slices.
Combining out kinetic and potential actions, we can write down our complete action for a single link as
[edit] Building code pieces
[edit] Path Object
A configuration of the electron paths in the system will be represented by a three-dimensional numpy array.
We will group array with other data and functions into a PathClass object. A single particle at a single time slice is sometimes called a bead. For this reason, the path array for an PathClass instance Path will be called Path.Beads. The path variables are indexed as
Path.beads[slice,ptcl,dim]
where dim is 0, 1, or 2, representing the x, y, or z, coordinate. Each slice is connected to the next, with the last slice also connected to the first.
Nota bene: Inside PathClass function definitions, you will access the same array through self.beads.
Through this tutorial you will primarily be working with the file PIMC.py
Test File: TestPath.py
We will start with a path object. This path object sets up our paths. We will be doing two electrons in this tutorial and consequently we will have two paths. These two paths will consist of a number of beads (representing a single coordinate of some particle at some point in imaginary time) that are attached as springs. Our paths will be represented by a three-dimensional array of the form
Path.beads[slice,ptcl,dim]
where the first dimension will be the slice in imaginary time, the second dimension will be the particle and the third the dimension (x,y, or z).
The general structure of PathClass looks like this:
class PathClass: def __init__(self,beads,tau,lam): self.tau=tau self.lam=lam self.beads=beads.copy() self.NumTimeSlices=len(beads) self.NumParticles=len(beads[0]) print "I have setup the path with a temperature of",\ 1.0/(tau*self.NumTimeSlices), "and",self.NumParticles,"particles" def SetCouplingConstant(self,c): self.c=c def SetPotential(self,externalPotentialFunction): self.VextHelper=externalPotentialFunction def Vee(self,R): # you will write this # using self.c return 0.0 def Vext(self,R): return self.VextHelper(R) def KineticAction(self,slice1,slice2): # you will fill this in return tot def PotentialAction(self,slice1,slice2): # you will fill this in return 0.0 def RelabelBeads(self): slicesToShift=random.randint(0,self.NumTimeSlices-1) l=range(slicesToShift,len(self.beads))+range(0,slicesToShift) self.beads=self.beads[l].copy() def KineticEnergy(self): # computes kinetic energy return KE/(self.NumTimeSlices+0.0) def PotentialEnergy(self): # computes potential energy return PE/(self.NumTimeSlices+0.0) def Energy(self): return self.PotentialEnergy()+self.KineticEnergy()
Let's take a quick look at the structure of our object. To begin with, there is an def __init__ function that is called whenever our class is initialized and sets up the class. It defines the number of particles, time slices, as well as the time step and lambda (lambda is a reserved keywork in Python, so we use lam intead). There is also def RelabelBeads which we will explain below. The rest of the class currently consists of pieces that you are going to fill in like the actions and the energies.
To create a new PathClass object instance, we can call
numParticles=2 numTimeSlices=5 tau=0.5 lam=0.5 Path=PathClass(pylab.load("TestPath.dat"),tau,lam) Path.SetPotential(ZeroFunction) Path.SetCouplingConstant(0.0) print "My slice 0 of the particle 1 (really the second particle) in the z dimension is ",Path.beads[0,1,2]
You should find that the slice is at -0.566223936211.
A Note about RelabelBeads:
Our path integrals are periodic in imaginary time and it doesn't matter which "slice" we call slice 0. To simplify our life, we are going to add a function to our class that relabels the slice labels so the one we call "slice 0" is selected randomly ( like so). This way, we can always move what is currently "slice 1", but we will still eventually move all of the slices. This saves us from worrying about wrapping around from the last slice to the first slice in all our functions.
There is nothing physical to this and it is simply a technical book-keeping apparatus. Toward that end, we have the function RelabelBeads in PathClass.
[edit] The Kinetic Action
Test File: TestKinetic.py
We need to evaluate both the Kinetic Action and the Potential Action. The Kinetic Action is of the form
Since the constant will not change when a path is moved, we need not evaluate it when decided whether to accept or reject a move.
Let's start out by adding a function to our PathClass of the form
def KineticAction(self,slice1,slice2): # Remember that a class variable named var is accessed through # self.var # calculate the kinetic action for the spring # between slice1 and slice2. Should sum up # ptcl1 and ptcl2
The indices of the slices are passed into this function as slice1 and slice2.
- you can get the dot product of two vectors by calling
numpy.dot - you can get
by calling self.lam * self.tau
To test our function, we can add to our initialization
tau=0.5 lam=0.5 Path=PathClass(ReadArray("TestPath.dat"),tau,lam) Path.SetPotential(ZeroFunction) Path.SetCouplingConstant(0.0) print "The value of the kinetic action is ", Path.KineticAction(1,2)
You should get a result of 0.567148772825.
[edit] Defining Vext
Test File: TestPotential.py
We will now define a function for the 
. Our PathClass object takes a function (which we are about to write) as a parameter. In python, functions are first-class objects. This means that you are allowed to pass around functions as arguments to functions which can then be used. We now are going to write a function that takes a single coordinate and returns the values of the harmonic potential
where we will set 
.
Later in the tutorial, we will generalize this to a double well quantum dot. Let's start by writing the function
def HarmonicOscillator(r1): #return the harmonic oscillator # potential for a single particle at position r1
Let us now test this function by calling
tau=0.5 lam=0.5 Path=PathClass(ReadArray("TestPath.dat"),tau,lam) Path.SetPotential(HarmonicOscillator) Path.SetCouplingConstant(0.0) print "The value of the external potential is", Path.Vext(numpy.array([0.1,0.3,0.1]))
You should get the answer 0.055.
[edit] Potential Action
Test File: TestPotentialAction.py
We are using the primitive approximation to the potential action. This means it is of the form
Fill in the potential action
def PotentialAction(self,slice1,slice2): # returns the potential action # (0.5*tau * (external potential+interaction potential) # for both slices. # make sure you use the functions # so if we change it the action changes automatically
that is currently part of our object. We will test the potential action by calling
tau=0.5 lam=0.5 Path=PathClass(ReadArray("TestPath.dat"),tau,lam) Path.SetPotential(HarmonicOscillator) Path.SetCouplingConstant(0.0) print "The value of the potential action between slice 1 and 2 (with a harmonic external potential is)", Path.PotentialAction(1,2)
The answer you expect is 2.46414551377
[edit] Your first PIMC Code
We are now going to go ahead and actually write our first PIMC code. This code will mainly be a big loop over two functions: a move you will write soon that implements a monte carlo step and a call to calculate the total energy of your system.
It will look like:
def PIMC(numSteps, Path, moveList, name='test'): observableSkip=10 EnergyTrace=[] numAccept = zeros((len(moveList)),int) PD = PathDump(name+'PathDump.h5') for steps in range(0,numSteps): for mi in range(0,len(moveList)): if (moveList[mi](Path)): numAccept[mi] += 1 if steps % observableSkip==0 and steps>1000: EnergyTrace.append(Path.Energy()) PairCorrelationFunction(Path,PairHistogram) CalcDensity(Path,DensityHistogram) PD.dump(Path) for mi in range(0,len(moveList)): print 'Accept ratio = %1.3f' % ((numAccept[mi]+0.0)/numSteps) print CalcStatistics.Stats(numpy.array(EnergyTrace)) pylab.plot(EnergyTrace) pylab.savefig(name+"Energy.png") PairHistogram.plotMeNorm(name+"PairCorrelation.png") DensityHistogram.plotMe(name+"Density.png") pylab.show()
which can be run by calling
PIMC(numSteps,Path,SingleSliceMove)
Note that in the above section of code, the energy is only evaluated every observableSkip steps (in this case 10). This is simply because evaluating the energy is computationally expensive, and it does not help to evaluate an observable much more often than every 
steps (recall that is the autocorrelation time).
Of course, all the work is in writing the move that implements the monte carlo step. Let's do that now
[edit] Our First Move: Single Slice Displacements
Test File: TestSingleSlice.py, TestSingleSlicePotential.py
We now would like to implement a single step of the monte carlo algorithm.
We will write a function of the form
def SingleSliceMove(Path): Path.RelabelBeads() #add your things here #make sure to remember if you reject the move, to restore the old location of the coordinates return accepted # return true if accepted, otherwise return false
Your move should move (or not if you reject) the appropriate bead and make sure that it updates the accepted and attempted number of moves.
The call to Path.ShiftBeads() relabels the time slices so the loop begins in a random starting place (this is done so we can always move slice 0)
This function must implement a single step of the monte carlo algorithm. Let's recall how the monte carlo algorithm works again:
- choose new coordinates for your system
- Evaluate the acceptance probability and accept if this ratio is greater then a random number uniformly chosen in (0,1)
Let's look at each of these pieces in turn:
[edit] Choosing new coordinates
Different moves will change the coordinates in different ways. For this move, we will change the coordinates of "slice 1" (we start counting from "slice 0") of a random particle (which you select). Without loss of generality, we can choose "slice 1" since we have randomly rotated the loops at the beginning of this move.
Then we will choose a new place for this coordinate by uniformly displacing it by an amount chosen uniformly at random. Remember, to store what and where you've moved things because you will have to move it back to its old location if your reject the move.
A useful function call to know is that numpy.random.random(3) will give you a numpy array of 3 elements that are between 0 and 1. Also, if you have a numpy array, you can multiply and add by a scalar as you would naturally expect.
[edit] Evaluating the acceptance probability
In metropolis markov chain monte carlo the probability of accepting is
The first term is the exponential of the actions. The actions consists of the kinetic and potential actions (part of which we won't implement till later. Just use it as 0 for now). You need to evaluate the difference between the old and new versions of the action. Although you could evaluate the total old action and the total new action this would be a bad idea and result in very slow code. Instead, just evaluate the difference in the pieces of the path that have changed.
Useful things to keep in mind
- it might be easier to work with the log of probabilities instead of the probabilities themselves
- Changing one bead changes two links in the kinetic action.
The second term,
is the probability you've selected the current new configuration over the probability you would return to the old configuration given that you started in the new configuration. Because we are simply moving something uniformly within a box, this ratio is going to be 1.0 for this move.
[edit] Putting it together
You should proceed to write this function by implementing these two pieces. Here are a couple things to keep note of:
- you will first evaluate the old action on what you are about to move, move it, and then evaluate the new action
- Remember if you reject to restore the pieces you've changed (which means you have to save it beforehand)
- Remember, when you move the "bead 1", you will have to compute the kinetic action both from bead 1 to 2 and from bead 0 to 1.
Once you have finished writing your function, we need to come up with a way to test it.
[edit] Kinetic Test:
In testing monte carlo codes, it's very important that we test things as incrementally as we can. Toward that end, let's start by testing things with the external (and interparticle) potentials to be set to 0. We can do this by calling
tau=0.5 lam=0.5 numTimeSlices=5 numParticles=2 Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam) Path.SetPotential(ZeroFunction) Path.SetCouplingConstant(0.0) print PIMC(numSteps,Path,[SingleSliceMove],"SingleSlice")
You need to call
./TestSingleSlice.py 5000
where 5000 is the number of steps. Throughout the rest of the tutorial you will be able to specify the number of steps to run by appending the number of steps desired to the command line. When debugging your code, you should first run with a few number of steps (<5000) and then once it appears to be working you should run it with more steps (~50000) to lower the error bars.
Now we will be working with two non-interacting particles. Let's now figure out what answer you should expect. Because there is no length scale in the problem for the polymer to "compare" itself too, it can't distinguish a quantum world from a classical world. Therefore, you expect the energy to reflect the classical energy of a particle at this temperature. By the equipartition theorem, this would be 3/2 KT per degree of freedom.
Compare your answer to this and verify that it agrees.
Don't continue until you can successfully get the kinetic energy. Also, at this point you should note what the acceptance ratio of your move is (in a little bit, we will write a better move that has a much better acceptance ratio).
[edit] External Potential Test:
Now, let's turn the external potential back on and see if we can reproduce the correct result for 2 quantum particles in a harmonic well. To do this we make the same calls with a different function sent to our class:
Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam) Path.SetPotential(HarmonicOscillator) Path.SetCouplingConstant(0.0) moveList = [SingleSliceMove] print PIMC(numSteps,Path,moveList,"SingleSlicePotential")
Because this is a harmonic oscillator, we should be able to calculate the energy analytically for any temperature. In order to do this, you can calculate all the eigenstates of the harmonic oscillator and then simply occupy them with the probability 
Add a function to PIMC.py
def HarmonicEnergy(tau,numTimeSlices,numParticles): # none return 0.0
that computes the harmonic oscillator answer.
Run your system and compare the analytic results with those calculated exactly. Note: you may notice that you are getting the incorrect answer! Let's try to understand why this is the case.
[edit] Visualizing the Paths
Your code should have produced a file: SingleSlicePotentialPathDump.h5. This file contains snapshots of your paths as a function of monte carlo time. We are going to visualize this now.
If you type
pathvis++ SingleSlicePathDump.h5 &an visualization window should open on your display. The display can be rotated by clicking and dragging the left mouse button in the white visualization area. The scroll wheel zooms and the right button can be used to translate the visualized area. The box is not physical, and is simply a guide to the eye.
The transparent green sphere you see is an isosurface of the external potential. On startup, it isosurface is defined by 
, i.e. the average potential energy per particle. The surface can be adjusted with the slider in the upper right of the window.
The display shows the two electron paths for a single instant in the Monte Carlo simulation. You can see later "snapshots" by using the slider on the bottom of the window. Each snapshot is 10 Monte Carlo steps apart. Note that paths move rather slowly the confining potential. This indicates that our present single-slice moves are not very efficient at moving the paths. This is reflected in the long autocorrelation time for the energy.
[edit] Our second move: Displace
To increase efficiency, we will add another type of move to our simulation. While the single-slice move does a reasonable job changing the shape of the path, it fails to move the center of mass very quickly. Fortunately, a very simple move exists which complements this move. The "displace" move rigidly translates an entire electron path by some random displacement vector. Please implement the move below in PIMC.py.
def DisplaceMove(Path): # First, compute the total potential action for all time slices. # This move won't change the kinetic action, so you don't need to compute it # Don't forget the link action from the last slice back to the first! oldV = 0.0 # Save a copy of the old path savePath = Path.beads.copy() # Now, create a random vector for displacement delta = 4.0*numpy.array([random.random()-0.5, random.random()-0.5, random.random()-0.5]) # move all the time slices # Compute the new potential action newV = 0.0 # Accept or reject based on the change in potential action # Remember to copy back savePath if you reject # Remember to return True if you accept and False if you reject
that shifts all the slices of the particle by some amount delta randomly selected from a box. Notice that when you do this rigid shift, none of the spring lengths will change (since the imaginary time slices don't change with respect to each other). Therefore, you just need to calculate the potential action. Remember that our loop is closed so there is a closed link between slice NumTimeSlices-1 and slice 0.
[edit] Potential Test, Take 2
Before we noticed when running our simulation that:
- the error bars were large
- the autocorrelation time that was calculated was large
- the autocorrelation time as estimated from looking at our energy trace was large
- the visualized paths moved very slowly.
We will now do the same simulation (with the same number of steps) but include a displace move. To test this move, run TestDisplaceMove.py
Notice, this time, that the error was much smaller (and gets the correct answer) and the autocorrelation time much better. We will visualize our system again by calling
pathvis++ DisplaceMove.h5 &
You should notice a major difference showing that between steps the paths move significantly further.
[edit] Our third move: Lowering the temperature and the BisectionMove
We now have a Monte Carlo code that works reasonably well at high temperatures. We are now going to lower our temperature by including 20 time slices. To test our system, run TestLowTemp.py
You should note the following things. To begin with, notice that
- the autocorrelation time is large
- the error bars are large.
This problem is due to the fact that we are moving only a single slice at a time. Since each "bead" in the path is connected to its neighbors through the kinetic action, each single-slice move cannot move the bead very far. (Imagine holding hands in a ring of 100 people, and each person, one at a time, takes a turn taking a step in a random direction.) To improve the efficiency, it is better to move a continuous section of time slices simultaneously. We are going to try to remedy that now. Our new move will move 7 slices (trivially 
slices at a time)
We have supplied this move. It is called Bisection and we will call it with the same calls instead replacing the call to PIMC with
PIMC(numSteps,Path,Bisection,"Bisection")
To appreciate the difference between these two moves, compare the error bars and autocorrelation time when the system is equilibrated with SingleSlice move as compared with Bisection. (Of course, this latter move is somewhat slower but you still normally win out significantly in the end)
Although the difference is certainly noticeable with the 20 slices that we have been using, it is especially important as we increase quantum effects (for example with 50 slices the SingleSliceMove becomes unusable)
[edit] The final steps
[edit] Observables
Beyond the energy, we would also like to know where our electrons are in our system as well as the correlations between them. This will become especially true when we have the double well potential. Some key properties we might care about include:
- the density along the x direction
- the pair correlation function between the electrons
We have supplied you with the function
def CalcDensity(Path,DensityHistogram):
that calculates the density in the x direction.
Using this as inspiration, you should write the function
def CalcPairCorrelation(Path,PairHistogram): #simply stores the current pair correlation in a histogram
Remember, that to add a value to a histogram you need to call
PairHistogram.add(myVal)
(Notice that in order to calculate these properties, we are going to want to sum over all the slices in imaginary time.)
Once this function is written, we have to initialize these histograms and plot them inside our PIMC function. You may have noticed the lines inside the PIMC function
DensityHistogram=Histogram(minVal,maxVal,numPoints)
and
Histogram.PlotMe("fileName.png")
that will initialize and plot our histogram.
You should now run your simulation and plot these observables. From the density plot, it should be clear that the particles fill out the harmonic potential effectively. You should be able to look at the pair correlation plot and notice that the two particles are independent.
[edit] Correlation: The Interaction Potential
So far, everything we have done has involved independent electrons. Of course, the whole point of doing quantum monte carlo is so you can calculate properties of correlated electrons. So far we have been construction our PathClass with the coupling constant of the interaction potential set to zero. Now we will write a function for the interaction potential of the form
where 
is the distance between the two particles
Depending on whether the coupling constant 
is positive or negative will determine whether or not the system is attractive or repulsive. Your function should call Path.c to access the coupling constant.
Let's make sure the function
def Vee(self,R): #compute coulumb potential #remembering to add the epsilon # for attractive potentials. # self.C accesses the C
is filled in. We can now run the code with the repulsive potential
Path.SetPotential(HarmonicOscillator) Path.SetCouplingConstant(1.0) PIMC(numSteps,Path,[DisplaceMove,BisectionMove],"Coupling")
How do the pair correlation functions of this function differ from the pair correlation when there was no attraction?
[edit] Adding the double well
At this point we are ready to add the second well to the system. In order to do this, we will write a potential for the double harmonic well
Use a=2 and s=0.1
Write a function
def DoubleWell(Path,r): # Return the double well potential # Remember that r is a 3-vector (stored as a numpy array)
Then we can call the code with
Path.SetCouplingConstant(0.0) print PIMC(numSteps,Path,[BisectionMove,DisplaceMove],"DoubleWellNoCoupling") Path.SetCouplingConstant(1.0) print PIMC(numSteps,Path,[BisectionMove,DisplaceMove],"DoubleWellRepulsive")
Now run your code and look at your density and pair correlation functions. You should notice that the particles (when repulsive) stay in opposing wells.
Now we will try turning off the interaction potential entirely (TestDoubleWellNone.py). What happens to the density profile and pair correlation function?
[edit] Manipulating a quantum dot
Congratulations! You have now put together a complete path integral code for examining two unpolarized electrons in a quantum dot. At this point, you might want to explore some properties of your quantum dot. For example you might examine the range of coupling constants an experimentalist might have to sample over in order to change between electrons staying in their respective dot in comparison with electrons staying in different dots. Beyond, this you could examine at what temperature, the coherence of these dots is effectively diminished.
