学术源码, 拓扑不变量

陈数Chern number的计算(定义法,附Python/Matlab代码)

其他几篇关于陈数的计算方法:

本篇采用的是定义法计算陈数,而不是用文献[1, 2]中高效的计算方法。定义法在数值计算上可能会比较耗时,但在理解上比较直接和方便:先是计算贝里联络 (Berry connection),然后再计算贝里曲率 (Berry curvature),积分得到陈数 (Chern number)。公式为[1]:

这里以文献[2]中的“量子反常霍尔模型”的哈密顿量为例子,用定义法在数值上计算陈数。

1. 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/3932
"""

import numpy as np
from math import * 
import time


def hamiltonian(kx, ky):  # 量子反常霍尔QAH模型(该参数对应的陈数为2)
    t1 = 1.0
    t2 = 1.0
    t3 = 0.5
    m = -1.0
    matrix = np.zeros((2, 2), dtype=complex)
    matrix[0, 1] = 2*t1*cos(kx)-1j*2*t1*cos(ky)
    matrix[1, 0] = 2*t1*cos(kx)+1j*2*t1*cos(ky)
    matrix[0, 0] = m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky)
    matrix[1, 1] = -(m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky))
    return matrix


def main():
    start_time = time.time()
    n = 100  # 积分密度
    delta = 1e-9  # 求导的偏离量
    chern_number = 0  # 陈数初始化
    for kx in np.arange(-pi, pi, 2*pi/n):
        for ky in np.arange(-pi, pi, 2*pi/n):
            H = hamiltonian(kx, ky)
            eigenvalue, eigenvector = np.linalg.eig(H)
            vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
            # print(np.argsort(np.real(eigenvalue))[0])  # 排序索引(从小到大)
            # print(eigenvalue)  # 排序前的本征值
            # print(np.sort(np.real(eigenvalue)))  # 排序后的本征值(从小到大)
           
            H_delta_kx = hamiltonian(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 = hamiltonian(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 = hamiltonian(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的波函数

            # 价带的波函数的贝里联络(berry connection) # 求导后内积
            A_x = np.dot(vector.transpose().conj(), (vector_delta_kx-vector)/delta)   # 贝里联络Ax(x分量)
            A_y = np.dot(vector.transpose().conj(), (vector_delta_ky-vector)/delta)   # 贝里联络Ay(y分量)
            
            A_x_delta_ky = np.dot(vector_delta_ky.transpose().conj(), (vector_delta_kx_ky-vector_delta_ky)/delta)  # 略偏离ky的贝里联络Ax
            A_y_delta_kx = np.dot(vector_delta_kx.transpose().conj(), (vector_delta_kx_ky-vector_delta_kx)/delta)  # 略偏离kx的贝里联络Ay

            # 贝里曲率(berry curvature)
            F = (A_y_delta_kx-A_y)/delta-(A_x_delta_ky-A_x)/delta

            # 陈数(chern number)
            chern_number = chern_number + F*(2*pi/n)**2
    chern_number = chern_number/(2*pi*1j)
    print('Chern number = ', chern_number)
    end_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)


if __name__ == '__main__':
    main()

计算结果:

Chern number = (2.0000002869141356-4.9548136005219324e-09j)

2. Matlab代码

该代码来源于读者提供的代码,经过一些修改。需要说明的是:这里对波函数连续性的问题只是对最简单情况的临时处理,只限于当前的模型和参数。正确的处理方式看本篇最后内容提供的链接。

% 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/3932

% 陈数定义法
clear;clc;
n=100; % 积分密度
delta=1e-9; % 求导的偏离量
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,2)=2*cos(kx)-1i*2*cos(ky);
H(2,1)=2*cos(kx)+1i*2*cos(ky);
H(1,1)=-1+2*0.5*sin(kx)+2*0.5*sin(ky)+2*cos(kx+ky);
H(2,2)=-(-1+2*0.5*sin(kx)+2*0.5*sin(ky)+2*cos(kx+ky));
end

以上matlab代码中,波函数符号会出现波动,因此要增加以下三段代码(原因是波函数整体符号变号,即相位为pi时,仍然是哈密顿量的波函数,但在做波函数相减时会完全不同):

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

此外,有时加了以上这三段还不够,波函数会时有以下这种形式(一个例子),仍然是不连续:

这时候取负号是不能解决波函数的不连续问题,这里临时用取模来处理(原理上是要通过选取规范进行变换,得到波函数(1, 0),下面会提到):

if abs(VV(1))<0.0001 ||abs(VV(2))<0.0001
      VV=abs(VV);
      Vkx=abs(Vkx);
      Vky=abs(Vky);
      Vkxky=abs(Vkxky);
 end

以上只考虑了2×2的哈密顿量的情况,如果维度再加上去,还会出现其他不连续的情况,在Python中也会遇到这个问题。

3. 总结

这里因为以下波函数都是哈密顿量的解:

所以规范 (gauge) 的选取在定义法中就显得很重要了,需要选定一个规范,使得局域的波函数连续。以下是文献[1]中截图:

在数值上可以写代码去搜索找个近似的规范 (前面提到的“取负号”、“取模”只是两种特殊情况) ,使得波函数连续,但在计算效率上又降低很多了。搜索连续规范的方法参考:在数值上寻找波函数近似的同一规范(二分查找)。或者采用固定某个波函数规范的方法 ,可参考这篇SSH模型中的Wilson loop(附Python代码)中的代码。这两种方法(二分查找和固定规范)都可以通过开源项目Guan进行函数的调用:https://py.guanjihuan.com

这里还是推荐使用高效法来计算(不要求波函数连续性),见参考文献[1, 2]或者参考这篇博文:陈数Chern number的计算(高效法,附Python/Matlab代码)。此外,还有方法:陈数Chern number的计算(Kubo公式,附Python代码)陈数Chern number的计算(Wilson loop方法,附Python代码)等。

4. Python代码(增加寻找连续规范的代码)

前面python代码如果取不同的n,在某些情况下会出现错误的结果。遇到这个问题是正常的,因为前面的波函数没有选取连续的规范。这里增加了寻找连续规范的代码,参考:在数值上寻找波函数近似的同一规范(二分查找,附Python代码)。或直接调用开源项目Guan:https://py.guanjihuan.com。此外也可以采用固定某个规范,参考这篇中的代码:SSH模型中的Wilson loop(附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/3932
"""

