Haldane Model笔记

有空好好整理下

六角蜂窝晶格

晶格规范与原子规范

图片

FIG.1: 色散关系，晶格规范(a)，原子规范(b)。Berry曲率，晶格规范(c)，原子规范(d)。[CODE]

边界态，Winding loops与二维系统Winding Number

图片

FIG.2: 理论模型，能带结构，winding loop以及winding number。分别对应Zigzag边界(a)(e)(i)(m)，Bearded边界(b)(f)(j)(n)，Armchair边界(c)(g)(k)(o)和twig边界(d)(h)(l)(p)。

Haldane Model

模型

图片

FIG.3: 色散关系及贝利曲率（对应陈数1，-1）。

图片

FIG.4: 能带结构。

相变与能带结构

图片
能带结构之前的图改下参数即可，画了没保存，相变交接不仅陈数会变为0，贝利曲率的正负也会出现突变。

手性边界态

图片

Modified Haldane Model

图片

FIG.5: 色散关系及能带结构。

核心代码

FIG.1

色散关系（晶格规范）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
clc
clear all

%% 参数
e_0=1;
a=sqrt(3)*e_0;
r1=a*[1,0]';
r2=a*[1/2,-sqrt(3)/2]';
t=1;

%% 循环次数
count_x=500;
count_y=500;

%% k_x,k_y 取值
k_xx=4*linspace(-1*pi,1*pi,count_x)/a;
k_yy=4*linspace(-1*pi,1*pi,count_y)/a;

%% 颜色
colorrd=[207,108,82]/255;
coloryl=[241,141,0]/255;
colorbl=[109,152,165]/255;

