用“杂质”产生1D的Anderson局域化
之前某视频的代码。（因为偷懒没有怎么加注释）
% 设置格点化的空间坐标
L = 400; % Interval Length
N = 2000; % No of points
x = linspace(-L/2,L/2-L/(N),N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
%***************************************************************
%设置 Hamiltonian
V=30; % 设置无磁性杂质对势能的扰动幅度
W = 30; % 设置磁性杂质对1/2自旋的作用强度
G = 1;%控制边界上的衰减
% 用U1把碰边的波尽可能消去
U1 = -heaviside(x-195)*(1i*G).*(x-195).*(x-195)-heaviside(-x-195)*(1i*G).*(x+195).*(x+195); %边界附近的虚势能，把跑到边界的粒子移出系统
%生成无序的“磁性杂质”的空间坐标和正负号
PM1 = heaviside(195-x).*heaviside(x+195).*(random('Discrete Uniform',2,[N,1])-1.5)*2;
Imp_m = W*floor(random('Discrete Uniform',200,[N,1])/196).*PM1;
%生成无序的“无磁性杂质”的空间坐标和正负号
PM1 = heaviside(195-x).*heaviside(x+195).*(random('Discrete Uniform',2,[N,1])-1.5)*2;
Imp_e = V*floor(random('Discrete Uniform',200,[N,1])/196).*PM1;
% 生成周期排布杂质的空间坐标和正负号
Imp_p = zeros(N,1);
pm = 1;
for j=1:N
if mod(j-1025,50)==0
Imp_p(j) = pm;
pm = -pm;
end
end
% Finite-difference representation of Laplacian and Hamiltonian
e = ones(N,1);
Lap = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2; % 周期性边界条件下的离散Laplacian
e2 = ones(2,1);
e3 = [1;1];
H1 = -(1/2)*kron(Lap,speye(2)) +...
kron(spdiags(Imp_m,0,N,N),spdiags([e3 e3],[-1 1],2,2))+...
kron(spdiags(U1,0,N,N),speye(2)); % Hamiltonian,磁性无序杂质
H0 = -(1/2)*kron(Lap,speye(2)) +...
kron(spdiags(Imp_e,0,N,N),spdiags(e3,0,2,2))+...
kron(spdiags(U1,0,N,N),speye(2)); % Hamiltonian，无磁性无序杂质
H00 = -(1/2)*kron(Lap,speye(2))+...
kron(spdiags(V*Imp_p,0,N,N),spdiags(e3,0,2,2))+...
kron(spdiags(U1,0,N,N),speye(2)) ; % Hamiltonian，无磁性周期性杂质
H000 = -(1/2)*kron(Lap,speye(2))+...
kron(spdiags(W*Imp_p,0,N,N),spdiags([e3 e3],[-1 1],2,2))+...
kron(spdiags(U1,0,N,N),speye(2)) ; % Hamiltonian，磁性周期性杂质
% Parameters for computing psi(x,T) at different times 0 < T < TF
NT = 1000; % No. of time steps
TF = 25; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
hbar = 1;
% Time displacement operator E=exp(-iHdT/hbar)
E = expm(-1i*full(H1)*dT/hbar); % time desplacement operator
E0 = expm(-1i*full(H0)*dT/hbar);
E00 = expm(-1i*full(H00)*dT/hbar);
E000 = expm(-1i*full(H000)*dT/hbar);
%***************************************************************
% 下面部分可以和上面部分分开运行。
% 仅修改波函数以及画图的相关参数时，不必重复运行上面的部分。
% Parameters for making intial wavefunction psi
ko1 = 5;
a = 2; % package width
xo=0;% peak center
% 离散化的空间坐标对最大动量有限制，不然会出现Bloch振荡这类现象
if ko1>(pi*N/L)
error('too large ko1');
end
S1 = [1;0];
S1 = S1/sqrt(S1'*S1);
% 设置波函数为 平面波*Gauss轮廓，及其归一化
psi1 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko1);
psi1 = psi1/sqrt(psi1'*psi1*dx); % normalize the psi(x,0)
psi1 = kron(psi1,S1);
psi2 = psi1;
psi3 = psi1;
psi4 = psi1;
%******************************************
% 画图的乱七八糟的参数
% y用于设置fft图的横坐标
D=floor(20*(L/2-2)/(2*pi));
y = 2*pi*(0:length(floor(N/2+N*(2/L)):N)-1)/(L/2-2);
N1=length(y);
H=0.3;
fig=figure;
fig.Position = [0 0 2000 1000];
dy = y(2)-y(1);
for t = 1:6*NT % time index for loop
% calculate probability density rho(x,T)
psi1 = E*psi1; % calculate new psi from old psi
rho1 = conj(psi1).*psi1; % rho(x,T)
psi2 = E0*psi2;
rho2 = conj(psi2).*psi2;
psi3 = E00*psi3;
rho3 = conj(psi3).*psi3;
psi4 = E000*psi4;
rho4 = conj(psi4).*psi4;
% 画图
subplot(414);
plot(x,rho1(1:2:end),'b',x,real(psi1(1:2:end)).^2,'g--',...
x,rho1(2:2:end),'r',x,real(psi1(2:2:end)).^2,'m--',...
x,rho1(1:2:end) + rho1(2:2:end),'k');
hold on
plot(x,H*(Imp_m/W*0.2),'Color',[0.5,0.7,0]);
plot(x,-H*(Imp_m/W*0.2),'Color',[0.5,0.3,1]);
hold off
text(-L*(0.4),0.4*H,'magnetic impurities, disorder');
axis([-L/2 L/2 0 H]); % set x,y axis parameters for plotting
xlabel('x', 'FontSize', 24);
subplot(413);
plot(x,rho2(1:2:end),'b',x,real(psi2(1:2:end)).^2,'g--',...
x,rho2(2:2:end),'r',x,real(psi2(2:2:end)).^2,'m--');
hold on
plot(x,H*(Imp_e/V*0.3),'c--');% 画势能
plot(x,-H*(Imp_e/V*0.3),'y--');% 画势能
hold off
text(-L*(0.4),0.4*H,'nonmagnetic impurities, disorder');
axis([-L/2 L/2 0 H]);
xlabel('x', 'FontSize', 24);
subplot(412);
plot(x,rho4(1:2:end),'b',x,real(psi4(1:2:end)).^2,'g--',...
x,rho4(2:2:end),'r',x,real(psi4(2:2:end)).^2,'m--',...
x,rho4(1:2:end) + rho4(2:2:end),'k');
text(-L*(0.4),0.4*H,'magnetic impurities, periodic');
hold on
plot(x,H*(Imp_p*0.2),'Color',[0.5,0.7,0]);
plot(x,-H*(Imp_p*0.2),'Color',[0.5,0.3,1]);
hold off
axis([-L/2 L/2 0 H]);
xlabel('x', 'FontSize', 24);
subplot(411);
plot(x,rho3(1:2:end),'b',x,real(psi3(1:2:end)).^2,'g--',...
x,rho3(2:2:end),'r',x,real(psi3(2:2:end)).^2,'m--');
hold on
plot(x,H*(Imp_p*0.3),'c--');% 画势能
plot(x,-H*(Imp_p*0.3),'y--');% 画势能
hold off
text(-L*(0.4),0.4*H,'nonmagnetic impurities, periodic');
axis([-L/2 L/2 0 H]);
xlabel('x', 'FontSize', 24);
drawnow
end
