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

Question

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
#was this move good?  
    R1 = r(e[0],e[1],e[2])
    Rp1 = r(ep[0],ep[1],ep[2])

    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
    #was this move good?  
    R1 = r(e[0],e[1],e[2])
    Rp1 = r(ep[0],ep[1],ep[2])

    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

Answer 1

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.

The Monte Carlo integration also has an error - both accepted and rejected moves should be included in the sum (this code only includes accepted moves). Otherwise individual energy values have the wrong weight in the final integral. And the final value of energy (E1) should be divided by the total number in the sum (10000), not just the number of acceptances.