matlab化石墨烯,锯齿石墨烯磁性hubbard matlab计算程序

主程序

Nx=3;

Ny=50;

n=Nx*Ny;

[x,y]=zigzagP(Ny,Nx);

t1=-2.7;

H=Hamiltonian_NN_graphene(x,y,t1);

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

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

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

x=x(n/3+1:2/3*n,1);

y=y(n/3+1:2/3*n,1);

m=length(HD); % 体系总的原子个数

Nupavg=0.5*eye(m,m); %将所有的取为1作初始值,自旋向上

Ndownavg=0.5*eye(m,m); %将所有的取为1作初始值，自旋向下

u=0; %费米能级

Ncc=20; %迭代次数

dk=0.01; %布里渊区步长

T=0.038; %温度

U=2; %Hubbard系数

k=0:dk:2*pi; %布里渊区路径

nk=length(k); % 布里渊区离散点个数

band_up=zeros(m,nk); % 记录最后自旋向上的能带

band_down=zeros(m,nk); % 记录最后自旋向下的能带

updensity=zeros(Ncc,m); %储存每一步迭代上自旋电子密度，观察是否收敛

downdensity=zeros(Ncc,m); % 储存每一步迭代下自旋电子密度，观察是否收敛

Nupavg(1:Ny:m,1:Ny:m)=0.65; %上边界上自旋电子密度

Ndownavg(1:Ny:m,1:Ny:m)=0.35; %上边界下自旋电子密度

Nupavg(Ny:Ny:m,Ny:Ny:m)=0.35; %下边界上自旋电子密度

Ndownavg(Ny:Ny:m,Ny:Ny:m)=0.65; %下边界下自旋电子密度

%%%%%%%%%%% 参数设置完成%%%%%%%%%%%

for nc=1:Ncc

Nupavg0=zeros(m,m); %初始化Nupavg的临时变量用于k空间的积分

Ndownavg0=zeros(m,m); %初始化Ndownavg的临时变量用于k空间的积分

for i=1:nk

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

[Adown,Edown]=eig(Hk+U*Nupavg-0.5*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量

[Aup,Eup]=eig(Hk+U*Ndownavg-0.5*U*eye(m,m)); %通过Ndownavg求解出Aup,Aup为up波函数系数向量

% 求电子密度的向量化写法

f_up=Fermi_funtion(diag(Eup),u,T)';

f_down=Fermi_funtion(diag(Edown),u,T)';

fermi_up = f_up(ones(m,1),:);

fermi_down = f_down(ones(m,1),:);

Nup = sum(Aup.*conj(Aup).*fermi_up,2);

Ndown = sum(Adown.*conj(Adown).*fermi_down,2);

Nupavg0=Nupavg0+diag(Nup)*dk;

Ndownavg0=Ndownavg0+diag(Ndown)*dk;

end

Nupavg=Nupavg0/(2*pi); %更新Nupavg

Ndownavg=Ndownavg0/(2*pi); %更新Ndownavg

updensity(nc,:)=diag(Nupavg);

downdensity(nc,:)=diag(Ndownavg);

if mod(nc,5)==0

fprintf('已完成第%d次迭代\n',nc)

end

end

%%%%%%%%%%%%%%%%迭代完成，开始计算能带%%%%%%%%%%%%%%

for i=1:nk

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

[Adown,Edown]=eig(Hk+U*Nupavg-0.5*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量

[Aup,Eup]=eig(Hk+U*Ndownavg-0.5*U*eye(m,m)); %通过Ndownavg求解出Aup,Aup为up波函数系数向量

band_up(:,i) = diag(Eup);

band_down(:,i) = diag(Edown);

end

plot(k,band_down,'b.')

hold on

plot(k,band_up,'r.')

set(gca,'YLim',[-2.7 2.7]);

set(gca,'XLim',[0 2*pi]);

生成体系座标程序

function [x,y]=zigzagP(N,number);

x=zeros(N*number,1);

y=zeros(N*number,1);

x(1)=sqrt(3)/2;

y(1)=0;

x(2)=0;

y(2)=-0.5;

for j=1:N-2

y(j+2)=y(j)-1.5;

if mod(j+2,4)==1||mod(j+2,4)==0

x(j+2)=sqrt(3)/2;

else

x(j+2)=0;

end

end

for j=1:number-1

for l=1:N

x(j*N+l)=x(l)+j*sqrt(3);

y(j*N+l)=y(l);

end

end

最近邻哈密顿量与前面程序相同

费米函数

function f=Fermi_funtion(E,u,T)

f=(exp((E-u)/T)+ones(length(E),1)).^-1;