学术源码, 拓扑不变量

Haldane模型中陈数的计算(附Python代码)

之前一篇关于Haldane模型的博文:Haldane模型哈密顿量与能带图(附Python代码)

在该博文中画出了Haldane模型条带(准一维)的能带,可以看到在体带隙中存在一个交叉的边缘态(可态密度验证属于边缘态)。该边缘态是受到拓扑保护的。本篇将给出Haldane模型中陈数的计算代码和计算结果。

方法仍然是采用这篇博文中的方法:陈数Chern number的计算(高效法,附Python/Matlab代码)

一、积分方法1:限制积分范围在六角晶格内

需要注意的是,在方格子中,布里渊区是方格子,因此kx和ky的积分范围都在[-pi, pi]。而在六角格子中,布里渊区是六角格子,在写程序的时候需要对六角格子的布里渊区积分。

六角格子实空间中坐标位置为:

对应倒空间中的坐标为:

以下是计算Haldane模型陈数的Python代码,积分时要用判断语句限制在布里渊区内。

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

import numpy as np
import matplotlib.pyplot as plt
from math import *   # 引入pi, cos等
import cmath
import time
import functools  # 使用偏函数functools.partial()


def hamiltonian(k1, k2, M, t1, t2, phi, a=1/sqrt(3)):  # Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    # 初始化为零矩阵
    h0 = np.zeros((2, 2))*(1+0j)   # 乘(1+0j)是为了把h0转为复数
    h1 = np.zeros((2, 2))*(1+0j)
    h2 = np.zeros((2, 2))*(1+0j)

    # 质量项(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()

    #次近邻项 # 对应陈数为-1
    h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    
    # # 次近邻项  # 对应陈数为1
    # h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    # h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))

    matrix = h0 + h1 + h2 + h2.transpose().conj()
    return matrix


def main():
    start_clock = time.perf_counter()
    delta = 0.005
    chern_number = 0  # 陈数初始化
    
    # 几个坐标中常出现的项
    a = 1/sqrt(3)
    aa1 = 4*sqrt(3)*pi/9/a
    aa2 = 2*sqrt(3)*pi/9/a
    bb1 = 2*pi/3/a

    hamiltonian0 = functools.partial(hamiltonian, M=2/3, t1=1, t2=1/3, phi=pi/4, a=a)   # 使用偏函数,固定一些参数

    for kx in np.arange(-aa1, aa1, delta):
        print(kx)
        for ky in np.arange(-bb1, bb1, delta):
            if (-aa2<=kx<=aa2) or (kx>aa2 and -(aa1-kx)*tan(pi/3)<=ky<=(aa1-kx)*tan(pi/3)) or (kx<-aa2 and  -(kx-(-aa1))*tan(pi/3)<=ky<=(kx-(-aa1))*tan(pi/3)):  # 限制在六角格子布里渊区内

                H = hamiltonian0(kx, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H)
                vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
            
                H_delta_kx = hamiltonian0(kx+delta, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
                vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]   # 略偏离kx的波函数

                H_delta_ky = hamiltonian0(kx, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
                vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离ky的波函数
                
                H_delta_kx_ky = hamiltonian0(kx+delta, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
                vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离kx和ky的波函数
                
                Ux = np.dot(np.conj(vector), vector_delta_kx)/abs(np.dot(np.conj(vector), vector_delta_kx))
                Uy = np.dot(np.conj(vector), vector_delta_ky)/abs(np.dot(np.conj(vector), vector_delta_ky))
                Ux_y = np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky))
                Uy_x = np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky))

                F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))

                # 陈数(chern number)
                chern_number = chern_number + F

    chern_number = chern_number/(2*pi*1j)
    print('Chern number = ', chern_number)
    end_clock = time.perf_counter()
    print('CPU执行时间(min)=', (end_clock-start_clock)/60)


if __name__ == '__main__':
    main()

计算结果为:

如果改变相位phi的符号,对应的计算结果为:

当选取的积分步长越小,计算结果越接近1。因为这里是用Python写的,所以计算效率不会太高;如果用Fortran写,计算所耗的时间会更短些。

量子自旋霍尔效应中比较简单的情况是这样的:通过引入自旋,让一个自旋对应的陈数为-1,一个自旋对应的陈数为1,这时整体的陈数为0,满足时间反演不变性,需要引入Z2不变量。可参考博文:Kane-Mele模型的哈密顿量和能带图(不考虑Rashba自旋轨道耦合的情况)

二、积分方法2:布里渊区平移成方形

示意图如上所示。代码为:

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

import numpy as np
import matplotlib.pyplot as plt
from math import *   # 引入pi, cos等
import cmath
import time
import functools  # 使用偏函数functools.partial()


def hamiltonian(k1, k2, M, t1, t2, phi, a=1/sqrt(3)):  # Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    # 初始化为零矩阵
    h0 = np.zeros((2, 2))*(1+0j)   # 乘(1+0j)是为了把h0转为复数
    h1 = np.zeros((2, 2))*(1+0j)
    h2 = np.zeros((2, 2))*(1+0j)

    # 质量项(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()

    #次近邻项 # 对应陈数为-1
    h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    
    # # 次近邻项  # 对应陈数为1
    # h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    # h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))

    matrix = h0 + h1 + h2 + h2.transpose().conj()
    return matrix


