pi flux计算程序

主计算程序 

Nx=18;  % 
Ny=4; % 体系宽度（y方向的长度）
[x,y]=zigzag_graphene(Nx,Ny);
n = 4*Nx*Ny;
t1=-2.7;
H=Hamiltonian_NN_graphene(x,y,t1);
H=Hamiltonian_NN_phase(Nx,4*Ny,H,t1);
H=Hamiltonian_defect(x,y,4*Ny,H,t1);
miu = sqrt(3)*t1;
H = H- miu*eye(n,n);
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.25*eye(m,m);  %将所有的<Nupi>取为1作初始值,自旋向上
Ndownavg=0.25*eye(m,m); %将所有的<Ndowni>取为1作初始值，自旋向下
u=0; %费米能级
Ncc=10; %迭代次数
dk=0.005; %布里渊区步长
T=0.038; %温度
U=1.5; %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.4; %上边界上自旋电子密度
Ndownavg(1:Ny:m,1:Ny:m)=0.1; %上边界下自旋电子密度




%Nupavg(Ny:Ny:m,Ny:Ny:m)=0.25; %下边界上自旋电子密度
%Ndownavg(Ny:Ny:m,Ny:Ny:m)=0.25; %下边界下自旋电子密度
%%%%%%%%%%% 参数设置完成%%%%%%%%%%%
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-1/4*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量
        [Aup,Eup]=eig(Hk+U*Ndownavg-1/4*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-1/4*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量
    [Aup,Eup]=eig(Hk+U*Ndownavg-1/4*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 4]);
set(gca,'XLim',[0 2*pi]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
local_density=abs(diag(Nupavg)-diag(Ndownavg));
Zs=0;  % 总的局域电子密度
sigma=0.2;  % 高斯函数的标准差
du=0.01; % x方向采样间隔
dv=0.01; % y方向采样间隔
u=min(x)-0.5:du:max(x)+0.5;
v=min(y)-0.5:dv:max(y)+0.5;
[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)

生成坐标构型的函数

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=[x2;x1];
    y2=[y2;y1+i*ones(4,1)*3];
end
x=x2;
y=y2;
n=length(x2);
for i=1:nx-1
    x=[x;x2+i*ones(n,1)*sqrt(3)];
    y=[y;y2];
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)<eps
            H(i,j)=t;
        end    
    end
end

生成pi flux的相位函数

function H=Hamiltonian_NN_phase(nx,ny,H,t)
for i=1:nx
    flag = mod(i,2);
    t1 = (-1)^flag*t;
    for j=1:ny
        k = (i-1)*ny+j;
        if mod(j,8)==0 && j~=ny
            H(k,k+1) = -t1;
        end
        if mod(j,8)==1 && j>1
            H(k,k-1) = -t1;
        end
        if mod(j,8)==2
            H(k,k+1)= t1;
        end
        if mod(j,8)==3
            H(k,k-1)=t1;
        end
        if mod(j,8)==4
            H(k,k+1)=t1;
        end
        if mod(j,8)==5
            H(k,k-1)=t1;
        end
        if mod(j,8)==6
            H(k,k+1)=-t1;
        end
        if mod(j,8)== 7
            H(k,k-1)=-t1;
        end
    end
end

生成缺陷的函数

function H=Hamiltonian_defect(x,y,ny,H,t)
N=length(x);
eps=0.01;
E = 1000000;
for i=1:ny:N
    if mod(i,2*ny)==1
        H(i,i)= E;
        for j=1:N
            if abs(sqrt((x(i)-x(j))^2+(y(i)-y(j))^2)-1)<eps
                H(i,j)=t;
            end
        end    
    end
end

费米函数

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