% Movie of the reconstruction of Si(111)7x7 starting with seven bulk atomic
% layers and allowing the top three to reconfigure to the 7x7
% reconstruction. The input data is generated by the program si7x7.m using
% atomic coordinates from the paper by K.D. Brommer et al., PRL, vol 68, pp
% 1355 - 1358 (1992). 

clear; close all;

si7x7 = load('Si7x7coords.dat'); 
% Arrange data in descending order of their vertical position 
si7x7 = sortrows(si7x7,3,'descend'); 

vid = VideoWriter('Si7x7.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

% Comment out one of the two below colormaps depending on your preference 

% % Customized colormap from dark grey through red, orange, and
% % yellow in order to highlight different layer heights 
% blackorange = customcolormap([0 0.35 0.88 1], ... 
%     {'#252525','#ff0000','#ff6a00','#ffff00'});
% blackorange = flip(blackorange);
% colormap(blackorange);

% Customized colormap from dark grey through the colors of the rainbow in 
% order to highlight different layer heights 
rainbow = customcolormap([0 1 2 3 4 5 6 7]/7, ... 
    {'#252525','#ff0000','#ff6a00','#ffff00','#00ff00','#00ffff','#0000ff','#7f00ff'});
rainbow = flip(rainbow);
colormap(rainbow);
colormap('jet'); 
c = jet; 

rt3 = sqrt(3);
rt24 = sqrt(24);

[x,y,z] = sphere(10);
r = 0.25;

xmin = 0;
xmax = 14*rt3;
ymin = -7;
ymax = 7;
zmin = -3.5;
zmax = 0.5; 

% Bulk atoms within the diamond reconstruction volume (7 x 7 x 3), composed 
% of three layers of 7x7, plus six further bulk atomic layers that remain
% untouched. The total number of atoms is thus 7 x 7 x 9 = 441. 
bulkCoords = zeros(size(441,1),4); % (x,y,z,index of associated (i.e., closest) surface atom)
indx = 1;
xshift = [0 0 -1 -1 -2 -2 -3 -3 -4]/rt3;
zshift = [1 4 5 8 9 12 13 16 17]/rt24;
for kk = 0:8 % z-layers 
    for jj = 0:6 % y-direction 
        for ii = 0:6 % x-direction 
            bulkCoords(indx,1) = ii*rt3/2 + jj*rt3/2 + xshift(kk+1);
            bulkCoords(indx,2) = (jj - ii)/2;
            bulkCoords(indx,3) = -zshift(kk+1);
            indx = indx+1;
        end
    end
end

for ii = 1:147 % Go through 147 surface atoms 
    % iadx = difference between x-coordinate of iith surface atom and all 
    % x-coordinates of the bulk atoms. Same for y- and z-coordinates 
    iadx = (si7x7(ii,1) - bulkCoords(:,1)); 
    iady = (si7x7(ii,2) - bulkCoords(:,2));
    iadz = (si7x7(ii,3) - bulkCoords(:,3)); 
    % Vector of the square of the interatomic distances between iith 
    % surface atom and all bulk atoms 
    iad = iadx.^2 + iady.^2 + iadz.^2; 
    [MMM,ind] = min(iad); % ind = index of bulk atom closest to surface atom ii
    
    while (bulkCoords(ind,4) ~= 0) % This bulk atom has already been selected 
        iad(ind) = iad(ind)+100; % Make this iad so large it will be ignored in next search 
        [MMM,ind] = min(iad); % Repeat search procedure to find new minimum
    end
    bulkCoords(ind,4) = ii;
end

recMaxHt = max(si7x7(:,3)); % Maximum height of atoms in reconstruction 
recMinHt = min(si7x7(:,3)); % Minimum height of atoms in reconstruction 

for tt = 0:0.005:1 % Loop over 201 frames between bulk surface and reconstructed surface 
    newplot
    tt % Monitor progress 
    tic % Monitor speed 
    hold off
    
    for jj = -2:2 % multiple diamond blocks along y-axis
        for kk = 0:2 % multiple diamond blocks along x-axis
            recShiftx = -3.5*rt3*mod(jj+1,2)+kk*7*rt3;
            recShifty = 3.5*jj;
            
            hold on 
            for ii = 1:size(bulkCoords,1)
                sindex = bulkCoords(ii,4); % Index of surface atom associated with iith bulk atom
                % Determine the coordinates of atom as it moves from the
                % bulk to the surface position
                if (bulkCoords(ii,4) ~= 0)
                    xposNow = bulkCoords(ii,1) + tt*(si7x7(sindex,1) ... 
                        - bulkCoords(ii,1)) + recShiftx;
                    yposNow = bulkCoords(ii,2) + tt*(si7x7(sindex,2) ... 
                        - bulkCoords(ii,2)) + recShifty;
                    zposNow = bulkCoords(ii,3) + tt*(si7x7(sindex,3) ... 
                        - bulkCoords(ii,3));
                    if (zposNow >= recMinHt)
                        % Index of colormap of surface atom at this point of 
                        % the development of the reconstruction
                        colNowIndex = ... 
                            1 + round(tt*255*(zposNow - recMinHt)/(recMaxHt - recMinHt));
                        colNow = c(colNowIndex,:); 
                    else
                        colNowIndex = 1; 
                        colNow = [25 25 25]/256; 
                    end
                else
                    xposNow = bulkCoords(ii,1) + recShiftx;
                    yposNow = bulkCoords(ii,2) + recShifty;
                    zposNow = bulkCoords(ii,3);
                    colNow = [25 25 25]/256; 
                end 
                % Initial positions of shown atoms avoids atoms popping in
                % and out of video 
                xposIn = bulkCoords(ii,1) + recShiftx;
                yposIn = bulkCoords(ii,2) + recShifty;
                zposIn = bulkCoords(ii,3);
                
                if (xposIn >= xmin-0.1) && (xposIn <= xmax+0.1) && ...
                        (yposIn >= ymin) && (yposIn <= ymax) && ...
                        (zposIn >= zmin) && (zposIn <= zmax)
                    surf(r*x+xposNow,r*y+yposNow,r*z+zposNow,...
                        'FaceAlpha',1,'FaceColor', ... 
                        colNow,'LineStyle',...
                        'none','FaceLighting','gouraud',...
                        'DiffuseStrength',1,'AmbientStrength',0.56);
                end
                
            end 
            % Fade in remaining surface atoms not picked by bulk atoms 
            for ii = 148:151 
                xposNow = si7x7(ii,1) + recShiftx;
                yposNow = si7x7(ii,2) + recShifty;
                zposNow = si7x7(ii,3);
                colNowIndex = ... % Index of colormap of surface atom
                    1 + round(255*(zposNow - recMinHt)/(recMaxHt - recMinHt));
                colNow = c(colNowIndex,:);
                if (xposNow >= xmin) && (xposNow <= xmax) && ...
                        (yposNow >= ymin) && (yposNow <= ymax) && ...
                        (zposNow >= zmin) && (zposNow <= zmax)
                    
                    surf(r*x+xposNow,r*y+yposNow,r*z+zposNow,...
                        'FaceAlpha',tt,'FaceColor', ...
                        colNow,'LineStyle',...
                        'none','FaceLighting','gouraud',...
                        'DiffuseStrength',1,'AmbientStrength',0.56);
                end
            end
        end
    end
            
    lp = [0.5 0 0.5];
    light('Position',lp,'Style','infinite');
    set(gca, 'Projection','perspective');
    view(41,50)
    axis off
    axis equal
    xlim([xmin-3.5*r xmax+3.5*r]);
    ylim([ymin-3.2*r ymax+3.2*r]);
    zlim([zmin-r zmax+r]);
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
    toc % Monitor speed 
end 
pauseMov(vid,77); % Pause for ca. 2.5 seconds to absorb the beauty 
% Output the movie as an mpg file
close(vid) 
