Python绘制Kagome晶格高对称方向能带结构遇问题求排查

我尝试用Python绘制Kagome晶格沿Γ-K-M-Γ高对称方向的能带结构，采用与石墨烯相同的晶格矢量和倒格矢参数，但结果不符合预期，推测高对称点位置存在问题。不过使用pythtb包以相同参数运行可得到正确结果，以下是我的代码、pythtb参考代码及对应结果，希望有人帮忙排查代码错误。

Kagome材料包含3个子晶格(A,B,C)，相关的晶格矢量、倒格矢、高对称点位置及子晶格连接矢量示意图已提供。

我的代码 #

import sympy as sp
import numpy as np
kx, ky = sp.symbols('kx ky')

def Pab(kx, ky):
    return 2*sp.cos( ((sp.sqrt(3) * kx) / 2) + ky / 2)

def Pbc(kx, ky):
    return 2*sp.cos(ky)

def Pac(kx, ky):
    return 2* sp.cos(((sp.sqrt(3) * kx) / 2 - ky / 2))

def CPab(kx, ky):
    return 2*sp.cos( ((sp.sqrt(3) * kx) / 2) + ky / 2)

def CPbc(kx, ky):
    return 2*sp.cos(ky)

def CPac(kx, ky):
    return 2* sp.cos(((sp.sqrt(3) * kx) / 2 - ky / 2))
Ea=0;Eb=0;Ec=0
tab=1;tbc=1;tca=1
H=sp.Matrix([[Ea , -tab* CPab(kx,ky), -tca* CPac(kx,ky)],
             [-tab* Pab(kx,ky), Eb, -tbc*CPbc(kx,ky)],[-tca*Pac(kx,ky), -tbc*Pbc(kx,ky),Ec]])
numeric_H = sp.lambdify((kx, ky), H, "numpy")

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
a=1
d1=(4*np.pi)/(3*np.sqrt(3)*a)

y1= (2*np.pi)/(3*np.sqrt(3)*a)
x1= (2*np.pi)/(3*a)
k1=[];k2=[];k3=[]
E1=[];E2=[];E3=[];En1=[];En2=[];En3=[];Ep1=[];Ep2=[];Ep3=[]
stepsizey = y1/500

#For Gamma to K path
kx=0
while kx<= x1:
    ky= kx/np.sqrt(3)
    k1.append(np.sqrt(kx**2+ky**2))
    numeric_matrix = numeric_H(kx, ky)
    eigenvalues = np.linalg.eigvals(numeric_matrix)
        
    E1.append(eigenvalues[0].real)
    En1.append(eigenvalues[1].real)
    Ep1.append(eigenvalues[2].real)
    kx+= stepsizey

#For K to M path
while ky >=0 :
    kx= x1
    k2.append(np.sqrt((kx-x1)**2+(ky-y1)**2))
    numeric_matrix = numeric_H(kx, ky)
    eigenvalues = np.linalg.eigvals(numeric_matrix)
        
    E2.append(eigenvalues[0].real)
    En2.append(eigenvalues[1].real)
    Ep2.append(eigenvalues[2].real)
    ky-= stepsizey
d2 =-( k2[0]-k2[len(k2)-1])


k2r =[]
for i in range (len(k2)):
    k2r.append(k2[i]+d1)
    i+=1

#For M to Gamma path

while kx>=0:
    ky=0
    k3.append(np.sqrt((x1-kx)**2))
    numeric_matrix = numeric_H(kx, ky)
    eigenvalues = np.linalg.eigvals(numeric_matrix)
        
    E3.append(eigenvalues[0].real)
    En3.append(eigenvalues[1].real)
    Ep3.append(eigenvalues[2].real)
    kx-= stepsizey
k3r=[]
for i in range (len(k3)):
    k3r.append(k3[i]+d1+d2)
#plotting

