关于石墨烯载流子密度图matlab仿真问题
请教一下各位专家，关于石墨烯载流子密度图matlab仿真问题。本人想复现文献中的图，但总是不对。
文献：
仿真的目标是文献support imformation中的公式S7，
仿真的目标图
根据文献所提供的引用中的公式推导：
本人推导，不知道对不对，请各位指教：
这是本人写的代码，对与τ的理解不够，所以随便设了个，所得到的函数图是一条直线与文献不符，我应该怎样修改代码
有没有相关的代码，提供一下，模拟试试。
有代码的话，模拟和实际相差很大吗。
你t_values的值首尾取的太小了，应该是不用乘以10的-15次方的，中间的步长可以乘以10的-15次方，所以画出来的图只是在10^-15这个数量级的一个范围内的一个图
你把代码发出来，不要截图，调试起来更方便
N_values改成for循环的形式计算，方便加断点调试查看参数值是否传递更新正确，
clc, clear, close all;
t_values = -350*10^-15:0.2*10^-15:700*10^-15; % 指定 t 的取值范围
N_values = zeros(size(t_values)); % 初始化 N_values 数组
hbar = 1.05457e-34;
V_F = 10e6;
hbar_w = 0.6;
D = hbar_w / (pi * hbar^2 * V_F^2);
pi_alphaf = 0.0023;
tau = 8.2e-17;
tau1 = 7e-15;
I = 1;
for i = 1:length(t_values)
t_temp = t_values(i); % 临时保存 t 的值
% 定义函数
f = @(x) exp(-x.^2/tau.^2 + (x-t_temp)/tau1 + ...
(pi_alphaf*I)/(D*hbar_w)*sqrt(pi)*tau*(erf((x-t_temp)/tau)-erf(0)));
% 计算函数值
N_values(i) = D - 2*pi_alphaf*I/(pi*hbar^2*V_F^2)*integral(f, -inf, t_temp);
end
plot(t_values, N_values)
xlabel('t')
ylabel('N')
title('关于 t 的函数曲线')
grid on
这样绘制出来的N_values还是定值，
随便篡改一下各参数的量纲(e-34这种数值太小了，对于integral函数来说计算精度可能精度没有那么高), N_values开始有变化：
说明程序语法和逻辑应该没错，需要对照paper检查一下参数的量纲是否有误。
这个文献有地址吗？
文献网盘
链接：https://pan.baidu.com/s/1E1K2BgHtyCxRF_QJFoxbJQ?pwd=edle
提取码：edle
文献地址发一下
matlab化石墨烯,锯齿石墨烯磁性hubbard matlab计算程序
主程序
Nx=3;
Ny=50;
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*eye(m,m); %将所有的取为1作初始值,自旋向上
Ndownavg=0.5*eye(m,m); %将所有的取为1作初始值，自旋向下
u=0; %费米能级
Ncc=20; %迭代次数
dk=0.01; %布里渊区步长
T=0.038; %温度
U=2; %Hubbard系数
k=0:dk:2*pi; %布里渊区路径
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:Ny:m)=0.65; %上边界上自旋电子密度
Ndownavg(1:Ny:m,1:Ny:m)=0.35; %上边界下自旋电子密度
Nupavg(Ny:Ny:m,Ny:Ny:m)=0.35; %下边界上自旋电子密度
Ndownavg(Ny:Ny:m,Ny:Ny:m)=0.65; %下边界下自旋电子密度
%%%%%%%%%%% 参数设置完成%%%%%%%%%%%
for nc=1:Ncc
Nupavg0=zeros(m,m); %初始化Nupavg的临时变量用于k空间的积分
Ndownavg0=zeros(m,m); %初始化Ndownavg的临时变量用于k空间的积分
for i=1:nk
Hk=HD+HDL*exp(-1i*k(i))+HDR*exp(1i*k(i));
[Adown,Edown]=eig(Hk+U*Nupavg-0.5*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量
[Aup,Eup]=eig(Hk+U*Ndownavg-0.5*U*eye(m,m)); %通过Ndownavg求解出Aup,Aup为up波函数系数向量
% 求电子密度的向量化写法
f_up=Fermi_funtion(diag(Eup),u,T)';
f_down=Fermi_funtion(diag(Edown),u,T)';
fermi_up = f_up(ones(m,1),:);
fermi_down = f_down(ones(m,1),:);
Nup = sum(Aup.*conj(Aup).*fermi_up,2);
Ndown = sum(Adown.*conj(Adown).*fermi_down,2);
Nupavg0=Nupavg0+diag(Nup)*dk;
Ndownavg0=Ndownavg0+diag(Ndown)*dk;
end
Nupavg=Nupavg0/(2*pi); %更新Nupavg
Ndownavg=Ndownavg0/(2*pi); %更新Ndownavg
updensity(nc,:)=diag(Nupavg);
downdensity(nc,:)=diag(Ndownavg);
if mod(nc,5)==0
fprintf('已完成第%d次迭代\n',nc)
end
end
%%%%%%%%%%%%%%%%迭代完成，开始计算能带%%%%%%%%%%%%%%
for i=1:nk
Hk=HD+HDL*exp(-1i*k(i))+HDR*exp(1i*k(i));
[Adown,Edown]=eig(Hk+U*Nupavg-0.5*U*eye(m,m)); %通过Nupavg求解出Adown,Adown为down波函数系数向量
[Aup,Eup]=eig(Hk+U*Ndownavg-0.5*U*eye(m,m)); %通过Ndownavg求解出Aup,Aup为up波函数系数向量
band_up(:,i) = diag(Eup);
band_down(:,i) = diag(Edown);
end
plot(k,band_down,'b.')
hold on
plot(k,band_up,'r.')
set(gca,'YLim',[-2.7 2.7]);
set(gca,'XLim',[0 2*pi]);
生成体系座标程序
function [x,y]=zigzagP(N,number);
x=zeros(N*number,1);
y=zeros(N*number,1);
x(1)=sqrt(3)/2;
y(1)=0;
x(2)=0;
y(2)=-0.5;
for j=1:N-2
y(j+2)=y(j)-1.5;
if mod(j+2,4)==1||mod(j+2,4)==0
x(j+2)=sqrt(3)/2;
else
x(j+2)=0;
end
end
for j=1:number-1
for l=1:N
x(j*N+l)=x(l)+j*sqrt(3);
y(j*N+l)=y(l);
end
end
最近邻哈密顿量与前面程序相同
费米函数
function f=Fermi_funtion(E,u,T)
f=(exp((E-u)/T)+ones(length(E),1)).^-1;
具体是哪篇文章啊，要从头看文章理解的基础上才好做。如果你还没解决，可以尝试联系论文作者。