# Demo entry 6350288

pimc

Submitted by anonymous on Mar 09, 2017 at 09:59
Language: Python 3. Code size: 5.8 kB.

```#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar  5 15:36:41 2017
Path Integral Monte-Carlo Code - PHYS 581 HW 3
author: Bora Basa
"""
#%%
import numpy
import pylab
#import CalcStatistics

def ZeroFunction(slice):
return 0.0

io = open(fname, 'r')
if (line[0] != '#'):
exit()
splitline = line.split()[1:]
dim = []
for i in splitline:
dim.append (int(i))
dimtuple = tuple(dim)

return data.reshape(dimtuple)
#%%
class PathClass:
self.tau=tau
self.lam=lam
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):
tot=0.0
lamt=self.lam*self.tau
for j in range(0,3):
for i in range(0,self.NumParticles):
def PotentialAction(self,slice1,slice2):
tot=0
for j in range(0,3):
for i in range(0,self.NumParticles):
slicesToShift=numpy.random.randint(0,self.NumTimeSlices-1)
def KineticEnergy(self):
tot = 0.0
norm = 1.0/(4.0*self.lam*self.tau*self.tau)
for tslice in range(self.NumTimeSlices):
tslicep1 = (tslice + 1) % self.NumTimeSlices
for ptcl in range(self.NumParticles):
tot = tot - norm*numpy.dot(delR,delR)
KE = 0.5*self.NumParticles/self.tau + tot/(self.NumTimeSlices)
return -KE
#return 0# KE/(self.NumTimeSlices+0.0)
#    def PotentialEnergy(self):
#        # computes potential energy
#        PE=self.PotentialAction(1,2)
#        return PE/(self.NumTimeSlices+0.0)

def PotentialEnergy(self):
'''The operator potential energy estimator.'''
PE = 0.0
for tslice in range(self.NumTimeSlices):
for ptcl in range(self.NumParticles):
PE = PE + self.Vext(R)
return PE/(self.NumTimeSlices)*0

def Energy(self):
return self.PotentialEnergy()+self.KineticEnergy()
def HarmonicOscillator(r1):
return numpy.dot(r1,r1)/2
def PIMC(numSteps, Path, moveList, name='test'):
observableSkip=10
EnergyTrace=[]
numAccept = numpy.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()
return numpy.mean(EnergyTrace)
def SingleSliceMove(Path):
move=(numpy.random.rand(3)-1/2)/4
ptcl=int(2*numpy.random.rand())
SoldV=Path.PotentialAction(0,1)
SoldK=Path.KineticAction(0,1)+Path.KineticAction(1,2)
Sold=SoldV+SoldK
SnewV=Path.PotentialAction(0,1)
SnewK=Path.KineticAction(0,1)+Path.KineticAction(1,2)
Snew=SnewV+SnewK
logPac=numpy.exp(-(Snew-Sold))
#print(logPac)
if abs(logPac)>numpy.random.rand():
return True
else:
return False

#
#    disp=numpy.random.rand(3)
#    coord_new=
#    #make sure to remember if you reject the move, to restore the old location of the coordinates
#    if Pnew>Pold*numpy.random.rand():
#        return 1.0
#    else :
#        return 0.0
#%%
tau=0.5
lam=0.5
numTimeSlices=20
numParticles=2
Path.SetPotential(HarmonicOscillator)
Path.SetCouplingConstant(0.0)
print("The value of the kinetic action is ", Path.KineticAction(1,2))
print("The value of the external potential is", Path.Vext(numpy.array([0.1,0.3,0.1])))
print("The value of the potential action between slice 1 and 2 (with a harmonic external potential is)", Path.PotentialAction(1,2))
Path=PathClass(numpy.zeros((numTimeSlices,numParticles,3),float),tau,lam)
Path.SetPotential(ZeroFunction)
Path.SetCouplingConstant(0.0)

numSteps=100000

print('<E>=',PIMC(numSteps,Path,[SingleSliceMove],"SingleSlice"))
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.