# Variational Monte Carlo: Variational energy is lower than ground state energy

I'm writing a VMC simulation for hydrogen and helium atoms, but in both my codes my variational energy for certain wavefunctions is not only statistically different from my expectation value, but it is also lower than my ground state energy.

For the hydrogen atom, I know that the ground state energy is $-0.5Ha$. If I use the trial wavefunction $e^{-1.2r}$ (where $r$ is the distance between my nucleus and electron) I'm getting an energy LOWER than my ground state energy (I've been told that my ground state energy is a lower bound for my variational energy). If I use the exact wavefunction ($e^{-r}$) I get $-0.5Ha$ exactly.

I'm quite confused by this - I'm guessing since my variational energy is lower than my ground state energy for some trial wavefunctions I must be doing something very wrong - but I do get the exact energy if I use the exact wavefuntion. Why might that be?

For reference here is my code. import math import random import numpy

b = numpy.float_(1.2)

def psiT(r): #my trial wavefunction
return (math.exp(-b*r))

def Elocal(r):
return -0.5*((b*b)-((2*b)/r))-(1/r)

def r(x,y,z):
return math.sqrt(math.pow(x,2)+math.pow(y,2)+math.pow(z,2))

ne = 3
a = 0
e = numpy.array([random.random() for i in range (ne)])

dr = 1.05
#equilibriate
for n in range (10000):
#move my electrons
dt = numpy.array([dr*random.uniform(-1,1) for i in range (ne)])

ep = e+dt
R1 = r(e,e,e)
Rp1 = r(ep,ep,ep)

psip = psiT(Rp1)
psi = psiT(R1)
W = math.pow(psip/psi,2)

if W >= random.random():
e = ep
a = a+1
#accumulate
El = 0
q = 0
srt = 0
for n in range (10000):
#move my electrons
dt = numpy.array([dr*random.uniform(-1,1) for i in range (ne)])

ep = e+dt
R1 = r(e,e,e)
Rp1 = r(ep,ep,ep)

psip = psiT(Rp1)
psi = psiT(R1)
W = math.pow(psip/psi,2)

if W >= random.random():
El = El+Elocal(Rp1)
q = q+1
#ignore srt and s, I use them for error analysis which I have not     included here
srt = srt+math.pow(Elocal(Rp1),2)
s = math.pow((Elocal(Rp1))-(El/q),2)
e = ep

print El/q
print q


The trial wavefunction, $\exp(-1.2r)$, does not respect the cusp condition - the derivatives of the wavefunction need to cancel the $1/r$ Coulomb term. Without a correct cusp condition, the local energy diverges at the origin, which is likely causing convergence problems with the value of the integral.
In general, the leading term should be $\exp(-Zr)$ (where $Z$ is the nuclear charge). To create a trial wavefunction for hydrogen with a correct cusp, try adding higher powers of $r$, $\exp(-r + a r^2)$, where $a$ is the variational parameter.