Unit cell to crystal
Matlab code
% Movie schematic from a single unit cell (of either NaCl, diamond, GaAs, or
% benzyl penicillin) to a crystal
clear all; close all;
prompt = 'Which unit cell?: diamond (d/D), GaAs (g/G), NaCl (s/S), or penicillin (p/P)? ';
str = input(prompt,'s');
if (str == 'P') || (str == 'p')
ucData = importdata('penicillinCoordinates.dat'); % (23 x 5) matrix
C = [1 1 0; 1 0 0; 1 0 0; 1 0 0; 1 0 0; 0 0 1; 0 0 1; 0.4 0.4 0.4; ...
0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ...
0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ...
0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4];
noAtoms = 23;
ucxStep = 1;
ucyStep = 0.5;
uczStep = 0.5;
ucMax = max(ucData(:,1));
ucMin = min(ucData(:,1));
ucSize = max(ucData(:,5)) + (ucMax - ucMin);
youChoose = 0.74; % How much smaller than u.c. the molecule is
sphereRes = 16;
azStep = 0.4; % Causes camera to swing azimuthally around the crystal as it
% is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'S') || (str == 's')
ucData = importdata('NaClCoordinates.dat'); % (8 x 5) matrix
% Fuschia for Na, green for Cl
C = [0.631 0.376 0.918; 0.631 0.376 0.918; 0.631 0.376 0.918; ...
0.631 0.376 0.918; 0.443 0.925 0.302; 0.443 0.925 0.302; ...
0.443 0.925 0.302; 0.443 0.925 0.302];
noAtoms = 8;
ucxStep = 1;
ucyStep = 1;
uczStep = 1;
ucSize = 5.64; % u.c. lattice constant of NaCl
youChoose = 1.00; % How much smaller than u.c. the molecule is
sphereRes = 35;
azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
% is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'D') || (str == 'd')
ucData = importdata('diamondCoordinates.dat'); % (8 x 5) matrix
C = [0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; ...
0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25];
noAtoms = 8;
ucxStep = 1;
ucyStep = 1;
uczStep = 1;
ucSize = 3.567; % u.c. lattice constant of diamond
youChoose = 1.00; % How much smaller than u.c. the molecule is
sphereRes = 16;
azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
% is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'G') || (str == 'g')
ucData = importdata('GaAsCoordinates.dat'); % (8 x 5) matrix
% Brown for Ga, purple for As
C = [0.395 0.215 0.215; 0.395 0.215 0.215; 0.395 0.215 0.215; 0.395 0.215 0.215; ...
0.773 0.502 0.949; 0.773 0.502 0.949; 0.773 0.502 0.949; 0.773 0.502 0.949];
noAtoms = 8;
ucxStep = 1;
ucyStep = 1;
uczStep = 1;
ucSize = 4.066; % u.c. lattice constant of GaAs
youChoose = 1.00; % How much smaller than u.c. the molecule is
sphereRes = 16;
azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
% is being formed. Rotates by (10 noUC + 3)*azStep degrees.
end
vid = VideoWriter('uc2crystal.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);
noUC = 4;
ii = 0;
az = -125;
pol = 20;
% Ensure uc within 1 x 1 x 1 volume of plotted space
ucMax = max(max(ucData(:,1:3)));
ucMin = min(min(ucData(:,1:3)));
scaleFac = youChoose/ucSize; % Reduces length of molecule along x-axis to 1.0*youChoose
[x,y,z] = sphere(140);
% Read unit-cell data
for i = 1:noAtoms
ii = ii+1;
rad = ucData(i,5)*scaleFac;
a = ucData(i,1)*scaleFac;
b = ucData(i,2)*scaleFac;
c = ucData(i,3)*scaleFac;
atomCol = C(i,:);
% Create concatenated surface image of the unit cell
HkB(ii) = surf(-noUC+rad*x+a,-noUC+rad*y+b,noUC+rad*z+c,...
'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
'none','FaceLighting','flat',... %'gouraud',...
'DiffuseStrength',1,'Visible','off');
hold on
end
set(HkB,'Visible','on')
lp = [-1 -0 0.];
light('Position',lp,'Style','infinite');
% Pan out in Ken Burns style
zoomOutFrames = 50;
for kenB = 0:zoomOutFrames
hold off
% frame limits increasing from centering on molecule to it being at the
% corner of a crystal
xlo = -noUC + min(ucData(:,1))/ucSize - max(ucData(:,5))/ucSize;
xhi = -noUC + max(ucData(:,1))/ucSize + max(ucData(:,5))/ucSize + 2*noUC*kenB/zoomOutFrames;
ylo = -noUC + min(ucData(:,2))/ucSize - max(ucData(:,5))/ucSize;
yhi = -noUC + max(ucData(:,2))/ucSize + max(ucData(:,5))/ucSize + 2*noUC*kenB/zoomOutFrames;
zlo = noUC + min(ucData(:,3))/ucSize - max(ucData(:,5))/ucSize - 2*noUC*kenB/zoomOutFrames;
zhi = noUC + max(ucData(:,3))/ucSize + max(ucData(:,5))/ucSize;
axis square
axis equal
axis off
set(gca,'View',[az,pol]);
set(gca, 'Projection','perspective');
xlim([xlo xhi]);
ylim([ylo yhi]);
zlim([zlo zhi]);
% Store the frame
frame = getframe(gcf);
writeVideo(vid,frame);
end
[x,y,z] = sphere(sphereRes); % Reduce resoution of spheres now they're tiny
hold off
% Generate first line of unit cells
for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
ii = 0;
for i = 1:noAtoms % Loop to draw atoms of one unit cell
ii = ii+1;
rad = ucData(i,5)*scaleFac;
a = ucData(i,1)*scaleFac;
b = ucData(i,2)*scaleFac;
c = ucData(i,3)*scaleFac;
atomCol = C(i,:);
H1d(ii) = surf(no1d+rad*x+a,-noUC+rad*y+b,noUC+rad*z+c,...
'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
'none','FaceLighting','flat',... %'gouraud',...
'DiffuseStrength',1,'Visible','off');
hold on
end
set(H1d,'Visible','on')
axis square
axis equal
axis off
set(gca, 'Projection','perspective');
xlim([xlo xhi]);
ylim([ylo yhi]);
zlim([zlo zhi]);
if (no1d == -noUC)
light('Position',lp,'Style','infinite');
end
for jj = 1:6
%az = az+azStep; % commented out if you don't want the line to
%rotate
set(gca,'View',[az,pol]);
% Store the frame
frame = getframe(gcf);
writeVideo(vid,frame);
end
end % End of drawing one line
hold off
% Generate first plane of unit cells
for no2d = -noUC:ucyStep:noUC % Loop to replicate molecule over a single plane
ii = 0;
for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
for i = 1:noAtoms % Loop to draw atoms of one unit cell
ii = ii+1;
rad = ucData(i,5)*scaleFac;
a = ucData(i,1)*scaleFac;
b = ucData(i,2)*scaleFac;
c = ucData(i,3)*scaleFac;
atomCol = C(i,:);
H2d(ii) = surf(no1d+rad*x+a,no2d+rad*y+b,noUC+rad*z+c,...
'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
'none','FaceLighting','flat',... %'gouraud',...
'DiffuseStrength',1,'Visible','off');
hold on
end % End of drawing one molecule
if (no1d == -noUC) && (no2d == -noUC)
light('Position',lp,'Style','infinite');
end
end % End of drawing one line
set(H2d,'Visible','on')
axis square
axis equal
axis off
set(gca, 'Projection','perspective');
xlim([xlo xhi]);
ylim([ylo yhi]);
zlim([zlo zhi]);
for jj = 1:6
az = az+azStep;
set(gca,'View',[az,pol]);
% Store the frame
frame = getframe(gcf);
writeVideo(vid,frame);
end
end % End of drawing one plane
hold off
newplot
% Generate multiple planes of unit cells
for no3d = noUC:-uczStep:-noUC % Loop to replicate planes
newplot
ii = 0;
for no2d = -noUC:ucyStep:noUC % Loop to replicate molecule over a single plane
for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
for i = 1:noAtoms % Loop to draw atoms of one unit cell
ii = ii+1;
rad = ucData(i,5)*scaleFac;
a = ucData(i,1)*scaleFac;
b = ucData(i,2)*scaleFac;
c = ucData(i,3)*scaleFac;
atomCol = C(i,:);
H3d(ii) = surf(no1d+rad*x+a,no2d+rad*y+b,no3d+rad*z+c,...
'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
'none','FaceLighting','flat',... %'gouraud',...
'DiffuseStrength',1,'Visible','off');
hold on
end % End of drawing one molecule
if (no1d == -noUC) && (no2d == -noUC) && (no3d == noUC)
light('Position',lp,'Style','infinite');
end
end % End of drawing one line
end % End of drawing one plane
set(H3d,'Visible','on')
set(0, 'DefaultFigureRenderer', 'opengl');
axis square
axis equal
axis off
set(gca, 'Projection','perspective');
xlim([xlo xhi]);
ylim([ylo yhi]);
zlim([zlo zhi]);
for jj = 1:6
az = az+azStep;
set(gca,'View',[az,pol]);
% Store the frame
drawnow
frame = getframe(gcf);
writeVideo(vid,frame);
end
clear H3d
end % End of drawing multiple planes
% Output the movie as an mpg file
close(vid);