Animation of the formation of the Si(111) 7x7 surface reconstruction
Matlab code 

% 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);

 

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

    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 = blackorange(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 = blackorange(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

end 

pauseMov(vid,77);

%Output the movie as an mpg file

close(vid)