Mathematica数值求解TB模型代码汇总

闲谈

本文介绍Mathematica数值求解TB模型，主要内容有

• 哈密顿量
• 正格子和倒格子
• 高对称路径能带
• 三维能带和费米面
• 贝里曲率和陈数
•      

哈密顿量

我们以Haldane模型为例，具体应用时请替换成自己的模型对应的哈密顿量。

k={kx,ky};
a=1;
aa1={0,a};aa2={-(Sqrt[3]/2)a,-(a/2)};aa3={Sqrt[3]/2 a,-(a/2)};
bb1=aa2-aa3;bb2=aa3-aa1;bb3=aa1-aa2;
{bb1,bb2,bb3}
{\[Sigma]0,\[Sigma]1,\[Sigma]2,\[Sigma]3}={PauliMatrix[0],PauliMatrix[1],PauliMatrix[2],PauliMatrix[3]};
\[Epsilon]=2*t2*Cos[\[Phi]]*(Cos[k . bb1]+Cos[k . bb2]+Cos[k . bb3]);
dx=t1*(Cos[k . aa1]+Cos[k . aa2]+Cos[k . aa3]);
dy=t1*(Sin[k . aa1]+Sin[k . aa2]+Sin[k . aa3]);
dz=M-2*t2*Sin[\[Phi]]*(Sin[k . bb1]+Sin[k . bb2]+Sin[k . bb3]);
H=\[Epsilon]*\[Sigma]0+dx*\[Sigma]1+dy*\[Sigma]2+dz*\[Sigma]3;
H//TraditionalForm

正格子和倒格子

已知正格子坐标，可以得到倒格子坐标，把它们画到一起方便写出高对称点坐标。

a=1;
a1={-Sqrt[3],0,0};a2={Sqrt[3]/2,-(3/2),0};a3={0,0,1};
{a1,a2,a3}(*正格子坐标*)
\[CapitalOmega]=a1 . Cross[a2,a3];
b1=2\[Pi] Cross[a2,a3]/\[CapitalOmega];b2=2\[Pi] Cross[a3,a1]/\[CapitalOmega];b3=2\[Pi] Cross[a1,a2]/\[CapitalOmega];
{b1,b2,b3}(*倒格子坐标*)
(*二维正格子和倒格子对应的第一布里渊区示意图*)
latticeVectorFig[x_,text_,color_]:=Graphics[{color,Arrow[{{0,0},x[[1;;2]]}],Style[Text[text,x[[1;;2]]*1.1],FontSize->20,color]}]
a1fig=latticeVectorFig[a1,"\!\(\*SubscriptBox[\(a\), \(1\)]\)",Green];a2fig=latticeVectorFig[a2,"\!\(\*SubscriptBox[\(a\), \(2\)]\)",Green];
b1fig=latticeVectorFig[b1,"\!\(\*SubscriptBox[\(b\), \(1\)]\)",Purple];b2fig=latticeVectorFig[b2,"\!\(\*SubscriptBox[\(b\), \(2\)]\)",Purple];
(*RegularPolygon[{x,y},{r,\[Theta]},n]*)
directLatticefig=Graphics[{EdgeForm[Green],FaceForm[None],RegularPolygon[{0,0},{1,30 Degree},6]},Axes->True];
reciprocalLatticefig=Graphics[{EdgeForm[Purple],FaceForm[None],RegularPolygon[{0,0},{(4\[Pi])/(3Sqrt[3]),0 Degree},6]},Axes->True];
Show[a1fig,a2fig,directLatticefig,b1fig,b2fig,reciprocalLatticefig,Axes->True,AxesLabel->{"x","y"}]

高对称路径能带

在《紧束缚模型的二次量子化表示和推导求解》中给出了沿着高对称路径画能带的代码，这里对其进行整理，只需要在HSPBand[{{A,B,num},...}]中指定高对称点坐标即可，用起来更加方便。

ktrans[x_]:=x[[1]]*b1+x[[2]]*b2+x[[3]]*b3;(*以笛卡尔坐标为基矢*)
kpoints[start_,end_,npoints_]:=Table[start+i/npoints (end-start),{i,0,npoints-1}];(*撒点*)
HSPBand[kpt2kpt_]:=Module[{kpath1,kpath2,kpath3,kk1,kk2,KK,energy,nbands,nkpts,band},
kpath1=Flatten[Table[kpoints[ktrans[kpt2kpt[[j,1]]],ktrans[kpt2kpt[[j,2]]],kpt2kpt[[j,3]]],{j,1,Length[kpt2kpt]}],1]//N;
AppendTo[kpath1,ktrans[Last[kpt2kpt][[2]]]];
kpath2=Table[Norm[kpath1[[i+1]]-kpath1[[i]]],{i,1,Length[kpath1]-1}];
kpath3=Join[{0},Accumulate[kpath2]];
kk1={"K","\[CapitalGamma]","M","K"};(*\[LightBulb] 高对称路径*)
kk2=Table[kpt2kpt[[j,3]],{j,1,Length[kpt2kpt]}];
PrependTo[kk2,1];
kk2=Accumulate[kk2];
kk2=Table[kpath3[[i]],{i,kk2}];
KK=Transpose[{kk2,kk1}];
(*求解本征值*)
energy=Table[Sort[Eigenvalues[Htb[kpath1[[i]]]]],{i,1,Length[kpath1]}];
(*处理数据*)
nbands=Length[Htb[kpath1[[1]]]];(*能带条数(矩阵维度)*)
nkpts=Length[kpath1];(*每条能带上的k点数*)
band=Table[Transpose[{kpath3,energy[[All,i]]}],{i,1,nbands}];
(*画图*)
ListLinePlot[Table[band[[i]][[All,1;;2]],{i,1,nbands}],
Axes->False,GridLines->{Transpose[KK][[1]],{0}},
GridLinesStyle->Directive[Orange, Dashed],
Frame->True,FrameTicks->{{Automatic,None},{KK,None}},
FrameStyle->Directive[Purple],
FrameLabel->{"Wave vector","Energy(eV)"},
AspectRatio->5/4,
PlotStyle->{Red,Blue},
PlotRange->{{KK[[1,1]],KK[[Length[kk1],1]]},{All,All}}]
]
(*\[LightBulb] 哈密顿量*)
Htb[{x_,y_,z_}]:=H/.{kx->x,ky->y,kz->z,M->0.1,t1->1,t2->0.2,\[Phi]->0.2*\[Pi]}
a=1;
a1={-Sqrt[3],0,0};a2={Sqrt[3]/2,-(3/2),0};a3={0,0,1};(*\[LightBulb] 正格子基矢*)
\[CapitalOmega]=a1 . Cross[a2,a3];
b1=2\[Pi] Cross[a2,a3]/\[CapitalOmega];b2=2\[Pi] Cross[a3,a1]/\[CapitalOmega];b3=2\[Pi] Cross[a1,a2]/\[CapitalOmega];
G0={0,0,0};K0={1/3,1/3,0};M0={1/2,0,0};(*\[LightBulb] 以倒格子晶格矢量为基矢的高对称点坐标*)
HSPBand[{{K0,G0,40},{G0,M0,40},{M0,K0,20}}]

三维能带和费米面

这里给出在直角坐标系下撒点得到的三维能带和费米面，当然，也可以在倒空间撒点。

(*\[LightBulb] 哈密顿量*)
hamiltonian[kx0_,ky0_]:=H/.{kx->kx0,ky->ky0,M->0.1,t1->1,t2->0.2,\[Phi]->0.2*\[Pi]}
n=101;(*\[LightBulb] 撒点数*)
kx0=Array[#&,n,{-\[Pi],\[Pi]}]//N;
ky0=Array[#&,n,{-\[Pi],\[Pi]}]//N;
eigensystem=Table[Transpose[Eigensystem[hamiltonian[kx0[[i]],ky0[[j]]]]//N//Chop]//Sort,{i,1,n},{j,1,n}];
vals=Table[eigensystem[[i,j]][[All,1]],{i,1,n},{j,1,n}];(*本征值*)
vecs=Table[eigensystem[[i,j]][[All,2]],{i,1,n},{j,1,n}];(*本征态*)
(*三维能带*)
a1=Table[{kx0[[i]],ky0[[j]],vals[[i,j,1]]},{i,1,n},{j,1,n}];
a2=Table[{kx0[[i]],ky0[[j]],vals[[i,j,2]]},{i,1,n},{j,1,n}];
ListPlot3D[{Flatten[a1,1],Flatten[a2,1]},Mesh->None,ColorFunction->"Rainbow"]
(*费米面*)
{ListContourPlot[Flatten[a1,1]],ListContourPlot[Flatten[a2,1]]}

贝里曲率和陈数

对于贝里曲率和陈数的计算，这里我们在倒空间撒点，避免了六角格子第一布里渊不好确定的麻烦。

下面的代码覆盖了完整的4个第一布里渊区，因此陈数等于-4，贝里曲率也很好的验证了结果。

(*数值求解*)
nx=50;(*\[LightBulb] kx方向格子数目*)
ny=50;(*\[LightBulb] ky方向格子数目*)
BZb1=2*b1;
BZb2=2*b2;
(*BZb1={-2\[Pi],0,0};
BZb2={0,-2\[Pi],0};*)
BZ[na_,nb_]:=na*BZb1+nb*BZb2;(*以倒格子BZa和BZb为单位定义函数*)
(*kx方向取整个布里渊区，ky方向取半个布里渊区*)
kk=Table[BZ[(i-1)/nx,(j-1)/ny]-BZ[1/2,1/2],{i,1,nx+1},{j,1,ny+1}];
(*撒点示意图*)
temp=Transpose[Transpose[Flatten[kk,1]][[1;;2]]];
tfig=ListPlot[temp,PlotStyle->Red];
Show[b1fig,b2fig,reciprocalLatticefig,tfig,Axes->True,AxesLabel->{"x","y"}]
(*\[LightBulb] 哈密顿量*)
Htb[{x_,y_,z_}]:=H/.{kx->x,ky->y,kz->z,M->0.1,t1->1,t2->0.2,\[Phi]->0.2*\[Pi]}
(*\[LightBulb] 需要计算陈数的能带*)
bnd=Table[i,{i,1,1}];
(*定义函数计算要求的能带对应的本征态*)
vectors[{kx_,ky_,kz_},nband_]:=Module[{eigensystem},
eigensystem=Transpose[Eigensystem[Htb[{kx,ky,kz}]//N]//Chop]//Sort;
Table[{k,eigensystem[[nband[[k]],2]]},{k,1,Length[nband]}]
]
(*创建矩阵索引*)
vec=ConstantArray[0,{nx+1,ny+1}];
(*计算Uab*)
Table[vec[[i,j]]=vectors[kk[[i,j]],bnd],{i,1,nx+1},{j,1,ny+1}];
(*贝里曲率F(k)和贝里联络A(k)*)
uab[vec1_,vec2_,nband_]:=Module[{k1,k2},
Det[Table[Conjugate[vec1[[k1,2]]] . vec2[[k2,2]],
{k1,1,Length[nband]},{k2,1,Length[nband]}]]
]
u12=Table[uab[vec[[i,j]],vec[[i+1,j]],bnd],{i,1,nx},{j,1,ny}];
u23=Table[uab[vec[[i+1,j]],vec[[i+1,j+1]],bnd],{i,1,nx},{j,1,ny}];
u34=Table[uab[vec[[i+1,j+1]],vec[[i,j+1]],bnd],{i,1,nx},{j,1,ny}];
u41=Table[uab[vec[[i,j+1]],vec[[i,j]],bnd],{i,1,nx},{j,1,ny}];
Fk=Table[Im[Log[u12[[i,j]]*u23[[i,j]]*u34[[i,j]]*u41[[i,j]]]],{i,1,nx},{j,1,ny}];
(*贝里曲率*)
phi=Fk;
BerryCurvature=Flatten[Table[Append[kk[[i,j]][[1;;2]],phi[[i,j]]],{i,1,nx},{j,1,ny}],1];
(*陈数*)
ChernNumber=Total[Flatten[Fk]]/(2\[Pi]);
Print["Chern number = " <> ToString[ChernNumber,TraditionalForm]]
(*三维图*)
ListPlot3D[BerryCurvature,PlotRange->Full,PlotLegends->Automatic,ColorFunction->Hue,AspectRatio->Automatic]
(*二维图*)
ListContourPlot[BerryCurvature,PlotRange->Full,PlotLegends->Automatic,ColorFunction->"Rainbow",AspectRatio->Automatic]
(*二维图*)
ListDensityPlot[BerryCurvature,PlotLegends->Automatic,PlotRange->Full,ColorFunction->"Rainbow",AspectRatio->Automatic]

对于    的计算，这里对结果做一些处理，使曲线平滑。

(*给出phi最小值在ky中的位置*)
position=Position[phi,Min[phi]][[1,1]];
(*处理数据*)
KY=Table[i,{i,1,ny}];
PHI=Table[Transpose[phi][[i]],{i,1,nc}];
PrependTo[PHI,KY];
PHI=Transpose[PHI];
(*图像有拐点,分两段处理*)
plotData=Flatten[{Table[PHI[[All,{1,i+1}]][[1;;position]],{i,1,nc}],Table[PHI[[All,{1,i+1}]][[position+1;;ny]],{i,1,nc}]},1];
(*画图*)
ListPlot[plotData,
Axes->False,AspectRatio->3/4,PlotStyle->Directive[Blue,PointSize[Tiny]],
Frame->True,FrameStyle->Directive[Purple],Joined->True,
FrameLabel->{Style["\!\(\*SubscriptBox[\(k\), \(y\)]\)/2\[Pi]",FontFamily->"Times New Roman",16],Style["\[Phi]/2\[Pi]",FontFamily->"Times New Roman",16]},
FrameTicks->{{{{-\[Pi],-0.5},{0,0},{\[Pi],0.5}},None},{{{0,0},{ny,0.5}},None}},
PlotRange->{{0,ny},{-\[Pi],\[Pi]}}]