python wannier90 基于wannier90的*_hr.dat文件选取截断hopping绘制能带图

我们知道wannier90可以根据选取TMDs的轨道信息生成详细的hopping energy *_hr.dat文件，选取所有的hopping绘制起来的时候比较简单，但是我们发现取几圈的近似hopping也可以将band表示出来，类似的思想有Pybinding的三带近似（DOI: 10.1103/PhysRevB.88.085433）和Fang的TMDs的紧束缚模型半经验参数近似（DOI: 10.1103/PhysRevB.92.205108）等等。

我们首先要知道wannier90.dat文件里的数据格式，我们需要的是类似于

-12   -6    0    1    1    0.000020   -0.000000
-12   -6    0    2    1    0.000000   -0.000000
-12   -6    0    3    1    0.000000    0.000000
-12   -6    0    4    1    0.000004    0.000000
-12   -6    0    5    1    0.000000    0.000000

这样的数据，这里前面三个代表directions[x,y,z] 中间两个代表[fromAtom, toAtom]也就是Rij，后面两个值代表了hopping，代表复数 a+bi 的[a, b]。

这里的画圈有点类似，也就是根据选择当前原子和目标原子的距离来简要进行阶段，例如下图展示了四层原子（包括[0，0]），可能第三圈的hopping贡献远大于第四圈，所以截取到第三圈即可绘制band。

这里给出 $MX_2$ 部分示例，选取M的 $d_(z^2) d_(xy) d_(x^2-y^2)$ 和总的S轨道计算的一个wannier90 *_hr.dat文件，任意选择截取1-10（[0,0]就算了几乎是onsite energy）圈的一个代码，其实取到第5圈就已经比较近似于全部圈层的结果了：

