Calculating integration using for loops

7 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2021-5-21
编辑： Luqman Saleem 2021-5-21
My question is about 2D square Brillouin zone integrations.
Problem Summary:
I am rederiving one of the results of this publication. The mathematical expressions for and are given in two different forms, Eq 10 and Eq 12. Expected results are . However, when I calculate using Eq.10 and Eq.12, I get confilicting results. I have attached my codes both for Eq10.m and Eq12.m.
Equation 10:
The mathematical expression for are
--- Eq. 10a
--- Eq. 10b
here
  • with as 2-by-2 matrices (given in code)
  • is th eigenvalue of matrix H
  • is eigenvalues of H
  • are constants
  • I use approximation
Equation 12:
The mathematical expression for are
--- Eq. 12a
--- Eq. 12b
here
  • are constants
My Results:
So far, for Eq.10, I get (which are expected results). For Eq.12, I get opposite results . I don't know what exactly I am missing here. Please help!
My Codes:
Eq10.m
% ============= Eq10.m =============
clear; clc;
% ===== inputs ====
NBZ = 100; %number of BZ points
eta = 1e-1; %for delta function approximation
EF = 100; %E_F
% ===== constants =====
me = 9.1e-31/(1.6e-19*(1e10)^2);
aR = 2.95;
muBx = 0.1;
m = 0.32*me;
hbar = 6.5911e-16;
Gamma = 1.9268e+13;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
sz = [1 0;0 -1];
% ===== k-points =====
kx1 = -pi; kx2 = pi;
ky1 = -pi; ky2 = pi;
kxs = linspace(kx1,kx2,NBZ+1); kxs(end) = [];
kys = linspace(ky1,ky2,NBZ+1); kys(end) = [];
dkx = kxs(2)-kxs(1);
dky = kys(2)-kys(1);
% ===== for loops to calculate \alpha =====
alpha_x = 0;
alpha_y = 0;
for ikx = 1:NBZ %loop over kx
kx = kxs(ikx);
for iky = 1:NBZ %loop over ky
ky = kys(iky);
% ===== constructing H matrix ====
k = sqrt(kx^2+ky^2);
HR = hbar^2*k^2/(2*m)*eye(2) + aR*(kx*sy - ky*sx);
Hzee = muBx*sx;
H = HR + Hzee;
% ===== diagonalizing H =====
[evecs,evals] = eig(H);
% ===== vx, vy matrix ====
vx = 1/hbar*[ 0, -aR*1i
aR*1i, 0];
vy = 1/hbar*[ 0, -aR
-aR, 0];
% ===== Sx, Sy matrix =====
Sx = hbar * (evecs' * sx * evecs);
Sy = hbar * (evecs' * sy * evecs);
% ===== X, Y, and J matrix =====
Xx = 0.5*(Sx*vx + vx*Sx);
Xy = 0.5*(Sx*vy + vy*Sx);
Yx = 0.5*(Sy*vx + vx*Sy);
Yy = 0.5*(Sy*vy + vy*Sy);
Jx = -e*vx;
Jy = -e*vy;
% ===== loop over n =====
for n = 1:2
Un = evecs(:,n); %nth eigenvector
En = evals(n,n); %nth eigenvalue
delta_fun = 1/(eta*sqrt(2*pi))*exp(-0.5*((En-EF)^2/(eta^2)));
Xx_k = Un'*Xx*Un; %value of X_x at given k
Xy_k = Un'*Xy*Un;
Yx_k = Un'*Yx*Un;
Yy_k = Un'*Yy*Un;
Jx_k = Un'*Jx*Un;
Jy_k = Un'*Jy*Un;
crossPx = Xx_k*Jy_k - Xy_k*Jx_k; %cross product for \alpha_x
crossPy = Yx_k*Jy_k - Yy_k*Jx_k; %cross product for \alpha_y
alpha_x = alpha_x + 1/(2*Gamma)*crossPx*delta_fun*dkx*dky;
alpha_y = alpha_y + 1/(2*Gamma)*crossPy*delta_fun*dkx*dky;
end
end
end
alpha_x = alpha_x/(1.6e-19); % 1.6e-19 factor comes from units (you can ignore it)
alpha_y = alpha_y/(1.6e-19);
[alpha_x, alpha_y]
%it gives: [0.0000 + 0.0000i -5.0997 - 0.0000i]
Eq12.m
% ============= Eq12.m =============
clear; clc;
% ===== inputs ====
NBZ = 100; %number of BZ points
EF = 100; %E_F
% ===== constants =====
me = 9.1e-31/(1.6e-19*(1e10)^2);
aR = 2.95;
muBx = 0.1;
m = 0.32*me;
hbar = 6.5911e-16;
Gamma = 1.9268e+13;
e = 1.6e-19;
sx = [0 1;1 0];
sy = [0 -1i;1i 0];
sz = [1 0;0 -1];
% ===== k-points =====
kx1 = -pi; kx2 = pi;
ky1 = -pi; ky2 = pi;
kxs = linspace(kx1,kx2,NBZ+1); kxs(end) = [];
kys = linspace(ky1,ky2,NBZ+1); kys(end) = [];
dkx = kxs(2)-kxs(1);
dky = kys(2)-kys(1);
% ===== for loops to X_x, X_y, Y_x, Y_y, J_x, J_y at each k and n=====
for ikx = 1:NBZ %loop over kx
kx = kxs(ikx);
for iky = 1:NBZ %loop over ky
ky = kys(iky);
% ===== constructing H matrix ====
k = sqrt(kx^2+ky^2);
HR = hbar^2*k^2/(2*m)*eye(2) + aR*(kx*sy - ky*sx);
Hzee = muBx*sx;
H = HR + Hzee;
% ===== diagonalizing H =====
[evecs,evals] = eig(H);
% ===== vx, vy matrix ====
vx = 1/hbar*[ 0, -aR*1i
aR*1i, 0];
vy = 1/hbar*[ 0, -aR
-aR, 0];
% ===== Sx, Sy matrix =====
Sx = hbar * (evecs' * sx * evecs);
Sy = hbar * (evecs' * sy * evecs);
% ===== X, Y, and J matrix =====
Xx = 0.5*(Sx*vx + vx*Sx);
Xy = 0.5*(Sx*vy + vy*Sx);
Yx = 0.5*(Sy*vx + vx*Sy);
Yy = 0.5*(Sy*vy + vy*Sy);
Jx = -e*vx;
Jy = -e*vy;
% ===== loop over n =====
for n = 1:2
Un = evecs(:,n); %nth eigenvector
En = evals(n,n); %nth eigenvalue
Xx_k(ikx,iky,n) = Un'*Xx*Un; %value of X_x at given k
Xy_k(ikx,iky,n) = Un'*Xy*Un;
Yx_k(ikx,iky,n) = Un'*Yx*Un;
Yy_k(ikx,iky,n) = Un'*Yy*Un;
Jx_k(ikx,iky,n) = Un'*Jx*Un;
Jy_k(ikx,iky,n) = Un'*Jy*Un;
E_k(ikx,iky,n) = En;
end
end
end
[X,Y] = meshgrid(kxs,kys);
% === finding curl for each n
for n=1:2
[curlX(:,:,n),~]=curl(X,Y,Xx_k(:,:,n),Xy_k(:,:,n));
[curlY(:,:,n),~]=curl(X,Y,Yx_k(:,:,n),Yy_k(:,:,n));
end
%% now evaluating integration over omega_x and omega_y
omega_x = 0;
omega_y = 0;
for n = 1:2
SumBZomega = zeros(1,2);
for ikx = 1:NBZ %over kx
kx = kxs(ikx);
for iky = 1:NBZ %over kx
ky = kys(iky);
En = E_k(ikx,iky,n); %En at each kx,ky, and n
if En<=EF %applying En<=EF condition
SumBZomega(1) = SumBZomega(1) + curlX(ikx,iky,n)*dkx*dky;
SumBZomega(2) = SumBZomega(2) + curlY(ikx,iky,n)*dkx*dky;
end
end
end
omega_x = omega_x + SumBZomega(1);
omega_y = omega_y + SumBZomega(2);
end
% ==== finally alpha_x, alpha_y ====
alpha_x = e/(2*Gamma*hbar*4*pi^2)*omega_x/(1.6e-19); %1.6e-19 coming from units
alpha_y = e/(2*Gamma*hbar*4*pi^2)*omega_y/(1.6e-19);
[alpha_x, alpha_y]
%it gives [0.4840 - 0.0000i -0.0000 + 0.0000i]

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by