计算科学 - 变分蒙特卡洛：变分能量低于基态能量 - 吾爱随笔录

我正在为氢和氦原子编写 VMC 模拟，但在我的两个代码中，某些波函数的变分能量不仅在统计上与我的期望值不同，而且还低于我的基态能量。

对于氢原子，我知道基态能量是 $-0.5Ha$ . 如果我使用试验波函数 $e^{-1.2r}$ （在哪里 $r$ 是我的原子核和电子之间的距离）我得到的能量低于我的基态能量（我被告知我的基态能量是我的变分能量的下限）。如果我使用确切的波函数（ $e^{-r}$ ) 我明白了 $-0.5Ha$ 确切地。

我对此感到非常困惑-我猜是因为对于某些试验波函数，我的变分能量低于我的基态能量，所以我必须做一些非常错误的事情-但是如果我使用确切的波函数，我确实得到了确切的能量。为什么会这样？

这里是我的代码供参考。导入数学导入随机导入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