import numpy as np
from math import *
import time
import cmath


def hamiltonian(kx, ky):  # 量子反常霍尔QAH模型(该参数对应的陈数为2)
    t1 = 1.0
    t2 = 1.0
    t3 = 0.5
    m = -1.0
    matrix = np.zeros((2, 2), dtype=complex)
    matrix[0, 1] = 2*t1*cos(kx)-1j*2*t1*cos(ky)
    matrix[1, 0] = 2*t1*cos(kx)+1j*2*t1*cos(ky)
    matrix[0, 0] = m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky)
    matrix[1, 1] = -(m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky))
    return matrix


def main():
    start_time = time.time()
    n = 20  # 积分密度
    delta = 1e-9  # 求导的偏离量
    chern_number = 0  # 陈数初始化
    for kx in np.arange(-pi, pi, 2*pi/n):
        for ky in np.arange(-pi, pi, 2*pi/n):
            H = hamiltonian(kx, ky)
            eigenvalue, eigenvector = np.linalg.eig(H)
            vector = eigenvector[:, np.argsort(np.real(eigenvalue))[0]]  # 价带波函数
           
            H_delta_kx = hamiltonian(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 = hamiltonian(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 = hamiltonian(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的波函数

            vector_delta_kx = find_vector_with_the_same_gauge(vector_delta_kx, vector)
            vector_delta_ky = find_vector_with_the_same_gauge(vector_delta_ky, vector)
            vector_delta_kx_ky  = find_vector_with_the_same_gauge(vector_delta_kx_ky, vector)

            # 价带的波函数的贝里联络(berry connection) # 求导后内积
            A_x = np.dot(vector.transpose().conj(), (vector_delta_kx-vector)/delta)   # 贝里联络Ax(x分量)
            A_y = np.dot(vector.transpose().conj(), (vector_delta_ky-vector)/delta)   # 贝里联络Ay(y分量)
            
            A_x_delta_ky = np.dot(vector_delta_ky.transpose().conj(), (vector_delta_kx_ky-vector_delta_ky)/delta)  # 略偏离ky的贝里联络Ax
            A_y_delta_kx = np.dot(vector_delta_kx.transpose().conj(), (vector_delta_kx_ky-vector_delta_kx)/delta)  # 略偏离kx的贝里联络Ay

            # 贝里曲率(berry curvature)
            F = (A_y_delta_kx-A_y)/delta-(A_x_delta_ky-A_x)/delta

            # 陈数(chern number)
            chern_number = chern_number + F*(2*pi/n)**2
    chern_number = chern_number/(2*pi*1j)
    print('Chern number = ', chern_number)
    end_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)


def find_vector_with_the_same_gauge(vector_1, vector_0):
    # 寻找近似的同一的规范
    phase_1_pre = 0
    phase_2_pre = pi
    n_test = 10001
    for i0 in range(n_test):
        test_1 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_1_pre) - vector_0))
        test_2 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_2_pre) - vector_0))
        if test_1 < 1e-8:
            phase = phase_1_pre
            # print('Done with i0=', i0)
            break
        if i0 == n_test-1:
            phase = phase_1_pre
            print('Gauge Not Found with i0=', i0)
        if test_1 < test_2:
            if i0 == 0:
                phase_1 = phase_1_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_1_pre
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
        else:
            if i0 == 0:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre 
        phase_1_pre = phase_1
        phase_2_pre = phase_2
    vector_1 = vector_1*cmath.exp(1j*phase)
    # print('二分查找找到的规范=', phase)   
    return vector_1


