I'm writing up a code to calculate the local energy for electrons in hydrogen like ions for a given wavefunction. My code is giving me weird results, which leads me to believe something is wrong.
Firstly, am I going about accepting/rejecting moves the right way? And am I doing anything weird in the accumulation phase which gives me the wrong answer?
Basically, I know from mathematica that my local energy for the trial wavefunction I'm using should be -0.48, however I keep getting a whole splash of results, with little consistency (eg 0.6, -0.38, 0.7 etc.). Am I doing anything weird?
EDIT: sorry about the hard to read code, I'm new to this coding thing and am not very sure I'm doing the right thing. Essentially, what I am trying to do is sample my local energy with the probability distribution of my trial wavefunction. I'm not sure if I'm accumulating the sampled energy correctly.
for reference here is my code
import math
import random
def psi(r): #my trial wavefunction
return (math.exp(-1.2*r))
#hamiltonian
def Ham(r): #the hamiltonian operator
return (1.44*psi(r)+(2/r)*(-1.2*psi(r)))
def Top(r): #numerator of the El expression
return (-0.5*psi(r)*Ham(r)*r*r)
def But(r): #denomenator of the El expression
return (psi(r)*psi(r)*r*r)
def Elocal(r): #defines the local energy
return (Top(r)/But(r))
def r(x,y,z):
return math.sqrt(math.pow(x,2)+math.pow(y,2)+math.pow(z,2))
x0 = random.uniform(0, 50)
y0 = random.uniform(0, 50)
z0 = random.uniform(0, 50)
#local energy
El = 0
#my movement weight
dr = 0.5
#begin equilibriate
n = 0
for n in range (1000):
dx = random.uniform(-1,1)*dr
dy = random.uniform(-1,1)*dr
dz = random.uniform(-1,1)*dr
x = x0+dx
y = y0+dy
z = z0+dz
dpsi = psi(r(x,y,z))
opsi = psi(r(x0,y0,z0))
W = (dpsi*dpsi)/(opsi*opsi)
if W >= random.random():
x0 = x
y0 = y
z0 = z
else:
x0 = x0
y0 = y0
z0 = z0
n = n+1
# begin accumulate
w = 0
N = 0
for N in range (1000):
dx = random.uniform(-1,1)*dr
dy = random.uniform(-1,1)*dr
dz = random.uniform(-1,1)*dr
x = x0+dx
y = y0+dy
z = z0+dz
dpsi = psi(r(x,y,z))
opsi = psi(r(x0,y0,z0))
W = (dpsi*dpsi)/(opsi*opsi)
if W >= random.random():
x0 = x
y0 = y
z0 = z
El = El+Elocal(r(x,y,z))
w = w+1
else:
x0 = x0
y0 = y0
z0 = z0
El = El
N = n+1
print El/w