# Variational Monte Carlo to calculate local energy of hydrogen like ions in python

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

• It is not clear to me what you want to achieve here: first I thought you want to minimize $E_{local}(x)$, but then you're talking about accumulations, so maybe you want to calculate the expectation value $\langle E \rangle=\int d^3x E(x) \left| \psi(x)\right|^2$. If you could edit your question to make that clear, we will be able to help. In both cases, your code contains several mistakes in the MCMC algorithm I think you're trying to implement. Also you might want to look into the numpy module to make calculations on vectors much more faster and readable. – Simeon Carstens Sep 10 '15 at 12:15
• Hi, Thanks for the reply. I'm trying to calculate local energy of an electron. Sorry if my question was worded poorly, it's probably due to the fact I'm not 100% sure of what I'm doing myself xD – istigatrice Sep 12 '15 at 2:20