if __name__ == '__main__':
    main()

参考文献:

[1] Chern Numbers in Discretized Brillouin Zone: Efficient Method of Computing (Spin) Hall Conductances

[2] Route towards Localization for Quantum Anomalous Hall Systems with Chern Number 2

[3] Berry phase effects on electronic properties

[4] Quantum anomalous Hall effect in graphene from Rashba and exchange effects

[5] Localization trajectory and Chern-Simons axion coupling for bilayer quantum anomalous Hall systems

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

9,657 次浏览

37 thoughts on “陈数Chern number的计算(定义法,附Python/Matlab代码)”

  1. 请问最后积分时为什么是乘以(1/n)^2而不是delta^2呢?是不是这样才能覆盖整个积分区域,那理想情况是不是将(1/n)和delta取相同的值

    1. 一个是积分的步长,一个是求导的步长。因为积分步长和循环次数有关,所以不能取太小,不然计算时间会很长。求导步长是越小越好,积分步长也是越小越好。两者可以取为相同的值,也可以取为不同的值,都只是对极限(极小)的近似。

  2. 博主,你好,我想请教一个问题,在一般的拓扑体系中,体系哈密顿量满足时间反演对称,那该体系对应的陈数必为零。那如果一个非平庸的拓扑体系哈密顿量,不满足时间反演对称,该体系的陈数有没有可能也是为零的?

    1. 不满足时间反演对称不一定是拓扑的。你指的非平庸的拓扑体系看指的是什么。如果是陈拓扑绝缘体,那么陈数就不为零。如果是Z2拓扑绝缘体,那么陈数为零,但这个是满足时间反演对称性的。如果是其他的拓扑体系,比如拓扑半金属、高阶拓扑绝缘体等,是由其他的不变量来描述。

  3. 请问一下我带入一个方格子哈密顿量算出来是陈数-1咋办,这是啥样的拓扑性质呢

  4. 您好,我有个问题想问一下:在SSH model中,我们先用周期性边界条件得到绕数,再用开边界得到边缘态,那么在二维的时候,是不是在两个方向上都得采用周期性边界条件,得到陈数,然后再对一个方向用开边界条件,得到边缘态?

    1. 不是周期性边界条件,而是本身考虑的就是二维的体系,包含了kx和ky,然后计算陈数。在开放边界条件下,会存在边缘态。

  5. 这个是只要是二阶矩阵哈密顿量就可以吗,我代入了一个正方格子模型,发现随着积分密度的变化,chern数有很大的变化

    1. 正方格子的chern数为零,值应该应该是在1e-8以下,数值变化很正常,理论上都是0。

  6. 这个是不是只要是二阶哈密顿量,带入H(1,1)H(1,2)H(2,1)H(2,2)就可以计算出chern数,我看到我带去运行的命令窗口H一直在闪

  7. 感谢博主的分享,我在重复博主分享的参考文献[1]的时候出现了一个问题,使用文献[1]中的(5)式计算的时候发现会出现非常大的值,而且求和为0,博主能帮我解决一下下吗?

  8. 请问一下,在python的方法中,main函数的作用是什么,能否把它改成其他函数。比如在最后用return语句

    1. def main()是主函数,其实也是普通的函数,可以用别的符号代替,例如def zhuhanshu()。因为在程序中,我直接运行main(),所以它就执行该函数的内容,该语句经常也写为

      if __name__ == '__main__':  # 如果是当前文件直接运行,执行main()函数中的内容;如果是import当前文件,则不执行。
          main()

  9. 你好,请问为何要设置偏离kx和ky的Ax,Ay呢?从公式上看delta的取值对于计算结果的大小有非常大的影响,在最后F的公式中就需要除以delta的平方,请问是什么原因?

    1. 没有平方。delta取小量就可以,影响不大。整个过程就是数值上的求导,把求导变成除法。

  10. 太感谢了!支持支持
    求得的陈数是一个近似的整数,请问有没有计算方法能得到一个确定的整数呢?

    1. 把积分的步长取小一点,就会接近于一个整数。步长越小,越接近,但计算的时间会增加很多,你可以计算试一下。所以取一个合理的步长就可以了,数值上的计算有一些误差也是正常的,误差很小时就可以忽略不计。

  11. 同学,你好!想问你两个问题。
    1. 为什么 matlab 里面中 C=C+F*(1/n)^2; 而python中 chern_number = chern_number + F*(2*pi/n)**2。为什么这里matlab中计算公式少了 2*pi 呢?
    2. 对于不含参的哈密顿量,你知道有什么方法可以计算chern number吗?

    1. (1) 两个代码的积分密度n表示的意思有点不一样。Python代码中步长是2*pi/n,Matlab中步长是1/n。也可以试着把代码改成一样的。
      (2)Chern number的计算一般都是在倒空间吧。如果是实空间的,要傅里叶变换到倒空间,然后就存在参数kx和ky。

    1. [vector,eigenvalue] = eig(H);
      这个语句。哈密顿量的本质态为vector和本质值为eigenvalue,eig(H)计算出来的结果是一一对应的。

  12. 非常感谢博主的分享,请问一下,如果是含三个参数,例如在三维空间的陈数能不能改成取特定值,在博主的程序的基础上加个循环取某些特定的值呢?

    1. 我没了解过,不知道三维的陈数是否存在?这个比较偏数学。你可以找一些文献看看,如果某维度可行的话,积分的话肯定是要在对应的维度积分。我目前只接触到二维的,即使是三维的Weyl半金属,也是固定一个参数,按切片来计算陈数的。(其他读者有了解的可以留言)

  13. 哈喽,同学你好,我想问一个问题,code里面的陈数的计算
    # 陈数(chern number)
    chern_number = chern_number + F*(2*pi/n)**2,
    中的 F*(2*pi/n)**2是怎么来的呀?

    1. 要对贝里曲率F积分。布里渊区正方形两个边长分别分成(2*pi/n)份,所以面积分成了(2*pi/n)**2份。

发表评论

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