SAD intensity differences of Friedel pairs
Matlab code
% Program to plot out in 3D reciprocal space the structure factors of
% benzylpenicillin across the sulfur edge and how the Friedel pairs differ
clear all; close all;
vid = VideoWriter('SADintensities.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',10);
dhklmin = 0.8; % In Angstroms
LL = 2.1*pi/dhklmin;
lp = [0.4 -0.7 0.7]; % Light projection angle
% Limits of real-space volume in Angstroms
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
axis equal
axis off
hold on
rgb = customcolormap([0 0.5 1], {'#0000ff','#ffffff','#ff0000'});
colormap(rgb);
set(gca,'View',[160,5]);
set(gca, 'Projection','perspective');
hold on
light('Position',lp,'Style','infinite');
% Lattice constants of a fictitious orthorhombic unit cell
a = 14; % u.c. size in x-direction in Angstroms
b = 8; % u.c. size in y-direction in Angstroms
c = 4; % u.c. size in z-direction in Angstroms
% (10 x 7) matrix: (# electrons + 9 form factor parameters) x 7 atom types
ffData = importdata('formFactorParameters_CHFClNOS.dat');
% Atomic coordinates = (5 x 23) matrix: x, y, z, Z, radius. Ignores
% hydrogen atoms
penData = importdata('penicillinCoordinates.dat');
% Coordinates of all atoms in terms of unit cell size
xa = penData(:,1)/a;
ya = penData(:,2)/b;
za = penData(:,3)/c;
zaSize = size(za);
phi = zeros(zaSize(1,1),1);
% Import f1 and f2 data for sulfur - this is not evenly spaced, so
% interpolation needed
sf12 = importdata('Sf1f2.dat');
xx = sf12(:,1);
f1 = sf12(:,2);
f2 = sf12(:,3);
xintp = 1000:0.25:8000; % Photon energy range in 4 eV steps
Sf1Intp = interp1(xx,f1,xintp,'pchip'); % Interpolated f1 data
Sf2Intp = interp1(xx,f2,xintp,'pchip'); % Interpolated f2 data
SfPrime = Sf1Intp - max(Sf1Intp); % f' for sulfur - negative values!!
hmax = floor(a/dhklmin); % largest Miller index h for given resolution dhklmin
kmax = floor(b/dhklmin); % ---"--- k
lmax = floor(c/dhklmin); % ---"--- l
Qmax = ceil(2*pi/dhklmin); % Largest Q-vector required for desired resolution
% Total number of electrons = (2 x 7) + (4 x 8) + (16 x 6) + (1 x 16) = 158
% Write direct (0,0) beam to file. Set dhkl to a large number so this comes
% first
fQzero = 158;
% Scan energy with different speeds
pt1 = (1:40:4801); % 10-eV steps
pt2 = (4840:12:5596); % 3-eV steps
pt3 = (5600:6160); % 0.25-eV steps across edge
pt4 = (6172:12:7000); % 3-eV steps
pt5 = (7121:120:28001); % 30-eV steps
pos1 = [0.02 0.1 0.64 0.9];
pos2 = [0.66 0.3 0.32 0.4];
for hv = [pt1 pt2 pt3 pt4 pt5]
subplot('Position',pos1); % Main plot
newplot
hold on
htext = text(1.07*LL,0,0,'h','FontSize',20);
ktext = text(0,1.07*LL,0,'k','FontSize',20);
ltext = text(0.0061*LL,0,1.07*LL,'l','FontSize',20);
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
axis equal
axis off
axis tight
rgb = colorbar;
set(rgb,'position',[.07 .23 .01 .2]); % Size and position of colorbar
rgb.Color = [0 0 0];
caxis([-2 2]);
set(rgb,'Ticks',[-2,-1,0,1,2],...
'TickLabels',{'-2','-1','0','1','2'},'FontSize',16);
set(gca,'View',[160,5]);
set(gca, 'Projection','perspective');
hold on
line([0,0], [0,0], [-LL,LL], 'LineWidth', 2, 'Color', 'k');
line([0,0], [-LL,LL], [0,0], 'LineWidth', 2, 'Color', 'k');
line([-LL,LL], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
sadInt = zeros((2*hmax+1)*(2*kmax+1)*(2*lmax+1),5);
bP = 0;
for h = -hmax:hmax
for k = -kmax:kmax
for l = -lmax:lmax
bP = bP+1;
% dhkl = 1/Y. Determine here Y
Y = ((h/a)^2 + (k/b)^2 + (l/c)^2)^0.5;
% Not the Q=0 forward scattered reflection or Q > Qmax
if (h^2 + k^2 + l^2 ~= 0) && (2*pi*Y <= Qmax)
dhkl = 1/Y; % Lattice spacing of (kl)-planes
Q = 2*pi/dhkl; % Associated Q-value in reciprocal Angstroms
% Now calculate form factors for each atom type according
% to standard of four Gaussians plus a constant
ffC = ffData(1,2)*exp(-ffData(1,3)*(Q/(4*pi))^2) + ...
ffData(1,4)*exp(-ffData(1,5)*(Q/(4*pi))^2) + ...
ffData(1,6)*exp(-ffData(1,7)*(Q/(4*pi))^2) + ...
ffData(1,8)*exp(-ffData(1,9)*(Q/(4*pi))^2) + ffData(1,10);
ffH = ffData(2,2)*exp(-ffData(2,3)*(Q/(4*pi))^2) + ...
ffData(2,4)*exp(-ffData(2,5)*(Q/(4*pi))^2) + ...
ffData(2,6)*exp(-ffData(2,7)*(Q/(4*pi))^2) + ...
ffData(2,8)*exp(-ffData(2,9)*(Q/(4*pi))^2) + ffData(2,10);
ffF = ffData(3,2)*exp(-ffData(3,3)*(Q/(4*pi))^2) + ...
ffData(3,4)*exp(-ffData(3,5)*(Q/(4*pi))^2) + ...
ffData(3,6)*exp(-ffData(3,7)*(Q/(4*pi))^2) + ...
ffData(3,8)*exp(-ffData(3,9)*(Q/(4*pi))^2) + ffData(3,10);
ffCl = ffData(4,2)*exp(-ffData(4,3)*(Q/(4*pi))^2) + ...
ffData(4,4)*exp(-ffData(4,5)*(Q/(4*pi))^2) + ...
ffData(4,6)*exp(-ffData(4,7)*(Q/(4*pi))^2) + ...
ffData(4,8)*exp(-ffData(4,9)*(Q/(4*pi))^2) + ffData(4,10);
ffN = ffData(5,2)*exp(-ffData(5,3)*(Q/(4*pi))^2) + ...
ffData(5,4)*exp(-ffData(5,5)*(Q/(4*pi))^2) + ...
ffData(5,6)*exp(-ffData(5,7)*(Q/(4*pi))^2) + ...
ffData(5,8)*exp(-ffData(5,9)*(Q/(4*pi))^2) + ffData(5,10);
ffO = ffData(6,2)*exp(-ffData(6,3)*(Q/(4*pi))^2) + ...
ffData(6,4)*exp(-ffData(6,5)*(Q/(4*pi))^2) + ...
ffData(6,6)*exp(-ffData(6,7)*(Q/(4*pi))^2) + ...
ffData(6,8)*exp(-ffData(6,9)*(Q/(4*pi))^2) + ffData(6,10);
ffS = ffData(7,2)*exp(-ffData(7,3)*(Q/(4*pi))^2) + ...
ffData(7,4)*exp(-ffData(7,5)*(Q/(4*pi))^2) + ...
ffData(7,6)*exp(-ffData(7,7)*(Q/(4*pi))^2) + ...
ffData(7,8)*exp(-ffData(7,9)*(Q/(4*pi))^2) + ffData(7,10);
ffScorr = ffS + SfPrime(hv); % f1 for present Q-value
% Now determine the phase of each atom
for atom = 1:zaSize(1,1)
phi(atom,1) = 2*pi*(h*xa(atom,1) + k*ya(atom,1) + l*za(atom,1));
end
% Now determine the real and imaginary components of the
% contributions of each atom type on the Argand diagram.
% First initialize components to zero
Creal = 0.0;
Cimag = 0.0;
Nreal = 0.0;
Nimag = 0.0;
Oreal = 0.0;
Oimag = 0.0;
Sreal = 0.0;
Simag = 0.0;
Srealbar = 0.0;
Simagbar = 0.0;
for Z = 1:zaSize(1,1)
if (penData(Z,4) == 6) % Carbon
Creal = Creal + ffC*cos(phi(Z,1));
Cimag = Cimag + ffC*sin(phi(Z,1));
elseif (penData(Z,4) == 7) % Nitrogen
Nreal = Nreal + ffN*cos(phi(Z,1));
Nimag = Nimag + ffN*sin(phi(Z,1));
elseif (penData(Z,4) == 8) % Oxygen
Oreal = Oreal + ffO*cos(phi(Z,1));
Oimag = Oimag + ffO*sin(phi(Z,1));
elseif (penData(Z,4) == 16) % Sulphur
Sf2 = Sf2Intp(hv);
Sreal = Sreal + ffScorr*cos(phi(Z,1))+Sf2*sin(phi(Z,1));
Simag = Simag + ffS*sin(phi(Z,1))+Sf2*cos(phi(Z,1));
Srealbar = Srealbar + ffScorr*cos(phi(Z,1))-Sf2*sin(phi(Z,1));
Simagbar = Simagbar + ffS*sin(phi(Z,1))-Sf2*cos(phi(Z,1));
end
end
% Now determine the real and imaginary components of the total
% structure factor plus its phase
FtotReal = Creal+Nreal+Oreal+Sreal;
FtotImag = Cimag+Nimag+Oimag+Simag;
FtotRealbar = Creal+Nreal+Oreal+Srealbar;
FtotImagbar = Cimag+Nimag+Oimag+Simagbar;
Ftot = (FtotReal^2 + FtotImag^2)^0.5; % Pythagoras
Ftotbar = (FtotRealbar^2 + FtotImagbar^2)^0.5; % Pythagoras
% phiTot = atan2(FtotImag,FtotReal);
sadInt(bP,1) = h*2*pi/a;
sadInt(bP,2) = k*2*pi/b;
sadInt(bP,3) = l*2*pi/c;
sadInt(bP,4) = Ftot^2;
sadInt(bP,5) = 2*(Ftot^2 - Ftotbar^2)/(Ftot^2 + Ftotbar^2);
else % Not (000) peak or Q > Qmax if statement
sadInt(bP,1) = NaN;
sadInt(bP,2) = NaN;
sadInt(bP,3) = NaN;
sadInt(bP,4) = NaN;
sadInt(bP,5) = NaN;
end
end % l-loop
end % k-loop
end % h-loop
scatter3(sadInt(:,1),sadInt(:,2),sadInt(:,3),2*sadInt(:,4).^0.5,sadInt(:,5),...
'MarkerFaceColor','Flat','MarkerEdgeColor','k');
hold off
subplot('Position',pos2); % f1f2vEnergy plot
hold off
newplot
plot(xintp,Sf1Intp, 'LineWidth', 2.0,'Color','blue')
hold on
plot(xintp,Sf2Intp, 'LineWidth', 2.0,'Color','red')
plot(xintp(hv),Sf1Intp(hv),'bo','MarkerFaceColor','blue')
plot(xintp(hv),Sf2Intp(hv),'ro','MarkerFaceColor','red')
set(gca,'linewidth',2);
xlabel('Photon energy [eV]');
ylabel('f_1, f_2');
set(gca,'FontName','Helvetica','fontsize',18);
set(gca,'TickLength',[0.016, 2]);
hText = text(7000, 15.1, 'f_1', 'FontSize',18, 'color', 'blue');
hText = text(7000, 2.5, 'f_2', 'FontSize',18, 'color', 'red');
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end % End scanning energy across S K-edge
close(vid);