#!/usr/bin python
# -*- encoding: utf-8 -*-
'''
@Author  :   Celeste Young
@File    :   基于fang2015提取MoS2轨道.py    
@Time    :   2023/4/3 14:47  
@E-mail  :   iamwxyoung@qq.com
@Tips    :   
'''import time
import matplotlib.pyplot as plt
import pandas as pd
# import pybinding as pb
import numpy as np
import functools
import threadingdef make_path(k0, k1, *ks, step=0.1):  # kpath k_points = [np.atleast_1d(k) for k in (k0, k1) + ks]if not all(k.shape == k_points[0].shape for k in k_points[:1]):raise RuntimeError("All k-points must have the same shape")k_paths = []point_indices = [0]for k_start, k_end in zip(k_points[:-1], k_points[1:]):num_steps = int(np.linalg.norm(k_end - k_start) // step)# k_path.shape == num_steps, k_space_dimensionsk_path = np.array([np.linspace(s, e, num_steps, endpoint=False)for s, e in zip(k_start, k_end)]).Tk_paths.append(k_path)point_indices.append(point_indices[-1] + num_steps)k_paths.append(k_points[-1])return np.vstack(k_paths)  # , point_indices)def isinDirect(v_, hvts):  # 判断是否在所选的directions里，上图的n层for vs in hvts:if ''.join(v_) == ''.join(list(map(str, vs))):return Truereturn Falsedef read_dat(*args, **kwargs):with open("data/wannier90_hr_MoS2.dat", "r") as f:lines = f.readlines()rij = [[[] for col in range(nb)] for row in range(nb)]tij = [[[] for col in range(nb)] for row in range(nb)]for line in lines:ll = line.split()if isinDirect(ll[:3], kwargs['hvts']):x, y, z, frAt, toAt = list(map(int, ll[:5]))t_real, t_image = list(map(float, ll[5:]))rij[frAt - 1][toAt - 1].append([x, y, z])tij[frAt - 1][toAt - 1].append([t_real + 1j * t_image])return rij, tijdef plot_rec(*args, **kwargs):  # 绘制布里渊区的积分路径fig = plt.figure(figsize=(5, 5))ax = plt.axes()ax.arrow(0, 0, a1[0], a1[1], length_includes_head=False, head_width=0.05, fc='b', ec='k')ax.arrow(0, 0, a2[0], a2[1], length_includes_head=False, head_width=0.05, fc='b', ec='k')ax.arrow(0, 0, b1[0], b1[1], length_includes_head=False, head_width=0.05, fc='r', ec='red')ax.arrow(0, 0, b2[0], b2[1], length_includes_head=False, head_width=0.05, fc='r', ec='red')ax.plot([0, Middle[0], K1[0], 0], [0, Middle[1], K1[1], 0])ax.scatter([0, Middle[0], K1[0], 0], [0, Middle[1], K1[1], 0])ax.set_xlim(-1, 4)ax.set_ylim(-2, 4)ax.grid()plt.show()def phase(*args, **kwargs):rij = kwargs['hr']wk = kwargs['wk']R_vec = 0for ii in range(3):R_vec += rij[ii] * basis_vector[ii]inner_product = np.dot(R_vec, wk)return np.exp(1j * inner_product)def Hamham(wak, tij, rij):h = np.zeros((nb, nb), dtype=complex)for ii in range(nb):for jj in range(nb):for vij in range(len(rij[ii][jj])):hr = rij[ii][jj][vij]h[ii, jj] = h[ii, jj] + tij[ii][jj][vij][0] * phase(hr=hr, wk=wak)return hdef get_sphere(*args, **kwargs):  #获取需要的层的directions，传入参数n=[1,11)sp_ = args[0]spDict = {}# sp = 0sp0 = [[0, 0, 0]]spDict['sp0'] = sp0# sp = 1 >> 3*2sp1 = [[1, 0, 0], [1, 1, 0], [0, 1, 0]]sp1 = sp1 + [list(-np.array(sp1[i])) for i in range(len(sp1))]spDict['sp1'] = sp1# sp = 2 >> 3*2sp2 = [[2, 1, 0], [1, 2, 0], [-1, 1, 0]]sp2 = sp2 + [list(-np.array(sp2[i])) for i in range(len(sp2))]spDict['sp2'] = sp2# sp = 3 >> 3*2sp3 = [[2, 0, 0], [2, 2, 0], [0, 2, 0]]sp3 = sp3 + [list(-np.array(sp3[i])) for i in range(len(sp3))]spDict['sp3'] = sp3# sp = 4 >> 6*2sp4 = [[3, 1, 0], [3, 2, 0], [2, 3, 0], [1, 3, 0], [-1, 2, 0], [-2, 1, 0]]sp4 = sp4 + [list(-np.array(sp4[i])) for i in range(len(sp4))]spDict['sp4'] = sp4# sp = 5 >> 3*2sp5 = [[3, 0, 0], [3, 3, 0], [0, 3, 0]]sp5 = sp5 + [list(-np.array(sp5[i])) for i in range(len(sp5))]spDict['sp5'] = sp5# sp = 6 >> 3*2sp6 = [[4, 2, 0], [2, 4, 0], [-2, 2, 0]]sp6 = sp6 + [list(-np.array(sp6[i])) for i in range(len(sp6))]spDict['sp6'] = sp6# sp = 7 >> 3*2sp7 = [[4, 1, 0], [4, 3, 0], [3, 4, 0], [1, 4, 0], [-1, 3, 0], [-3, 1, 0]]sp7 = sp7 + [list(-np.array(sp7[i])) for i in range(len(sp7))]spDict['sp7'] = sp7# sp = 8 >> 3*2sp8 = [[4, 0, 0], [4, 4, 0], [0, 4, 0]]sp8 = sp8 + [list(-np.array(sp8[i])) for i in range(len(sp8))]spDict['sp8'] = sp8sp9 = [[5, 2, 0], [5, 3, 0], [2, 5, 0], [3, 5, 0], [-2, 3, 0], [-3, 2, 0]]sp9 = sp9 + [list(-np.array(sp9[i])) for i in range(len(sp9))]spDict['sp9'] = sp9sp10 = [[5, 1, 0], [5, 4, 0], [1, 5, 0], [4, 5, 0], [-1, 4, 0], [-4, 1, 0]]sp10 = sp10 + [list(-np.array(sp10[i])) for i in range(len(sp10))]spDict['sp10'] = sp10if sp_ == 0:return sp0tmpl = []for i in range(sp_):tmpl += [j for j in spDict['sp%d' % i]]return tmpldef selectAllOrbital():with open("data/wannier90_hr_MoS2.dat", "r") as f:lines = f.readlines()rij = [[[] for col in range(nb)] for row in range(nb)]tij = [[[] for col in range(nb)] for row in range(nb)]for line in lines:ll = line.split()if len(ll) == 7:x, y, z, frAt, toAt = list(map(int, ll[:5]))t_real, t_image = list(map(float, ll[5:]))rij[frAt - 1][toAt - 1].append([x, y, z])tij[frAt - 1][toAt - 1].append([t_real + 1j * t_image])hamilton = functools.partial(Hamham, tij=tij, rij=rij)idx = 0result = np.zeros([len(path), 11])for kxy in range(len(path)):k = np.r_[path[kxy], [0]]w, t = np.linalg.eig(hamilton(k))w = list(w)w.sort()result[idx, :] = np.real(w)  # 将本征值进行保存idx += 1xk = [0, rt3, rt3 + 1, rt3 + 3]kk = np.linspace(0, 4.7, num=len(path))  # 作为x轴，使其和本征值矩阵每一列的y的值个数相同plt.figure(figsize=(4, 5))plt.plot(kk, result, c="r")plt.xticks(xk, ["Γ", "M", "K", "Γ"])plt.ylabel("Energy(eV)", fontsize=14)plt.axvline(rt3, color='gray', linestyle='--')plt.axvline(rt3 + 1, color='gray', linestyle='--')plt.title('all orbital')plt.savefig('all orbital.png')print('Saved all orbital figure!')def selectOrbital(orbi):print('started %d/%d' % (orbi, 10))hvts = get_sphere(orbi)result = np.zeros([len(path), 11])  # 解的矩阵Rij, Tij = read_dat(hvts=hvts)hamilton = functools.partial(Hamham, tij=Tij, rij=Rij)idx = 0for kxy in range(len(path)):k = np.r_[path[kxy], [0]]w, t = np.linalg.eig(hamilton(k))w = list(w)w.sort()result[idx, :] = np.real(w)  # 将本征值进行保存idx += 1xk = [0, rt3, rt3 + 1, rt3 + 3]kk = np.linspace(0, 4.7, num=len(path))  # 作为x轴，使其和本征值矩阵每一列的y的值个数相同plt.figure(figsize=(4, 5))plt.plot(kk, result, c="r")plt.xticks(xk, ["Γ", "M", "K", "Γ"])plt.ylabel("Energy(eV)", fontsize=14)plt.axvline(rt3, color='gray', linestyle='--')plt.axvline(rt3 + 1, color='gray', linestyle='--')plt.title('%d orbital' % orbi)plt.savefig('%d orbital.png' % orbi)print('finish %d/%d' % (orbi, 10))if __name__ == '__main__':# start heretime1 = time.time()nb = 11  # number of bandsa = 3.160  # 3.18c = 12.29  # distance of layersdxx = 3.13  # distance of orbital X-Xdxm = 2.41  # distance of orbital X-M# constant checkedrt3 = 3 ** 0.5pi = np.pia1 = a * np.array([rt3 / 2, -1 / 2, 0])a2 = a * np.array([0, 1, 0])a3 = a * np.array([0, 0, 5])basis_vector = np.array([a1, a2, a3])b1 = 2 * pi / a * np.array([1 / rt3, 1])b2 = 4 * pi / a / rt3 * np.array([1, 0])Gamma = np.array([0, 0])Middle = 1 / 2 * b1K1 = 1 / 3 * (2 * b1 - b2)K2 = -1 / 3 * (2 * b1 - b2)# plot_rec()path = make_path(Gamma, Middle, K1, Gamma, step=0.01)# selectAllOrbital()  # 这里取消注释会绘制总的圈层的bandfor orb in range(1, 11):selectOrbital(orbi=orb)time2 = time.time()print('>> Finished, use time %d s' % (time2 - time1))

附上我的hr.dat的结果图: