拓扑不变量, 学术

陈数Chern number的计算(多条能带的高效法,附Python代码)

本篇主要参考:

高效法和Wilson loop方法在公式是接近的,只是高效法在内积时考虑了归一化系数,所以计算结果中虚部可以接近于零。参考多条能带的Wilson loop方法,自然而然地可以把高效法推广到多条能带的情况。自己目前没去找相关的文献和公式证明,只是对公式做了猜想,并在数值上做了计算验证,不保证公式或算法的正确性,但大概率是正确的,仅作参考。本篇的思路主要来源于Song M.R.同学以及共同讨论。

这里采用和这篇“陈数Chern number的计算(多条能带的Wilson loop方法,附Python代码)”中相同的例子,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/25107
"""

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


def hamiltonian(kx, ky, Ny, B):
    h00 = np.zeros((Ny, Ny), dtype=complex)
    h01 = np.zeros((Ny, Ny), dtype=complex)
    t = 1
    for iy in range(Ny-1):
        h00[iy, iy+1] = t
        h00[iy+1, iy] = t
    h00[Ny-1, 0] = t*cmath.exp(1j*ky)
    h00[0, Ny-1] = t*cmath.exp(-1j*ky)
    for iy in range(Ny):
        h01[iy, iy] = t*cmath.exp(-2*np.pi*1j*B*iy)
    matrix = h00 + h01*cmath.exp(1j*kx) + h01.transpose().conj()*cmath.exp(-1j*kx)
    return matrix


def main():
    Ny = 20

    H_k = functools.partial(hamiltonian, Ny=Ny, B=1/Ny)

    chern_number = calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(H_k, index_of_bands=range(int(Ny/2)-1))
    print('价带:', chern_number)
    print()

    chern_number = calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(H_k, index_of_bands=range(int(Ny/2)+2))
    print('价带(包含两个交叉能带):', chern_number)
    print()

    chern_number = calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(H_k, index_of_bands=range(Ny))
    print('所有能带:', chern_number)

    # 函数可通过Guan软件包调用。安装方法:pip install --upgrade guan
    # import guan
    # chern_number = guan.calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0)


def calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0): 
    delta = 2*math.pi/precision
    chern_number = 0
    for kx in np.arange(-math.pi, math.pi, delta):
        if print_show == 1:
            print(kx)
        for ky in np.arange(-math.pi, math.pi, delta):
            H = hamiltonian_function(kx, ky)
            eigenvalue, vector = np.linalg.eigh(H) 
            H_delta_kx = hamiltonian_function(kx+delta, ky) 
            eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) 
            H_delta_ky = hamiltonian_function(kx, ky+delta)
            eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) 
            H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta)
            eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky)
            dim = len(index_of_bands)
            det_value = 1
            # first dot product
            dot_matrix = np.zeros((dim , dim), dtype=complex)
            i0 = 0
            for dim1 in index_of_bands:
                j0 = 0
                for dim2 in index_of_bands:
                    dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2])
                    j0 += 1
                i0 += 1
            dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
            det_value = det_value*dot_matrix
            # second dot product
            dot_matrix = np.zeros((dim , dim), dtype=complex)
            i0 = 0
            for dim1 in index_of_bands:
                j0 = 0
                for dim2 in index_of_bands:
                    dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2])
                    j0 += 1
                i0 += 1
            dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
            det_value = det_value*dot_matrix
            # third dot product
            dot_matrix = np.zeros((dim , dim), dtype=complex)
            i0 = 0
            for dim1 in index_of_bands:
                j0 = 0
                for dim2 in index_of_bands:
                    dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2])
                    j0 += 1
                i0 += 1
            dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
            det_value = det_value*dot_matrix
            # fourth dot product
            dot_matrix = np.zeros((dim , dim), dtype=complex)
            i0 = 0
            for dim1 in index_of_bands:
                j0 = 0
                for dim2 in index_of_bands:
                    dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2])
                    j0 += 1
                i0 += 1
            dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
            det_value= det_value*dot_matrix
            chern_number += cmath.log(det_value)
    chern_number = chern_number/(2*math.pi*1j)
    return chern_number


if __name__ == '__main__':
    main()

运行结果:

价带: (-9.000000000000043-1.6279778224388367e-15j)

价带(包含两个交叉能带): (8.000000000000052-2.4257029167245346e-15j)

所有能带: (-5.694283434593287e-15-3.8760359718102057e-13j)

可以看出结果和之前的相同,但虚部会小很多。

3,798 次浏览

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

11 thoughts on “陈数Chern number的计算(多条能带的高效法,附Python代码)”

  1. 关老师,这个高效法是不是无法应用到石墨烯简并能带计算chern数。完全按您的haldan模型中的定义,在简并的情况下,算出的数也是极其古怪。即使只有石墨烯跃迁项再加上一个打开能隙的项,也算不出零。将您的程序改成Fortran后,单带的情况计算结果一样,加上简并,求行列式,结果就错了。

  2. 关老师,已经将代码的布里渊区的积分范围更改了,但是计算出来的数值还是不正确。价带: (10.502001226218756-1.9712108027775957e-11j)价带(包含两个交叉能带): (-2.8720927329182986e-15-8.164510371939192e-11j)所有能带: (-2.8720927329182986e-15-8.164510371939192e-11j)。您能帮我看看嘛/跪求,呜呜呜~

    import numpy as np
    import math
    from math import *
    
    import cmath
    import functools
    
    
    def hamiltonian(kx, ky, Ny, B):
    
        t2=0.1;
    
        gamma=0.446;
        m=0.35;
        matrix = np.zeros((4, 4))*(1+0j)
        matrix[0, 0] = m+np.complex(0, gamma)+2*t2*np.sin(kx)
        matrix[0, 1] =2*np.cos(kx/2)*cmath.exp(1j*ky/2/sqrt(3))
        matrix[0, 2] =-4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[0, 3] =cmath.exp(-1j*ky/sqrt(3))
        
        matrix[1, 0] = 2*np.cos(kx/2)*cmath.exp(-1j*ky/2/sqrt(3))
        matrix[1, 1] = -m+np.complex(0, gamma)-2*t2*np.sin(kx)
        matrix[1, 2] = cmath.exp(1j*ky/sqrt(3))
        matrix[1, 3] =4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        
        matrix[2, 0] = -4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[2, 1] = cmath.exp(-1j*ky/sqrt(3))
        matrix[2, 2] = m-np.complex(0, gamma)+2*t2*np.sin(kx)
        matrix[2, 3] =2*np.cos(kx/2)*cmath.exp(1j*ky/2/sqrt(3))
        
        matrix[3, 0] = cmath.exp(1j*ky/sqrt(3))
        matrix[3, 1] = 4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[3, 2] = 2*np.cos(kx/2)*cmath.exp(-1j*ky/2/sqrt(3))
        matrix[3, 3] =-m-np.complex(0, gamma)-2*t2*np.sin(kx)
        return matrix
    
    
    
    def main():
        Ny = 4
    
        H_k = functools.partial(hamiltonian, Ny=Ny, B=1/Ny)
    
        chern_number = calculate_chern_number(H_k, index_of_bands=range(int(Ny/2)-1))
        print('价带:', chern_number)
        print()
    
        chern_number = calculate_chern_number(H_k, index_of_bands=range(int(Ny/2)+2))
        print('价带(包含两个交叉能带):', chern_number)
        print()
    
        chern_number = calculate_chern_number(H_k, index_of_bands=range(Ny))
        print('所有能带:', chern_number)
    
      
    
    def calculate_chern_number(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0):   
        delta = 0.005  
        chern_number = 0 
        
        a = 1/sqrt(3)
        bb1 = 2*sqrt(3)*pi/3/a
        bb2 = 2*pi/3/a 
        for kx in np.arange(0 , bb1, delta):
            print(kx)
            for ky in np.arange(0, 2*bb2, delta):
                
                H = hamiltonian_function(kx, ky)
                eigenvalue, vector = np.linalg.eigh(H) 
                
                H_delta_kx = hamiltonian_function(kx+delta, ky) 
                eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx)
                
                H_delta_ky = hamiltonian_function(kx, ky+delta)
                eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky)
                
                H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta)
                eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky)
                dim = len(index_of_bands)
                det_value = 1
                # first dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # second dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # third dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # fourth dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value= det_value*dot_matrix
                chern_number += cmath.log(det_value)
        chern_number = chern_number/(2*math.pi*1j)
        return chern_number
    
    
    if __name__ == '__main__':
        main()

    1. 你确定你这个哈密顿量的布里渊区范围是这个吗?即使同样是蜂窝格子,数值范围可能也是有区别的。建议检查下布里渊区的大小,然后也画个能带来验证一下。

    2. 我也是用这个方法计算石墨烯模型出现问题,即使是只有跃迁项加上自旋后(二重简并),也算的不对。请问您的问题解决了吗?

  3. 关老师,我将您上面计算方晶格多能带的高效法改成计算具有一个原胞4个格点的Haldane模型(哈密顿量是四阶),计算出来的数值不正确
    价带: (7.4477953185066434+8.817204591573345e-16j)
    价带(包含两个交叉能带):(-8.780877525850066e-16-4.423798166951144e-13j)
    所有能带:(-8.780877525850066e-16-4.423798166951144e-13j)。您能帮我看看是代码哪里有问题嘛?

    import numpy as np
    import math
    from math import *
    import cmath
    import functools
    
    
    def hamiltonian(kx, ky, Ny, B):
    
        t2=0.1;
    
        gamma=0;
        m=0.5;
        matrix = np.zeros((4, 4))*(1+0j)
        matrix[0, 0] = m+np.complex(0, gamma)+2*t2*np.sin(kx)
        matrix[0, 1] =2*np.cos(kx/2)*cmath.exp(1j*ky/2/sqrt(3))
        matrix[0, 2] =-4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[0, 3] =cmath.exp(-1j*ky/sqrt(3))
        
        matrix[1, 0] = 2*np.cos(kx/2)*cmath.exp(-1j*ky/2/sqrt(3))
        matrix[1, 1] = -m+np.complex(0, gamma)-2*t2*np.sin(kx)
        matrix[1, 2] = cmath.exp(1j*ky/sqrt(3))
        matrix[1, 3] =4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        
        matrix[2, 0] = -4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[2, 1] = cmath.exp(-1j*ky/sqrt(3))
        matrix[2, 2] = m-np.complex(0, gamma)+2*t2*np.sin(kx)
        matrix[2, 3] =2*np.cos(kx/2)*cmath.exp(1j*ky/2/sqrt(3))
        
        matrix[3, 0] = cmath.exp(1j*ky/sqrt(3))
        matrix[3, 1] = 4*t2*np.sin(kx/2)*np.cos(sqrt(3)*ky/2)
        matrix[3, 2] = 2*np.cos(kx/2)*cmath.exp(-1j*ky/2/sqrt(3))
        matrix[3, 3] =-m-np.complex(0, gamma)-2*t2*np.sin(kx)
        return matrix
    
    
    
    def main():
        Ny = 4
    
        H_k = functools.partial(hamiltonian, Ny=Ny, B=1/Ny)
    
        chern_number = calculate_chern_number(H_k, index_of_bands=range(int(Ny/2)-1))
        print('价带:', chern_number)
        print()
        chern_number = calculate_chern_number(H_k, index_of_bands=range(int(Ny/2)+2))
        print('价带(包含两个交叉能带):', chern_number)
        print()
        chern_number = calculate_chern_number(H_k, index_of_bands=range(Ny))
        print('所有能带:', chern_number)
    
    def calculate_chern_number(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0): 
        delta = 2*math.pi/precision
        chern_number = 0
        for kx in np.arange(-math.pi, math.pi, delta):
            if print_show == 1:
                print(kx)
            for ky in np.arange(-math.pi, math.pi, delta):
                H = hamiltonian_function(kx, ky)
                eigenvalue, vector = np.linalg.eigh(H) 
                H_delta_kx = hamiltonian_function(kx+delta, ky) 
                eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) 
                H_delta_ky = hamiltonian_function(kx, ky+delta)
                eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) 
                H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta)
                eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky)
                dim = len(index_of_bands)
                det_value = 1
                # first dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # second dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # third dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value = det_value*dot_matrix
                # fourth dot product
                dot_matrix = np.zeros((dim , dim), dtype=complex)
                i0 = 0
                for dim1 in index_of_bands:
                    j0 = 0
                    for dim2 in index_of_bands:
                        dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2])
                        j0 += 1
                    i0 += 1
                dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
                det_value= det_value*dot_matrix
                chern_number += cmath.log(det_value)
        chern_number = chern_number/(2*math.pi*1j)
        return chern_number
    
    
    if __name__ == '__main__':
        main()
    

  4. 关老师,我是改变模型的磁场强度来算其陈数,但是算出来的陈数有些是0,但有很多都是小数,例如-0.2260,-0.2594,0.3255等,请问这可能是什么原因呢?我的kx ky步长是0.005

发表评论

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

Captcha Code