def main():
    start_clock = time.perf_counter()
    delta = 0.005
    chern_number = 0  # 陈数初始化
    
    # 常出现的项
    a = 1/sqrt(3)
    bb1 = 2*sqrt(3)*pi/3/a
    bb2 = 2*pi/3/a

    hamiltonian0 = functools.partial(hamiltonian, M=2/3, t1=1, t2=1/3, phi=pi/4, a=a)   # 使用偏函数,固定一些参数

    for kx in np.arange(0 , bb1, delta):
        print(kx)
        for ky in np.arange(0, 2*bb2, delta):
            H = hamiltonian0(kx, ky) 
            eigenvalue, eigenvector = np.linalg.eig(H)
            vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
        
            H_delta_kx = hamiltonian0(kx+delta, ky) 
            eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
            vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]   # 略偏离kx的波函数

            H_delta_ky = hamiltonian0(kx, ky+delta)  
            eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
            vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离ky的波函数
            
            H_delta_kx_ky = hamiltonian0(kx+delta, ky+delta)  
            eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
            vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离kx和ky的波函数
            
            Ux = np.dot(np.conj(vector), vector_delta_kx)/abs(np.dot(np.conj(vector), vector_delta_kx))
            Uy = np.dot(np.conj(vector), vector_delta_ky)/abs(np.dot(np.conj(vector), vector_delta_ky))
            Ux_y = np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky))
            Uy_x = np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky))

            F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))
            # 陈数(chern number)
            chern_number = chern_number + F

    chern_number = chern_number/(2*pi*1j)
    print('Chern number = ', chern_number)
    end_clock = time.perf_counter()
    print('CPU执行时间(min)=', (end_clock-start_clock)/60)


if __name__ == '__main__':
    main()

计算结果为:

三、积分方法3:坐标变换

此外,还有一种计算六角格子的方法,是直接对六角格子的基矢做积分。这个计算方法是由知乎网友“青春银河”提供:

“用基矢生成一组倒矢之后,两个倒矢的长度分别是k1与k2的周期,每一组kx与ky的值可以用每一组k1与k2的值分别乘以两个单位倒矢,再进行矢量加法得到。这么做的好处是不同的模型只需要改哈密顿量以及基矢,计算陈数部分全部都统一了。”

这种方法其实挺不好理解的,没有前面两种直观。编程时需要注意:for循环是以k1和k2来取点,变换到kx和ky,取点的个数保持不变,但步长发生了变化。做了变换后,积分的个数和步长不是很统一。

根据文献中给的基矢(位置刚好对调),我画了示意图如下:

其中,基矢和基矢方向上某个长度为

\overrightarrow{k_{10}}=(\frac{2\sqrt{3}\pi}{3a}, \frac{2\pi}{3a})

\overrightarrow{k_{20}}=(-\frac{2\sqrt{3}\pi}{3a}, \frac{2\pi}{3a})

k_{1}\overrightarrow{k_{10}}=(k_{1}\frac{2\sqrt{3}\pi}{3a}, k_{1}\frac{2\pi}{3a})

k_{2}\overrightarrow{k_{20}}=(-k_{2}\frac{2\sqrt{3}\pi}{3a}, k_{2}\frac{2\pi}{3a})

因此,k_xk_y用基矢方向上的长度来表示为

k_x=(k_{1}-k_{2})\frac{2\sqrt{3}\pi}{3a}

k_y=(k_{1}+k_{2})\frac{2\pi}{3a}

需要注意的是坐标变换后,步长也要跟着一起变化,要保证步长与积分个数的乘积刚好为布里渊区面积,即S(布里渊区)=∫dkx*dky。

坐标变换前:

积分步长为:delta

积分个数:1/delta

坐标变换后:

令 bb1 = 2*sqrt(3)*pi/3/a 以及 bb2 = 2*pi/3/a

积分步长为:delta_x = delta*bb1*constant;delta_y= delta*bb2*constant

积分个数仍然是:1/delta

如果不包括积分函数项时,这里直接积分的结果是:bb1*bb2*constant*constant

