。/._CalcStatistics.py 000755 000765 000024 00000000357 13216751510 015504 0 ustar 00bkclark staff 000000 000000 Mac OS X 2 � � ATTR � � S � S com.dropbox.attributes x��VJ)�/Hʯ�O��I�L���ON�Q�R�V�ML����%����RK� %O�䈪D7�"GK���Ȝ��@[[���Z Ä T CalcStatistics.py 000755 000765 000024 00000002636 13216751510 015134 0 ustar 00bkclark staff 000000 000000 #!/bin/env python import numpy import pylab from math import * def Mean (myArray): return sum(myArray)/len(myArray) def Var (myArray): mean = Mean(myArray) mean2 = Mean(myArray**2) N = len(myArray)+0.0 return N/(N-1) * (mean2 - mean**2) def NaiveStandardError(myArray): # calculate the standard error by calling the standard deviation function return sqrt(Var(myArray))/sqrt(len(myArray)+0.0) def BinData(data, n): nbins = len(data)/n binned = numpy.zeros(nbins,float) # fill in the binned array here for i in range(0,nbins): binned[i] = numpy.average(data[n*i:(i+1)*n+1]) return binned def BinningMethod(data): errList = [] for n in range(1,len(data)/20): errList.append(NaiveStandardError(BinData(data,n))) return errList def C(data, i, mean, var): delta = data - mean N = len(data) return numpy.sum(delta[0:N-i]*delta[i:N])/((N-i)*var) def Kappa(data): done = False N = len(data) i = 1 mean = Mean (data) var = Var (data) csum = 1.0 while ((not done) and i
\ u001b [0m in \ u001b [0; 36m
\ u001b [0; 34m（）\ u001b [0m \ n \ u001b [1; 32m 1 \ u001b [0m \ u001b [0; 32mimport \ u001b [0m \ u001b [0m \ u001b [0mnumpy \ u001b \ u001b [0m \ u001b [0m \ u00001b \ u \ u001b \ \ u \ u \ u \ u \ 0; 34mmmmnmmm \ u;U001B [0m \ u001b [0m \ n \ u001b [1; 32m 2 \ u001b [0m \ u001b [0; 32mimport \ u001b [0m \ u001b [0mpylab \ u001b [0mpylab \ u001b \ u001b [0m \ u001b \ u001b [0m \ u001b [0; 34mm \ u001b \ u001b [0m \ u001b [0m \ u001b [0m \ usim \ us.34mm \ us.34mm \ us.34mm \ us.34mm \ us.34mm \ us.34mm \ us.34mm \ u [0m \ u.34mm \ u [0m \ u im][0m \ n \ u001b [0; 32m ----> 3 \ u001b [0; 31m \ u001b [0; 32mimport \ u001b [0m \ u001b [0mcalcstatistics \ u001b \ u001b [0m \ u001b [0m \ u001b [0; u001b [0; 34m \ u00001b \ u00001b [0m \ u00001b\ u001b [0m \ n \ u001b [0m \ u001b [1; 32m 4 \ u001b [0m \ u001b [0; 34m \ u001b [0m \ u001b [0m \ u001b [0m \ n \ n \ n \ n \ n \ u001b [1; 32m 5; 32m 5 \ u001b [0m \ u001b [0m \ u001b [0m \ u001b [0m \ u001b [0m \ u001b [[0; 32mdef \ u001b [0m \ u001b [0mzerofunction \ u001b [0m \ u001b [0; 34m（\ u001b [0m \ u001b [0mslice \ u001b [0mslice \ u001b [0m \ u001b [0m \ u001b [0m \ u001b [0; 34m）\ u001b [0m \ u0000b [0m \ u00001b [0m \ u00001b [0m \ u.00m [0m \ us [0m \ [0m [0m \ [0m \ [0m \ [0m \ [0m \ [0m \ [0M [0M [0M [0M [0M [0M]; 34m：\ u001b [0m \ u001b [0; 34m \ u001b [0m \ u001b [0m \ n“，”，“，” \ u001b [0; 31mimporterror \ u001b [0m：no Module nos calctatistics calcstistics} calcstatistics}]，calcstatistics}，}]，source，“ source”源。[“＃pimc.py \ n”，“导入numpy \ n”，“ import pylab \ n”，“导入calcstatistics \ n”，“从直方图导入 *\ n”，“ \ n”，“ def zerofunction（slice slice））：\ n“，”返回0.0 \ n“，” \ n“，” def readarray（fname）：\ n“，” io = open（fname，'r'r'）\ n“），” line = io.readline（）\ n“，” if（line [0]！='＃'）：\ n“”，“ print'ReadArray！'\ n“，” exit（）\ n“，” splitline = line.split（）[1：] \ n“，” dim = [] \ n“，” for splitline：\ n“，” dim.Append（int（i））\ n“，” dimtuple = tuple（dim）\ n“，” \ n“，” data = numpy.loadtxt（fname）\ n“，” return data.Reshape（dimtuple）“]}，{“ cell_type”：“ code”，“ execution_count”：null，“ metadata”：{“ collapsed”：true}，“ outputs”：[]：[]，“ source”，“ source”：[]}，{“ cell_type”：“ markdown”，“ metadata”：{}，“ source”：[“在下面，我们在文件\“ data/testpath.dat \”“]}中保存一个用于测试的路径，{” cell_type“：“代码”，“ execution_count”：null，“ metadata”：{“ collapsed”：false}，“ outputs”：[]，“ source”：[]}，{“ cell_type”：“ markdown”：“ markdown”，“ metadata“：{}，“ source”：[“ \ n”，“ ###路径对象\ n”，“ \ n”，“系统中电子路径的配置将由三维numpy数组表示。 We will group arry with other data and functions into a PathClass object. A single particle at a 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\n", "\n", "```python\n", "Path.beads[slice,ptcl,dim]\n", "```\n", "\n", "where `dim` is 0,1,2, representing the x,y, or z, coordinate. Each slice is connected to the next, with the last slice also connected to the first. \n", "\n", "*Note:* Inside `PathClass` function definitions, you will access the sam array through `self.beads.`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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\n", "```python\n", "Path.beads(slice,ptcl,dim) \n", "```\n", "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).\n", "\n", "The general structure of PathClass looks like this:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class PathClass:\n", " def __init__(self,beads,tau,lam):\n", " self.tau=tau\n", " self.lam=lam\n", " self.beads=beads.copy()\n", " self.NumTimeSlices=len(beads)\n", " self.NumParticles=len(beads[0])\n", " print(self)\n", " def __str__(self):\n", " rep = \"I have setup the path with a temperature of %6.4f and %d particles.\" % \\\n", " (1.0/(self.tau*self.NumTimeSlices),self.NumParticles)\n", " return rep\n", " def SetCouplingConstant(self,c):\n", " self.c=c\n", " def SetPotential(self,externalPotentialFunction):\n", " self.VextHelper=externalPotentialFunction\n", " def Vee(self,R):\n", " # you will write this\n", " # using self.c\n", " return 0.0\n", " def Vext(self,R):\n", " return self.VextHelper(R)\n", " def KineticAction(self,slice1,slice2):\n", " # you will fill this in\n", " tot=0.0\n", " return tot\n", " def PotentialAction(self,slice1,slice2):\n", " # you will fill this in\n", " return 0.0\n", " def RelabelBeads(self):\n", " slicesToShift=random.randint(0,self.NumTimeSlices-1)\n", " l=range(slicesToShift,len(self.beads))+range(0,slicesToShift)\n", " self.beads=self.beads[l].copy()\n", " def KineticEnergy(self):\n", " # computes kinetic energy\n", " KE=0.0\n", " return KE/(self.NumTimeSlices+0.0)\n", " def PotentialEnergy(self):\n", " # computes potential energy\n", " PE=0.0\n", " return PE/(self.NumTimeSlices+0.0)\n", " def Energy(self):\n", " return self.PotentialEnergy()+self.KineticEnergy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 keyword in Python, so we use `lam` instead). There is also a `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. \n", "\n", "To create a new PathClass object instance we can call\n", "\n", "#### TestPath.py " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#TestPath.py \n", "from PIMC import *\n", "numParticles=2\n", "numTimeSlices=5\n", "tau=0.5\n", "lam=0.5\n", "Path=PathClass(ReadArray(\"data/TestPath.dat\"),tau,lam)\n", "Path.SetPotential(ZeroFunction)\n", "Path.SetCouplingConstant(0.0)\n", "print( \"My slice 0 of the particle 1 (really the second particle) in the z dimension is %14.12f. \"%Path.beads[0,1,2] )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should find that the slice is at -0.566223936211." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*A note about RelabelBeads:*\n", "\n", "Our path integrals are periodic in imaginary time and it doesn't matter which \"slice\" we call slice 0. To simpliy 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. This way, we can always move what is currently \"slice 1, but we will still eventually move all 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.`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The Kinetic Action\n", "\n", " \n", "We need to evaluate both the Kinetic Action and the Potential Action. The Kinetic Action is of the form \n", "\n", "$$\\sum_i -\\ln\\left[\\langle R(t) | \\exp[-\\tau \\lambda \\nabla^2]R(t+\\tau)\\rangle \\right] = \\textrm{const}+\\sum_{t,i} \\frac{|r_{i,t}-r_{i,t+\\tau}|^2}{4\\lambda\\tau}$$\n", "\n", "Since the constant will not change when a path is moved, we need not evaluate it when we decide whether to accept or reject a move.\n", "\n", "Let's start out by modifying the function to our `PathClass` of the form `def KineticAction(self,slice1,slice2):\n", "`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The indices of the slices are passed into this function as slice1 and slice2.\n", "\n", "1. You can get the dot product of two vectors by calling `numpy.dot`\n", "2. You can get $\\lambda \\tau$ by calling `self.lam * self.tau`\n", "\n", "To test our function, we can add to our initialization" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#TestPath.py \n", "tau=0.5\n", "lam=0.5\n", "Path=PathClass(ReadArray(\"data/TestPath.dat\"),tau,lam)\n", "Path.SetPotential(ZeroFunction)\n", "Path.SetCouplingConstant(0.0)\n", "print \"The value of the kinetic action is \", Path.KineticAction(1,2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should get a result of 0.567148772825." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Defining Vext\n", "\n", "We will now define a function for the $V_\\textrm{ext}$. Our `PathClass` object takes a function (which we are about to write) as a paramter. In python, functions are first-class objects (https://en.wikipedia.org/wiki/First-class_function). This means that you are allowed to pass around functions as arguments to functions which can then be used. We are now going to write a function that takes a single coordinate and returns the values of the harmonic potential\n", "\n", "$$V_\\textrm{ext}(r)=\\frac{1}{2}m\\omega^2|r|^2$$\n", "\n", "where we will set $m=\\omega=1.$ Later in this tutorial, we will generalize this to a double well quantum dot. Let's start by writing the function.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#top of PIMC.py\n", "def HarmonicOscillator(pos,mass=1.0,omega=1.0):\n", " pot = 0.0\n", " return pot\n", "#return the harmonic oscillator\n", "# potential for a single particle at position r1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us now test this function by using" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#TestPath.py\n", "tau=0.5\n", "lam=0.5\n", "Path=PathClass(ReadArray(\"data/TestPath.dat\"),tau,lam)\n", "Path.SetPotential(HarmonicOscillator)\n", "Path.SetCouplingConstant(0.0)\n", "print \"The value of the external potential is\", Path.Vext(numpy.array([0.1,0.3,0.1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should get the answer 0.055." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Potential Action\n", "\n", "We are using the primitive approximation to the potential. This means that it is of the form $\\tau V.$ Fill in the potential action `PotentialAction(self,slice1,slice2):` that is currently part of our object. We will test the potential action by calling\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#TestPath.py\n", "tau=0.5\n", "lam=0.5\n", "Path=PathClass(ReadArray(\"data/TestPath.dat\"),tau,lam)\n", "Path.SetPotential(HarmonicOscillator)\n", "Path.SetCouplingConstant(0.0)\n", "print \"The value of the potential action between slice 1 and 2 (with a harmonic external potential is)\", Path.PotentialAction(1,2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The answer you expect is 2.46414551377" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Your First PIMC Code\n", "\n", "We are now going to go ahead and actually write our first PIMC code. This code will mainly\n", "be a big loop over two functions: a move you will write soon that implements a monte carlo\n", "step and a call to calculate the total energy of your system.\n", "\n", "It will look like:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#PIMC.py\n", "from numpy import *\n", "def PIMC(numSteps, Path, moveList, name='test'):\n", " observableSkip=10\n", " EnergyTrace=[]\n", " numAccept = zeros((len(moveList)),int)\n", "# PD = PathDump(name+'PathDump.h5')\n", " for steps in range(0,numSteps):\n", " for mi in range(0,len(moveList)):\n", " if (moveList[mi](Path)):\n", " numAccept[mi] += 1\n", " if steps % observableSkip==0 and steps>1000:\n", " EnergyTrace.append(Path.Energy())\n", "# PairCorrelationFunction(Path,PairHistogram)\n", "# CalcDensity(Path,DensityHistogram)\n", "# PD.dump(Path)\n", " for mi in range(0,len(moveList)):\n", " print 'Accept ratio = %1.3f' % ((numAccept[mi]+0.0)/numSteps)\n", " print CalcStatistics.Stats(numpy.array(EnergyTrace))\n", " pylab.plot(EnergyTrace)\n", " pylab.savefig(name+\"Energy.png\")\n", "# PairHistogram.plotMeNorm(name+\"PairCorrelation.png\")\n", "# DensityHistogram.plotMe(name+\"Density.png\")\n", " pylab.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You will also need to add a `SingleSliceMove` which you will fill in below." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#PIMC.py\n", "def SingleSliceMove(Path,sigma=0.1):\n", " Path.RelabelBeads()\n", " #add your things here\n", " #make sure to remember if you reject the move, to restore the old location of the coordinates\n", " return 0.0 #accepted # return true if accepted, otherwise return false" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This can all be run by calling:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tau=0.1\n", "lam=0.5\n", "numTimeSlices=5\n", "numParticles=2\n", "Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam)\n", "Path.SetPotential(ZeroFunction)\n", "Path.SetCouplingConstant(0.0)\n", "NumSteps=2000\n", "print PIMC(NumSteps,Path,[SingleSliceMove])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that in the above section of code, the energy is only evaluated every `observableSkip`\n", "steps (in this case 10). This is simply because evaluating the energy is computationally\n", "expensive, and it does not help to evaluate an observable much more often than every\n", "steps (recall that is the autocorrelation time).\n", "\n", "Of course, all the work is in writing the move that implements the monte carlo step. Let's\n", "do that now" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Our First Move: Single Slice Displacements\n", "\n", "We now would like to implement a single step of the Monte Carlo algorithm.\n", "\n", "We will write a function of the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Your move should move (or not if you reject) the appropriate bead and make sure that it\n", "updates the accepted and attempted number of moves.\n", "\n", "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 1).\n", "\n", "This function must implement a single step of the Monte Carlo algorithm. Let's recall how a Monte Carlo algorithm works again:\n", "\n", "1. Choose new coordinates for your system\n", "\n", "2. Evaluate the acceptance probability and accept if the ratio is greater then a random number uniformly chosen in $(0,1)$\n", "\n", "Let's look at each of these pieces in turn:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Choosing new coordinates \n", "\n", "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.\n", "\n", "Then we will choose a new place fo rthis 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 you reject the move.\n", "\n", "A useful function call to know if 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.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Evaluating the acceptance probability\n", "\n", "In metropolis markov chain Monte Carlo, the probability of accepting is \n", "\n", "$$\\exp[-S_\\textrm{new}-S_\\textrm{old}] \\frac{T(\\textrm{new}\\rightarrow \\textrm{old})}{T(\\textrm{old}\\rightarrow \\textrm{new}]}$$\n", "\n", "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 dfference 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 a very slow code. instead just evalaute the difference in the pieces of the path that you have changed.\n", "\n", "Useful things to keep in mind:\n", "\n", "1. It might be easier to work with log of probabilities instead of the probablities themsleves.\n", "2. Change one bead changes two links in the kinetic action.\n", "\n", "The second term,\n", "$$\\frac{T(\\textrm{new}\\rightarrow \\textrm{old})}{T(\\textrm{old}\\rightarrow \\textrm{new}]}$$\n", "\n", "is the probability that you've selected the current new configuration over the probability that you would return to the old configuration given that you started in the new configuration. Because we are simply moving something uniformly with a box, this ration is going to be 1.0 for this move.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Putting it together\n", "\n", "You should proceed to write this function by implementing these two pieces. Here are a copule things to keep note of:\n", "\n", "1. you will first evaluate the old action on what you are about to move, move it, and then evaluate the new action.\n", "\n", "2. Remember that if you reject to resotre the pieces you've changed (which means you have to save it beforehand).\n", "\n", "3. 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.\n", "\n", "Once you have finished writing your function, we need to come up with a way to test it. To go about this, let's go ahead and implement the kinetic energy. Once we have that, then we can go ahead and test a non-interacting particle in a box. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Kinetic Test\n", "\n", "In testing monte carlo codes, it's very important that we test things as incrementallly 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\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "#TestPath.py\n", "\n", "tau=0.5\n", "lam=0.5\n", "numTimeSlices=5\n", "numParticles=2\n", "Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam)\n", "Path.SetPotential(ZeroFunction)\n", "Path.SetCouplingConstant(0.0)\n", "numSteps=5\n", "print PIMC(numSteps,Path,[SingleSliceMove],\"SingleSlice\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When debugging your code, you should first run with a few number of steps (<5000) and\n", "then once it appears to be working you should run it with more steps (~50000) to lower the\n", "error bars.\n", "\n", "Now we will be working with two non-interacting particles. Let's now figure out what answer\n", "you should expect. Because there is no length scale in the problem for the polymer to\n", "\"compare\" itself too, it can't distinguish a quantum world from a classical world.\n", "Therefore, you expect the energy to reflect the classical energy of a particle at this\n", "temperature. By the equipartition theorem, this would be $3/2 KT$ per degree of freedom.\n", "\n", "Compare your answer to this and verify that it agrees.\n", "\n", "Don't continue until you can succesfully 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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### External Potential Test:\n", "\n", "Now, let's turn the external potential back on and see if we can reproduce the correct\n", "result for 2 quantum particles in a harmonic well. To do this we make the same calls with a\n", "different function sent to our class:\n", "TestSingleSlicePotential\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#PIMC.py\n", "\n", "Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam)\n", "Path.SetPotential(HarmonicOscillator)\n", "Path.SetCouplingConstant(0.0)\n", "moveList = [SingleSliceMove]\n", "print PIMC(numSteps,Path,moveList,\"SingleSlicePotential\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because this is a harmonic oscillator, we should be able to calculate the energy\n", "analytically for any temperature. In order to do this, you can calculate all the\n", "eigenstates of the harmonic oscillator and then simply occupy them with the probability\n", "Add a function to your code like" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def HarmonicEnergy(tau,numTimeSlices,numParticles):\n", " # none\n", " return 0.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "that computes the harmonic oscillator answer.\n", "\n", "Run your system and compare the analytic results with those calculated exactly. Note: you\n", "may notice that you are getting the incorrect answer! Let's try to understand why this is\n", "the case." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add the following function to your PathClass" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ " @staticmethod\n", " def draw_beads_3d(ax,beads):\n", " \"\"\" draw all beads in 3D\n", " Inputs:\n", " ax: matplotlib.Axes3D object\n", " beads: 3D numpy array of shape (nslice,nptcl,ndim)\n", " Output:\n", " ptcls: a list of pairs of plot objects. There is ony entry for each particle. Each entry has two items: line representing the particle and text labeling the particle.\n", " Effect:\n", " draw all particles on ax \"\"\"\n", "\n", " nslice,nptcl,ndim = beads.shape\n", " com = beads.mean(axis=0) # center of mass of each particle, used to label the particles only\n", "\n", " ptcls = []\n", " for iptcl in range(nptcl):\n", " mypos = beads[:,iptcl,:] # all time slices for particle iptcl\n", " pos = np.insert(mypos,0,mypos[-1],axis=0) # close beads\n", "\n", " line = ax.plot(pos[:,0],pos[:,1],pos[:,2],marker='o') # draw particle\n", " text = ax.text(com[iptcl,0],com[iptcl,1],com[iptcl,2],'ptcl %d' % iptcl,fontsize=20) # label particle\n", " ptcls.append( (line,text) )\n", " return ptcls" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and use to plot things" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D # enable 3D\n", "fig = plt.figure()\n", "ax = fig.add_subplot(1,1,1,projection='3d',aspect=1)\n", "ax.set_xlabel('x',fontsize=16)\n", "ax.set_ylabel('y',fontsize=16)\n", "ax.set_zlabel('z',fontsize=16)\n", "\n", "Path.draw_beads_3d(ax,Path.beads)\n", "\n", "fig.tight_layout()\n", "#fig.savefig('figures/test_beads.png',dpi=300)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you visualize your paths enough, you will see that the system is not very ergodic. It's not moving very far from your initial configuration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Our second move: Displace\n", "\n", "To increase efficiency, we will add another type of move to our simulation. While the\n", "single-slice move does a reasonable job changing the shape of the path, it fails to move\n", "the center of mass very quickly. Fortunately, a very simple move exists which complements\n", "this move. The \"displace\" move rigidly translates an entire electron path by some random\n", "displacement vector. Please implement the move below" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#PIMC.py \n", "def DisplaceMove(Path):\n", " # First, compute the total potential action for all time slices.\n", " # This move won't change the kinetic action, so you don't need to compute it\n", " # Don't forget the link action from the last slice back to the first!\n", " oldV = 0.0\n", " # Save a copy of the old path\n", " savePath = Path.beads.copy()\n", " # Now, create a random vector for displacement\n", " delta = 4.0*numpy.array([random.random()-0.5, random.random()-0.5, random.random()-0.5])\n", " # move all the time slices\n", " # Compute the new potential action\n", " newV = 0.0\n", " # Accept or reject based on the change in potential action\n", " # Remember to copy back savePath if you reject\n", " # Remember to return True if you accept and False if you reject" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "that shifts all the slices of the particle by some amount delta randomly selected from a\n", "box. Notice that when you do this rigid shift, none of the spring lengths will change\n", "(since the imaginary time slices don't change with respect to each other). Therefore, you\n", "just need to calculate the potential action. Remember that our loop is closed so there is a\n", "closed link between slice NumTimeSlices-1 and slice 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Potential Test: Take 2\n", "\n", "Before we noticed when running our simulation that:\n", "1. the error bars were large\n", "2. the autocorrelation time that was calculated was large\n", "3. the autocorrelation time as estimated from looking at our energy trace was large\n", "4. the visualized paths moved very slowly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we noticed when running our simulation that:\n", "\n", "1. the error bars were large\n", "2. the autocorrelation time that was calculated was large\n", "3. the autocorrelation time as estimated from looking at our energy trace was large\n", "4. the visualized paths moved very slowly.\n", "\n", "We will now do the same simulation (with the same number of steps) but include a displace\n", "move. To test this move, run the simulatino again. Make sure you add the displace move to your move list!\n", "\n", "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\n", "\n", "Go ahead and visualize your paths again. \n", "\n", "You should notice a major difference showing that between steps the paths move significantly further." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Our third move: Lowering the temperature and the BisectionMove\n", "\n", "We now have a Monte Carlo code that works reasonably well at high temperatures. We are now\n", "going to lower our temperature by including 20 time slices. To test our system, run\n", "your simulation but now include 20 time slices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "You should note the following things. To begin with, notice that\n", "1. the autocorrelation time is large\n", "2. the error bars are large.\n", "\n", "This problem is due to the fact that we are moving only a single slice at a time. Since\n", "each \"bead\" in the path is connected to its neighbors through the kinetic action, each\n", "single-slice move cannot move the bead very far. (Imagine holding hands in a ring of 100\n", "people, and each person, one at a time, takes a turn taking a step in a random direction.)\n", "To improve the efficiency, it is better to move a continuous section of time slices\n", "simultaneously. We are going to try to remedy that now. Write a bisection move that moves 7 slices at a time. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To appreciate the difference between these two moves, compare the error bars and\n", "autocorrelation time when the system is equilibrated with SingleSlice move as compared with\n", "Bisection. (Of course, this latter move is somewhat slower but you still normally win out\n", "significantly in the end)\n", "Although the difference is certainly noticeable with the 20 slices that we have been using,\n", "it is especially important as we increase quantum effects (for example with 50 slices the\n", "SingleSliceMove becomes unusable)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The final steps\n", "\n", "#### Observables\n", "\n", "Beyond the energy, we would also like to know where our electrons are in our system as well\n", "as the correlations between them. This will become especially true when we have the double\n", "well potential. Some key properties we might care about include:\n", "1. the density along the x direction\n", "2. the pair correlation function between the electrons\n", "\n", "\n", "Use the following function:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#top of PIMC.py\n", "def CalcDensity(self,DensityHistogram):\n", " for slice in range(0,self.NumTimeSlices):\n", " dist=numpy.sqrt(dot(self.beads[slice,0],self.beads[slice,0]))\n", " DensityHistogram.add(self.beads[slice,0,0])\n", " dist=numpy.sqrt(dot(self.beads[slice,1],self.beads[slice,1]))\n", " DensityHistogram.add(self.beads[slice,1,0])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using this as inspiration, you should write the function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def CalcPairCorrelation(Path,PairHistogram):\n", "#simply stores the current pair correlation in a histogram" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember, that to add a value to a histogram you need to call" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "PairHistogram.add(myVal)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Notice that in order to calculate these properties, we are going to want to sum over all\n", "the slices in imaginary time.)\n", "Once this function is written, we have to initialize these histograms and plot them inside\n", "our PIMC function. You may have noticed the lines inside the PIMC function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "DensityHistogram=Histogram(minVal,maxVal,numPoints)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Histogram.PlotMe(\"fileName.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "that will initialize and plot our histogram.\n", "You should now run your simulation and plot these observables. From the density plot, it\n", "should be clear that the particles fill out the harmonic potential effectively. You should\n", "be able to look at the pair correlation plot and notice that the two particles are\n", "independent." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correlation: The Interaction Potential\n", "\n", "So far, everything we have done has involved independent electrons. Of course, the whole\n", "point of doing quantum monte carlo is so you can calculate properties of correlated\n", "electrons. So far we have been construction our PathClass with the coupling constant of the\n", "interaction potential set to zero. Now we will write a function for the interaction\n", "potential of the form\n", "\n", "$$V_{ee}(r)=\\frac{c}{r_{12}}$$\n", "\n", "where $r_{12}$ is the distance between the two particles. Depending on whether the coupling\n", "constant $c$ is positive or negative will determine whether or not the system is attractive\n", "or repulsive. Your function should call Path.c to access the coupling constant. Let's make\n", "sure the function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Vee(self,R):\n", " #compute coulumb potential\n", " #remembering to add the epsilon\n", " # for attractive potentials.\n", " # self.C accesses the C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "is filled in. We can now run the code with the repulsive potential. Make sure you modify PotentialAction to call this!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Path.SetPotential(HarmonicOscillator)\n", "Path.SetCouplingConstant(1.0)\n", "PIMC(numSteps,Path,[DisplaceMove,BisectionMove],\"Coupling\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do the pair correlation functions of this function differ from the pair correlation\n", "when there was no attraction?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding the double well\n", "\n", "At this point we are ready to add the second well to the system. In order to do this, we\n", "will write a potential for the double harmonic well\n", "\n", "$$V(r)= s((x-a)^2+y^2+z^2))((x+a)^2+y^2+z^2))$$\n", "\n", "Use $a=2$ and $s=0.1$\n", "\n", "Write a function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def DoubleWell(Path,r):\n", "# Return the double well potential\n", "# Remember that r is a 3-vector (stored as a numpy array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we can call the code with " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Path.SetCouplingConstant(0.0)\n", "print PIMC(numSteps,Path,[BisectionMove,DisplaceMove],\"DoubleWellNoCoupling\")\n", "Path.SetCouplingConstant(1.0)\n", "print PIMC(numSteps,Path,[BisectionMove,DisplaceMove],\"DoubleWellRepulsive\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now run your code and look at your density and pair correlation functions. You should\n", "notice that the particles (when repulsive) stay in opposing wells.\n", "Now we will try turning off the interaction potential entirely (TestDoubleWellNone.py).\n", "What happens to the density profile and pair correlation function?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Manipulating a quantum dot\n", "\n", "Congratulations! You have now put together a complete path integral code for examining two\n", "unpolarized electrons in a quantum dot. At this point, you might want to explore some\n", "properties of your quantum dot. For example you might examine the range of coupling\n", "constants an experimentalist might have to sample over in order to change between electrons\n", "staying in their respective dot in comparison with electrons staying in different dots.\n", "Beyond, this you could examine at what temperature, the coherence of these dots is\n", "effectively diminished." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 } ./._data 000755 000765 000024 00000000357 13154331532 013050 0 ustar 00bkclark staff 000000 000000 Mac OS X 2 � � ATTR � � S � S com.dropbox.attributes x��VJ)�/Hʯ�O��I�L���ON�Q�R�V�ML����%����RK� ��(����d��"�b���@_Ǥ@[[���Z �6 � data/ 000755 000765 000024 00000000000 13154331532 012543 5 ustar 00bkclark staff 000000 000000 data/._TestPath.dat 000644 000765 000024 00000000357 13131767615 015225 0 ustar 00bkclark staff 000000 000000 Mac OS X 2 � � ATTR � � S � S com.dropbox.attributes x��VJ)�/Hʯ�O��I�L���ON�Q�R�V�ML����%����RK� %��`����"� ��tSw� �,��r[[���Z Ǝ � data/TestPath.dat 000644 000765 000024 00000001302 13131767615 014777 0 ustar 00bkclark staff 000000 000000 # 5 2 3 -3.5178018924153720e-01 6.2875679856876687e-01 4.2611778784508614e-01 1.1981425277030739e+00 2.5681512468888439e+00 -5.6622393621149902e-01 -1.1562430517111322e-01 3.0660854553761441e-01 5.5583028021140923e-01 1.3850326160173481e+00 2.8440658233256326e+00 -8.0243729576306200e-01 2.5309785160086284e-01 3.2006380503135867e-01 4.5944141171511987e-01 1.5814316112731144e+00 2.2311004424598333e+00 -8.8858694738666821e-01 -1.7950804157158107e-02 6.3631711955333259e-01 4.9708738798390195e-01 1.3167148487412281e+00 2.2582335413248185e+00 -5.3700637074127122e-01 -2.8268013770187550e-01 5.2561455853800843e-01 5.0711287143310191e-01 9.6393011514981985e-01 2.4074020605693587e+00 -1.2504473280922879e-01