{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Quantum MonteCarlo (Feynmann's Path Integral)

\n", "\n", "\"\"\" From \"COMPUTATIONHL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\" by RH Landau, MJ Paez, and CC Bordeianu (deceased). Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, C Bordeianu, Univ Bucharest, 2020. Please respect copyright & acknowledge our work.\"\"\"\n", " \n", "Compute the harmonic oscillator ground state using Feynmann's path integral " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "ename": "ImportError", "evalue": "No module named vpython", "output_type": "error", "traceback": [ "\u001b[1;31m\u001b[0m", "\u001b[1;31mImportError\u001b[0mTraceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mrandom\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0mvpython\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[1;33m*\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 7\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mImportError\u001b[0m: No module named vpython" ] } ], "source": [ "# QuantumMC.ipynb Feynman Path integral for the harmonic oscillator \n", "\n", "%matplotlib notebook\n", "\n", "import random\n", "from vpython import *\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "N = 101; M = 101; xscale = 10.\n", "path = np.zeros((M), float); prob =np.zeros([M], float) # Initialize # Initialize \n", "\n", "trajec = canvas(width = 300, height = 500, title = 'Spacetime Trajectories')\n", "trplot = curve( color = color.magenta, display = trajec, radius=0.8)\n", "\n", "def trjaxs(): # plot axis for trajectories\n", " trax = curve(pos = [( - 97, - 100,0), (100, - 100,0)], color = color.cyan, canvas = trajec)# x ax\n", " curve(pos = [(0, - 100,0), (0, 100,0)], color = color.cyan, canvas = trajec) # y label(pos = (0, 110), text = 't', box = 0, display = trajec)\n", " label(pos =vector (0, - 110,0), text = '0', box = 0, canvas = trajec)\n", " label(pos =vector (60, - 110,0), text = 'x', box = 0, canvas = trajec) \n", "\n", "wvgraph = canvas(x = 340, y = 150, width = 500, height = 300, title = 'Ground State Probability')\n", "wvplot = curve(x = range(0, 100), canvas = wvgraph) # for probability\n", "wvfax = curve(color = color.cyan)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def wvfaxs(): # axis for probability\n", " wvfax = curve(pos = [(-600,-155,0), (800,-155,0)], canvas=wvgraph, color=color.cyan)\n", " curve(pos = [(0,-150,0), (0,400,0)], display=wvgraph, color=color.cyan) \n", " label(pos = vector(-80,450,0), text='Probability', box = 0, canvas = wvgraph)\n", " label(pos = vector(600,-220,0), text='x', box=0, canvas=wvgraph)\n", " label(pos = vector(0,-220,0), text='0', box=0, canvas=wvgraph) \n", "\n", "trjaxs()\n", "wvfaxs() # plot axes \n", "\n", "def energy(path): # HO energy\n", " sums = 0. \n", " for i in range(0, N-2): sums += (path[i+1] - path[i])*(path[i+1] - path[i]) \n", " sums += path[i+1]*path[i+1]; \n", " return sums \n", "\n", "def plotpath(path): # plot trajectory in xy scale\n", " for j in range (0, N): \n", " trplot.append(pos=vector(20*path[j], 2*j - 100,0))\n", " \n", "def plotwvf(prob): # plot prob\n", " for i in range (0, 100):\n", " wvplot.color = color.yellow\n", " wvplot.append(pos=vector( 8*i - 400 , 4.0*prob[i] - 150,0)) # for centered figure\n", " \n", "oldE = energy(path) # find E of path\n", "\n", "#while True: # pick random element\n", "for i in range (0,1500):\n", " rate(50) # slows the paintings\n", " element = int(N*random() ) # Metropolis algorithm\n", " change = 2.0*(random() - 0.5) \n", " path[element] += change # Change path\n", " newE = energy(path); # Find new E\n", " if newE > oldE and np.exp( - newE + oldE)<= random():\n", " path[element] -= change # Reject\n", " trplot.clear() #erase previous trajctory for animation\n", " plotpath(path)\n", " trplot.visible=True #make visible new trajectory for animation\n", " elem = int(path[element]*16 + 50) # if path = 0, elem = 50\n", " # elem = m *path[element] + b is the linear transformation\n", " # if path = - 3.0, elem = 2 if path = 3.0, elem = 98 = > b = 50, m = 16 linear trnsf.\n", " # this way x = 0 correspond to prob[50]\n", " if elem < 0: \n", " elem = 0 # negative case not allowed\n", " if elem > 100: \n", " elem = 100 # if exceed max\n", " prob[elem] += 1 # increase probability for that x value\n", " plotwvf(prob) # plot prob\n", " oldE = newE " ] } ], "metadata": { "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.13" } }, "nbformat": 4, "nbformat_minor": 2 }