当constant=1时,直接积分的面积是红色长方形的四份之一,是黄色布里渊区面积的二分之一。因此constant要取\sqrt{2},或者一个constant取2,一个constant取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/5133
"""

import numpy as np
import matplotlib.pyplot as plt
from math import *   # 引入pi, cos等
import cmath
import time
import functools  # 使用偏函数functools.partial()


def hamiltonian(k1, k2, M, t1, t2, phi, a=1/sqrt(3)):  # Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    # 初始化为零矩阵
    h0 = np.zeros((2, 2))*(1+0j)   # 乘(1+0j)是为了把h0转为复数
    h1 = np.zeros((2, 2))*(1+0j)
    h2 = np.zeros((2, 2))*(1+0j)

    # 质量项(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()

    #次近邻项 # 对应陈数为-1
    h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    
    # # 次近邻项  # 对应陈数为1
    # h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    # h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))

    matrix = h0 + h1 + h2 + h2.transpose().conj()
    return matrix


def main():
    start_clock = time.perf_counter()
    delta = 0.002
    chern_number = 0  # 陈数初始化
    
    # 常出现的项
    a = 1/sqrt(3)
    bb1 = 2*sqrt(3)*pi/3/a
    bb2 = 2*pi/3/a

    hamiltonian0 = functools.partial(hamiltonian, M=2/3, t1=1, t2=1/3, phi=pi/4, a=a)   # 使用偏函数,固定一些参数

    for k1 in np.arange(0 , 1, delta):
        print(k1)
        for k2 in np.arange(0, 1, delta):
            # 坐标变换
            kx = (k1-k2)*bb1
            ky = (k1+k2)*bb2

            # 这里乘2或除以2是为了保证“步长与积分个数的乘积刚好为布里渊区面积”
            delta_x = delta*bb1*2    
            delta_y = delta*bb2*2/2

            H = hamiltonian0(kx, ky) 
            eigenvalue, eigenvector = np.linalg.eig(H)
            vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
        
            H_delta_kx = hamiltonian0(kx+delta_x, ky) 
            eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
            vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]   # 略偏离kx的波函数

            H_delta_ky = hamiltonian0(kx, ky+delta_y)  
            eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
            vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离ky的波函数
            
            H_delta_kx_ky = hamiltonian0(kx+delta_x, ky+delta_y)  
            eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
            vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离kx和ky的波函数
            
            Ux = np.dot(np.conj(vector), vector_delta_kx)/abs(np.dot(np.conj(vector), vector_delta_kx))
            Uy = np.dot(np.conj(vector), vector_delta_ky)/abs(np.dot(np.conj(vector), vector_delta_ky))
            Ux_y = np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky))
            Uy_x = np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky))

            F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))
            # 陈数(chern number)
            chern_number = chern_number + F

    chern_number = chern_number/(2*pi*1j)
    print('Chern number = ', chern_number)
    end_clock = time.perf_counter()
    print('CPU执行时间(min)=', (end_clock-start_clock)/60)


if __name__ == '__main__':
    main()

计算结果为:

这里计算耗时会更少一些,可能的原因有:(1)没有if判断,节省计算时间;(2)经过变换后,可能是数值上恰好更容易收敛。

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

3,724 次浏览

39 thoughts on “Haldane模型中陈数的计算(附Python代码)”

  1. 你好,请问使用坐标变换方法,“这么做的好处是不同的模型只需要改哈密顿量以及基矢,计算陈数部分全部都统一了”这句话怎么理解?如何改哈密顿量以及基失?

    1. 这句话不是我说的,大概可能是这样的意思:一般来说布里渊区是由基矢来决定,所以使用基矢量时,不管是什么模型,都可以用相同的方式积分,而不用适配于不同哈密顿量中不同形状的布里渊区。

      1. 对于三角格子而言,布里渊区都是正六边形的。是不是所有的三角格子都可以这样积分,不需要管哈密顿量的具体形式(不同配置的三角格子哈密顿量不同,但布里渊区相似)

  2. 您好,我想问一个如果考虑磁子之间的相互作用,怎么把哈密顿量对角化来画色散关系呢

    1. 不清楚,没考虑过,得看哈密顿量的具体形式吧。能带一般体现的是单体的性质,多体作用可能需要处理后才可以并入。

      1. 关老师,我现在计算石墨烯条带的陈数,因为考虑了边界效应,感觉计算陈数非常复杂,我没有理清楚思路,你能否给我介绍一下计算思路呢,万分感谢

        1. 条带是没法算陈数的,因为陈数是bulk的性质,需要有kx和ky,维度为偶数维。如果需算了解石墨烯条带的拓扑性质,可以直接算二维石墨烯的陈数。石墨烯不考虑其他项时,陈数应该为零。
          还有一种是算实空间的陈数,那是另外一回事了。

  3. 博主这里只考虑了能带非简并情况,简并情况的代码如下,哈密顿量是一个kagome格子的:

    %chen数高效法,布里渊区平移成矩形,简并陈数计算,算多条能带简并时,改的是函数M1、M2、u,有几个
    %band_index就有几个u函数
    function []=degen_chen_number_main2()
    clear;clc
    global delta band_index1 band_index2   band_index3   a
    tic
    delta = 0.01;
    chern_number =0;  %陈数初始化
    band_index1=1;%第几条能带索引?本征值最小的能带到本征值最大的能带,最大为矩阵阶数
    band_index2=2;%简并带索引
    band_index3=3;%简并带索引
    %算某一能带的陈数可能在index中设置
    % 几个坐标中常出现的项
    %matlab的dot函数算复数似乎不是很好
    a=1;
    aa1=pi/a;
    bb1=2*sqrt(3)*pi/(3*a);
    ky0=0:delta:bb1;
    kx0=0:delta:aa1;
    temp1=length(kx0);
    temp2=length(ky0);
    for i=1:temp1
        i/temp1
        for j=1:temp2
            
                kx=kx0(i);ky=ky0(j);
                % 陈数(chern number)
                chern_number=chern_number + F(kx,ky);
            
        end
    end
     
    chern_number=chern_number/(2*pi*1j);
    disp(['Chern number=',num2str(chern_number)])
    toc
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%
    
    function [oput]=F(kx,ky)
    global delta
    oput=log(U1(kx,ky)*U2(kx+delta,ky)*(U1(kx, ky +delta)^(-1))*(U2(kx, ky)^(-1)));
    end
    
    
    function [op1]=U1(kx,ky)
    op1=det(M1(kx,ky))/abs(det(M1(kx,ky)));
    end
    
    function [op2]=U2(kx,ky)
    op2=det(M2(kx,ky))/abs(det(M2(kx,ky)));
    end
    
    function [o1]=M1(kx,ky)
    global delta
    o1=zeros(3);
    o1(1,1)=u1(kx,ky)'*u1(kx + delta, ky);%u的排列和它的位置是有关的
    o1(1,2)=u1(kx,ky)'*u2(kx + delta, ky);
    o1(1,3)=u1(kx,ky)'*u3(kx + delta, ky);
    o1(2,1)=u2(kx,ky)'*u1(kx + delta, ky);
    o1(2,2)=u2(kx,ky)'*u2(kx + delta, ky);
    o1(2,3)=u2(kx,ky)'*u3(kx + delta, ky);
    o1(3,1)=u3(kx,ky)'*u1(kx + delta, ky);
    o1(3,2)=u3(kx,ky)'*u2(kx + delta, ky);
    o1(3,3)=u3(kx,ky)'*u3(kx + delta, ky);
    
    end
    
    
    function [o2]=M2(kx,ky)
    global delta
    o2=zeros(3);
    o2(1,1)=u1(kx,ky)'*u1(kx , ky+ delta);
    o2(1,2)=u1(kx,ky)'*u2(kx , ky+ delta);
    o2(1,3)=u1(kx,ky)'*u3(kx , ky+ delta);
    o2(2,1)=u2(kx,ky)'*u1(kx , ky+ delta);
    o2(2,2)=u2(kx,ky)'*u2(kx , ky+ delta);
    o2(2,3)=u2(kx,ky)'*u3(kx , ky+ delta);
    o2(3,1)=u3(kx,ky)'*u1(kx , ky+ delta);
    o2(3,2)=u3(kx,ky)'*u2(kx , ky+ delta);
    o2(3,3)=u3(kx,ky)'*u3(kx , ky+ delta);
    
    end
    
    
    function [output]=u1(kx,ky)
    global band_index1
    getH=hamiltonian0(kx,ky);
    [vector,value]=eig(getH);
    [~,index]=sort(real(diag(value)));
    output=vector(:,index(band_index1));%列向量
    end
    
    function [output]=u2(kx,ky)
    global band_index2
    getH=hamiltonian0(kx,ky);
    [vector,value]=eig(getH);
    [~,index]=sort(real(diag(value)));
    output=vector(:,index(band_index2));%列向量
    end
    
    function [output]=u3(kx,ky)
    global band_index3
    getH=hamiltonian0(kx,ky);
    [vector,value]=eig(getH);
    [~,index]=sort(real(diag(value)));
    output=vector(:,index(band_index3));%列向量
    end
    
    function [matrix]=hamiltonian0(kx, ky) %晶格常数a一般等于1
    global a
    phi=0;
    t2=0;
    t1=1;
     
    a1=[-0.5,-sqrt(3)/2]*a;
    a2=[1,0]*a;
    a3=[-0.5,sqrt(3)/2]*a;
    h1=zeros(3);
    h2=zeros(3);
    
            k=[kx,ky];
            ezhen=exp(1i*phi/3);
            efu=exp(-1i*phi/3);
            c1=cos(k*a1');
            c2=cos(k*a2');
            c3=cos(k*a3');
            
            c32=cos(k*(a3-a2)');
            c21=cos(k*(a2-a1)');
            c31=cos(k*(a3-a1)');
            h1(1,2)=c1*efu;
            h1(2,1)=c1*ezhen;
            h1(1,3)=c3*ezhen;
            h1(3,1)=c3*efu;
            h1(2,3)=c2*efu;
            h1(3,2)=c2*ezhen;
            
            h1=h1*(-2*t1);%最终,最近临
            
            h2(1,2)=c32;
            h2(1,3)=c21;
            h2(2,1)=c32;
            h2(2,3)=c31;
            h2(3,1)=c21;
            h2(3,2)=c31;
            h2=h2*(-2*t2);%最终,次近临
            
           matrix=h1+h2;
     
    end

    1. 嗯,感谢提供代码。问下:简并情况遇到的是什么问题?代码中需要修改和注意的是哪个地方?

          1. 如果能够区分开简并能带的波函数,然后就可以各自算陈数。在数值计算过程中可能会存在一些问题,比如每个k点得到的简并能带的波函数完全不同,需要额外的代码进行变换处理。

            此外,在数值上还有一种方法是加个极小的微扰使得能带退简并,这种方法比较简单,缺点是可能会略微破坏原本的对称性。如果对称性对整个体系的性质没太大影响,可以这么做。

            应该有更好的算法来计算简并的情况,有需要可以查些文献看看。我之后可能也会考虑下。

    2. 您好,我有一个问题,就是对于二阶简并陈数,是否可以先用这个程序算出F_xy,F_zw, F_xz,…..代入二阶陈数表达式中?

  4. 我把积分方法1的python代码改成了matlab的,但是算的陈数是几百,不知道错在哪里了,哪位大神可以帮忙看下错在哪?

    function []=chen_number_main()
        tic
        delta = 0.005;
        chern_number = 0;  %陈数初始化
        
        % 几个坐标中常出现的项
        a = 1/sqrt(3);
        aa1 = 4*sqrt(3)*pi/9/a;
        aa2 = 2*sqrt(3)*pi/9/a;
        bb1 = 2*pi/3/a;
    
       
        for kx =-aa1:delta:aa1%in np.arange(-aa1, aa1, delta):
            kx
            for ky=-bb1:delta:bb1%in np.arange(-bb1, bb1, delta):
                if (-aa2<=kxaa2 && -(aa1-kx)*tan(pi/3)<=ky<=(aa1-kx)*tan(pi/3))...
                    || (kx<-aa2 && -(kx-(-aa1))*tan(pi/3)<=ky<=(kx-(-aa1))*tan(pi/3))  % 限制在六角格子布里渊区内
    
                    H=hamiltonian0(kx, ky); 
                    [eigenvector,eigenvalue]=eig(H);
                    [~,index]=sort(real(eigenvalue));
                    vector = eigenvector(:,index(1));%价带波函数,本征值最小的点对应的波函数列向量
                   
                    
                    H_delta_kx=hamiltonian0(kx+delta, ky); 
                    [eigenvector, eigenvalue]=eig(H_delta_kx);
                    [~,index]=sort(real(eigenvalue));
                    vector_delta_kx=eigenvector(:,index(1));   % 略偏离kx的波函数
    
                    H_delta_ky=hamiltonian0(kx, ky+delta);  
                    [eigenvector, eigenvalue]=eig(H_delta_ky);
                    [~,index]=sort(real(eigenvalue));
                    vector_delta_ky=eigenvector(:,index(1));  % 略偏离ky的波函数
                    
                    H_delta_kx_ky=hamiltonian0(kx+delta, ky+delta);  
                   [eigenvector, eigenvalue]=eig(H_delta_kx_ky);
                   [~,index]=sort(real(eigenvalue));
                   vector_delta_kx_ky=eigenvector(:,index(1));  %略偏离kx和ky的波函数
                    
                    Ux=dot(conj(vector),vector_delta_kx)./abs(dot(conj(vector),...
                                                                             vector_delta_kx));
                    Uy=dot(conj(vector),vector_delta_ky)./abs(dot(conj(vector),...
                                                                             vector_delta_ky));
                    Ux_y=dot(conj(vector_delta_ky),vector_delta_kx_ky)./abs(...
                    dot(conj(vector_delta_ky),vector_delta_kx_ky));
                    Uy_x=dot(conj(vector_delta_kx),vector_delta_kx_ky)./abs(...
                    dot(conj(vector_delta_kx),vector_delta_kx_ky));
    
                    F = log(Ux.*Uy_x.*(1/Ux_y).*(1/Uy));
    
                    % 陈数(chern number)
                    chern_number=chern_number + F;
                end
            end
        end
        
        chern_number=chern_number/(2*pi*1j);
        disp(['Chern number=',num2str(chern_number)])
        toc
    end
    
    
    function [out]=hamiltonian0(kx,ky)
    out=hamiltonian(kx, ky, 2/3, 1, 1/3, pi/4, 1/sqrt(3));
    end
    
    
    function [matrix]=hamiltonian(k1, k2, M, t1, t2, phi, a) %Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
      
        h0 = zeros(2,2);  
        h2 = zeros(2,2);
    
       % 质量项(mass term), 用于打开带隙
        h0(1,1) = M;
        h0(2,2)= -M;
    
       %最近邻项
        h1(2,1) = t1*(exp(1j*k2*a)+...
        exp(1j*(sqrt(3)/2)*k1*a-(1j/2)*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-(1j/2)*k2*a));
        h1(1,2) = conj(h1(2,1));
    
        %次近邻项 # 对应陈数为-1
        h2(1,1) = t2*exp(-1j*phi)*(exp(1j*sqrt(3)*k1*a)+...
        exp(-1j*(sqrt(3)/2)*k1*a+1j*3/2*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-1j*(3/2)*k2*a));
        h2(2,2) = t2*exp(1j*phi)*(exp(1j*sqrt(3)*k1*a)+...
        exp(-1j*(sqrt(3)/2)*k1*a+1j*(3/2)*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-1j*(3/2)*k2*a));
        
        % 次近邻项  # 对应陈数为1
       % h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
        %h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    
        matrix = h0 + h1 + h2 + h2';
        
    end 

    1. 如果是一一对应的,应该不会有太大问题。你可以关闭循环,以某个值代替,然后在某个语句前后都分别输出数值结果,一步步对比验证,找到导致数值不一样的语句,就是出错的地方。

      1. 修改好的matlab代码,可以运行,速度较python快了很多,方法1:

        function []=chen_number_main()
        clear;clc
        tic
        delta = 0.005;
        chern_number = 0;  %陈数初始化
        %算某一能带的陈数可能在index中设置
        % 几个坐标中常出现的项
        %matlab的dot函数算复数似乎不是很好
        a = 1/sqrt(3);
        aa1 = 4*sqrt(3)*pi/9/a;
        aa2 = 2*sqrt(3)*pi/9/a;
        bb1 = 2*pi/3/a;
        ky0=-bb1:delta:bb1;
        kx0=-aa1:delta:aa1;
        temp1=length(kx0);
        temp2=length(ky0);
        for i=1:temp1
            kx0(i)
            for j=1:temp2
                if abs(kx0(i))<=aa2+tand(30)*(bb1-abs(ky0(j)))  % 限制在六角格子布里渊区内,只适用于正六角格子,新的限制算法,ref:https://blog.csdn.net/cuixiping/article/details/11716805?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-0&spm=1001.2101.3001.4242
                    kx=kx0(i);ky=ky0(j);
                    H=hamiltonian0(kx, ky);
                    [eigenvector,eigenvalue]=eig(H);
                    [~,index]=sort(real(diag(eigenvalue)));
                    vector = eigenvector(:,index(1));%价带波函数,本征值最小的点对应的波函数列向量
                    
                    
                    H_delta_kx=hamiltonian0(kx+delta, ky);
                    [eigenvector, eigenvalue]=eig(H_delta_kx);
                    [~,index]=sort(real(diag(eigenvalue)));
                    vector_delta_kx=eigenvector(:,index(1));   % 略偏离kx的波函数
                    
                    H_delta_ky=hamiltonian0(kx, ky+delta);
                    [eigenvector, eigenvalue]=eig(H_delta_ky);
                    [~,index]=sort(real(diag(eigenvalue)));
                    vector_delta_ky=eigenvector(:,index(1));  % 略偏离ky的波函数
                    
                    H_delta_kx_ky=hamiltonian0(kx+delta, ky+delta);
                    [eigenvector, eigenvalue]=eig(H_delta_kx_ky);
                    [~,index]=sort(real(diag(eigenvalue)));
                    vector_delta_kx_ky=eigenvector(:,index(1));  %略偏离kx和ky的波函数
                    
                    
                    Ux=(vector'*vector_delta_kx)./abs((vector'*vector_delta_kx));
                    Uy=(vector'*vector_delta_ky)./abs((vector'*vector_delta_ky));
                    Ux_y=(vector_delta_ky'*vector_delta_kx_ky)./abs((vector_delta_ky'*vector_delta_kx_ky));
                    Uy_x=(vector_delta_kx'*vector_delta_kx_ky)./abs((vector_delta_kx'*vector_delta_kx_ky));
                    
                    
                    F = log(Ux.*Uy_x.*(1/Ux_y).*(1/Uy));
                    
                    % 陈数(chern number)
                    chern_number=chern_number + F;
                end
            end
        end
        
        chern_number=chern_number/(2*pi*1j);
        disp(['Chern number=',num2str(chern_number)])
        toc
        end
        
        
        function [out]=hamiltonian0(kx,ky)
        out=hamiltonian(kx, ky, 2/3, 1, 1/3, pi/4, 1/sqrt(3));
        end
        
        
        function [matrix]=hamiltonian(k1, k2, M, t1, t2, phi, a) %Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
        
        h0 = zeros(2,2);
        h2 = zeros(2,2);
        
        % 质量项(mass term), 用于打开带隙
        h0(1,1) = M;
        h0(2,2)= -M;
        
        %最近邻项
        h1(2,1) = t1*(exp(1j*k2*a)+...
            exp(1j*(sqrt(3)/2)*k1*a-(1j/2)*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-(1j/2)*k2*a));
        h1(1,2) = conj(h1(2,1));
        
        %次近邻项 # 对应陈数为-1
        h2(1,1) = t2*exp(-1j*phi)*(exp(1j*sqrt(3)*k1*a)+...
            exp(-1j*(sqrt(3)/2)*k1*a+1j*3/2*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-1j*(3/2)*k2*a));
        h2(2,2) = t2*exp(1j*phi)*(exp(1j*sqrt(3)*k1*a)+...
            exp(-1j*(sqrt(3)/2)*k1*a+1j*(3/2)*k2*a)+exp(-1j*(sqrt(3)/2)*k1*a-1j*(3/2)*k2*a));
        
        % 次近邻项  # 对应陈数为1
        % h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
        %h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
        
        matrix = h0 + h1 + h2 + h2';
        
        end 

        1. 嗯,在有一些情况下(比如可能是循环),Matlab会比Python快很多。我之所以不怎么爱用Matlab,是因为不是开源的,而且软件安装比较麻烦,体积太大。如果追求运行效率,我一般会用Fortran或者C语言。

              1. 我又弄了新的更简洁的限制在六边形布里渊区内的代码,如下:
                function [out]=shape_f(xx,yy,xp,yp)
                %xx、yy是输入的要判断的点,
                %xp 、yp是多边形的顶点(顺时针或逆时针,首尾相连)
                %xp、yp是向量
                %如果在多边形内或边上,out为1,否则为0
                out=double(inpolygon(xx,yy,xp,yp));

                end

                1. 嗯,都是可以的。我这里的代码也没有做过多的优化,只是作为例子实现对应的功能。此外,这个对运行速度的影响应该不是很大。

  5. # 这里乘2或除以2是为了保证“步长与积分个数的乘积刚好为布里渊区面积”
    delta_x = delta*bb1*2
    delta_y = delta*bb2*2/2

    这个地方的delta_x,y根据文中说法不应该是在bb1,2后面成以sqrt(2)吗?
    所以是不是手误打错了?

    1. 不算是错的,两种都可以。constant取\sqrt{2},或者一个constant取2,一个constant取1。只要坐标变换后直接积分的结果是布里渊区的面积,步长的选取不影响结果。尤其是当步长足够小时,在数值上没有太大区别。(在博文中也添加这个评论)

  6. 博主你好

    不同于这篇文章的陈数,请问博主接触过高阶拓扑绝缘体的拓扑不变量(Bulk Polarization)吗?

    例如计算出kagome模型的陈数是变化的,但是其Bulk Polarization就是 (1/3, 1/3) ,近期一直尝试计算,可是一直一直出错。。。:(

    https://zhuanlan.zhihu.com/p/175953375

    1. 我目前还没怎么接触高阶拓扑绝缘体的拓扑不变量,之后可能也会考虑。你可以多看些文献,确保对该不变量的理解是正确的,或者用简单体系做一下测试,找找程序的问题。

  7. 您好,感谢分享。我觉得你这里程序中的“a = 1/sqrt(3)”是错误的,应该是”a=1″,如果a是 carbon-carbon distance的话。可以参考review “The electronic properties of graphene”。

    1. 不影响的,a可以取任何相对的值,对应的布里渊区大小会有变化。在博文中我是取次近邻原子间的距离为1,所以a = 1/sqrt(3)。如果a取为1也是可以的。

  8. 我用你的“陈数Chern number的计算(定义法)(附Python、Matlab代码)”博文里面提到的matlab代码,并且结合你在“Haldane模型哈密顿量与能带图(附Python代码)”博文开头提到的Haldane模型的哈密顿量(含泡利矩阵那个),选取原子间距为1时,得到陈数为2.8690,原子间距为1/sqrt(3)时得到陈数为0.6313(m=2/3,t1=1,t2=2/3,phi=0.25pi),一直想不通哪里出了问题,回到你现在这篇博文的代码里发现你这里的哈密顿量好像有一点点不一样…

    matlab代码(a=1):

    %陈数定义法
    clear;clc;
    n=100;%积分密度
    delta=1e-5;%求导的偏离量
    C=0;%陈数初始化
    for kx=-pi:(1/n):pi
        for ky=-pi:(1/n):pi
            VV=get_vector(HH(kx,ky));
            Vkx=get_vector(HH(kx+delta,ky));%略偏离kx的波函数
            Vky=get_vector(HH(kx,ky+delta));%略偏离ky的波函数
            Vkxky=get_vector(HH(kx+delta,ky+delta));%略偏离kx,ky的波函数
            if  sum((abs(Vkx-VV)))>0.01   %为了波函数的连续性(这里的不连续只遇到符号问题,所以可以直接这么处理)
                Vkx=-Vkx;
            end
            
            if  sum((abs(Vky-VV)))>0.01
                Vky=-Vky;
            end
            
            if  sum(abs(Vkxky-VV))>0.01
                Vkxky=-Vkxky;
            end
            %价带的波函数的berry connection,求导后内积
            Ax=VV'*(Vkx-VV)/delta;%Berry connection Ax
            Ay=VV'*(Vky-VV)/delta;%Berry connection Ay
            Ax_delta_ky=Vky'*(Vkxky-Vky)/delta;%略偏离ky的berry connection Ax
            Ay_delta_kx=Vkx'*(Vkxky-Vkx)/delta;%略偏离kx的berry connection Ay
            %berry curvature
            F=((Ay_delta_kx-Ay)-(Ax_delta_ky-Ax))/delta;
            %chern number
            C=C+F*(1/n)^2;  
        end
    end
    C=C/(2*pi*1i)
    
    
    function vector_new = get_vector(H)
    [vector,eigenvalue] = eig(H);
    [eigenvalue, index]=sort(diag(eigenvalue), 'descend');
    vector_new = vector(:, index(2));
    end
    
    function H=HH(kx,ky)
    H(1,1)=-2*1/3*cos(0.25*pi).*(cos(sqrt(3)*1*kx)+cos(-0.5*sqrt(3)*1*kx+3/2*1*ky)+cos(0.5*sqrt(3)*1*kx+3/2*1*ky))+2/3+2*1/3*sin(0.25*pi).*(sin(sqrt(3)*1*kx)+sin(-0.5*sqrt(3)*1*kx+3/2*1*ky)+sin(-0.5*sqrt(3)*1*kx-3/2*1*ky));
    H(1,2)=-1*(cos(1*ky)+cos(-0.5*sqrt(3)*1*kx-0.5*1*ky)+cos(0.5*sqrt(3)*1*kx-0.5*1*ky))-1i*(1*(sin(1*ky)+sin(-0.5*sqrt(3)*1*kx-0.5*1*ky)+sin(0.5*sqrt(3)*1*kx-0.5*1*ky)));
    H(2,1)=-1*(cos(1*ky)+cos(-0.5*sqrt(3)*1*kx-0.5*1*ky)+cos(0.5*sqrt(3)*1*kx-0.5*1*ky))+1i*(1*(sin(1*ky)+sin(-0.5*sqrt(3)*1*kx-0.5*1*ky)+sin(0.5*sqrt(3)*1*kx-0.5*1*ky)));
    H(2,2)=-2*1/3*cos(0.25*pi).*(cos(sqrt(3)*1*kx)+cos(-0.5*sqrt(3)*1*kx+3/2*1*ky)+cos(0.5*sqrt(3)*1*kx+3/2*1*ky))-(2/3+2*1/3*sin(0.25*pi).*(sin(sqrt(3)*1*kx)+sin(-0.5*sqrt(3)*1*kx+3/2*1*ky)+sin(-0.5*sqrt(3)*1*kx-3/2*1*ky)));
    end


    然后我就去按照你有关python博文里面搞了一下python,把我上面matlab代码里面的哈密顿量化简后替代现在这篇博文里面代码的哈密顿量,得到
    chern number=1.0000000000000133+2.363979843223985e-13j

    python代码:

    import numpy as np
    import matplotlib.pyplot as plt
    from math import *   # 引入pi, cos等
    import cmath
    import time
    import functools  # 使用偏函数functools.partial()
    
    
    def hamiltonian(k1, k2, M, t1, t2, phi, a=1):  # Haldane哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
        # 初始化为零矩阵
        h0 = np.zeros((2, 2))*(1+0j)   # 乘(1+0j)是为了把h0转为复数
        h1 = np.zeros((2, 2))*(1+0j)
        h2 = np.zeros((2, 2))*(1+0j)
    
        # 质量项(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*(cmath.exp(1j*k2*a)+2*cmath.cos(0.5*sqrt(3)*k1*a)*cmath.exp(-1j*0.5*k2*a))
        h1[0, 1] = t1*(cmath.exp(-1j*k2*a)+2*cmath.cos(0.5*sqrt(3)*k1*a)*cmath.exp(1j*0.5*k2*a))
        #次近邻项 # 对应陈数为-1
    #     h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    #     h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
        h2[0, 0] = -2*t2*(-cmath.cos(phi+sqrt(3)*k1*a)-2*cmath.cos(3/2*k2*a)*cmath.cos(phi-0.5*sqrt(3)*k1*a))
        h2[1, 1] = 2*t2*(cmath.cos(phi-sqrt(3)*k1*a)+2*cmath.cos(3/2*k2*a)*cmath.cos(phi+0.5*sqrt(3)*k1*a))
        
        # # 次近邻项  # 对应陈数为1
        # h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
        # h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    
        matrix = h0 + h1 + h2
        return matrix
    
    
    def main():
        start_clock = time.perf_counter()
        delta = 0.002
        chern_number = 0  # 陈数初始化
        
        # 常出现的项
        a = 1
        bb1 = 2*sqrt(3)*pi/3/a
        bb2 = 2*pi/3/a
    
        hamiltonian0 = functools.partial(hamiltonian, M=2/3, t1=1, t2=1/3, phi=pi/4, a=a)   # 使用偏函数,固定一些参数
    
        for k1 in np.arange(0 , 1, delta):
            print(k1)
            for k2 in np.arange(0, 1, delta):
                # 坐标变换
                kx = (k1-k2)*bb1
                ky = (k1+k2)*bb2
    
                # 这里乘2或除以2是为了保证“步长与积分个数的乘积刚好为布里渊区面积”
                delta_x = delta*bb1*2    
                delta_y = delta*bb2*2/2
    
                H = hamiltonian0(kx, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H)
                vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
            
                H_delta_kx = hamiltonian0(kx+delta_x, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
                vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]   # 略偏离kx的波函数
    
                H_delta_ky = hamiltonian0(kx, ky+delta_y)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
                vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离ky的波函数
                
                H_delta_kx_ky = hamiltonian0(kx+delta_x, ky+delta_y)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
                vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 略偏离kx和ky的波函数
                
                Ux = np.dot(np.conj(vector), vector_delta_kx)/abs(np.dot(np.conj(vector), vector_delta_kx))
                Uy = np.dot(np.conj(vector), vector_delta_ky)/abs(np.dot(np.conj(vector), vector_delta_ky))
                Ux_y = np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_ky), vector_delta_kx_ky))
                Uy_x = np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky)/abs(np.dot(np.conj(vector_delta_kx), vector_delta_kx_ky))
    
                F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))
                # 陈数(chern number)
                chern_number = chern_number + F
    
        chern_number = chern_number/(2*pi*1j)
        print('Chern number = ', chern_number)
        end_clock = time.perf_counter()
        print('CPU执行时间(min)=', (end_clock-start_clock)/60)
    
    
    if __name__ == '__main__':
        main()


    嗯…可能那篇博文给的matlab代码算法上不适合算蜂窝格子吧…
    对我帮助很大,非常感谢!!!

    1. 你前面matlab代码遇到的问题是因为你积分没有限制在布里渊区内。六角格子的布里渊区还是六角格子,已经不能在[-pi, pi]*[-pi, pi]区域积分了。有两种处理办法,一个是用条件语句把积分限制在布里渊区,还有一个办法是坐标变换。本篇列出的就是这两种处理办法。

      此外,也尽量少用定义的方法吧。定义法可以作为入门,比较直观,但在计算中经常会遇到波函数不连续的问题。这里是2乘2的矩阵,波函数还不会遇到麻烦,而对于矩阵比较大的矩阵时,求解时波函数的连续性会比较难处理。

发表评论

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