学术源码, 模型和能带

石墨烯哈密顿量与能带图(附Python代码)

石墨烯示意图:

该图片来源于Hideo Aoki Mildred S. Dresselhaus的“Physics of Graphene”

只考虑最近邻跃迁t=1。以下是画二维石墨烯能带的代码:

"""
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/408
"""

import numpy as np
import matplotlib.pyplot as plt
from math import *
import cmath 
import functools  


def hamiltonian(k1, k2, M, t1, a=1/sqrt(3)):  # graphene哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    # 初始化为零矩阵
    h0 = np.zeros((2, 2), dtype=complex)
    h1 = np.zeros((2, 2), dtype=complex)

    # 质量项(mass term),用于打开带隙
    h0[0, 0] = M
    h0[1, 1] = -M

    # 最近邻项
    h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
    h1[0, 1] = h1[1, 0].conj()

    # # 最近邻项也可写成这种形式
    # h1[1, 0] = t1+t1*cmath.exp(1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)+t1*cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)
    # h1[0, 1] = h1[1, 0].conj()

    matrix = h0 + h1
    return matrix


def main():
    hamiltonian0 = functools.partial(hamiltonian, M=0, t1=1, a=1/sqrt(3))  # 使用偏函数,固定一些参数
    k1 = np.linspace(-2*pi, 2*pi, 500)
    k2 = np.linspace(-2*pi, 2*pi, 500)
    plot_bands_two_dimension(k1, k2, hamiltonian0)


def plot_bands_two_dimension(k1, k2, hamiltonian): 
    from matplotlib import cm
    dim = hamiltonian(0, 0).shape[0]
    dim1 = k1.shape[0]
    dim2 = k2.shape[0]
    eigenvalue_k = np.zeros((dim2, dim1, dim))
    i0 = 0
    for k10 in k1:
        j0 = 0
        for k20 in k2:
            matrix0 = hamiltonian(k10, k20)
            eigenvalue, eigenvector = np.linalg.eig(matrix0)
            eigenvalue_k[j0, i0, :] = np.sort(np.real(eigenvalue[:]))
            j0 += 1
        i0 += 1
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    k1, k2 = np.meshgrid(k1, k2)
    for dim0 in range(dim):
        ax.plot_surface(k1, k2, eigenvalue_k[:, :, dim0], rcount=200, ccount=200, cmap=cm.coolwarm, linewidth=0, antialiased=False)  
    plt.show()


if __name__ == '__main__':
    main()

能带图为:

下面考虑准一维的情况( 这里只考虑Zigzag边的石墨烯条带)。

在我自己写的代码里,原子的编号方式如下图所示。按既定规则,把元胞内的跃迁矩阵h00写和元胞间的跃迁矩阵h01都写出来,包括所有的最近邻跃迁, 然后整体傅里叶变换即可。关于离散傅里叶变换可参考博文:离散格子的傅里叶变换和反傅里叶变换 。这里是以元胞为单元(红色虚线圈了两个元胞作为例子),在一个方向做傅里叶变换转到倒空间(k空间),求本征值后就可以找到k与能量E的关系,即能带。

此图像的alt属性为空;文件名为image-16.png

以下是代码:

"""
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/408
"""

import numpy as np
import matplotlib.pyplot as plt
from math import *  
import cmath  
import functools 


def hamiltonian(k, N, M, t1):  # graphene哈密顿量(N是条带的宽度参数)
    # 初始化为零矩阵
    h00 = np.zeros((4*N, 4*N), dtype=complex)
    h01 = np.zeros((4*N, 4*N), dtype=complex)

    # 原胞内的跃迁h00
    for i in range(N):
        h00[i*4+0, i*4+0] = M
        h00[i*4+1, i*4+1] = -M
        h00[i*4+2, i*4+2] = M
        h00[i*4+3, i*4+3] = -M

        # 最近邻
        h00[i*4+0, i*4+1] = t1
        h00[i*4+1, i*4+0] = t1
        h00[i*4+1, i*4+2] = t1
        h00[i*4+2, i*4+1] = t1
        h00[i*4+2, i*4+3] = t1
        h00[i*4+3, i*4+2] = t1
    for i in range(N-1):
        # 最近邻
        h00[i*4+3, (i+1)*4+0] = t1
        h00[(i+1)*4+0, i*4+3] = t1

    # 原胞间的跃迁h01
    for i in range(N):
        # 最近邻
        h01[i*4+1, i*4+0] = t1
        h01[i*4+2, i*4+3] = t1

    matrix = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)
    return matrix


