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