%% 计算
bar = waitbar(0,'色散关系计算中...'); %显示进度条
for x = 1:count_x
    per=x/count_x;
    str=['色散关系计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1:count_y
        k_x=k_xx(x);
        k_y=k_yy(y);
        k=[k_x k_y];
        % Hamiltonian
        H=zeros(2);
        H(1,2)=t+t*exp(-1i*k*r1)+t*exp(-1i*k*r2);
        H(2,1)=conj(H(1,2));
        E(x,y,:)=eig(H);
    end
end

%% 绘图
% 设置图窗
f1=figure(1);
position0=get(f1,'position');
set(f1,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

%结果
nexttile
pcolor(k_xx/(pi/a),k_yy/(pi/a),E(:,:,1)')
hold on
pcolor(k_xx/(pi/a),k_yy/(pi/a),E(:,:,2)')
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
clim([0,max(max(max(E)))])
load('colorData.mat');
axis([min(k_xx/(pi/a)),max(k_xx/(pi/a)),min(k_yy/(pi/a)),max(k_yy)/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[min(k_xx/(pi/a)):1:max(k_xx/(pi/a))], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[min(k_yy/(pi/a)):1:max(k_yy/(pi/a))], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
colormap(Purples6);

%% 第一布里渊区

% 计算倒格子的基矢量
b1 = 2 * pi * [r2(2), -r2(1)] / (r1(1)*r2(2) - r1(2)*r2(1));
b2 = 2 * pi * [-r1(2), r1(1)] / (r1(1)*r2(2) - r1(2)*r2(1));

% 生成倒格子的点
point_N = 5; % 在每个方向生成的点的数量
points = [];
for i = -point_N:point_N
    for j = -point_N:point_N
        point = i * b1 + j * b2;
        points = [points; point];
    end
end

% 计算到原点的距离
distances = sqrt(sum((points - [0,0]).^2, 2));

% 找到距离最小的6个点
[~, indices] = mink(distances, 7);
nearest_points = points(indices, :);

% 计算第一布里渊区的顶点
K = convhull(nearest_points);

% 绘制第一布里渊区
plot(nearest_points(K, 1)/(pi/a), nearest_points(K, 2)/(pi/a),':','Color','k');
K1=[K(1),K(3),K(5),K(1)];
K2=[K(2),K(4),K(6),K(2)];
plot(nearest_points(K1, 1)/(pi/a), nearest_points(K1, 2)/(pi/a),':','Color','k');
plot(nearest_points(K2, 1)/(pi/a), nearest_points(K2, 2)/(pi/a),':','Color','k');

% 绘制倒格子点
plot(points(:, 1)/(pi/a), points(:, 2)/(pi/a), '.','Color','black');

% 正六边形
n = 6;
radius = 4*sqrt(3)/3/a;
theta = linspace(0, 2*pi, n+1);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x, y, '--','LineWidth',1.5,'Color',colorrd);

close(bar)%关闭进度条

色散关系（原子规范）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
clc
clear all

%% 参数
e_0=1;
a=sqrt(3)*e_0;
e1=e_0*[0,1]';
e2=e_0*[sqrt(3)/2,-0.5]';
e3=e_0*[-sqrt(3)/2,-0.5]';
a1=a*[1,0]';
a2=a*[1/2,-sqrt(3)/2]';
t=1;
%% 循环次数
count_x=500;
count_y=500;

%% k_x,k_y 取值
k_xx=4*linspace(-1*pi,1*pi,count_x)/a;
k_yy=4*linspace(-1*pi,1*pi,count_y)/a;

%% 颜色
colorrd=[207,108,82]/255;
coloryl=[241,141,0]/255;
colorbl=[109,152,165]/255;

%% 计算
bar = waitbar(0,'色散关系计算中...'); %显示进度条
for x = 1:count_x
    per=x/count_x;
    str=['色散关系计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1:count_y
        k_x=k_xx(x);
        k_y=k_yy(y);
        k=[k_x k_y];
        % Hamiltonian
        H=zeros(2);
        H(1,2)=t*exp(-1i*k*e1)+t*exp(-1i*k*e2)+t*exp(-1i*k*e3);
        H(2,1)=conj(H(1,2));
        E(x,y,:)=eig(H);
    end
end

%% 绘图
% 设置图窗
f1=figure(1);
position0=get(f1,'position');
set(f1,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

%结果
nexttile
pcolor(k_xx/(pi/a),k_yy/(pi/a),E(:,:,1)')
hold on
pcolor(k_xx/(pi/a),k_yy/(pi/a),E(:,:,2)')
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
clim([0,max(max(max(E)))])
load('colorData.mat');
axis([min(k_xx/(pi/a)),max(k_xx/(pi/a)),min(k_yy/(pi/a)),max(k_yy)/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[min(k_xx/(pi/a)):1:max(k_xx/(pi/a))], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[min(k_yy/(pi/a)):1:max(k_yy/(pi/a))], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
colormap(Purples6);

%% 第一布里渊区

% 计算倒格子的基矢量
b1 = 2 * pi * [a2(2), -a2(1)] / (a1(1)*a2(2) - a1(2)*a2(1));
b2 = 2 * pi * [-a1(2), a1(1)] / (a1(1)*a2(2) - a1(2)*a2(1));

% 生成倒格子的点
point_N = 5; % 在每个方向生成的点的数量
points = [];
for i = -point_N:point_N
    for j = -point_N:point_N
        point = i * b1 + j * b2;
        points = [points; point];
    end
end

% 计算到原点的距离
distances = sqrt(sum((points - [0,0]).^2, 2));

% 找到距离最小的6个点
[~, indices] = mink(distances, 7);
nearest_points = points(indices, :);

% 计算第一布里渊区的顶点
K = convhull(nearest_points);

% 绘制第一布里渊区
plot(nearest_points(K, 1)/(pi/a), nearest_points(K, 2)/(pi/a),':','Color','k');
K1=[K(1),K(3),K(5),K(1)];
K2=[K(2),K(4),K(6),K(2)];
plot(nearest_points(K1, 1)/(pi/a), nearest_points(K1, 2)/(pi/a),':','Color','k');
plot(nearest_points(K2, 1)/(pi/a), nearest_points(K2, 2)/(pi/a),':','Color','k');

% 绘制倒格子点
plot(points(:, 1)/(pi/a), points(:, 2)/(pi/a), '.','Color','black');

% 正六边形
n = 6;
radius = 4*sqrt(3)/3/a;
theta = linspace(0, 2*pi, n+1);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x, y, '--','LineWidth',1.5,'Color',colorrd);

close(bar)%关闭进度条

贝里曲率（晶格规范）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
clc
clear all

%% 参数
e_0=1;
a=sqrt(3)*e_0;
a1=a*[1,0]';
a2=a*[1/2,-sqrt(3)/2]';
t=1;

% 颜色
colorrd = [207,108,82]/255;
coloryl = [241,141,0]/255;
colorbl = [109,152,165]/255;

% 循环次数
count_x = 200;
count_y = 200;

% 区间
[x_min,x_max]=deal(4*-pi/a,4*pi/a);
[y_min,y_max]=deal(4*-pi/a,4*pi/a);

%处理数据
delta_k_x = (x_max-x_min)/count_x;
delta_k_y = (y_max-y_min)/count_y;
k_x_valuex = (x_min:delta_k_x:x_max);
k_y_valuey = (y_min:delta_k_y:y_max);
% n = 1;

%% 计算
bar = waitbar(0,'贝里曲率计算中...'); %显示进度条
for x = 1: count_x+1
    per=x/(4*count_x);
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    k_x = k_x_valuex(x);
    for y = 1: count_y+1
        k_y = k_y_valuey(y);
        k = [k_x, k_y];
        % Hamiltonian
        H(1,2,x,y)=t+t*exp(-1i*k*a1)+t*exp(-1i*k*a2);
        H(2,1,x,y)=conj(H(1,2));
    end
end
for x = 1: count_x
    per=(count_x+x)/(4*count_x);;
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1: count_y
        H_k_x_delta_k_x(:,:,x,y) = H(:,:,x+1,y);
        H_k_y_delta_k_y(:,:,x,y) = H(:,:,x,y+1);
        H__k_x_k_y_delta_k_x_k_y(:,:,x,y) = H(:,:,x+1,y+1);
    end
end
H(:,:,:,count_y+1)=[];
H(:,:,count_x+1,:)=[];

for n = 1:2
for x = 1: count_x
    per=((n+1)*count_x+x)/(4*count_x);;
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1: count_y
        [vecs, Energy(x,y,:)] = eigen(H(:,:,x,y));
        [vecs_k_x_delta_k_x, Energy_k_x_delta_k_x(x,y,:)] = eigen(H_k_x_delta_k_x(:,:,x,y));
        [vecs_k_y_delta_k_y, Energy_k_y_delta_k_y(x,y,:)] = eigen(H_k_y_delta_k_y(:,:,x,y));
        [vecs__k_x_k_y_delta_k_x_k_y, Energy__k_x_k_y_delta_k_x_k_y(x,y,:)] = eigen(H__k_x_k_y_delta_k_x_k_y(:,:,x,y));

        % 归一化
        U_x(n) = dot(vecs(:,n), vecs_k_x_delta_k_x(:,n)) / abs(dot(vecs(:,n), vecs_k_x_delta_k_x(:,n)));
        U_y(n) = dot(vecs(:,n), vecs_k_y_delta_k_y(:,n)) / abs(dot(vecs(:,n), vecs_k_y_delta_k_y(:,n)));
        U_x_y(n) = dot(vecs_k_y_delta_k_y(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)) / abs(dot(vecs_k_y_delta_k_y(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)));
        U_y_x(n) = dot(vecs_k_x_delta_k_x(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)) / abs(dot(vecs_k_x_delta_k_x(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)));

        F_tilde(x,y,n) = log(U_x(n) * U_y_x(n) / (U_x_y(n) * U_y(n)));
        Berry_Curvature(x,y,n) = real(1i * F_tilde(x,y,n) / (delta_k_x * delta_k_y));
        Chern_number_point(x,y,n)=-F_tilde(x,y,n)/(2*pi*1i); % 从Field_strength求陈数
    end
end
Chern_number(n)=sum(Chern_number_point(:,:,n),"all")
end


%% 绘图
% 设置图窗
f1=figure(1);
position0=get(f1,'position');
set(f1,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

% 绘制曲率
nexttile
surf(k_x_valuex(1:count_x)/(pi/a), k_y_valuey(1:count_y)/(pi/a), Berry_Curvature(:,:,1))
hold on
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
load('colorData.mat');
axis([min(k_x_valuex(1:count_x)/(pi/a)),max(k_x_valuex(1:count_x)/(pi/a)),min(k_y_valuey(1:count_y)/(pi/a)),max(k_y_valuey(1:count_y))/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[x_min/(pi/a):1:x_max/(pi/a)], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[y_min/(pi/a):1:y_max/(pi/a)], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
view(-44.5,85);
load('colorData.mat');
colormap(Purples6);

f2=figure(2);
position0=get(f2,'position');
set(f2,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

nexttile
surf(k_x_valuex(1:count_x)/(pi/a), k_y_valuey(1:count_y)/(pi/a), Berry_Curvature(:,:,2))
hold on
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
load('colorData.mat');
axis([min(k_x_valuex(1:count_x)/(pi/a)),max(k_x_valuex(1:count_x)/(pi/a)),min(k_y_valuey(1:count_y)/(pi/a)),max(k_y_valuey(1:count_y))/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[x_min/(pi/a):1:x_max/(pi/a)], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[y_min/(pi/a):1:y_max/(pi/a)], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
view(-44.5,85);
colormap(Purples6);

close(bar)%关闭进度条
%% 求本征值并排序
function [vecs, Energy] = eigen(H)
    [v, E] = eig(H);
    [Energy, index] = sort(real(diag(E)));
    vecs = v(:, index);
end

贝里曲率（原子规范）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
clc
clear all

%% 参数
e_0=1;
a=sqrt(3)*e_0;
e1=e_0*[0,1]';
e2=e_0*[sqrt(3)/2,-0.5]';
e3=e_0*[-sqrt(3)/2,-0.5]';
t=1;

% 颜色
colorrd = [207,108,82]/255;
coloryl = [241,141,0]/255;
colorbl = [109,152,165]/255;

% 循环次数
count_x = 200;
count_y = 200;

% 区间
[x_min,x_max]=deal(2*-pi,2*pi);
[y_min,y_max]=deal(2*-pi,2*pi);

%处理数据
delta_k_x = (x_max-x_min)/count_x;
delta_k_y = (y_max-y_min)/count_y;
k_x_valuex = (x_min:delta_k_x:x_max);
k_y_valuey = (y_min:delta_k_y:y_max);
% n = 1;

%% 计算
bar = waitbar(0,'贝里曲率计算中...'); %显示进度条
for x = 1: count_x+1
    per=x/(4*count_x);
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    k_x = k_x_valuex(x);
    for y = 1: count_y+1
        k_y = k_y_valuey(y);
        k = [k_x, k_y];
        % Hamiltonian
        H(1,2,x,y)=t*exp(-1i*k*e1)+t*exp(-1i*k*e2)+t*exp(-1i*k*e3);
        H(2,1,x,y)=conj(H(1,2));
    end
end
for x = 1: count_x
    per=(count_x+x)/(4*count_x);;
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1: count_y
        H_k_x_delta_k_x(:,:,x,y) = H(:,:,x+1,y);
        H_k_y_delta_k_y(:,:,x,y) = H(:,:,x,y+1);
        H__k_x_k_y_delta_k_x_k_y(:,:,x,y) = H(:,:,x+1,y+1);
    end
end
H(:,:,:,count_y+1)=[];
H(:,:,count_x+1,:)=[];

for n = 1:2
for x = 1: count_x
    per=((n+1)*count_x+x)/(4*count_x);;
    str=['贝里曲率计算中...',num2str(100*per),'%']; % 百分比形式显示
    waitbar(per,bar,str)
    for y = 1: count_y
        [vecs, Energy(x,y,:)] = eigen(H(:,:,x,y));
        [vecs_k_x_delta_k_x, Energy_k_x_delta_k_x(x,y,:)] = eigen(H_k_x_delta_k_x(:,:,x,y));
        [vecs_k_y_delta_k_y, Energy_k_y_delta_k_y(x,y,:)] = eigen(H_k_y_delta_k_y(:,:,x,y));
        [vecs__k_x_k_y_delta_k_x_k_y, Energy__k_x_k_y_delta_k_x_k_y(x,y,:)] = eigen(H__k_x_k_y_delta_k_x_k_y(:,:,x,y));

        % 归一化
        U_x(n) = dot(vecs(:,n), vecs_k_x_delta_k_x(:,n)) / abs(dot(vecs(:,n), vecs_k_x_delta_k_x(:,n)));
        U_y(n) = dot(vecs(:,n), vecs_k_y_delta_k_y(:,n)) / abs(dot(vecs(:,n), vecs_k_y_delta_k_y(:,n)));
        U_x_y(n) = dot(vecs_k_y_delta_k_y(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)) / abs(dot(vecs_k_y_delta_k_y(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)));
        U_y_x(n) = dot(vecs_k_x_delta_k_x(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)) / abs(dot(vecs_k_x_delta_k_x(:,n), vecs__k_x_k_y_delta_k_x_k_y(:,n)));

        F_tilde(x,y,n) = log(U_x(n) * U_y_x(n) / (U_x_y(n) * U_y(n)));
        Berry_Curvature(x,y,n) = real(1i * F_tilde(x,y,n) / (delta_k_x * delta_k_y));
        Chern_number_point(x,y,n)=-F_tilde(x,y,n)/(2*pi*1i); % 从Field_strength求陈数
    end
end
Chern_number(n)=sum(Chern_number_point(:,:,n),"all")
end


%% 绘图
% 设置图窗
f1=figure(1);
position0=get(f1,'position');
set(f1,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

% 绘制曲率
nexttile
surf(k_x_valuex(1:count_x)/(pi/a), k_y_valuey(1:count_y)/(pi/a), Berry_Curvature(:,:,1))
hold on
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
load('colorData.mat');
axis([min(k_x_valuex(1:count_x)/(pi/a)),max(k_x_valuex(1:count_x)/(pi/a)),min(k_y_valuey(1:count_y)/(pi/a)),max(k_y_valuey(1:count_y))/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[x_min/(pi/a):1:x_max/(pi/a)], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[y_min/(pi/a):1:y_max/(pi/a)], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
view(-44.5,85);
load('colorData.mat');
colormap(Purples6);

f2=figure(2);
position0=get(f2,'position');
set(f2,'position',position0+[0*position0(3),0*position0(4),-1*position0(3)+1.05*position0(4),0*position0(4)]);

nexttile
surf(k_x_valuex(1:count_x)/(pi/a), k_y_valuey(1:count_y)/(pi/a), Berry_Curvature(:,:,2))
hold on
shading interp
grid off
set(gca,'linewidth',1)
set(gca,'boxstyle','full','box','on')
set(gca,'xticklabel','','yticklabel','','zticklabel','','xtick','','ytick','','ztick','')
xlabel('$k_{x}\ (\pi/a)$','Interpreter','latex','FontSize',22);
ylabel('$k_{y}\ (\pi/a)$','interpreter','latex','FontSize',22,'FontName','Arial')
load('colorData.mat');
axis([min(k_x_valuex(1:count_x)/(pi/a)),max(k_x_valuex(1:count_x)/(pi/a)),min(k_y_valuey(1:count_y)/(pi/a)),max(k_y_valuey(1:count_y))/(pi/a)]*1.05)
xtickangle(0);
set(gca,'xminortick','on');%小刻度打开
set(gca,'yminortick','on');%小刻度打开
set(gca,'ticklength',[0.015 0.025]);%刻度长度
set(gca, 'xtick',[x_min/(pi/a):1:x_max/(pi/a)], 'XTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
set(gca, 'ytick',[y_min/(pi/a):1:y_max/(pi/a)], 'yTickLabel', {'$-4$', '', '', '', '$0$', '', '', '','$4$'}, 'TickLabelInterpreter', 'latex', 'FontSize', 18)
view(-44.5,85);
colormap(Purples6);

close(bar)%关闭进度条
%% 求本征值并排序
function [vecs, Energy] = eigen(H)
    [v, E] = eig(H);
    [Energy, index] = sort(real(diag(E)));
    vecs = v(:, index);
end