plt.scatter(k1,E1,color='r', s=1)
plt.scatter(k2r,E2,color='r', s=1)
plt.scatter(k3r,E3,color='r', s=1)
plt.scatter(k1,En1,color='r', s=1)
plt.scatter(k2r,En2,color='r', s=1)
plt.scatter(k3r,En3,color='r', s=1)
plt.scatter(k1,Ep1,color='r', s=1)
plt.scatter(k2r,Ep2,color='r', s=1)
plt.scatter(k3r,Ep3,color='r', s=1)
plt.axvline(x=d1, color='b', linestyle='--', label='K')
plt.axvline(x=0, color='b', linestyle='--', label='K')
plt.axvline(x=d1+d2, color='b', linestyle='--', label='K')
plt.axvline(x=k3r[len(k3r)-1], color='b', linestyle='--', label='K')
custom_tick_positions = [0, d1, d1+d2, k3r[len(k3r)-1]]
custom_labels = ['$\Gamma$', 'K', 'M','$\Gamma$' ]
plt.xticks(custom_tick_positions, custom_labels)
plt.show()

我的代码运行结果 #

绘制出的能带结构形态异常，不符合Kagome晶格的特征。

Pythtb参考代码 #

参考自《Topological properties of flat bands in generalized Kagome lattice materials》, Daniela Pinto Dias

#!/usr/bin/env python

from __future__ import print_function  # python3 style print

# Tight-binding 2D Kagome Lattice

from pythtb import *  # import TB model class
import matplotlib.pyplot as plt

# Set model parameters
E = 0.  # on-site energy
t = -1.0  # spin-independent first-neighbor hop


def set_model(E, t):
    # Set up Kagome model
    lat = [[np.sqrt(3.0) / 2.0, -0.5], [np.sqrt(3.0) / 2.0, 0.5]]
    orb = [[0., 0.], [1. / 2., 0.], [0., 1. / 2.]]
    model = tb_model(2, 2, lat, orb, nspin=2)
    model.set_onsite([E, E, E])

    # Spin-independent first-neighbor hops
    # Atom A hopping terms
    model.set_hop(t, 1, 0, [0, 0])
    model.set_hop(t, 1, 0, [1, 0])

    # Atom B hopping terms
    model.set_hop(t, 2, 1, [0, 0])
    model.set_hop(t, 2, 1, [-1, 1])

    # Atom C hopping terms
    model.set_hop(t, 0, 2, [0, 0])
    model.set_hop(t, 0, 2, [0, -1])
    return model, orb


# Print tight-binding model
my_model, orb = set_model(E, t)
my_model.display()

# Generate list of k-points following a segmented path in the BZ
# List of nodes (high-symmetry points) that will be connected
path = [ [0., 0.], [2. / 3., 1. / 3.], [0.5,0.5],[0.,0.]]
# Labels of the nodes
label = (r'$\Gamma$', r'$K$', r'$M$', r'$\Gamma$')

# Total number of interpolated k-points along the path
nk = 121

# Call function k_path to construct the actual path
(k_vec, k_dist, k_node) = my_model.k_path(path, nk)
# Inputs:
# path, nk: see above
# my_model: the pythtb model
# Outputs:
# k_vec: list of interpolated k-points
# k_dist: horizontal axis position of each k-point in the list
# k_node: horizontal axis position of each original node

print('---------------------------------------')
print('starting calculation')
print('---------------------------------------')
print('Calculating bands...')

# Obtain eigenvalues to be plotted
evals = my_model.solve_all(k_vec)

# Figure for band structure
fig, ax = plt.subplots()

# Specify horizontal axis details
# Set range of horizontal axis
ax.set_xlim(k_node[0], k_node[-1])
# Put tick marks and labels at node positions
ax.set_xticks(k_node)
ax.set_xticklabels(label, size=22)
ax.tick_params(axis="y", labelsize=22)
# Add vertical lines at node positions
for n in range(len(k_node)):
    ax.axvline(x=k_node[n], linewidth=0.5, color='k')
# Put labels
ax.set_xlabel("Path in k-space", size=22)
ax.set_ylabel("Band energy", size=22)

# Plot bands
for i in range(2 * len(orb)):
    ax.plot(k_dist, evals[i], color='r')
plt.tight_layout()
plt.show()

Pythtb代码运行结果 #

得到了正确的Kagome晶格能带结构，呈现出该体系特征性的平带与色散带。

内容的提问来源于stack exchange，提问作者Sinchan