用“雜質(zhì)”產(chǎn)生1D的Anderson局域化
之前某視頻的代碼。（因為偷懶沒有怎么加注釋）
% 設(shè)置格點化的空間坐標(biāo)
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
%***************************************************************
%設(shè)置 Hamiltonian
V=30; % 設(shè)置無磁性雜質(zhì)對勢能的擾動幅度
W = 30; % 設(shè)置磁性雜質(zhì)對1/2自旋的作用強度
G = 1;%控制邊界上的衰減
% 用U1把碰邊的波盡可能消去
U1 = -heaviside(x-195)*(1i*G).*(x-195).*(x-195)-heaviside(-x-195)*(1i*G).*(x+195).*(x+195); %邊界附近的虛勢能，把跑到邊界的粒子移出系統(tǒng)
%生成無序的“磁性雜質(zhì)”的空間坐標(biāo)和正負(fù)號
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;
%生成無序的“無磁性雜質(zhì)”的空間坐標(biāo)和正負(fù)號
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;
% 生成周期排布雜質(zhì)的空間坐標(biāo)和正負(fù)號
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,磁性無序雜質(zhì)
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，無磁性無序雜質(zhì)
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，無磁性周期性雜質(zhì)
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，磁性周期性雜質(zhì)
% 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);?
%***************************************************************
% 下面部分可以和上面部分分開運行。
% 僅修改波函數(shù)以及畫圖的相關(guān)參數(shù)時，不必重復(fù)運行上面的部分。
% Parameters for making intial wavefunction psi
ko1 = 5;?
a = 2; % package width
xo=0;% peak center
% 離散化的空間坐標(biāo)對最大動量有限制，不然會出現(xiàn)Bloch振蕩這類現(xiàn)象
if ko1>(pi*N/L)
? ? error('too large ko1');
end
S1 = [1;0];
S1 = S1/sqrt(S1'*S1);
% 設(shè)置波函數(shù)為 平面波*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;
%******************************************
% 畫圖的亂七八糟的參數(shù)
% y用于設(shè)置fft圖的橫坐標(biāo)
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
