首先解析推导出表面态开边界的哈密顿量，这一步也可以数值完成，但是需要写比较多的代码，不过有很多现成的工具可以使用
def ham(t, kx):
h01 = np.array([[0, 0], [t, 0]])
h00 = np.zeros((2, 2), dtype=complex)
h00[0, 1] = t + t * cmath.exp(-1j * kx)
h00[1, 0] = h00[0, 1].conjugate()
return h00, h01
def full_ham(t, n_site, kx):
# 大哈密顿矩阵
h01 = np.array([[0, 0], [t, 0]])
h00 = np.zeros((2, 2), dtype=complex)
h00[0, 1] = t + t * cmath.exp(-1j * kx)
h00[1, 0] = h00[0, 1].conjugate()
h00 = np.kron(np.eye(n_site), h00)
h01 = np.kron(np.eye(n_site, k=1), h01) + np.kron(np.eye(n_site, k=-1), h01.T.conjugate())
h1 = h00 + h01
return h1一维情况比较容易求解，这里重点看表面态的迭代解法：
def Green_Iteration(dim, E, h00, h01):
n_iterm = 500
err = 1e-10
epsilon_i = epsilon_s = h00
alpha_i = h01
beta_i = h01.T.conjugate()
for k in range(n_iterm):
mat = np.linalg.inv(E * np.eye(dim) - epsilon_i)
epsilon_s = epsilon_s + np.dot(np.dot(alpha_i, mat), beta_i)
epsilon_i = epsilon_i + np.dot(np.dot(alpha_i, mat), beta_i) + np.dot(np.dot(beta_i, mat), alpha_i)
alpha_i = np.dot(np.dot(alpha_i, mat), alpha_i)
beta_i = np.dot(np.dot(beta_i, mat), beta_i)
if np.sum(abs(alpha_i)) < err:
break
GL = np.trace(np.linalg.inv(E * np.eye(dim) - epsilon_s))
GB = np.trace(np.linalg.inv(E * np.eye(dim) - epsilon_i))
return GL, GB计算表明格林函数：
def Green_Surface():
star = time.perf_counter()
tij = 1
eta = 1e-4
n_step = 100
E_min = -0.5
E_max = 0.5
nE = 100
kx_list = np.linspace(-np.pi, np.pi, n_step)
E_list = np.linspace(E_min, E_max, nE)
GS = np.zeros((nE, n_step), dtype=complex)
GB = np.zeros((nE, n_step), dtype=complex)
for i in range(nE):
E = E_list[i] + 1j * eta
for j in range(nE):
h00, h01 = ham(tij, kx_list[j])
GS[i, j], GB[i, j] = Green_Iteration(2, E, h00, h01)
A = -1 * GS.imag / np.pi
AB = -1 * GB.imag / np.pi
end = time.perf_counter()
print("CPU time:", end - star)
Draw_Surface(kx_list, E_list, A, AB)画图：
def Draw_Surface(x, y, A, AB):
plt.subplot(121)
plt.pcolormesh(x, y, A, shading='auto', cmap="magma")
plt.colorbar()
plt.xlabel("Kx")
plt.ylabel("Energy")
plt.title("Surface DOS")
plt.subplot(122)
plt.pcolormesh(x, y, AB, shading='auto', cmap="magma")
plt.colorbar()
plt.yticks([])
plt.xlabel("Kx")
plt.title("Bulk DOS")
plt.savefig("Zigzag_Surface.jpg", dpi=300)
plt.show()
