SU(4)自洽计算实空间磁矩

这张写法得优势在于有坐标信息，因此可以实空间可视化。

Matlab 

Nx=6;  % 
Ny=8;  % 体系宽度（y方向的长度）
[x,y]=zigzag_graphene(Nx,Ny);
n = 4*Nx*Ny;
t1=-1;
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);           % 体系总的原子个数
% N1avg=0.25*eye(m,m);        
% N2avg=0.25*eye(m,m);        
% N3avg=0.25*eye(m,m);       
% N4avg=0.25*eye(m,m);
u= 0.056;                    % 费米能级
Ncc=6;                  % 迭代次数
dk=0.005;               % 布里渊区步长
T=0.001;                % 温度
U=0.5;                  % Hubbard系数
k=0:dk:2*pi;            % 布里渊区路径
nk=length(k);           % 布里渊区离散点个数
% 
% n1 =  unifrnd (0,1);
% n2 =  unifrnd (0,1);
% n3 =  unifrnd (0,1);
% n4 =  unifrnd (0,1);
% sumN = n1+n2+n3+n4;
% n1 = n1/sumN;
% n2 = n2/sumN;
% n3 = n3/sumN;
% n4 = n4/sumN;
% N1avg(1:Ny:m,1:Ny:m) = n1;  %上边界自旋1电子密度
% N2avg(1:Ny:m,1:Ny:m) = n2;  %上边界自旋2电子密度
% N3avg(1:Ny:m,1:Ny:m) = n3;  %上边界自旋3电子密度
% N4avg(1:Ny:m,1:Ny:m) = n4;  %上边界自旋4电子密度

% N1avg(Ny:Ny:m, Ny:Ny:m)= 0.1;   %下边界自旋1电子密度
% N2avg(Ny:Ny:m,Ny:Ny:m) = 0.4;   %下边界自旋2电子密度
% N3avg(Ny:Ny:m,Ny:Ny:m) = 0.4;   %下边界自旋3电子密度
% N4avg(Ny:Ny:m,Ny:Ny:m) = 0.1;   %下边界自旋3电子密度

load('N4avg1-2.mat')
load('N3avg1-2.mat')
load('N2avg1-2.mat')
load('N1avg1-2.mat')

% N1avg(1,1)=0;
% N2avg(1,1)=0;
% N3avg(1,1)=0;
% N4avg(1,1)=0;

band_1=zeros(m,nk); % 记录最后自旋向上的能带
band_2=zeros(m,nk); % 记录最后自旋向下的能带
band_3=zeros(m,nk); 
band_4=zeros(m,nk);


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


%%%%%%%%%%% 参数设置完成%%%%%%%%%%%
for nc=1:Ncc
    N1avg0=zeros(m,m); %初始化N1avg的临时变量用于k空间的积分
    N2avg0=zeros(m,m); %初始化N2avg的临时变量用于k空间的积分
    N3avg0=zeros(m,m); 
    N4avg0=zeros(m,m);
    m1avg = (N1avg-N2avg)/2;
    m2avg = (N1avg+N2avg-2*N3avg)/(2*sqrt(3));
    m3avg = (N1avg+N2avg+N3avg-3*N4avg)/(2*sqrt(6));
    C = 8/5*U*(m1avg.^2+m2avg.^2+m3avg.^2);
    U1 = m1avg/2+m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
    U2 = -m1avg/2+m2avg/(2*sqrt(3))+ m3avg/(2*sqrt(6));
    U3 =  -2*m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
    U4 =  -3*m3avg/(2*sqrt(6));
    for i=1:nk
        Hk=HD+HDL*exp(-1i*k(i))+HDR*exp(1i*k(i));
        [A1,E1] = eig(Hk-16/5*U*U1+C);
        [A2,E2] = eig(Hk-16/5*U*U2+C);
        [A3,E3] = eig(Hk-16/5*U*U3+C);
        [A4,E4] = eig(Hk-16/5*U*U4+C);
        f1 = Fermi_funtion(diag(E1),u,T)';
        f2 = Fermi_funtion(diag(E2),u,T)';
        f3 = Fermi_funtion(diag(E3),u,T)';
        f4 = Fermi_funtion(diag(E4),u,T)';
        fermi_1 = f1(ones(m,1),:);
        fermi_2 = f2(ones(m,1),:);
        fermi_3 = f3(ones(m,1),:);
        fermi_4 = f4(ones(m,1),:);
        N1 = sum(A1.*conj(A1).*fermi_1,2);
        N2 = sum(A2.*conj(A2).*fermi_2,2);
        N3 = sum(A3.*conj(A3).*fermi_3,2);
        N4 = sum(A4.*conj(A4).*fermi_4,2);
        N1avg0=N1avg0+diag(N1)*dk;
        N2avg0=N2avg0+diag(N2)*dk;
        N3avg0=N3avg0+diag(N3)*dk;
        N4avg0=N4avg0+diag(N4)*dk;
    end
    N1avg=N1avg0/(2*pi); %更新Nupavg
    N2avg=N2avg0/(2*pi); %更新Ndownavg
    N3avg=N3avg0/(2*pi); %更新Ndownavg
    N4avg=N4avg0/(2*pi); 
    density1(nc,:)=diag(N1avg);
    density2(nc,:)=diag(N2avg);
    density3(nc,:)=diag(N3avg);
    density4(nc,:)=diag(N4avg);
    if mod(nc,1)==0
        fprintf('已完成第%d次迭代\n',nc)
    end
