我正在编写一个代码来计算给定波函数的氢中电子的局部能量，例如离子。我的代码给了我奇怪的结果，这让我相信出了点问题。

首先，我是否以正确的方式接受/拒绝移动？我是否在积累阶段做了什么奇怪的事情，这给了我错误的答案？

基本上，我从mathematica 知道，我正在使用的试验波函数的局部能量应该是-0.48，但是我不断得到一个完整的结果，几乎没有一致性（例如0.6、-0.38、0.7 等）。我在做什么奇怪的事吗？

编辑：对难以阅读的代码感到抱歉，我对这个编码很陌生，我不太确定我在做正确的事情。本质上，我想做的是用我的试验波函数的概率分布对我的局部能量进行采样。我不确定我是否正确地累积了采样能量。

供参考这里是我的代码

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