平面波法公式

$\det|\{\frac{\hbar^2}{2m}(\mathbf{k}+\mathbf{K_n})^2-E\}\delta_{KK'}+\langle\mathbf{k+K_n'}|V|\mathbf{k+K_n}\rangle|=0$

其中 $\mathbf{k}$ 为第一布里渊区内的点， $\mathbf{K_n}$ 为整数倍倒格矢。求解该方程需要选择合适的平面波截断维数以及势能的傅里叶级数展开数，以保证在计算效率。

势能处理

简单立方元胞内，势能

$V(x,y,z)=\sum_{i=1}^8\frac{A}{\sqrt{(x-x_i)^2+(y-y_i)^2}+\epsilon}$

其中

$\begin{aligned} R_1&=(0,0,0) ,R_2=(a,0,0)\\ R_3&=(0,a,0) ,R_4=(0,0,a)\\ R_5&=(a,a,0) ,R_6=(a,0,a)\\ R_7&=(0,a,a) ,R_8=(a,a,a)\\ \end{aligned}$

对势能进行傅里叶级数展开并离散化得到

$\begin{aligned} c_{\mathbf{k}}&=\int_0^a V(x,y,z)e^{-i\frac{2\pi}{a}(k_1x+k_2y+k_3z)}dxdydz\\ &=\frac{a^3}{N^3}\sum_{n_i=0}^{N-1}V(\frac{an_1}{N},\frac{an_2}{N},\frac{an_3}{N})e^{-i\frac{2\pi}{a}(k_1\frac{an_1}{N}+k_2\frac{an_2}{N}+k_3\frac{an_3}{N})}\\ &=\frac{a^3}{N^3}\sum_{n_i=0}^{N-1}V(\frac{an_1}{N},\frac{an_2}{N},\frac{an_3}{N})e^{-i\frac{2\pi}{N}(n_1k_1+n_2k_2+n_3k_3)} \end{aligned}$

由此可利用FFT处理势能。

具体代码如下

自由哈密顿量与相互作用哈密顿量的定义

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as inte
import scipy.fft as scifft
import scipy.linalg as lin
def H_3D(k1,k2,k3,a,N):
    H_3d=np.zeros([2*N-1,2*N-1,2*N-1,2*N-1,2*N-1,2*N-1])
    V_3d=np.zeros([2*N-1,2*N-1,2*N-1,2*N-1,2*N-1,2*N-1])
    G=2*np.pi/a
    for i in range(0,2*N-1):
        for j in range(0,2*N-1):
            for k in range(0,2*N-1):
                H_3d[i,j,k,i,j,k]=(k3+(i-N+1)*G)**2+(k2+(j-N+1)*G)**2+(k1+(k-N+1)*G)**2
            #H_gra[i,j,i,j]=(k1+np.abs(i-N-1))**2*b1@b1+(k2+np.abs(j-N-1))**2*b2@b2

    return 0.05*H_3d.reshape((2*N-1)*(2*N-1)*(2*N-1),(2*N-1)*(2*N-1)*(2*N-1))

def HV3d(V_3df,N):
    VN=len(V_3df)
    
    V_3d=np.zeros([2*N-1,2*N-1,2*N-1,2*N-1,2*N-1,2*N-1])
    for i in range(2*N-1):
        for j in range(2*N-1):
            for k in range(2*N-1):
                for l in range(2*N-1):
                    for m in range(2*N-1):
                        for n in range(2*N-1):
                            if (np.abs(i-l)>VN-1)or(np.abs(j-m)>VN-1)or(np.abs(k-n)>VN-1):
                                continue
                            else:
                                V_3d[i,j,k,l,m,n]=V_3df[np.abs(i-l),np.abs(j-m),np.abs(k-n)]
    return V_3d.reshape((2*N-1)*(2*N-1)*(2*N-1),(2*N-1)*(2*N-1)*(2*N-1))

势能

def V3dcubic(N,a):
    epsilon=0.05*a
    x=np.linspace(0,1,N)
    V=np.zeros([N,N,N])
    for i in range(N):
        for j in range(N):
            for k in range(N):
                V[i,j,k]=1/(epsilon+np.sqrt(x[i]**2+x[j]**2+x[k]**2))+1/(epsilon+np.sqrt((x[i]-a)**2+x[j]**2+x[k]**2))+1/(epsilon+np.sqrt((x[i]-a)**2+(x[j]-a)**2+x[k]**2))+1/(epsilon+np.sqrt((x[i]-a)**2+(x[j]-a)**2+(x[k]-a)**2))+1/(epsilon+np.sqrt((x[i]-a)**2+(x[j])**2+(x[k]-a)**2))+1/(epsilon+np.sqrt((x[i])**2+(x[j])**2+(x[k]-a)**2))+1/(epsilon+np.sqrt((x[i])**2+(x[j]-a)**2+(x[k])**2))+1/(epsilon+np.sqrt((x[i])**2+(x[j]-a)**2+(x[k]-a)**2))
    return V

由于势函数关于中心对称，故只需要傅里叶展开系数的实部

a=1
n=15
l=2*n+1
V3d=V3dcubic(l,a)
V_3df=scifft.fftn(V3d)/(l**3)
V_3df=np.round(V_3df[:n,:n,:n],5).real

沿着高对称性路径扫描

a=1
n=50
ln=2*n+1
V3d=V3dcubic(l,a)
V_3df=scifft.fftn(V3d)/(l**3)
V_3df=V_3df[:20,:20,:20]
N=5
E1=[]
E2=[]
E3=[]
E4=[]
e1=np.array([np.pi,0,0])
e2=np.array([0,np.pi,0])
e3=np.array([0,0,np.pi])
lam=np.arange(0,ln,3)/ln
#Gamma to X
for i in range(len(lam)):
    k=lam[i]*e1
    w,v=lin.eigh(H_3D(k[0],k[1],k[2],a,N)+HV3d(V_3df.real,N))
    E1.append(w)
#X to M
for i in range(len(lam)):
    k=lam[i]*e2+e1
    w,v=lin.eigh(H_3D(k[0],k[1],k[2],a,N)+HV3d(V_3df.real,N))
    E2.append(w)
#M to R
for i in range(len(lam)):
    k=lam[i]*e3+e2+e1
    w,v=lin.eigh(H_3D(k[0],k[1],k[2],a,N)+HV3d(V_3df.real,N))
    E3.append(w)
#R to Gamma
for i in range(len(lam)):
    k=lam[::-1][i]*(e3+e2+e1)
    w,v=lin.eigh(H_3D(k[0],k[1],k[2],a,N)+HV3d(V_3df.real,N))
    E4.append(w)