Create a function to choose spatial dimensions of twisted atomic layers

Dear all

I am analyzing a situation on which an atomic layer is rotated with respect to other one by a given angle. These atomic sets are given, initially, by two small in-plane 2d space data, given by the files "Magnetic_AA_NCS.txt" and "Non_Magnetic_AA_NCS.txt". As you can see in the code below, I extend the spatial dimensions of them, then I perform a rotation with respect to the center of the sample, and then I try to keep only points along the x- and y-th spatial directions within a range with respect to the rotational center:

clc;
clear all;
close all force;
%% AA bilayer (atomic positions)
% Magnetic atoms (8 Cr atoms. First, second, and third columns represent
% the spatial coordinates along the orthogonal x-, y-, and z-th basis. The
% coordinates are given in Angstroms)
AA_magnetic_coordinates=readmatrix(strcat(pwd,'\Different-Bilayer-Stackings\Magnetic_AA_NCS.txt')); % Angstroms
% Non-magnetic atoms (24 I atoms. First, second, and third columns represent
% the spatial coordinates along the orthogonal x-, y-, and z-th basis. The
% coordinates are given in Angstroms)
AA_non_magnetic_coordinates=readmatrix(strcat(pwd,'\Different-Bilayer-Stackings\Non_Magnetic_AA_NCS.txt')); % Angstroms
%% Additionally, let's upload relevant information for each of these two stackings
% First line: three columns, which gives the size, in Angstroms, along the
% x-th, y-th, and z-th spatial directions
% Second line: first column: 1st nearest intralayer distance (Angstroms); second
% column: 2nd nearest intralayer distance (Angstroms)
% Third line: one element, with the lattice parameter a0 (Angstroms)
info_AA=dlmread(strcat(pwd,'\Different-Bilayer-Stackings\Relevant_Information_AA_NCS.txt')); % Angstroms
%% Now, let's expand the atomic positions along the x- and y-th spatial directions
% Considering that what we have in the AA is a sort of (0,0,0) reference that 
% can be replicated in space if needed, let's add to the system magnetic 
% and non-magnetic atoms. For that, we have to expand the in-plane atomic 
% space through the information in the first and second column of the first 
% line of the variable info_AA
counter_Cr=0;
counter_I=0;
for dx=[0:1:110]
    for dy=[0:1:200]
        for i=1:length(AA_magnetic_coordinates(:,1))
            counter_Cr=counter_Cr+1;
            AA_magnetic_supercell(counter_Cr,1)=AA_magnetic_coordinates(i,1)+dx*info_AA(1,1); % Angstroms
            AA_magnetic_supercell(counter_Cr,2)=AA_magnetic_coordinates(i,2)+dy*info_AA(1,2); % Angstroms
            AA_magnetic_supercell(counter_Cr,3)=AA_magnetic_coordinates(i,3); % Angstroms
        end
        for i=1:length(AA_non_magnetic_coordinates(:,1))
            counter_I=counter_I+1;
            AA_non_magnetic_supercell(counter_I,1)=AA_non_magnetic_coordinates(i,1)+dx*info_AA(1,1); % Angstroms
            AA_non_magnetic_supercell(counter_I,2)=AA_non_magnetic_coordinates(i,2)+dy*info_AA(1,2); % Angstroms
            AA_non_magnetic_supercell(counter_I,3)=AA_non_magnetic_coordinates(i,3); % Angstroms
        end
    end
