模型和能带, 学术

Kagome晶格中的平带(附Mathematica、Python代码)

Kagome lattice(笼目结构晶格)示意图[1]:

Kagome晶格在倒空间的哈密顿量[1]:

H=\sum\limits_{ \bm{k} }\Phi_{k}^{\dagger}H_{ \bm{k} }\Phi_{ \bm{k} }

其中,\Phi_{k}=(c_{1\bm{k}},  c_{2\bm{k}},  c_{3\bm{k}}  )^{T},以及

H_{\bm{k}}=-t\begin{pmatrix}0   & 2 \cos (\vec{k_1}\cdot\vec{a_1})   & 2\cos (\vec{k_2}\cdot\vec{a_2})\\2\cos (\vec{k_1}\cdot\vec{a_1})  &  0  &  2\cos (\vec{k_3}\cdot\vec{a_3})  \\ 2\cos (\vec{k_2}\cdot\vec{a_2})   & \ 2\cos (\vec{k_3}\cdot\vec{a_3})    &  0\end{pmatrix}

Mathematica符号计算:

Clear["`*"]
k1a1 = kx;
k2a2 = kx/2 + ky*Sqrt[3]/2;
k3a3 = -kx/2 + ky*Sqrt[3]/2;
H = -2*t*({{0, Cos[k1a1], Cos[k2a2]}, {Cos[k1a1], 0, Cos[k3a3]}, {Cos[k2a2], Cos[k3a3], 0}});
MatrixForm[H]
eigenvalue = MatrixForm[Eigenvalues[H]]

计算结果:

Python数值计算(可以用到开源软件包Guan: https://py.guanjihuan.com):

"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/16730
"""

import numpy as np
from math import *

def hamiltonian(kx, ky):  # kagome lattice
    k1_dot_a1 = kx
    k2_dot_a2 = kx/2+ky*sqrt(3)/2
    k3_dot_a3 = -kx/2+ky*sqrt(3)/2
    h = np.zeros((3, 3), dtype=complex)
    h[0, 1] = 2*cos(k1_dot_a1)
    h[0, 2] = 2*cos(k2_dot_a2)
    h[1, 2] = 2*cos(k3_dot_a3)
    h = h + h.transpose().conj()
    t = 1
    h = -t*h
    return h

def main():
    kx_array = np.linspace(-pi ,pi, 500)
    ky_array = np.linspace(-pi ,pi, 500)
    eigenvalue_array = calculate_eigenvalue_with_two_parameters(kx_array, ky_array, hamiltonian)
    plot_3d_surface(kx_array, ky_array, eigenvalue_array, xlabel='kx', ylabel='ky', zlabel='E', rcount=200, ccount=200)
    
    # import guan
    # eigenvalue_array = guan.calculate_eigenvalue_with_two_parameters(kx_array, ky_array, hamiltonian)
    # guan.plot_3d_surface(kx_array, ky_array, eigenvalue_array, xlabel='kx', ylabel='ky', zlabel='E', rcount=200, ccount=200)

def calculate_eigenvalue_with_two_parameters(x_array, y_array, hamiltonian_function, print_show=0, print_show_more=0):  
    dim_x = np.array(x_array).shape[0]
    dim_y = np.array(y_array).shape[0]
    if np.array(hamiltonian_function(0,0)).shape==():
        eigenvalue_array = np.zeros((dim_y, dim_x, 1))
        i0 = 0
        for y0 in y_array:
            j0 = 0
            for x0 in x_array:
                hamiltonian = hamiltonian_function(x0, y0)
                eigenvalue_array[i0, j0, 0] = np.real(hamiltonian)
                j0 += 1
            i0 += 1
    else:
        dim = np.array(hamiltonian_function(0, 0)).shape[0]
        eigenvalue_array = np.zeros((dim_y, dim_x, dim))
        i0 = 0
        for y0 in y_array:
            j0 = 0
            if print_show==1:
                print(y0)
            for x0 in x_array:
                if print_show_more==1:
                    print(x0)
                hamiltonian = hamiltonian_function(x0, y0)
                eigenvalue, eigenvector = np.linalg.eigh(hamiltonian)
                eigenvalue_array[i0, j0, :] = eigenvalue
                j0 += 1
            i0 += 1
    return eigenvalue_array

def plot_3d_surface(x_array, y_array, matrix, xlabel='x', ylabel='y', zlabel='z', title='', fontsize=20, labelsize=15, show=1, save=0, filename='a', file_format='.jpg', dpi=300, z_min=None, z_max=None, rcount=100, ccount=100): 
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator
    matrix = np.array(matrix)
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    plt.subplots_adjust(bottom=0.1, right=0.65) 
    x_array, y_array = np.meshgrid(x_array, y_array)
    if len(matrix.shape) == 2:
        surf = ax.plot_surface(x_array, y_array, matrix, rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    elif len(matrix.shape) == 3:
        for i0 in range(matrix.shape[2]):
            surf = ax.plot_surface(x_array, y_array, matrix[:,:,i0], rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman')
    ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.set_zlabel(zlabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.zaxis.set_major_locator(LinearLocator(5)) 
    ax.zaxis.set_major_formatter('{x:.2f}')  
    if z_min!=None or z_max!=None:
        if z_min==None:
            z_min=matrix.min()
        if z_max==None:
            z_max=matrix.max()
        ax.set_zlim(z_min, z_max)
    ax.tick_params(labelsize=labelsize) 
    labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels()
    [label.set_fontname('Times New Roman') for label in labels] 
    cax = plt.axes([0.8, 0.1, 0.05, 0.8]) 
    cbar = fig.colorbar(surf, cax=cax)  
    cbar.ax.tick_params(labelsize=labelsize)
    for l in cbar.ax.yaxis.get_ticklabels():
        l.set_family('Times New Roman')
    if save == 1:
        plt.savefig(filename+file_format, dpi=dpi) 
    if show == 1:
        plt.show()
    plt.close('all')

if __name__ == '__main__':
    main()

计算结果:

参考资料:

[1] Topological insulator on the kagome lattice

附:

6,273 次浏览

【说明:本站主要是个人的一些笔记和代码分享,内容可能会不定期修改。为了使全网显示的始终是最新版本,这里的文章未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com

13 thoughts on “Kagome晶格中的平带(附Mathematica、Python代码)”

  1. 博主,你好!请问你有做过磁子相关内容吗?在写出磁子哈密顿量后,根据自旋算符的表达式,也可以按“跃迁能”、“在位能”来构造哈密顿量矩阵吗

    1. 我没做过,但听说是可以写成类似的本征方程的形式。矩阵中每一项代表什么物理意义,具体你可以看看文献

  2. 博主你好,这个基矢的数目是等于原胞内原子个数吗?然后c如果是湮灭算符的话,c1k就是空间点1的湮灭算符吗

      1. 好的,谢谢博主。那这个序号的选取因人而异,写出来的哈密顿量可能不同,这是没有影响的吧

        1. 嗯,哈密顿量形式不一样,相当于取了不同的基矢,但物理结果是不影响的

            1. 如果选取不同的基的话,那不同的哈密顿量区别是不是在于他们可能有相同的矩阵元,但是在矩阵中的位置和相位不同

              1. 嗯,如果只是调整基矢顺序,那么只是矩阵元位置不一样。如果是其他基矢,例如平面波等非紧束缚基,那么矩阵会完全不同。

发表评论

您的邮箱地址不会被公开。 必填项已用 * 标注

Captcha Code