扭角进阶30和0度为对称轴

分为对称轴30和0度

%%%Nie-Wei,Wang 2023-12-11 超胞构建摩尔角，消除晶格失配平
%%%角度是算出来的，不能用于计算
%%限定0-30容易丢数据
clear all

a1=[1;0]; % a1,a2是基矢
a2=[-1/2;sqrt(3)/2];



%%%%%%%%%%%以30°为轴
angle_60=a1+a2;   %60
angle_30=2*a1+a2; %30度
angle_0=a1;       %0度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%以0°为轴
angle_60=2*a1+a2;   %60
angle_30=a1; %30度
angle_0=-a2;       %0度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num=0;      %计数
beta=120.0; %%bete是a1和a2的夹角
N=20;       %遍历
Rz=[cosd(beta) -sind(beta); sind(beta) cosd(beta)];  %%120绕z轴旋转矩阵
for m=-N:N
    for n=-N:N
        for p=-N:N
            for q=-N:N
                b1=m*a1+n*a2;
                b2=p*a1+q*a2;
                theta=acosd(dot(b1,b2)/(norm(b1)*norm(b2)));
                %%%%%%%%%%%%%%%
                b1_angle_max=cross([b1;0],[angle_30;0]);  %%% b1 0-30
                b1_angle_min=cross([b1;0],[angle_0;0]);  %%%%
                b2_angle_max=cross([b2;0],[angle_60;0]);  %%% b2 30-60
                b2_angle_min=cross([b2;0],[angle_30;0]);  %%%%
                b2b1=cross([b2;0],[b1;0]);   %%%%            b2在b1的顺时针方向
                type1=cross([b1;0],[angle_30;0]);
                type2=cross([b2;0],[angle_30;0]);
                theta1=acosd(dot(b1,angle_30)/(norm(b1)*norm(angle_30)))*(type1(3)/norm(type1));
                theta2=acosd(dot(b2,angle_30)/(norm(b2)*norm(angle_30)))*(type2(3)/norm(type2));
                %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
               if  abs(norm(b1)-norm(b2))<0.5 && b1_angle_max(3)>0 && b1_angle_min(3)<0 && b2_angle_max(3)>0 && b2_angle_min(3)<0 && b2b1(3)<0  ...
                       && abs(theta1+theta2)<0.1 %&& abs(norm(b2-b1)-2*norm(a2))<0.01 %（差为整数，为1个b2 abs(norm(b2-b1)-norm(a2))<0.01）
%%%%abs(theta1+theta2)<0.1 设置30/0度为旋转对称轴
                   
                                num=num+1;
                                a(num,:)=[m,n,p,q,theta,4*norm(b1)^2,norm(b1)^2];
                                %%%%%%  旋转矩阵   总角度  b1_30夹角  b2_30夹角   原子数
                           
               end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            end  
        end
    end
end
 

% 按照第九列进行排序
% 提取第九列和第十二列
% 按照第五列的值进行排序
data=a;
sortedData = sortrows(data, 5);
% 保留第五列小数点5位
data(:, 5) = round(data(:, 5), 1);

% 按照第五列的值进行排序
[~, idx] = sort(data(:, 6));
sortedData = data(idx, :);

% % 找到相同值中第六列最小的那一行
% [~, uniqueIdx, ~] = unique(sortedData(:, 6), 'stable');
% resultData = sortedData(uniqueIdx, :);
% 
% % 显示结果
% disp(resultData);

30°结果  %%%%%% 旋转矩阵  夹角  原子数  超胞的晶格常数的平方

%%%%%% 旋转矩阵 总角度 b1_30夹角 b2_30夹角 原子数

#root7 30±8.2°
root13 30±2.2°