end
clear dx dy counter_Cr counter_I
%% Now, let's characterize the sorted coordinates along the z-th spatial direction
sorted_z=unique([vertcat(unique(AA_magnetic_coordinates(:,3)),unique(AA_non_magnetic_coordinates(:,3)))]); % Angstroms
%% Let's rearrange the data set so that the current midpoint goes to (0,0)
shift=[(min(AA_magnetic_supercell(:,1))+max(AA_magnetic_supercell(:,1)))/2 (min(AA_magnetic_supercell(:,2))+max(AA_magnetic_supercell(:,2)))/2];
% Global magnetic system
AA_magnetic_supercell(:,1)=AA_magnetic_supercell(:,1)-shift(1); % Angstroms
AA_magnetic_supercell(:,2)=AA_magnetic_supercell(:,2)-shift(2); % Angstroms
% Global non_magnetic system
AA_non_magnetic_supercell(:,1)=AA_non_magnetic_supercell(:,1)-shift(1); % Angstroms
AA_non_magnetic_supercell(:,2)=AA_non_magnetic_supercell(:,2)-shift(2); % Angstroms
clear shift
%% Now, we are prepared to rotate one monolayer with respect to the other
% Let's define an arbitrary twisting angle, theta
theta=2; % º
% Now, we have to apply the rotation, around (0,0) to all the magnetic and 
% non-magnetic atoms of the top monolayer
for i=1:length(AA_magnetic_supercell(:,1))
    if AA_magnetic_supercell(i,3)==sorted_z(5)
        % Let's apply the rotation
        AA_magnetic_supercell(i,1)=AA_magnetic_supercell(i,1).*cosd(theta)-AA_magnetic_supercell(i,2).*sind(theta); % Angstroms
        AA_magnetic_supercell(i,2)=AA_magnetic_supercell(i,1).*sind(theta)+AA_magnetic_supercell(i,2).*cosd(theta); % Angstroms
    end
end
for i=1:length(AA_non_magnetic_supercell(:,1))
    if AA_non_magnetic_supercell(i,3)>sorted_z(3)
        % Let's apply the rotation
        AA_non_magnetic_supercell(i,1)=AA_non_magnetic_supercell(i,1).*cosd(theta)-AA_non_magnetic_supercell(i,2).*sind(theta); % Angstroms
        AA_non_magnetic_supercell(i,2)=AA_non_magnetic_supercell(i,1).*sind(theta)+AA_non_magnetic_supercell(i,2).*cosd(theta); % Angstroms
    end
end
clear theta
%% Now, let's separate the magnetic atoms of each monolayer
counter_Cr1=0;
counter_Cr2=0;
for i=1:length(AA_magnetic_supercell(:,1))
    if AA_magnetic_supercell(i,3)==sorted_z(2)
        counter_Cr1=counter_Cr1+1;
        Cr1_atoms(counter_Cr1,:)=AA_magnetic_supercell(i,1:3); % Angstroms
    elseif AA_magnetic_supercell(i,3)==sorted_z(5)
        counter_Cr2=counter_Cr2+1;
        Cr2_atoms(counter_Cr2,:)=AA_magnetic_supercell(i,1:3); % Angstroms  
    end
end
clear counter_Cr1 counter_Cr2
%% Now, let's plot the magnetic atoms
u1=figure(1)
scatter3(Cr1_atoms(:,1),Cr1_atoms(:,2),Cr1_atoms(:,3),'MarkerEdgeColor','b','MarkerFaceColor','b');
hold on
scatter3(Cr2_atoms(:,1),Cr2_atoms(:,2),Cr2_atoms(:,3),'MarkerEdgeColor','r','MarkerFaceColor','r');
plot((min(AA_magnetic_supercell(:,1))+max(AA_magnetic_supercell(:,1)))/2,(min(AA_magnetic_supercell(:,2))+max(AA_magnetic_supercell(:,2)))/2,'x','HandleVisibility','off');
l1=legend('Cr(1)','Cr(2)','FontSize',14,'Location','northeastoutside');
set(l1,'interpreter','latex','FontSize',14);
xlim([min(AA_magnetic_supercell(:,1)) max(AA_magnetic_supercell(:,1))]);
ylim([min(AA_magnetic_supercell(:,2)) max(AA_magnetic_supercell(:,2))]);
view([0 90]);
box on;
grid off;
% pbaspect([1 3.3868/5.9900 1]);
clear u1 l1 Cr1_atoms Cr2_atoms AA_magnetic_coordinates AA_non_magnetic_coordinates
%% Now, let's separate the magnetic atoms of each monolayer
counter_I1=0;
counter_I2=0;
for i=1:length(AA_non_magnetic_supercell(:,1))
    if AA_non_magnetic_supercell(i,3)==sorted_z(1) || AA_non_magnetic_supercell(i,3)==sorted_z(3)
        counter_I1=counter_I1+1;
        I1_atoms(counter_I1,:)=AA_non_magnetic_supercell(i,1:3); % Angstroms
    elseif AA_non_magnetic_supercell(i,3)==sorted_z(4) || AA_non_magnetic_supercell(i,3)==sorted_z(5)
        counter_I2=counter_I2+1;
        I2_atoms(counter_I2,:)=AA_non_magnetic_supercell(i,1:3); % Angstroms  
    end
