磷烯计算程序

生成坐标程序

function [x,y,z]=zigzag_phos(nx,ny)
x1=2.164*sin(98.15*pi/360);
y1=2.164*cos(98.15*pi/360);
z1=2.14;zy=sqrt(2.207^2-2.14^2);
x=zeros(1,4*ny*nx);
y=zeros(1,4*ny*nx);
z=zeros(1,4*ny*nx);
for m=1:nx
    for n=1:ny
        x(1+(n-1)*4+(m-1)*ny*4)=x1+(m-1)*2*x1;
        y(1+(n-1)*4+(m-1)*ny*4)=-(n-1)*(2*y1+2*zy);
        z(1+(n-1)*4+(m-1)*ny*4)=z1;
        x(2+(n-1)*4+(m-1)*ny*4)=(m-1)*2*x1;
        y(2+(n-1)*4+(m-1)*ny*4)=-y1-(n-1)*(2*y1+2*zy);
        z(2+(n-1)*4+(m-1)*ny*4)=z1;
        x(3+(n-1)*4+(m-1)*ny*4)=(m-1)*2*x1;
        y(3+(n-1)*4+(m-1)*ny*4)=-zy-y1-(n-1)*(2*y1+2*zy);
        z(3+(n-1)*4+(m-1)*ny*4)=0;
        x(4+(n-1)*4+(m-1)*ny*4)=x1+(m-1)*2*x1;
        y(4+(n-1)*4+(m-1)*ny*4)=-zy-2*y1-(n-1)*(2*y1+2*zy);
        z(4+(n-1)*4+(m-1)*ny*4)=0;
    end
end

生成动能哈密顿量程序

function H=Hamiltonian_NN_phos(x,y,z,ny)
t1=-1.22;t2=3.665;t3=-0.205;t4=-0.105;t5=-0.055;
d3=sqrt((x(5)-x(ny*4+2))^2+(y(5)-y(ny*4+2))^2+(z(5)-z(ny*4+2))^2);
d4=sqrt((x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2);
d5=sqrt((x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2);
N=length(x);
H=zeros(N,N);
for m=1:N
    for n=1:N
        d=sqrt((x(m)-x(n))^2+(y(m)-y(n))^2+(z(m)-z(n))^2);
        if abs(d-2.164)<0.001
            H(m,n)=t1;
        end
         if abs(d-2.207)<0.001
            H(m,n)=t2;
         end
         if abs(d-d3)<0.001
            H(m,n)=t3;
         end
         if abs(d-d4)<0.001
            H(m,n)=t4;
         end  
         if abs(d-d5)<0.001
            H(m,n)=t5;
         end  
    end
end

Hubbard计算程序

nx=6;
ny=8;
n=4*nx*ny;
Nx=nx;
Ny=4*ny;
[x,y,z]=zigzag_phos(nx,ny);
%plot3(x,y,z,'.','MarkerSize',20);
H=Hamiltonian_NN_phos(x,y,z,ny);
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(1,n/3+1:2/3*n);
y=y(1,n/3+1:2/3*n);
m=length(HD); % 体系总的原子个数
Nupavg=0.5*ones(m,1);  %将所有的<Nupi>取为1作初始值,自旋向上
Ndownavg=0.5*ones(m,1); %将所有的<Ndowni>取为1作初始值，自旋向下
u=0; %费米能级
Ncc=20; %迭代次数
dk=0.005*2*pi*3/Nx; %布里渊区步长
T=0.026; %温度
U=2; %Hubbard系数
k=-pi*3/Nx:dk:pi*3/Nx; %布里渊区路径
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)=0.8; %上边界上自旋电子密度
Ndownavg(1:Ny:m,1)=0.2; %上边界下自旋电子密度
Nupavg(Ny:Ny:m,1)=0.2; %下边界上自旋电子密度
Ndownavg(Ny:Ny:m,1)=0.8; %下边界下自旋电子密度
% Nupavg(1:2*Ny:m,1)=0.8; %上边界前半段上自旋电子密度
% Ndownavg(1:2*Ny:m,1)=0.2; %上边界前半段下自旋电子密度
% Nupavg(Ny+1:2*Ny:m,1)=0.2; %上边界后半段上自旋电子密度
% Ndownavg(Ny+1:2*Ny:m,1)=0.8; %上边界后半段下自旋电子密度
% Nupavg(Ny:2*Ny:m,1)=0.8; %下边界前半段上自旋电子密度
% Ndownavg(Ny:2*Ny:m,1)=0.2; %下边界前半段下自旋电子密度
% Nupavg(2*Ny:2*Ny:m,1)=0.2; %下边界后半段上自旋电子密度
% Ndownavg(2*Ny:2*Ny:m,1)=0.8; %下边界后半段下自旋电子密度
Nupavg=diag(Nupavg);
Ndownavg=diag(Ndownavg);
%%%%%%%%%%% 参数设置完成%%%%%%%%%%%
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)*Nx/3)+HDR*exp(1i*k(i)*Nx/3);
        [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*3/Nx); %更新Nupavg
    Ndownavg=Ndownavg0/(2*pi*3/Nx); %更新Ndownavg
    updensity(nc,:)=diag(Nupavg);
    downdensity(nc,:)=diag(Ndownavg);
    if mod(nc,5)==0
        fprintf('已完成第%d次迭代\n',nc)
    end
end
%%%%%%%%%%%%%%%%迭代完成，开始计算能带%%%%%%%%%%%%%%
Energy_up=0;
Energy_down=0;
for i=1:nk
    Hk=HD+HDL*exp(-1i*k(i)*Nx/3)+HDR*exp(1i*k(i)*Nx/3);
    [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);
    Energy_up=Energy_up+sum(diag(Eup).*f_up);
    Energy_down=Energy_down+sum(diag(Edown).*f_down);
    band_up(:,i) = diag(Eup);
    band_down(:,i) = diag(Edown);
end
Energy_up=Energy_up/nk/m;
Energy_down=Energy_down/nk/m;
plot(k,band_down,'b.')
hold on
plot(k,band_up,'r.')
set(gca,'YLim',[-2.7 2.7]);
set(gca,'XLim',[-pi*3/Nx pi*3/Nx]);

%%%%%%%%%%%%%下面画出自旋极化密度%%%%%%%%%%%%%%%%
figure(2)
local_density=diag(Nupavg)-diag(Ndownavg);
Zs=0;  % 总的局域电子密度
sigma=0.5;  % 高斯函数的标准差
du=0.01; % x方向采样间隔
dv=0.01; % y方向采样间隔
dx=2;  % x方向画图边界
dy=2;  % y方向画图边界
u=min(x)-dx:du:max(x)+dx;
v=min(y)-dy:dv:max(y)+dy;
[X,Y]=meshgrid(u,v);
 for i=1:m
     Z=real(local_density(i))*exp(-((X-x(i))/sigma).^2-((Y-y(i))/sigma).^2);
     Zs=Z+Zs;     
 end
mesh(X,Y,Zs)
set(gca,'YLim',[min(y)-dy max(y)+dy]);
set(gca,'XLim',[min(x)-dx max(x)+dx]);

猜你喜欢