我们知道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。
这里给出部分示例，选取M的
和总的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 threading
def 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_dimensions
k_path = np.array([np.linspace(s, e, num_steps, endpoint=False)
for s, e in zip(k_start, k_end)]).T
k_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 True
return False
def 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, tij
def 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 = 0
for 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 h
def get_sphere(*args, **kwargs): #获取需要的层的directions，传入参数n=[1,11)
sp_ = args[0]
spDict = {}
# sp = 0
sp0 = [[0, 0, 0]]
spDict['sp0'] = sp0
# sp = 1 >> 3*2
sp1 = [[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*2
sp2 = [[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*2
sp3 = [[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*2
sp4 = [[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*2
sp5 = [[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*2
sp6 = [[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*2
sp7 = [[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*2
sp8 = [[4, 0, 0], [4, 4, 0], [0, 4, 0]]
sp8 = sp8 + [list(-np.array(sp8[i])) for i in range(len(sp8))]
spDict['sp8'] = sp8
sp9 = [[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'] = sp9
sp10 = [[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'] = sp10
if sp_ == 0:
return sp0
tmpl = []
for i in range(sp_):
tmpl += [j for j in spDict['sp%d' % i]]
return tmpl
def 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 = 0
result = 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 += 1
xk = [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 = 0
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 += 1
xk = [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 here
time1 = time.time()
nb = 11 # number of bands
a = 3.160 # 3.18
c = 12.29 # distance of layers
dxx = 3.13 # distance of orbital X-X
dxm = 2.41 # distance of orbital X-M
# constant checked
rt3 = 3 ** 0.5
pi = np.pi
a1 = 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 * b1
K1 = 1 / 3 * (2 * b1 - b2)
K2 = -1 / 3 * (2 * b1 - b2)
# plot_rec()
path = make_path(Gamma, Middle, K1, Gamma, step=0.01)
# selectAllOrbital() # 这里取消注释会绘制总的圈层的band
for orb in range(1, 11):
selectOrbital(orbi=orb)
time2 = time.time()
print('>> Finished, use time %d s' % (time2 - time1))
附上我的hr.dat的结果图:
