matlab计算石墨烯能带结构的算法是什么

这个程序是初步优化后的matlab版本，主要思路是先生成体系的格点坐标，再运用坐标生成体系的哈密顿量，然后进行对角化计算能带，能带的计算使用一维体系超原胞的处理方法。可以进一步优化

主程序

nx=3%

ny=100% 体系宽度(y方向的长度)

[x,y]=zigzag_graphene(nx,ny)

%plot(x,y,'.','MarkerSize',20)

t1=-2.7

t2=0.0038/3/sqrt(3)

H=Hamiltonian_NN_graphene(x,y,t1)

Hsp=Hamiltonian_Haldane(x,y,sqrt(3),t2)

H=H+Hsp

N=length(H)

HDL=H(N/3+1:N*2/3,1:N/3)

HD=H(N/3+1:N*2/3,N/3+1:N*2/3)

HDR=H(N/3+1:N*2/3,N*2/3+1:N)

n = length(HD)

dk = 0.01

kx=0:dk:2*pi% k空间路径

Ek=band_calculate(kx,HD,HDL,HDR)

plot(kx,Ek,'.')

set(gca,'YLim',[-0.5 0.5])%X轴的数据显示范围

坐标生成函数

function [x,y]=zigzag_graphene(nx,ny)

x1=zeros(4,1)

y1=zeros(4,1)

x1(1,1)=sqrt(3)/2

x1(2,1)=0

x1(3,1)=0

x1(4,1)=sqrt(3)/2

y1(1,1)=0

y1(2,1)=0.5

y1(3,1)=1.5

y1(4,1)=2

x2=x1

y2=y1

for i=1:ny-1

x2=[x2x1]

y2=[y2y1+i*ones(4,1)*3]

end

x=x2

y=y2

n=length(x2)

for i=1:nx-1

x=[xx2+i*ones(n,1)*sqrt(3)]

y=[yy2]

end

最近邻相互作用哈密顿量的生成

function H=Hamiltonian_NN_graphene(x,y,t)

%t=-2.7

N=length(x)

H=zeros(N,N)

eps=0.01

for i=1:N

for j=1:N

if abs(sqrt((x(i)-x(j))^2+(y(i)-y(j))^2)-1)

Haldane模型哈密顿量

function H=Hamiltonian_Haldane(x,y,a,t2)

N=length(x)

H=zeros(N,N)

for l=1:N

for j=1:N

if x(l)>x(j)&&y(l)==y(j)&&mod(j,2)==1&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=-1i*t2

end

if x(l)x(j)&&y(l)>y(j)&&mod(j,2)==1&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=1i*t2

end

if x(l)y(j)&&mod(j,2)==1&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=-1i*t2

end

if x(l)>x(j)&&y(l)x(j)&&y(l)==y(j)&&mod(j,2)==0&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=1i*t2

end

if x(l)x(j)&&y(l)>y(j)&&mod(j,2)==0&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=-1i*t2

end

if x(l)y(j)&&mod(j,2)==0&&abs(sqrt((x(j)-x(l))^2+(y(j)-y(l))^2)-a)<0.001

H(j,l)=1i*t2

end

if x(l)>x(j)&&y(l)

能带计算函数

function Ek=band_calculate(kx,HD,HDL,HDR)

dN = length(kx)

n = length(HD)

Ek = zeros(n,dN)

for i = 1:dN

Hk=HDL*exp(-1i*kx(i))+HD+HDR*exp(1i*kx(i))

[~,E]=eig(Hk)

Ek(:,i) = diag(E)

end

用第一原理计算软件开展的工作，分析结果主要是从以下三个方面进行定性/定量的讨论： 1、电荷密度图（charge density）； 2、能带结构（Energy Band Structure）； 3、态密度（Density of States，简称DOS）。电荷密度图是以图的形式出现在文章中，非常直观，因此对于一般的入门级研究人员来讲不会有任何的疑问。唯一需要注意的就是这种分析的种种衍生形式，比如差分电荷密图（def-ormation charge density）和二次差分图（difference charge density）等等，加自旋极化的工作还可能有自旋极化电荷密度图（spin-polarized charge density）。所谓“差分”是指原子组成体系（团簇）之后电荷的重新分布，“二次”是指同一个体系化学成分或者几何构型改变之后电荷的重新分布，因此通过这种差分图可以很直观地看出体系中个原子的成键情况。通过电荷聚集（accumulation）/损失（depletion）的具体空间分布，看成键的极性强弱；通过某格点附近的电荷分布形状判断成键的轨道（这个主要是对d轨道的分析，对于s或者p轨道的形状分析我还没有见过）。分析总电荷密度图的方法类似，不过相对而言，这种图所携带的信息量较小。能带结构分析现在在各个领域的第一原理计算工作中用得非常普遍了。但是因为能带这个概念本身的抽象性，对于能带的分析是让初学者最感头痛的地方。关于能带理论本身，我在这篇文章中不想涉及，这里只考虑已得到的能带，如何能从里面看出有用的信息。首先当然可以看出这个体系是金属、半导体还是绝缘体。判断的标准是看费米能级和导带（也即在高对称点附近近似成开口向上的抛物线形状的能带）是否相交，若相交，则为金属，否则为半导体或者绝缘体。对于本征半导体，还可以看出是直接能隙还是间接能隙：如果导带的最低点和价带的最高点在同一个k点处，则为直接能隙，否则为间接能隙。在具体工作中，情况要复杂得多，而且各种领域中感兴趣的方面彼此相差很大，分析不可能像上述分析一样直观和普适。不过仍然可以总结出一些经验性的规律来。主要有以下几点： 1） 因为目前的计算大多采用超单胞（supercell）的形式，在一个单胞里有几十个原子以及上百个电子，所以得到的能带图往往在远低于费米能级处非常平坦，也非常密集。原则上讲，这个区域的能带并不具备多大的解说/阅读价值。因此，不要被这种现象吓住，一般的工作中，我们主要关心的还是费米能级附近的能带形状。 2） 能带的宽窄在能带的分析中占据很重要的位置。能带越宽，也即在能带图中的起伏越大，说明处于这个带中的电子有效质量越小、非局域（non-local）的程度越大、组成这条能带的原子轨道扩展性越强。如果形状近似于抛物线形状，一般而言会被冠以类sp带（sp-like band）之名。反之，一条比较窄的能带表明对应于这条能带的本征态主要是由局域于某个格点的原子轨道组成，这条带上的电子局域性非常强，有效质量相对较大。 3） 如果体系为掺杂的非本征半导体，注意与本征半导体的能带结构图进行对比，一般而 言在能隙处会出现一条新的、比较窄的能带。这就是通常所谓的杂质态(doping state)，或者按照掺杂半导体的类型称为受主态或者施主态。 4） 关于自旋极化的能带，一般是画出两幅图：majority spin和minority spin。经典的说，分别代表自旋向上和自旋向下的轨道所组成的能带结构。注意它们在费米能级处的差异。如果费米能级与majority spin的能带图相交而处于minority spin的能隙中，则此体系具有明显的自旋极化现象，而该体系也可称之为半金属（half metal）。因为majority spin与费米能级相交的能带主要由杂质原子轨道组成，所以也可以此为出发点讨论杂质的磁性特征。 5） 做界面问题时，衬底材料的能带图显得非常重要，各高对称点之间有可能出现不同的情况。具体地说，在某两点之间，费米能级与能带相交；而在另外的k的区间上，费米能级正好处在导带和价带之间。这样，衬底材料就呈现出各项异性：对于前者，呈现金属性，而对于后者，呈现绝缘性。因此，有的工作是通过某种材料的能带图而选择不同的面作为生长面。具体的分析应该结合试验结果给出。（如果我没记错的话，物理所薛其坤研究员曾经分析过$\beta$-Fe的(100)和(111)面对应的能带。有兴趣的读者可进一步查阅资料。）原则上讲，态密度可以作为能带结构的一个可视化结果。很多分析和能带的分析结果可以一一对应，很多术语也和能带分析相通。但是因为它更直观，因此在结果讨论中用得比能带分析更广泛一些。简要总结分析要点如下： 1） 在整个能量区间之内分布较为平均、没有局域尖峰的DOS，对应的是类sp带，表明电子的非局域化性质很强。相反，对于一般的过渡金属而言，d轨道的DOS一般是一个很大的尖峰，说明d电子相对比较局域，相应的能带也比较窄。 2） 从DOS图也可分析能隙特性：若费米能级处于DOS值为零的区间中，说明该体系是半导体或绝缘体；若有分波DOS跨过费米能级，则该体系是金属。此外，可以画出分波（PDOS）和局域（LDOS）两种态密度，更加细致的研究在各点处的分波成键情况。 3） 从DOS图中还可引入“赝能隙”（pseudogap）的概念。也即在费米能级两侧分别有两个尖峰。而两个尖峰之间的DOS并不为零。赝能隙直接反映了该体系成键的共价性的强弱：越宽，说明共价性越强。如果分析的是局域态密度（LDOS），那么赝能隙反映的则是相邻两个原子成键的强弱：赝能隙越宽，说明两个原子成键越强。上述分析的理论基础可从紧束缚理论出发得到解释：实际上，可以认为赝能隙的宽度直接和Hamiltonian矩阵的非对角元相关，彼此间成单调递增的函数关系。 4） 对于自旋极化的体系，与能带分析类似，也应该将majority spin和minority spin分别画出，若费米能级与majority的DOS相交而处于minority的DOS的能隙之中，可以说明该体系的自旋极化。

%我想用matlab对函数dy/dx=（x-y）/（1-y-x）进行模拟，请问该怎样做，

%解微分方程！！

%归一化：

%令y=y（1）

%x=t=y(2)

%dy(1)/dt=(y(2)-y(1))/(1-y(1)-y(2))

%dy(2)/dt=1

%函数文件

founction dy=fun1(t,y)

dy=zeros(2,1)

dy=[(y(2)-y(1))/(1-y(1)-y(2))1]

%以上保存为fun1.m文件

%以下是脚本程序

clear

ts=0:0.01:?%时间范围

y0=[??]%函数边界条件y(0),x(0)

[t,y]=ode45('fun1',ts,y0)

plot(y(:,1),y(:,2))

---