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

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

    data = numpy.loadtxt(fname)
    return data.reshape(dimtuple)
#%%
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):
        tot=0.0
        lamt=self.lam*self.tau
        for j in range(0,3):
            for i in range(0,self.NumParticles):
                tot=tot+abs(self.beads[slice1,i,j]-self.beads[slice2,i,j])**2/(4*lamt) 
        return tot
    def PotentialAction(self,slice1,slice2):
        tot=0
        for j in range(0,3):
            for i in range(0,self.NumParticles):
                tot=tot+abs(self.Vext(self.beads[slice1,i,j])+self.Vext(self.beads[slice2,i,j]))
        return tot*self.tau/2
    def RelabelBeads(self):
        slicesToShift=numpy.random.randint(0,self.NumTimeSlices-1)
        l=list(range(slicesToShift,len(self.beads)))+list(range(0,slicesToShift))
        self.beads=self.beads[l].copy()
    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):
                delR = self.beads[tslicep1,ptcl,:] - self.beads[tslice,ptcl,:]
                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):
                R = self.beads[tslice,ptcl,:]
                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):
    Path.RelabelBeads()
    move=(numpy.random.rand(3)-1/2)/4
    ptcl=int(2*numpy.random.rand())
    bead_old=numpy.copy(Path.beads[1,ptcl,:])
    SoldV=Path.PotentialAction(0,1)
    SoldK=Path.KineticAction(0,1)+Path.KineticAction(1,2)
    Sold=SoldV+SoldK
    Path.beads[1,ptcl,:]=Path.beads[1,ptcl,:]+move
    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:
        Path.beads[1,ptcl,:]=numpy.copy(bead_old)
        return False
        
    
#    
#    disp=numpy.random.rand(3)
#    coord_new=
#    #add your things here
#    #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=PathClass(ReadArray("TestPath.dat"),tau,lam)
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.

Delete this entry (admin only).