end
clear counter_I1 counter_I2
%% Now, let's plot the non-magnetic atoms
u1=figure(2)
scatter3(I1_atoms(:,1),I1_atoms(:,2),I1_atoms(:,3),'MarkerEdgeColor','b','MarkerFaceColor','b');
hold on
scatter3(I2_atoms(:,1),I2_atoms(:,2),I2_atoms(:,3),'MarkerEdgeColor','r','MarkerFaceColor','r');
plot((min(AA_non_magnetic_supercell(:,1))+max(AA_non_magnetic_supercell(:,1)))/2,(min(AA_non_magnetic_supercell(:,2))+max(AA_non_magnetic_supercell(:,2)))/2,'x','HandleVisibility','off');
l1=legend('I(1)','I(2)','FontSize',14,'Location','northeastoutside');
set(l1,'interpreter','latex','FontSize',14);
xlim([min(AA_non_magnetic_supercell(:,1)) max(AA_non_magnetic_supercell(:,1))]);
ylim([min(AA_non_magnetic_supercell(:,2)) max(AA_non_magnetic_supercell(:,2))]);
view([0 90]);
box on;
grid off;
% pbaspect([1 3.3868/5.9900 1]);
clear u1 l1 I1_atoms I2_atoms AA_non_magnetic_coordinates
%% Now, let's keep only the atoms in the desired interval
% Magnetic atoms
counter=0;
for i=1:length(AA_magnetic_supercell(:,1))
    if abs(AA_magnetic_supercell(i,1))<=500 && abs(AA_magnetic_supercell(i,2))<=500
        counter=counter+1;
        magnetic_atoms(counter,1:3)=AA_magnetic_supercell(i,:); % Angstroms
        if AA_magnetic_supercell(i,3)==sorted_z(2)
            magnetic_atoms(counter,4)=0; % Untwisted monolayer
        else
            magnetic_atoms(counter,4)=1; % Twisted monolayer
        end
    end
end
clear counter AA_magnetic_supercell
% Non-magnetic atoms
counter=0;
for i=1:length(AA_non_magnetic_supercell(:,1))
    if abs(AA_non_magnetic_supercell(i,1))<=550 && abs(AA_non_magnetic_supercell(i,2))<=550
        counter=counter+1;
        non_magnetic_atoms(counter,1:3)=AA_non_magnetic_supercell(i,:); % Angstroms
    end
