沿着高对称线计算磁化率

首先计算出布里渊区规则区域的磁化率

kx=load('KX.grd');      
ky=load('KY.grd'); 
Ed =load('BAND.B9.grd');
Ep = load('BAND.B8.grd');
mu = 0;
T = 0.01;             
delta = 0.001;             
nx = length(kx);
ny = length(ky);
chi = zeros(nx,ny);
for i = 1:nx
    for j = 1:ny
        for k=1:nx
            for l=1:ny
                index_kq_x = mod(i+k,nx)+1;
                index_kq_y = mod(j+l,ny)+1;
                Ek = Ed(k,l);
                Ekq = Ed(index_kq_x,index_kq_y);
                %chi(i,j) = chi(i,j) + (Fermi_funtion(Ek,mu,T)-Fermi_funtion(Ekq,mu,T))/(Ek-Ekq+1i*delta);
                chi(i,j) = chi(i,j) + delta_function(Ek)*delta_function(Ekq);
            end
        end
    end
end

然后截取我们需要的磁化率，这里有如下技巧可以让我们顺利拿到所需数据

 

 

chi = real(chi)/nx/ny;
chi_line_Gammar_M = chi(1:nx/2,1);
chi_line = diag(chi);
chi_line_M_Gammar = flip(chi_line(1:nx/2));
chi_line_plot = [chi_line_Gammar_M;chi_line_M_Gammar];

plot([25 25], [0 25],'-b');
hold on
plot([50 50], [0 25],'-k');
plot([60 60], [0 25],'-b');
plot([66 66], [0 25],'-k');
plot([79 79], [0 25],'-b');
plot(chi_line_plot);

其中delta函数用洛伦兹分布近似

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

function delta =  delta_function(x)
a = 0.01;
delta = 1/pi*a/(a^2+x^2);
end