def main():
    hamiltonian0 = functools.partial(hamiltonian, N=40, M=0, t1=1)
    k = np.linspace(-pi, pi, 300)
    plot_bands_one_dimension(k, hamiltonian0)


def plot_bands_one_dimension(k, hamiltonian):
    dim = hamiltonian(0).shape[0]
    dim_k = k.shape[0]
    eigenvalue_k = np.zeros((dim_k, dim))
    i0 = 0
    for k0 in k:
        matrix0 = hamiltonian(k0)
        eigenvalue, eigenvector = np.linalg.eig(matrix0)
        eigenvalue_k[i0, :] = np.sort(np.real(eigenvalue[:]))
        i0 += 1
    for dim0 in range(dim):
        plt.plot(k, eigenvalue_k[:, dim0], '-k')
    plt.show()


if __name__ == '__main__':
    main()

这是画出的能带图:

【未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com

7,194 次浏览

22 thoughts on “石墨烯哈密顿量与能带图(附Python代码)”

  1. clear
    close all
    clc
    N=4;
    t=1;
    M=0;
    h00=zeros(4*N);
    h01=zeros(4*N);
    for i=0:N-1
        %包内跃迁
        h00(i*4+1,i*4+1)=M;
        h00(i*4+2,i*4+2)=-M;
        h00(i*4+3,i*4+3)=M;
        h00(i*4+4,i*4+4)=-M;
        %最近邻
        h00(i*4+1,i*4+2)=t;
        h00(i*4+2,i*4+1)=t;
        h00(i*4+2,i*4+3)=t;
        h00(i*4+3,i*4+2)=t;
        h00(i*4+3,i*4+4)=t;
        h00(i*4+4,i*4+3)=t;
        
        for i=0:N-2
            h00(i*4+4,(i+1)*4+1)=t;
            h00((i+1)*4+1,i*4+4)=t;   
        end
        for i=0:N-1%-------%原胞间的跃迁 最近邻
            h01(i*4+2,i*4+1)=t;
            h01(i*4+3,i*4+4)=t;
            
        end    
    end
    kx=linspace(-pi,pi,100);
    for i=1:100
       hh=h00+h01*exp(1i*kx(i))+h01*exp(-1i*kx(i));
       e(:,i)=sort(eig(hh));
    end
    plot(kx,e,'b')

    1. 不知道是哪里出了问题,运行的结果在接近-pi和pi的时候图像不太对

      1. h01少了一撇。这句改成:hh=h00+h01*exp(1i*kx(i))+h01’*exp(-1i*kx(i));
        在Matlab中,一撇等于转置加共轭,即Python中的.transpose().conj()。

  2. 请问下博主,在zigzag模型中,哈密顿量中如果有kx,ky两个,那应该怎么写

    1. zigzag石墨烯条带是准一维体系,只有一个方向有k,为无限长。还有一个方向是有限宽度,因此不可能同时有kx和ky。

  3. 想问一下博主,为什么第一个程序里面不用写 原胞间的跃迁项h01矩阵,而第二个程序里面需要构造h01*cmath.exp(1j*k) 呢?

    1. 第一个程序的哈密顿量已经是倒空间的形式了,因此不需要傅里叶变换,或者说是已经做完傅里叶变换。第二个程序是从实空间出发,写出条带的元胞,然后是按准一维条带进行傅里叶变换,得到倒空间的哈密顿量。参考:离散格子的傅里叶变换和反傅里叶变换

  4. 博主你好,我问一个简单的问题,为什么zigzag边界的计算和你之前那个准一维的计算过程差不多?这个结构胞间和胞内的跃迁都是有斜线的,不需要考虑矢量问题吗?

    1. 如果是以原子进行傅里叶变换,傅里叶变换产生的相位是需要考虑在x方向的投影。但如果是以元胞进行傅里叶变换,类似于一维原子链,其他为内部自由度。折线结构已经体现在hopping上了。

      1. 所以在写hamiltonian 的h01[4 i +0, 4 i +1]这一项的时候,应该是: t(1+exp(i k A); 而不应该是: t(1+exp(i k A* \sqrt{3}/2);
        想问一下博主是这样么;A对应的应该是超胞的晶格常数,A=\sqrt{3} a;
        a是石墨烯单元的两原子间距;
        实际上在楼主代码里面A默认是1;而并非a默认是1;
        请问博主大哥哥,我的理解对么

        1. 跟同学讨论过了,hamiltonian 的h01[4 i +0, 4 i +1]这一项,应该是: t(1+exp(i k A); 而不应该是: t(1+exp(i k A* \sqrt{3}/2);

          1. 括号中没有1吧。代码中是这个:matrix = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)。

  5. 不好意思,我问一个很trivial的问题,您这里要计算zigzag的边缘态,但是您画的红色的条带里面难道不是armchair型的么?所谓原胞内原子数是4N是怎么理解呢?zigzag边缘不就考虑相邻两个原子之间的hopping就可以了么?

    1. 嗯,这里我图标明的不是太完整。示意图水平方向是无限长的,垂直方向是有限宽度的,红色圈出来的不是条带,是表示元胞,边缘是zigzag型的。

      因为4个原子重复一次,所以在写程序时,我就按4的倍数来写循环了。

      有连接的地方都要写上hopping。

  6. 请问,我类推到armchair纳米带的方法,出来的图像大相径庭呢?

    % 原胞内的跃迁h0
      for i=0:1:N-1
        h0(i*4+1,i*4+1) = M;
        h0(i*4+2,i*4+2) = -M;  %2*cos(kx)-1i*2*cos(ky)
        h0(i*4+3,i*4+3) = M;
        h0(i*4+4,i*4+4) = -M;
        
        % 最近邻
        h0(i*4+1,i*4+2) = t1;
        h0(i*4+2,i*4+1) = t1;
        h0(i*4+1,i*4+1) = t1;
        h0(i*4+3,i*4+1) = t1;
        h0(i*4+2,i*4+4) = t1;
        h0(i*4+4,i*4+2) = t1;
      end
        
      for i=0:1:N-2
        %最近邻
        h0(i*4+4,(i+1)*4+2) = t1;
        h0((i+1)*4+2,i*4+4) = t1;
        h0(i*4+3,(i+1)*4+1) = t1;
        h0((i+1)*4+1,i*4+3) = t1;    % airchair原胞内,上下有两组最近邻
      end
      
      % 原胞间的跃迁h1
      for i=0:1:N-1
        h1(i*4+4,i*4+3) = t1;  % airchair原胞内间,左右有一组最近邻
      end
      H_matrix = h0+h1*exp(-1j*k)+conj(h1')*exp(1j*k); 

    所以请教博主。

    1. armchair的能带和zigzag的能带本身就是不一样的。我也不确定你写的对不对。你可以找armchair石墨烯的一些文章或者综述,应该有挺多的,看是否能重复对应的图。

  7. 博主可以讲一下Armchair型边界的思路吗?应该选取什么样形状做元胞。

    1. 画出来后,找条带的最小重复单元。怎么选取都可以,只要是最小单元,同一个体系计算出的能带应该都是一样的。

  8. 想问下博主,这个质量项代表什么呀,科研新手,有点不太明白,谢谢啦
    # 质量项(mass term), 用于打开带隙
    h0[0, 0] = M
    h0[1, 1] = -M

    1. 其实就是交错势(staggered potential),一般写成$m*\sigma_z$。至于为什么称为质量项(mass term),我也不是很清楚,但文献中都是这么使用的。我猜可能是因为这里的m会使得能带打开带隙,而对于有带隙的体系,大多数时候电子是存在有效质量的,所以称为质量项,不是很严谨。不一定是这样的,有需要可以找找看有没有文献中提到这个,我目前是没看到。

发表评论

您的电子邮箱地址不会被公开。