本篇主要参考：
- 陈数Chern number的计算（高效法，附Python/Matlab代码）
- 陈数Chern number的计算（Wilson loop方法，附Python代码）
- 陈数Chern number的计算（多条能带的Wilson loop方法，附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)可以看出结果和之前的相同，但虚部会小很多。
【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】
关老师，这个高效法是不是无法应用到石墨烯简并能带计算chern数。完全按您的haldan模型中的定义，在简并的情况下，算出的数也是极其古怪。即使只有石墨烯跃迁项再加上一个打开能隙的项，也算不出零。将您的程序改成Fortran后，单带的情况计算结果一样，加上简并，求行列式，结果就错了。
应该是可以计算的。可以检查积分/求和时布里渊区是否选对，参考这篇：Haldane模型中陈数的计算（附Python代码）。
关老师，已经将代码的布里渊区的积分范围更改了，但是计算出来的数值还是不正确。价带： (10.502001226218756-1.9712108027775957e-11j)价带（包含两个交叉能带）： (-2.8720927329182986e-15-8.164510371939192e-11j)所有能带： (-2.8720927329182986e-15-8.164510371939192e-11j)。您能帮我看看嘛/跪求，呜呜呜~
你确定你这个哈密顿量的布里渊区范围是这个吗？即使同样是蜂窝格子，数值范围可能也是有区别的。建议检查下布里渊区的大小，然后也画个能带来验证一下。
我也是用这个方法计算石墨烯模型出现问题，即使是只有跃迁项加上自旋后（二重简并），也算的不对。请问您的问题解决了吗？
关老师，我将您上面计算方晶格多能带的高效法改成计算具有一个原胞4个格点的Haldane模型（哈密顿量是四阶），计算出来的数值不正确
价带： (7.4477953185066434+8.817204591573345e-16j)
价带（包含两个交叉能带）：(-8.780877525850066e-16-4.423798166951144e-13j)
所有能带：(-8.780877525850066e-16-4.423798166951144e-13j)。您能帮我看看是代码哪里有问题嘛？
布里渊区不再是方形的了，可能是积分区域的问题。可以参考这篇：Haldane模型中陈数的计算（附Python代码）。
关老师，请问处理简并或交叉能带的思路是什么呢？
公式类似于这个：https://topocondmat.org/w4_haldane/ComputingChern.html。需要求多带波函数内积后的矩阵行列式。严格的证明我也没推导过。
关老师，我是改变模型的磁场强度来算其陈数，但是算出来的陈数有些是0，但有很多都是小数，例如-0.2260，-0.2594，0.3255等，请问这可能是什么原因呢？我的kx ky步长是0.005
有简并或者交叉的能带要一起算，不能分开算。