今天看到有些同学在求Chern number的计算方法，正好前段时间我有花时间做过相关的工作，本着科研共享的精神，来分享一下我的方法。
由于我不是专业做这个方向的研究人员，所以我就不能做更多理论方面的解释了，直接给大家一步步说明我的算法就好了。0.0
准备工具：COMSOL +MATLAB
参考文献：PRL 100,013905 (2008) + J.phys.Soc.Jpn. vol,74 no,6 (2005)
步骤
利用COMSOL求解光子晶体结构的本征频率：把第一布里渊区等分成n1*n2个点，点命名为N(ki,kj)（i=1:n1,j=1:n2）好了，于是可以得到n1*n2个本征频率解，保存他们的电场或磁场分布(按介质的不同保存，比如说是Si和Air就分别保存两个区域的场分布)。特别注意的是由于求解器的问题有可能有些结构在算本征频率的时候，同一条能带的解的位置会不一样。例如我想求4条能带的Chern number，但是N14和N20（只是举例）的解并不是正好跟能带位置对应的地方（即解的位置不等于能带位置的话会影响后面的计算）。所以建议每次求解完之后简单的在COMSOL中把这些解放在图中画一下，仔细比对一下（不用把n1*n2个点都画出来，只用画简约布里渊区边界的那些点就足够精确了（从效率上来说，当然时间够有兴趣可以导出来在其他工具中画出三维的（k1,k2,freq）图来确认））。
得到了n1*n2个本征频率下的电场/磁场强度分布后导出保存为EH1.txt+EH2.txt(两种介质)文件：（注意每个介质的实部虚部的导出，我是比较麻烦每次分实部虚部导出两个文件（总共4个文件）再相加的）这里多说一句，对于那些对COMSOL不是很熟悉的人，首先COMSOL导出的电场数据我假设是4条能带11*11=121个点即会有484个不同的场分布，每个场分布中假设有10000个点（这个跟网格划分有关），则我导出的数据为10000*484的矩阵，加上COMSOL的位置（x,y）信息则导出数据为10000*486的矩阵。
在得到EH1.txt+EH2.txt文件之后下一步就是对这两个10000*486的矩阵进行插值处理了：首先在进行插值之前要把各个能带之间的数据分割开来，在这个案例中能带数为4，则对3：end列的数据每隔4个为一个能带的场分布（即3,7,11，...为第一能带，4,8,12，...为第二能带）命名为Eh1,Eh2,Eh3,Eh4（Air）+Eh5,Eh6,Eh7,Eh8(Si)。（Eh为10000*121的矩阵对应了121个不同的点N（ki,kj））。
以求第一个能带的chern number为例。首先对Eh1,Eh5进行插值，griddata命令即可。插值以后的场分布分别称为eH1和eH5。
计算Chern number
理论上的东西大家可以参考上面推荐的两个文献资料，计算过程简单来说是，
：
做面积分，
同理，
：
做面积分，
：
做面积分，
：
做面积分，
然后
，
最后
Chern number 。
很多不够准确的地方大家见谅，如果有什么错误后续会改，希望大家科研顺利。合作共享是科研的出路！
有问题的同学可以去我某乎上的文章找一些答案，网址
