-
Notifications
You must be signed in to change notification settings - Fork 0
/
kDoF_tri.m
170 lines (127 loc) · 6.34 KB
/
kDoF_tri.m
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
classdef kDoF_tri < handle
%kDoF creates and stores the k-space Degrees of Freedom for a trilayer
%updated 11/22/2019, change how
properties
layer_list % collection of Layer.m objects (input)
max_k % cutoff in momentum space (input)
num_k % total number of k-points (output)
k_list % list of k-points (output)
grid_search; % maximum size of the searching square (input)
end
methods
function obj = kDoF_tri(Layers,max_k,grid_search)
% set input variables (see properties for more info)
obj.layer_list = Layers;
obj.max_k = max_k;
obj.grid_search = grid_search;
end
function [] = gen_dof(obj)
% void function to create the k-point list
clear obj.k_list;
clear obj.num_k;
grid_s = obj.grid_search;
k_idx = 0;
[~,ind] = min(abs([obj.layer_list(1).theta, obj.layer_list(3).theta]));
kcut = norm(obj.layer_list(2*(ind-1)+1).G(:,1)-obj.layer_list(2).G(:,1));
for i = 1:3
K(:,i) = 1/3 * (2*obj.layer_list(i).G(:,1) + obj.layer_list(i).G(:,2));
end
for l = 1:3
l_1 = l; % the layer we are in
l_2 = mod(l,3)+1;
l_3 = mod(l+1,3)+1;
kcut1 = obj.max_k*norm(K(:,2))*sqrt(3);
b1 = obj.layer_list(l_1).G;
b2 = obj.layer_list(l_2).G;
b3 = obj.layer_list(l_3).G;
% kcut2 = 4*norm(b1(:,2));
kcut2 = kcut1; % constraining the norm of monolayer G vectors
shift_vec=K(:,l_1)-K(:,2);
grid_s1 = grid_s;
grid_s2 = grid_s;
% if l == 3
% grid_s1 = 0;
% elseif l == 2
% grid_s2 = 0;
% else
% grid_s1 = 0;
% grid_s2 = 0;
% end
for n2_1 = -grid_s1:grid_s1
for n2_2 = -grid_s1:grid_s1
k2 = b2(:,1)*n2_1 + b2(:,2)*n2_2;
cond1 = norm(k2)<kcut2;
if cond1
for n3_1 = -grid_s2:grid_s2
for n3_2 = -grid_s2:grid_s2
k3 = b3(:,1)*n3_1 + b3(:,2)*n3_2;
k1 = b1(:,1)*(n3_1+n2_1) + b1(:,2)*(n3_2+n2_2);
k1p = (b2(:,1)*n3_1+b2(:,2)*n3_2)+(b3(:,1)*n2_1+b3(:,2)*n2_2);
cond2 = norm(k2-k3)<kcut1;
% cond2 = 1;
cond3 = norm(k3)<kcut2;
% cond3 = 1;
if cond2 && cond3
% mapping back to be near the origin of the given layer
% and the origin is the K point of layer 2 (shifting by shift_vec)
k_here = k2+k3-k1+shift_vec;
k_in(1) = k_here(1);
k_in(2) = k_here(2);
k_in(3) = l_1;
l1_index = 4 + (l_1 - 1)*2; % location of layer 1 in the data set
l2_index = 4 + (l_2 - 1)*2; % location of layer 2 in the data set
l3_index = 4 + (l_3 - 1)*2; % location of layer 3 in the data set
k_in(l1_index) = 0;
k_in(l1_index+1) = 0;
k_in(l2_index) = n2_1; % save grid index for this lattice
k_in(l2_index+1) = n2_2;
k_in(l3_index) = n3_1; % save grid index for this lattice
k_in(l3_index+1) = n3_2;
k_idx = k_idx+1;
k_data(k_idx,:) = k_in; % push back the new DoF
end
end
end
end
end
end
end
obj.num_k = size(k_data, 1);
% get rid of repeated indices
% disp(tar_dofs)
k_vec = k_data(:, 1:3);
disp(length(k_vec))
idx = 1;
idx_list = ones([length(k_vec), 1]);
tol = 1e-4;
for a = 1:2
for i = 1:size(k_vec,1)
k_here = k_vec(i,:);
if idx_list(i)
if a == 1
jbeg = 1;
jend = i-1;
else
jbeg = i+1;
jend = size(k_vec,1);
end
for j = jbeg:jend
if abs(k_vec(j,1) - k_here(1)) < tol && ...
abs(k_vec(j,2) - k_here(2)) < tol && abs(k_vec(j,3) - k_here(3)) < tol
if ~min( k_data(j,3:end) == [k_data(j,3) 0 0 0 0 0 0])
idx_list(j) = 0;
idx = idx + 1;
end
end
end
end
end
end
idx_list = boolean(idx_list);
disp(size(k_data))
k_data = k_data(idx_list, :);
% disp(size(k_data))
obj.k_list = k_data;
end
end
end