    
 %% Sorting first perovskite block  
    idxy = mod(atom_pos.posRefineM(:,8),1) == 0 & mod(atom_pos.posRefineM(:,9),1) == 0; % find index (0,0)
    up = atom_pos.posRefineM(idxy,8);
    vp = atom_pos.posRefineM(idxy,9);
    
    Ca_list = zeros(length(up),6);% create a list containing coordinates of Ca atoms
    Ca_list(:,1) = atom_pos.posRefineM(idxy,3);
    Ca_list(:,2) = atom_pos.posRefineM(idxy,4);    
    imcomp_cell = [];
    
    %iteratively find Ca atoms with index (0,0.2) and (0,0.4)
    for ii = 1:length(up)
        try
            id1 = atom_pos.posRefineM(:,8) == up(ii,1)+0.2 & atom_pos.posRefineM(:,9) == (vp(ii,1));% find (0,0.2)
            Ca_list(ii,3) = atom_pos.posRefineM(id1,3);
            Ca_list(ii,4) = atom_pos.posRefineM(id1,4);

            id2 = atom_pos.posRefineM(:,8) == up(ii,1)+0.4 & atom_pos.posRefineM(:,9) == (vp(ii,1));% find (0,0.4)
            Ca_list(ii,5) = atom_pos.posRefineM(id2,3);
            Ca_list(ii,6) = atom_pos.posRefineM(id2,4);        

        catch
            imcomp_cell = [imcomp_cell ii];% find incomplete unit cells
        end

    end
    Ca_list(imcomp_cell,:) = [];%delete incomplete unit cells
    relative_shift = Ca_list(:,3:4)-(Ca_list(:,1:2)+Ca_list(:,5:6))/2;
    
    
    clearvars up vp idxy id1 id2 ii imcomp_cell
  %% Sorting second perovskite block    
    idxy = mod(atom_pos.posRefineM(:,8),1) == 0.5 & mod(atom_pos.posRefineM(:,9),1) == 0.5; % find index (0.5,0.5)
    up = atom_pos.posRefineM(idxy,8);
    vp = atom_pos.posRefineM(idxy,9);
    Ca_list2 = zeros(length(up),6);
    Ca_list2(:,1) = atom_pos.posRefineM(idxy,3);
    Ca_list2(:,2) = atom_pos.posRefineM(idxy,4);    
    imcomp_cell = [];
    
    
    for ii = 1:length(up)
        try
            id1 = atom_pos.posRefineM(:,8) == up(ii,1)+0.2 & atom_pos.posRefineM(:,9) == (vp(ii,1));% find (0.7,0.5)
            Ca_list2(ii,3) = atom_pos.posRefineM(id1,3);
            Ca_list2(ii,4) = atom_pos.posRefineM(id1,4);

            id2 = atom_pos.posRefineM(:,8) == up(ii,1)+0.4 & atom_pos.posRefineM(:,9) == (vp(ii,1));% find (0.9,0.5)
            Ca_list2(ii,5) = atom_pos.posRefineM(id2,3);
            Ca_list2(ii,6) = atom_pos.posRefineM(id2,4);        
            
            
        catch
            imcomp_cell = [imcomp_cell ii];
        end

    end
    
    Ca_list2(imcomp_cell,:) = [];    
    
    %% Plot sorting results 
    figure(1);imagesc(atom_pos.image);axis image off;colormap(gray);...
        hold on;
        scatter(Ca_list(:,2),Ca_list(:,1),'r.');
        scatter(Ca_list(:,4),Ca_list(:,3),'g.');
        scatter(Ca_list(:,6),Ca_list(:,5),'b.');
        scatter(Ca_list2(:,2),Ca_list2(:,1),'rx');
        scatter(Ca_list2(:,4),Ca_list2(:,3),'gx');
        scatter(Ca_list2(:,6),Ca_list2(:,5),'bx');
        
    %% Calculate shift and plot vectors
    relative_shift = [relative_shift; (Ca_list2(:,3:4)-(Ca_list2(:,1:2)+Ca_list2(:,5:6))/2)];
    Ca_list = [Ca_list;Ca_list2];             
    figure(2);imagesc(atom_pos.image);axis image off;colormap(gray);...
        hold on;
        quiver(Ca_list(:,4),Ca_list(:,3),relative_shift(:,2),relative_shift(:,1),'Color','y','LineWidth',2,'AutoScaleFactor',0.5);