计算科学 - 具有中心电位的薛定谔方程的基态原点发生了什么 - 吾爱随笔录

具有中心电位的薛定谔方程的基态原点发生了什么

计算科学数值分析计算物理学量子力学

2021-12-04 06:00:14

我有代码尝试在 1-D 中使用 numerov 方法来构造波函数，其中存在中心势（或多或少地考虑氢）的薛定谔方程。它不关心状态（N = 1，..），因为我的目标是理解理论和软件之间的联系。以下是一些能量的一些结果，无论它们是否是特征值，我不确定：

[6.054999999999915, 6.214999999999912, 6.534999999999905]

和一张图。我试图实现的算法在第 6 节 PDF中有说明，但我的潜力是不同的。我的问题是，中心势在原点有 1/r，我将如何处理？（我通过添加任意 0.00001 来消除 0 在代码中做了一个回避）。

在这里看到的另一个问题中有一个更通用的代码，但我没有遵循从方程到 3d 矩阵的转变。

这是我的代码：

import numpy as np
import scipy
import matplotlib.pyplot as plt 

h=np.float64(1)
m=60
c=m/2
estep = 1.0/(100)
delta = 0.0002
x=np.linspace(-1000,1000,num=m)
x=(x - x[m/2]+0.00001 )
v =  1.0/(x**2)
sig=None
E = []

def get_left(wf,k):
    for i in range(2,int(c)):
        wf[i] = 2*(( 1- 1/12* h**2* k[i-2]**2)*wf[i-2] - (1 + 1/12 *h**2* k[i-1]**2)*wf[i-1])/(1 + 1/12* h**2* k[i]**2)
    wf[0:c] = wf[0:c]/(np.sqrt(np.sum(wf[0:c]**2)))
    return wf

def get_right(wf,k):
    mi = m-1
    for i in range(2,int(c)):
        wf[mi-i] = 2*(( 1- 1/12* h**2* k[mi-i+2]**2)*wf[mi-i+2] - (1 + 1/12 *h**2* k[mi-i+1]**2)*wf[mi-i+1])/(1 + 1/12* h**2*k[mi-i]**2)
    wf[c:mi] = wf[c:mi]/(np.sqrt(np.sum(wf[c:mi]**2)))
    return wf

def root_search(epsilonu, epsilonl):
    en = None
    it =0
    while it < 10000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + (epsilonu+epsilonl)/2)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        erri = diffl[m/2]- diffr[m/2]
        if erri < delta:
            plt.plot(x, (phil+phir)**2 )
            plt.savefig( str('sho'+'.png' ))
            return (epsilonu+epsilonl)/2
        if erri < 0 : 
            epsilonl = (epsilonu+epsilonl)/2
        if erri > 0:
            epsilonu =epsilonl 
            epsilonl = (epsilonu+epsilonl)/2
        it+=1
    return None

def start():
    epsilonl=0.01
    eps = epsilonl
    epsilonu=None
    iteration=0
    sig = None
    while iteration < 100000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + eps)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        err = diffl[m/2]- diffr[m/2]
        ssig = 1 if err>0 else -1
        if sig:
            if sig != ssig:
                epsilonu= eps 
                E.append(root_search(epsilonu, epsilonl))
                print(E)
        sig = ssig
        epsilonl = eps
        eps = eps+estep
        iteration+=1

if __name__ == "__main__":
    start()

理论：我试图求解一维薛定谔方程，

(- \frac{h^{2}}{2 m} \frac{d^{2}}{d x^{2}} + V (x)) ϕ (x) = ϵ ϕ (x)

$\left( -\frac{h^2}{2m}\frac{d^2}{dx^2} + V\left(x\right)\right) \phi\left(x\right) = \epsilon \phi\left( x\right)$ 随着v作为核心势（库仑力，逆量子关系），最终的等式看起来像（具有正确选择的单位给q ^ 2 = 1，h-1 ...，hartree的能量）

\frac{d^{2}}{d x^{2}} ϕ (x) + 2 (ϵ - \frac{1}{x^{2}}) ϕ (x) = 0

$\frac{d^2}{dx^2}\phi\left(x\right) + 2 \left(\epsilon -\frac{1}{x^2} \right) \phi\left(x\right) =0$

1个回答

正如@Kirill 提到的，氢的微分方程略有不同。变量分离后，你最终得到

[- \frac{ℏ^{2}}{2 μ} (\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial R (r)}{\partial r}) - \frac{l (l + 1) R (r)}{r^{2}}) + V (r) R (r)] = E R (r),

$\left[ - \frac{\hbar^2}{2\mu} \left({1 \over r^2}{\partial \over \partial r}\left(r^2 {\partial R(r)\over \partial r}\right) - {l(l+1)R(r)\over r^2} \right) + V(r)R(r) \right]= E R(r),$

您可以将等式重写为

[- \frac{1}{2} \frac{d^{2}}{d y^{2}} + \frac{1}{2} \frac{l (l + 1)}{y^{2}} - \frac{1}{y}] u_{l} = W u_{l}

$\left[ -\frac{1}{2} \frac{d^2}{dy^2} + \frac{1}{2} \frac{l(l+1)}{y^2} - \frac{1}{y}\right] u_l = W u_l$

下一个代码（基于我之前的帖子）计算这个方程的特征值/特征向量。看到我为分母添加了一个小值，因此它们不会在零处消失。

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

l = 1
k = 0
T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 0.5*l*(l + 1)/(x**2 + 1e-6) - 1/(np.abs(x) + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))
print(np.round([-1/(2*n**2) for n in range(k + l + 1, k + l + 7)], 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvectors.png", dpi=600) 
plt.show()

这给出了解决方案：

[-0.125106 -0.0556   -0.031272 -0.020012 -0.013896 -0.010209]
[-0.125    -0.055556 -0.03125  -0.02     -0.013889 -0.010204]

对于您的排斥潜力（与拉普拉斯算子不同的符号），我使用了相同的代码

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 1/(x**2 + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvecs.png", dpi=600)  
plt.show()

有答案

[ 0.000408  0.001207  0.002405  0.004001  0.005997  0.008392]

虽然，我不确定这意味着什么。因为它不会让电子在原点将其推开，并且域应该是真正无界的。它类似于在左侧具有圆形边缘的方形电位。对于更实际的情况，您可以使用非均匀有限差分并以指数方式增加它们或吸收边界条件。

其它你可能感兴趣的问题

上一篇将高斯散度定理应用于FVM后，有限差分法和有限体积法有什么不同吗？下一篇了解纽曼 BC

具有中心电位的薛定谔方程的基态 原点发生了什么

具有中心电位的薛定谔方程的基态原点发生了什么