计算科学 - 光学晶格中Wannier函数的数值计算 - 吾爱随笔录

光学晶格中Wannier函数的数值计算

计算科学 Python 计算物理学

2021-12-14 02:47:41

我正在研究一些光学晶格带结构（例如这里）。我对设置特征值方程没有任何问题：

H_{j j^{'}} c_{j^{'}} = E c_{j^{'}}

$H_{jj'}c_{j'}=Ec_{j'}$ 在哪里

H

$H$ 是方程化不同傅里叶分量的三对角矩阵。布洛赫波由下式给出

ψ = e^{i q x} u_{n}^{q} = e^{i q x} \sum_{j} a_{c} e^{2 i k x j}

$\psi =e^{iqx} u_n^q = e^{iqx} \sum_j a_c e^{2ikxj}$

Wannier 函数是

W = \int d q u_{n}^{q} e^{- i q x}

$\mathcal{W}=\int dq \ u_n^q \ e^{-iqx}$

所以我的伪代码是：

对于每个 $x$ ：
对于每个 $q$ ：
找到特征向量 $c_j$ . 乘以 $c_j$ 通过其傅里叶分量 $\exp(2ik j)$
将这些项相加并乘以 $\exp(iq_i x_q)$
重复 $q$ ，然后添加所有项
对所有人重复 $x$

我的python代码如下。当我绘制x,|w|^2时，我没有得到任何类似于高斯近似的东西。此外，我认为这可以矢量化，但我正在努力让它以循环形式工作。

x = np.linspace(-2,2,101)
lmax = 10
l=np.arange(-lmax,lmax+1)
V0 = 5
wavelength = np.pi
k_lattice = 2*np.pi/wavelength
qx = np.linspace(-1,1,101) #This is qx/k_lattice
wave_dict = {q:None for q in qx}

for q in qx:
    diags = [(q+2*k)**2 for k in l]
    Hmat = np.diag(diags)
    Hmat = np.add(Hmat,-V0/4 *(np.diag(np.ones(len(l)-1),1)+np.diag(np.ones(len(l)-1),-1)),casting='unsafe')
    evals, evecs = LA.eigh(Hmat)

    coefs = evecs[:,0][:,None]
    planewaves = coefs*np.exp(1j*(2*k_lattice)*np.outer(l,x)) # exp(2ik_l x)
    psi = planewaves.sum(axis=0)
    wave_dict[q] = psi

w = np.zeros(x.shape,dtype = 'complex128')
for q in qx:

    w+= wave_dict[q]* np.exp(1j*x*q*k_lattice) 

w/=len(qx)    
plt.plot(x,np.abs(w)**2)
plt.plot(x,np.sin(k_lattice*x)**2)
plt.show()

1个回答

正如我在评论中所说，您不应该采用具有最高能量的特征向量，因为最高频带通常对应于非常离域的波函数，因此相应的 Wannier 函数肯定不会是高斯函数。

您还应该明确写下您的选择 $k$ 作为一个参数，因为它可以帮助防止错误。我认为除了最后一个特征向量的选择之外，你犯了两个主要错误：

傅里叶级数发展写道 $u_n^q(x) = \sum_j c_j e^{i 2 \pi x j/a}$ ，在哪里 $a$ 是你的格子的周期。但因为你已经采取 $k = 1$ ，和 $V(x) = \frac{V_0}{2} \cos(2 k x)$ ，这里格子的周期是 $a=\pi$ . 所以你的傅里叶级数应该是 $u_n^q(x) = \sum_j c_j e^{i 2 x j}$ ，并不是 $u_n^q(x) = \sum_j c_j e^{i 2 \pi x j}$ 正如你所写的。这就是为什么如果你不想被混淆，明确地命名你的变量是特别重要的。
Wannier 波函数应该在势的最小值处取。但是看看你选择的矩阵 Hmat，你似乎已经采取了 $V(x) = \frac{V_0}{2} \cos(2 k x)$ . 你可以看到 $V(x)$ 最大在_ $x=0$ ，并且您正在尝试查看 Wannier 函数 $X = 0$ ，这是没有意义的。相反，如果你想在 $X=0$ ，你应该采取 $V(x) = -\frac{V_0}{2} \cos(2 k x)$ 相反（切换非对角项的符号）。

通过这些更正，代码应该可以更好地工作。我还建议一些改进：

使用数组。它会让你的生活更轻松。例如，在您的代码中，您正在为每个位置对角化整个矩阵 $x$ 并且对于每个 $q$ ，而每次对角化一次就足够了 $q$ 仅（值不依赖于 $x$ ）。
使用 scipy.linalg.eigh 来最小化你的矩阵。它适用于 Hermitian/实对称矩阵，并为您提供实特征值，从最低到最高排序。它更有效，也可以省去您自己订购特征值/特征向量的麻烦。
试着用你知道的东西来测试你的算法。例如，在考虑计算 Wannier 函数之前，请检查您的程序是否为能带提供了正确的形式。完成此操作后，尝试获取 Wannier 函数并将它们与谐波陷阱的高斯近似值进行比较。

您可以在下面找到具有更正错误和一些改进的代码（我只是对循环进行了部分矢量化，因此它仍然可以稍微优化，但至少您没有将同一个矩阵对角化 100 次）：


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

Ei = []
phii = []
w = []

k_light = 1 #wavevector of the light beam creating the lattice
a_lattice = np.pi/k_light #period of the lattice
lmax = 20
l=np.arange(-lmax,lmax+1)
V0 = 5 #V0/E_recoil
V0 *= k_light**2 #"true" V0
x= np.linspace(-a_lattice/2*1.5,a_lattice/2*1.5,4000)
dx = x[1]-x[0]
qx = np.linspace(-k_light,k_light,100, endpoint=False)
for q in qx:
    u = 0
    diags = [(q+2*k*k_light)**2 for k in l]
    Hmat = np.diag(diags)
    Hmat += -V0/4 *(np.diag(np.ones(len(l)-1),1)+np.diag(np.ones(len(l)-1),-1))
    evals, evecs = LA.eigh(Hmat)
    Ei.append(evals)
    phii.append(evecs)
phii = np.array(phii)
Ei = np.array(Ei)

for xi in x:
    b = 0
    a = np.exp(1j*2*np.pi*xi*l/a_lattice)
    for p in range(len(qx)):
        b += np.sum(a*phii[p, :, 0])*np.exp(1j*qx[p]*xi)    
    w.append(b/(len(qx)))
w = np.array(w)    
plt.plot(x,abs(w)**2/np.sum(abs(w)**2)/dx, label='Wannier wavefunction')

X_harmonic_sq = 1/(k_light*V0**0.5) 

plt.plot(x, np.exp(-x**2/X_harmonic_sq)/(np.sqrt(np.pi*X_harmonic_sq)), label='Harmonic approximation')

plt.legend()

#plt.plot(qx, Ei[:, 0])
#plt.plot(qx, Ei[:, 1])
#plt.plot(qx, Ei[:, 2])

Wannier 波函数图 $k=1$ ，和 $V_0=5$ ，与谐波近似相比：

如果您还有其他问题，请告诉我。

其它你可能感兴趣的问题

上一篇（稀疏）2D numpy 数组的每行/列的快速非零索引下一篇使用 FFT 精确实现卷积的问题