SU(N)平均场Hubbard模型求解代码

nk = 512;                        
k = linspace(0,2*pi,nk);
t = 1;          % 最近邻hooping
U = 1;          % Hubbard系数
mu = 0;         % 化学势1/4填充取值
N = 40;         % N/2是对应的y方向原子个数
m = 4*N;        % 总的原子的个数
T = 0.001;      % 温度
N1avg = 0.25*eye(m);
N2avg = 0.25*eye(m);
N3avg = 0.25*eye(m);
N4avg = 0.25*eye(m);

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

N1avg(2,2) = 0.1;
N2avg(2,2) = 0.3;
N3avg(2,2) = 0.3;
N4avg(2,2) = 0.3;

Ncc = 5;        % 自洽计算迭代次数
band1=zeros(m,nk);
band2=zeros(m,nk);
band3=zeros(m,nk);
band4=zeros(m,nk);

% delta = 0.5;
% Hp = -delta + 2*delta*rand(m,1);
% Hp = diag(Hp);

for ncc = 1:Ncc
    N1avg0=zeros(m,m);                   
    N2avg0=zeros(m,m);                  
    N3avg0=zeros(m,m); 
    N4avg0=zeros(m,m);
    E1avg = 0;
    E2avg = 0;
    E3avg = 0;
    E4avg = 0;
    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 =  Hamiltonian_H0(k(i),t,N,mu);
        [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),0,T)';
        f2 = Fermi_funtion(diag(E2),0,T)';
        f3 = Fermi_funtion(diag(E3),0,T)';
        f4 = Fermi_funtion(diag(E4),0,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);
        N2avg0 = N2avg0+diag(N2);
        N3avg0 = N3avg0+diag(N3);
        N4avg0 = N4avg0+diag(N4);
        E1avg = E1avg + sum(f1'.*diag(E1));
        E2avg = E2avg + sum(f1'.*diag(E2));
        E3avg = E3avg + sum(f1'.*diag(E3));
        E4avg = E4avg + sum(f1'.*diag(E4));
    end
    N1avg = N1avg0/nk;
    N2avg = N2avg0/nk;
    N3avg = N3avg0/nk;
    N4avg = N4avg0/nk;
    fprintf('已完成第%d次迭代\n',ncc)
end
fprintf("自洽计算完成\n")
density = sum(diag(N1avg)+diag(N2avg)+diag(N3avg)+diag(N4avg));
Eavg = (E1avg + E2avg + E3avg + E4avg)/nk/m;
doping = 1-density/(m-1);
fprintf('总的电子密度%.4f/总的格点数%d\n',density,m)
fprintf('平均电子能量%.4f\n',Eavg);
fprintf('掺杂率%.4f\n',doping);


mu_up = zeros(m,nk);
mu_down = zeros(m,nk);

for i = 1:nk
    Hk =  Hamiltonian_H0(k(i),t,N,mu);
    [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);
    mu_up(:,i) = A1(:,N);
    mu_down(:,i) = A2(:,N);
    band1(:,i) = diag(E1);
    band2(:,i) = diag(E2);
    band3(:,i) = diag(E3);
    band4(:,i) = diag(E4);
end
plot(k,band1,'r');
hold on
plot(k,band2,'b');
plot(k,band3,'g');
plot(k,band4,'y');
set(gca,'YLim',[-0.5 0.5]);
set(gca,'XLim',[0 2*pi]);

S1 = (N1avg - N2avg) / 2;
S2 = (N1avg + N2avg - 2 * N3avg ) / 2 / sqrt(3);
S3 = (N1avg + N2avg + N3avg - 3 * N4avg) / 2 / sqrt(6); 
S1 = diag(S1);
S2 = diag(S2);
S3 = diag(S3);

plot(S1(2:32),'-r.')
hold on
plot(S2(2:32),'-g.')
plot(S3(2:32),'-b.')