end
%%%%%%%%%%%%%%%%迭代完成，开始计算能带%%%%%%%%%%%%%%
m1avg = (N1avg-N2avg)/2;
m2avg = (N1avg+N2avg-2*N3avg)/(2*sqrt(3));
m3avg = (N1avg+N2avg+N3avg-3*N4avg)/(2*sqrt(6));
C = 8/5*U*(m1avg.^2+m2avg.^2+m3avg.^2);
U1 = m1avg/2+m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
U2 = -m1avg/2+m2avg/(2*sqrt(3))+ m3avg/(2*sqrt(6));
U3 =  -2*m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
U4 =  -3*m3avg/(2*sqrt(6));

E1_total = 0;
E2_total = 0;
E3_total = 0;
E4_total = 0;

for i=1:nk
    Hk=HD+HDL*exp(-1i*k(i))+HDR*exp(1i*k(i));
    [A1,E1] = eig(Hk-16/5*U1+C);
    [A2,E2] = eig(Hk-16/5*U2+C);
    [A3,E3] = eig(Hk-16/5*U3+C);
    [A4,E4] = eig(Hk-16/5*U4+C);
    band_1(:,i) = diag(E1);
    band_2(:,i) = diag(E2);
    band_3(:,i) = diag(E3);
    band_4(:,i) = diag(E4);
    f1 = Fermi_funtion(diag(E1),u,T);
    f2 = Fermi_funtion(diag(E2),u,T);
    f3 = Fermi_funtion(diag(E3),u,T);
    f4 = Fermi_funtion(diag(E4),u,T);
    E1_total = E1_total + sum(f1.*diag(E1));
    E2_total = E2_total + sum(f2.*diag(E2));
    E3_total = E3_total + sum(f3.*diag(E3));
    E4_total = E4_total + sum(f4.*diag(E4));
end
E1_total = E1_total/nk/m;
E2_total = E2_total/nk/m;
E3_total = E3_total/nk/m;
E4_total = E4_total/nk/m;
E_total = E1_total + E2_total + E3_total + E4_total;
density = diag(N1avg)+diag(N2avg)+diag(N3avg)+diag(N4avg);
density_total = sum(density);
r = 1-density_total/(m-1);
fprintf('掺杂浓度%.4f\n',r)
fprintf('总能量%.4f\n',E_total)

% band_1 = band_1/abs(t1);
% band_2 = band_2/abs(t1);
% band_3 = band_3/abs(t1);
% band_4 = band_4/abs(t1);
plot(k,band_1,'r')
hold on
plot(k,band_2,'b')
plot(k,u*k./k);

kx = pi;
Hk=HD+HDL*exp(-1i*kx)+HDR*exp(1i*kx);
[~,Eh] = eig(Hk-16/5*U1+C);
[~,Ed] = eig(Hk-16/5*U2+C);
Eh = diag(Eh);
Ed = diag(Ed);
gap = abs(Eh(m/4)-Ed(m/4));
fprintf('gap的值是%d\n',gap);



% hold ondensity
% plot(k,band_3,'y')
% hold on
% plot(k,band_4,'g')
% plot(k,band_1,'color',[0.5 0.5 0.5])
% hold on 
% plot(k,band_2,'color',[0.5 0.5 0.5])
% plot(k,band_3,'color',[0.5 0.5 0.5])
% plot(k,band_4,'color',[0.5 0.5 0.5])
set(gca,'YLim',[-0.2 0.2]);
set(gca,'XLim',[0 2*pi]);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(2)
% local_density=diag(N1avg)-diag(N2avg);
%local_density = diag(N1avg)+diag(N2avg)-2*diag(N3avg);
% %local_density = diag(N1avg)+diag(N2avg)+diag(N3avg)-3*diag(N4avg);
% local_density = density;
% 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 H=Hamiltonian_NN_graphene(x,y,t)
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

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 = 100000;
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 [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 f=Fermi_funtion(E,u,T)
f=(exp((E-u)/T)+ones(length(E),1)).^-1;