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

2024年6月10日 202次阅读

主程序

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);  %将所有的<Nupi>取为1作初始值,自旋向上
Ndownavg=0.5*eye(m,m); %将所有的<Ndowni>取为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;

 

点赞