子函数与前面程序相同,使用了超胞计算的技术
Nx=30;
Ny=32;
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*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.038; %温度
U=2; %Hubbard系数
k=0:dk:2*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.8