end
clear counter AA_non_magnetic_supercell
%% Now, let's shift the origin
shift=[abs(min(magnetic_atoms(:,1))) abs(min(magnetic_atoms(:,2)))]; % Angstroms
% Magnetic atoms
magnetic_atoms(:,1:2)=magnetic_atoms(:,1:2)+shift; % Angstroms
% Non-magnetic atoms
non_magnetic_atoms(:,1:2)=non_magnetic_atoms(:,1:2)+shift; % Angstroms
clear shift
%% Now, let's create the reduced magnetic atoms-based matrix
reduced_magnetic_atoms(:,1)=[0:1:(length(magnetic_atoms(:,1))-1)];
reduced_magnetic_atoms(:,2)=magnetic_atoms(:,1)./info_AA(1,1);
reduced_magnetic_atoms(:,3)=magnetic_atoms(:,2)./info_AA(1,2);
reduced_magnetic_atoms(:,4)=magnetic_atoms(:,3)./info_AA(1,3);
reduced_magnetic_atoms(:,5)=magnetic_atoms(:,4);
reduced_magnetic_atoms(:,6:7)=0;
%% Now, let's save the atomic positions
% Magnetic atoms
writematrix(magnetic_atoms(:,1:3),'Twisted_Monolayer_Magnetic_Atoms_Twisting_2_Degrees.txt','delimiter','space');
writematrix(reduced_magnetic_atoms,'Twisted_Monolayer_Reduced_Magnetic_Atoms_Twisting_2_Degrees.txt','delimiter','space');
% Non-magnetic atoms
writematrix(non_magnetic_atoms,'Twisted_Monolayer_Non_Magnetic_Atoms_Twisting_2_Degrees.txt','delimiter','space');
%% Now, let's plot the magnetic atoms
u1=figure(3)
plot(magnetic_atoms(:,1),magnetic_atoms(:,2),'MarkerEdgeColor','b','MarkerFaceColor','b');
xlim([min(magnetic_atoms(:,1)) max(magnetic_atoms(:,1))]);
ylim([min(magnetic_atoms(:,2)) max(magnetic_atoms(:,2))]);
box on;
clear u1
%% Now, let's plot the non-magnetic atoms
u1=figure(4)
plot(non_magnetic_atoms(:,1),non_magnetic_atoms(:,2),'MarkerEdgeColor','b','MarkerFaceColor','b');
xlim([min(non_magnetic_atoms(:,1)) max(non_magnetic_atoms(:,1))]);
ylim([min(non_magnetic_atoms(:,2)) max(non_magnetic_atoms(:,2))]);
box on;

In this case, I tried to guarantee that there were magnetic atoms belonging to both layers within a range of 100x100 nm2, with the center (50,50) nm corresponding with respect to which the rotation was performed. For the non-magnetic atoms, I saved the atoms within a range [-5,105]x[-5,105] nm2. So, basically, far enough from the borders of the layers (it does not need to that much as in the figure below):

As you can see, I kind of guess the needed dx and dy to fulfill my requirements. Could be a good strategy to create a function where I only provide the desired dimensions, in nm, along x- and y-th directions and the twisting angle? It is not needed to make the plots that I included in my code, that was just to check I had something reasonable according to my requirements.

Answer 1

[Ayush]

Hi Richard,

I am assuming that in the context of your problem, dx and dy are expansion factors that determine how many times you need to replicate or expand your original atomic lattice in the x- and y-directions to ensure that the rotated area fully covers your desired dimensions.

I have written an example function calculateExpansion to get the dx and dy based on desired_x, desired_y and rotation_center.

% Example usage
desired_x = 100;  % in nm
desired_y = 100;  % in nm
angle_deg = 30;   % degrees
rotation_center = [50, 50];  % center of rotation

[dx, dy] = calculateExpansion(desired_x, desired_y, angle_deg, rotation_center);
fprintf('Required expansion: dx = %d, dy = %d\n', dx, dy);

Required expansion: dx = 2, dy = 2

function [dx, dy] = calculateExpansion(desired_x, desired_y, angle_deg, rotation_center)
    % Convert angle from degrees to radians
    angle_rad = deg2rad(angle_deg);
    
    % Calculate the cosine and sine of the angle
    cos_angle = cos(angle_rad);
    sin_angle = sin(angle_rad);
    
    % Rotation center
    cx = rotation_center(1);
    cy = rotation_center(2);
    
    % Define the corners of the rectangle
    corners = [
        0, 0;
        desired_x, 0;
        desired_x, desired_y;
        0, desired_y
    ];
    
    % Translate corners to rotation center
    translated_corners = corners - [cx, cy];
    
    % Perform the rotation
    rotation_matrix = [
        cos_angle, -sin_angle;
        sin_angle, cos_angle
    ];
    rotated_corners = (rotation_matrix * translated_corners')';
    
    % Translate corners back
    final_corners = rotated_corners + [cx, cy];
    
    % Calculate the bounding box of the rotated rectangle
    min_x = min(final_corners(:, 1));
    max_x = max(final_corners(:, 1));
    min_y = min(final_corners(:, 2));
    max_y = max(final_corners(:, 2));
    
    % Calculate the needed expansion in x and y directions
    new_x = max_x - min_x;
    new_y = max_y - min_y;
    
    dx = ceil(new_x / desired_x);
    dy = ceil(new_y / desired_y);
end

Let me know if this helps!