拓扑不变量, 学术

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

定义法计算陈数参考这篇:陈数Chern number的计算(定义法,附Python/Matlab代码)

此外还有:

本篇采用的是参考文献[1]中的方法,这里截图一部分出来:

以参考文献[2]中量子反常霍尔模型为例子。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/4179
"""

import numpy as np
import matplotlib.pyplot as plt
from math import *   # 引入pi, cos等
import cmath
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))*(1+0j)
    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 = 2*pi/n
    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的波函数
            
            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_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)


if __name__ == '__main__':
    main()

计算结果为:

Chern number = (1.9999999999999998-1.7272599290109276e-16j)

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/4179

% 陈数高效法
clear;clc;
n=1000 % 积分密度
delta=2*pi/n;
C=0; 
for kx=-pi:(2*pi/n):pi
    for ky=-pi:(2*pi/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的波函数
        
        Ux = VV'*Vkx/abs(VV'*Vkx);
        Uy = VV'*Vky/abs(VV'*Vky);
        Ux_y = Vky'*Vkxky/abs(Vky'*Vkxky);
        Uy_x = Vkx'*Vkxky/abs(Vkx'*Vkxky);
        
        % berry curvature
        F=log(Ux*Uy_x*(1/Ux_y)*(1/Uy));
        % chern number
        C=C+F;  
    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

补充说明1:高效法中的偏移量和积分步长要保持一致才行。

补充说明2:公式定义中的\tilde{F}本身已经包含了步长,所以最后只要直接求和就可以了,在文献[1]中的Eq.12到Eq.13之间有提到,见下面的截图。所以,如果要得到贝里曲率,\tilde{F}还需要除以两个方向的步长。可以大概这么去理解:通常定义的贝里联络A需要对k求导,但在文献[1]中定义的\tilde{A}中只有偏移量,没有除以步长。如果对\tilde{A}做一个泰勒展开,和A刚好就差一个步长。

参考资料:

[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

9,902 次浏览

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

评论说明:
(1)在保留浏览器缓存的前提下,目前支持72小时自主修改或删除个人评论。如果自己无法修改或删除评论,可再次评论或联系我。如有发现广告留言,请勿点击链接,博主会不定期删除。
(2)评论支持Latex公式。把latexpage作为标签放在任何位置,评论中的公式可正常编译,示例:
$Latex formula$  [latexpage]

56 thoughts on “陈数Chern number的计算(高效法,附Python/Matlab代码)”

  1. 老师您好我目前有关于C4算拓扑体极化的代码(其中数据由comsol导出),我想在这基础上算有关C6的拓扑体极化,需要在C4基础上改哪些内容

    1. 布里渊区不一样,而且可能需要考虑具体某个极化的方向吧,不同极化方向选择的k路径不一样。可以查这方面的文献看看,我目前没做过这个

  2. 您好,请问这个计算方法对于布里渊区不为长方形的情形有效吗?比如对菱形区域的布里渊区?是否要乘上一个系数?

      1. 老师您好!为什么您的Matlab程序最后求陈数积分的时候的时候不需要乘以步长

        1. 公式定义中的\tilde{F}本身已经包含了步长,所以最后只要直接求和就可以了,在文献[1]中的Eq.12到Eq.13之间有提到。所以如果要得到贝里曲率,\tilde{F}还需要除以两个方向的步长。可以大概这么去理解:通常定义的贝里联络A需要对k求导,但在定义的\tilde{A}中只有偏移量,没有除以步长。如果对\tilde{A}做一个泰勒展开,和A刚好就差一个步长。

  3. 您好,我用了这些代码之后发现有些带算出来的结果和我取得n有关,请问是什么原因,比如我取n=40,就会有些带的陈数是80

  4. 关博士您好,请问能带的陈数一定是整数吗?我计算过程中遇到用价带陈数相加作为体系的陈数的情况,这时各个价带的陈数可以不是整数吗?

    1. 如果不存在能带简并或交叉以及电子半填充的情况,单条能带携带的陈数应该都要为整数。

  5. 请问大佬,这四种方法有可以用来算那种能带存在简并的体系吗(6*6的矩阵只有三个带),我自己试了一下,发现算出来的陈数对不上,请问有啥方法解决吗,跪求!!!

        1. 非常感谢。我也在您的另一篇博客(haldan模型)中看到评论中有人推广到简并情况,我最终的目的是计算简并二阶陈数,c2=1/(4pi∧2)sum(F_xyF_zw+F_xzF_yw+F_xwF_zy)请问可以先用一阶的方法算出F_xy, F_zw, F_xz等六项再带入计算二阶陈数吗?

          1. 我不了解你这个二阶陈数的概念,要确保F_xy等项是被良好定义的,是不是和这里定义的一样。另外如果是简并的,波函数可能会出现问题,因为波函数的任意线性组合都是简并能带的波函数。我不确定是不是可以直接用,你也可以测试下看算得对不对。

            1. 我有做过测试计算4维狄拉克模型(二重简并)的second chern number。但是很奇怪,计算结果随着delta_k 的取值有很大变化(四个纬度delta_k一致,由于计算量的关系 ,我最小试到delta_k=2pi/50),偶尔才能与理论结果一致。如果您愿意并且有时间的话,我可以加您微信一起讨论一下可能的原因吗?非常感谢

  6. 关老师你好,请问一下这个方法适用于平移对称被破坏的实空间哈密顿量吗?

    1. 这个是在倒空间计算的,是二维的体系。如果没有空间平移对称性,可以考虑周期性边界条件,认为是无限大的,也可以考虑实空间计算陈数的方法,我目前还没算过,你可以查查文献。

  7. vector_new = vector(:, index(2));
    请问这里的index(2)是要随着哈密顿的维度而改变吗,假如我要算的是8X8的哈密顿量那就是index(8)吧

    1. 这是取价带,因为这里只有两条带,所以取第2条。如果哈密顿量是8x8的,能带可以取第5、6、7、8条,前提是要在排序后。

      1. 请问老师如果我是8X8的无自旋超导哈密顿量,一个原胞四个电子,四个空穴;电子掺杂浓度为2.12;doping0.12;那么还是取价带的四条能带算陈数吗,,

  8. 请问老师在matlab代码里的偏函数HH中,我新输入一个8X8的哈密顿,算出来陈数是-1.8514+0000i,这代表是拓扑平庸还是非平庸呢

    1. 可能哪里有问题吧,比如存在简并什么的。一般来说,算出来应该为整数。

  9. 暂且称呼您为关老师吧!
    关老师您好,我用您的高效法的matlab代码,设置n=500,算出的C=2.002638...,我看您用python设置n=100算出的C都比我精确,我算的是有什么问题吗?

    1. 我也不确定。如果逻辑都一样的话应该没有太大问题。之所以结果不一样,有可能是因为Matlab和Python计算特征向量时有不同的结果,恰巧一种算法得到的特征向量在后续的计算中更收敛。

  10. 高效法和定义法计算chern number,积分范围都是选取的(-pi,pi),不用对应格子的第一布里渊区吗?

    1. 这里的模型类似于方格子,第一布里渊区布的范围刚好是kx,ky为(-pi,pi)。如果是其他形式的,需要对布里渊区积分。

  11. 哈密顿量二次量后,矩阵就是所要表达的哈密顿量,直接代入就可?

      1. 明白了,谢谢。不过我用文献中方法计算chern number结果为0,正负1,正负2,但是用高效法计算结果一直为0(不论是否赋值)

  12. 博主在取离散kx,ky时,用的代码是 for kx in np.linspace(-pi, pi, n).

    这里用np.linspace给出的kx包括左右边界,会导致重复求和,因此得到的Chern number也精度也不高(例题主代码中n=1000,C=1.997...)。改成for kx in np.arange(-pi, pi, 2*pi/n) 的话,n很小就能得到正确结果(n=4得到C=2)。

    1. 感谢!用np.linspace是会存在多取点的问题,之前没关注到这个,代码已更新。

    1. HH(kx,ky)返回的是哈密顿量矩阵。如果你已知的是每个(kx,ky)对应的波函数,那么直接替换get_vector()就可以了。

  13. 您好博主,我想问一下,是否有类似的方法算出贝里联络,现在需要算贝里联络(还是2*2的矩阵),最后积分,请问有类似避开随机相位的方法么?

  14. 您好,请问下参考文献1公式(7)后面有一段话“ The link variables are well defined as long as N_mu neq 0, which can always be assumed to be the case (one can avoid a singularity N_mu = 0 by an infinitesimal shift of the lattice).”这里“an infinitesimal shift of the lattice”该如何理解?计算中又是怎么操作的了?

    1. 要求是分母不为零,这时候U才可以定义。如果在某些取值下,分母出现为零的情况,这时候可以稍微加一些微小量,例如1e-6,使得在数值上避开这个奇异点。目前我是没遇到这个情况,你如果有遇到这个情况,可以试试看这种处理方式。

  15. 您好,感谢分享,对我帮助很大。还想请教一个问题:我在MATLAB下使用这种高效法计算一个多能带体系的时候,发现本征矢矩阵中出现了很多零项,导致最后有些能带陈数为NaN。请问这可能是什么原因了?

    1. 本征矢为零?我也不清楚是什么原因导致的。可以试着看错误结果在每一行代码中的中间数据,找找什么是问题。

      1. 感谢回复。我思考了下,应该是我哈密顿量的结构问题,波戈留波夫哈密顿量类似于超导中电子空穴对的哈密顿量。在对角有两个互为共轭的block矩阵,逆对角线上是0,最后导致本征矢对应部分也是0。这种情况可能不能直接用高效法计算。

发表回复

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