首先计算出布里渊区规则区域的磁化率
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