Au nanocrystals :: Calculation of diffraction pattern

• What are we going to calculate?
  • positions and intensities of diffractions of Au nanocrystals = structure Au/fcc
  • as we calculate powder diffraction pattern $\Rightarrow$ positions $\sim$ distances of diffraction rings from the origin
• Approximations:
  • kinematic diffraction theory (for small nanocrystals, multiple scattering can be neglected even in SAED)
  • atom scattering factors from XRD (i.e. we calculate XRD pattern, but for small nanocrystals it is very similar to SAED)

0. Initialization

import numpy as np
import matplotlib.pyplot as plt

from math import sqrt,pi,sin,cos,exp,asin

# Define parameters for plotting
# (global setting valid for the whole notebook
# (optimized parameters for export of graphs to PPTX
plt.rcParams['figure.figsize'] = (6,4)
plt.rcParams['axes.titlesize'] = 14
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['font.size'] = 13

1. Input parameters

• This section contains all important input parameters, defined by user:
  • Crystal structure: unit cell dimensions + atom positions
  • Basic information about diffraction experiment: wavelenght and maximum diffraction angle
• Crystal structure should not be changed
  • if we changed atoms, we would also have to change atomic scattering factors below
  • this is possible, but impractical: this program is a simple demo, not a real diffraction software
• Wavelength and maximumm diffraction angle can be adjusted for XRD or SAED
  • both parameters (wavelength, max_theta) are connected (and should be changed reasonably!)
  • even if we change these parameters, the final diffraction pattern $I = f(q)$ will be the same
  • this is due to the fact that diffraction vector $q = 4\pi \sin\theta/\lambda$ is independent of $\lambda$

# Unit cell of Au
# (a = b = c = 4.08A, alpha = beta = gamma = 90deg
a = 4.08

# Atom positions in Au unit cell
# (given in fractional coordinates
# (just four atoms per unit cell; the other atoms genrated by translations
Au1 = [0.0, 0.0, 0.0]
Au2 = [0.5, 0.5, 0.0]
Au3 = [0.5, 0.0, 0.5]
Au4 = [0.0, 0.5, 0.5]

# Wavelenght during diffraction experiment in Angstroms
# (for XRD, typical wavelenghts are CuKa = 1.54A, MoKa = 0.71A
# (for ED, typical wavelenghts range from ~0.04A (for 100kV) and ~0.02A (for 300kV)
# (technical note: lambda = reserved Python keyword => lamda conventional replacement
lamda = 0.71

# Typical maximal diffraction angle
# (for XRD: maximum diffraction angle = 90deg
# (for SAED: maximum diffractin angles are usualy below 3deg
theta_max = 90

2. Atomic scattering factor f(Au)

• Atomic scattering factor $f$ describes intensity of radiation scattered by one atom
• Atomic scattering factor is a decreasing function of scattering/diffraction angle of $\sin\theta/\lambda$
• Atomic scattering factor is usually calculated as: $f(\sin\theta/\lambda) = \sum_{i=1}^{4} a_i \exp(-b_i \sin\theta/\lambda) + c$
• The constants for calculation ($a_i, b_i, c$) are tabulated in International Tables for Crystallography

# Constants for calculation of fAu
Au_scattering_coefficients = {
    'A': [16.882, 18.591, 25.558,  5.860],
    'B': [0.461,   8.622,  1.483, 36.396],
    'C': 12.066 }

# Function for calculation of fAu
# (the function employs constants defined above
# (x = sin(theta)/lambda, theta = diffraction angle, lambda = wavelength
def fAu(x):
    # initialization: set fAu = 0 (separately for np.array or single number)
    if type(x)==np.array: fAu = np.zeros(len(x))
    else: fAu = 0
    # read the coefficients (into one-letter-variables for more readable calculation)
    (A,B,C) = A,B,C = Au_scattering_coefficients.values()
    # calculate fAu, using tabulated coefficients, which we read above
    for i in range(4): fAu = fAu + A[i] * np.exp(-B[i] * x**2)
    fAu = fAu + C 
    return fAu

Check: atomic scattering factor should be a decreasing function of sin(theta)/lambda

# (a) Calculate typical maximum values of sin(theta)/lambda for diffraction experiment
# (constants lambda and theta_max are given above, in the input section
sinTL_max = sin(theta_max)/lamda
print('%.2f' % sinTL_max)

1.26

# (b) Prepare X-variable = range of values sin(theta)/lambda
# (we will use even higher values of sin(theta)/lambda - to include exceptional cases
X = np.linspace(0,2,201, endpoint=True)

# (c) Plot atomic scattering factor a function of sin(theta)/lambda
plt.plot(X,fAu(X))
plt.title('Atomic scattering factor f(Au)')
plt.xlabel('sin(theta)/lambda [1/A]')
plt.ylabel('Scattering intensity')
plt.ylim(0,80)
plt.grid()
plt.show()

3. Structure factor: F(hkl) = A(hkl) + iB(hkl)

• Structure factor $F(hkl)$ describes scattering by one unit cell
• For crystalline samples, it holds: $I(hkl) \propto |F(hkl)|^2$, where...
  • $I(hkl)$ is the intensity of diffracted radiation for given interplanar distance $d(hkl) \sim \theta \sim q$
  • $F(hkl) = A(hkl) + iB(hkl)$ is a (complex) structure factor, which depends on atom positions in the unit cell
    • $F(hkl) = \sum_j^N f_j \exp[2 \pi i(hx_j + ky_j + lz_j)]$
    • $F(hkl) = \sum_j^N f_j \cos[2 \pi(hx_j + ky_j + lz_j)] + i \sum_j^N f_j \sin[2 \pi(hx_j + ky_j + lz_j)]$
    • where $f_j$ and ($x_j$, $y_j$, $z_j$) are atomic scattering factor and positions of $j$-th atom in the unit cell, respectively
    • and $N$ is number of atoms in the unit cell (and $i$ is imaginary unit)
• Brief practical conclusions:
  • once we calculate $F(hkl)$, we can calculate $I(hkl) = f(q)$ which is the diffraction pattern that we want
  • derivation of $F(hkl)$: EM lectures and crystallography textbooks
  • calculation of $F(hkl)$: demonstrated below

# Supplementary functions
# (these functions are needed in many following calculations
def d(h,k,l)    : return sqrt(a**2 / ( h**2 + k**2 + l**2 ))
def theta(h,k,l): return asin( lamda / (2*d(h,k,l)) ) 
def sinTL(h,k,l): return 1 / (2*d(h,k,l))

# Ahkl = real part of structure factor Fhkl
def Ahkl(h,k,l):
    Ahkl = \
        fAu(sinTL(h,k,l)) * cos(2 * pi * (h * Au1[0] + k * Au1[1] + l * Au1[2])) + \
        fAu(sinTL(h,k,l)) * cos(2 * pi * (h * Au2[0] + k * Au2[1] + l * Au2[2])) + \
        fAu(sinTL(h,k,l)) * cos(2 * pi * (h * Au3[0] + k * Au3[1] + l * Au3[2])) + \
        fAu(sinTL(h,k,l)) * cos(2 * pi * (h * Au4[0] + k * Au4[1] + l * Au4[2]))
    return(Ahkl)

# Bhkl = imaginary part of structure factor Fhkl
def Bhkl(h,k,l):
    Bhkl = \
        fAu(sinTL(h,k,l)) * sin(2 * pi * (h * Au1[0] + k * Au1[1] + l * Au1[2])) + \
        fAu(sinTL(h,k,l)) * sin(2 * pi * (h * Au2[0] + k * Au2[1] + l * Au2[2])) + \
        fAu(sinTL(h,k,l)) * sin(2 * pi * (h * Au3[0] + k * Au3[1] + l * Au3[2])) + \
        fAu(sinTL(h,k,l)) * sin(2 * pi * (h * Au4[0] + k * Au4[1] + l * Au4[2]))
    return(Bhkl)

# ahsFhkl = absolute value of structure factor Fhkl
def absFhkl(h,k,l):
    Fhkl = sqrt(Ahkl(h,k,l)**2 + Bhkl(h,k,l)**2)
    return(Fhkl)

Check: Fhkl values should be non-zero only if h,k,l are all odd or all even

• note: zero is rearded as even number
• reason: symmetry of Au crystal => fcc structure => F unit cell
• this is called systematic extinctions; justification: crystallography textbooks

# Max = maximum value of diffraction indexes h,k,l
# Eps = minimal value, below which absFhkl is regarded as zero
Max = 3
Eps = 1e-5

# Calculate non-zero Fhkl values
# (...check: non-zero values only if h,k,l are all odd or all even
# (...note: in cubic system: h00,0k0,00l are symmetrically equivalent
for h in range(0,Max+1):
    for k in range(h,Max+1):
        for l in range(k,Max+1):
            # Skip F(0,0,0); due to division by zero
            if not(h==0 and k==0 and l==0):
                # Skip F(h,k,l) = 0; systematic extinctions
                absFhkl_value = absFhkl(h,k,l)
                if not(absFhkl_value < Eps): 
                    print('F(%d,%d,%d) = %d' % (h,k,l,round(absFhkl_value)))

F(0,0,2) = 254
F(0,2,2) = 224
F(1,1,1) = 265
F(1,1,3) = 209
F(1,3,3) = 181
F(2,2,2) = 204
F(3,3,3) = 163

4. Corrections to F(hkl): Multiplicity, Temperature factor, LP factor

• In general, intensity of radiation diffracted by a crystal: $I(hkl) \propto |F(hkl)|^2$ as discussed in previous section
• For real-life calculations, we need several corrections: $I(hkl) \propto |F(hkl)|^2 \cdot T(hkl)^2 \cdot M(hkl) \cdot \mathrm{LP}(hkl)$
  • $T(hkl)$ = temperature factor, thermal motion of atoms (decreases value of structure factor)
  • $M(hkl)$ = multiplicity, correction for symmetrically equivalent reflections (their positions overlap)
  • $\mathrm{LP}(hkl)$ = Lorentz + Polarization factor (decrease in intensity due to exp.geometry and polarization)
• Technical notes:
  • $T(hkl)$ is correction of $F(hkl)$ $\Rightarrow$ we must square it together with Fhkl (see below)
  • $M(hkl)$ and $\mathrm{LP}(hkl)$ are corrections of $I(hkl)$ $\Rightarrow$ their values are not squared with Fhkl

4.1. Calculation of T(hkl) = temperature factor

• Simple isotropic temperature factor (from diffraction theory):
  • $T(hkl) = T = \exp(-B_{iso} \sin\theta^2/\lambda^2)$
• For atomic scattering factor:
  • $f = f_0 T = f_0 \cdot \exp(-B_{iso} \sin\theta^2/\lambda^2)$
• For structure factor:
  • $F(hkl) = \sum f_j T_j exp[2\pi i(hx_j + ky_j + lz_j)]$ $= F_0(hkl) \exp(-B_{iso} \sin\theta^2/\lambda^2)$
• For intensity:
  • $I(hkl) \propto |F(hkl)|^2 = |F(hkl)^2|$ $= |F_0(hkl)^2| \cdot \exp[-B_{iso} \sin\theta^2/\lambda^2]^2 = |F_0(hkl)^2| \cdot T(hkl)^2$
• Mathematical and technical notes:
  • $T(hkl)^2 = [\exp(-B_{iso} \sin\theta^2/\lambda^2)]^2 = \exp(-2B_{iso} \sin\theta^2/\lambda^2)$ ...because exp(x)*exp(x) = exp(2*x)
  • $B_{iso} = 8\pi^2 U_{iso}$ ...typical values of U(iso) = 0.005-0.02 A$^2$ $\Rightarrow$ typical value of B(iso) $\approx$ 1.2 A$^2$

# Thkl = temperature factor
# (temperature factor says, how the intensity is decreased due to thermal motion
# (the higher the diffraction angle, the higher the decrease of Fhkl due to Thkl

# Biso = constant describing isotropic thermal motion of atoms
# (atomic scattering factor: f = f(0) * exp[-Biso*sinTL**2]
# (Biso = 8*pi**2*Uiso; Uiso from various CIFs = 0.005-0.02 => Biso ~ 1.2
Biso = 0.8

# Thkl = average isotropic thermal factor calculated with above-defined Biso
# (rough approximation: the same constant for all atoms, isotropic vibrations...
# (we use the original temperature factor; when calculating Ihkl, this must be squared
def Thkl(h,k,l):
    return( exp(-Biso * sinTL(h,k,l)**2) )

Check: Thkl should decrease with increasing diffraction angle

# Numerical verification: calculation of T(hkl) for several reflections

print('%3s %8s %10s %10s' % ('hkl', 'Th[deg]','T(hkl)','T(hkl)**2'))
for h in range(1,10):
    print('%d00 %8.2f %10.3f %10.3f' % (h, theta(h,0,0)*180/pi, Thkl(h,0,0), Thkl(h,0,0)**2))

hkl  Th[deg]     T(hkl)  T(hkl)**2
100     4.99      0.988      0.976
200    10.02      0.953      0.908
300    15.13      0.898      0.806
400    20.37      0.825      0.681
500    25.79      0.741      0.548
600    31.47      0.649      0.421
700    37.52      0.555      0.308
800    44.11      0.464      0.215
900    51.54      0.378      0.143

# Graphical verification: calculation of T[sin(theta)/lambda] for the whole range

# Prepare X,Y-data for plotting
# X = sinTL = sin(theta)/lambda; range (0,2) as justified in section 2 above
# Y = T[sin(theta)/lambda] = like Thkl function defined above, but continuous
X = np.linspace(0,2,201, endpoint=True)
Y = np.exp(-Biso * X**2)
Y = Y/np.max(Y)

# Plot T[sin] a function of sin(theta)/lambda
plt.plot(X,Y)
plt.title('Temperature factor T = f[sin(theta)/lambda)]')
plt.xlabel('sin(theta)/lambda [1/A]')
plt.ylabel('Normalized temp.factor')
plt.grid()
plt.show()

4.2. Calculation of M(hkl) = multiplicity

# Mhkl = multiplicity
# (multiplicity says, how many symmetrically equivalent reflections overlap
# (example: (1,0,0) = (-1,0,0) = (0,1,0) = (0,-1,0) = (0,0,1) = (0,0,-1) => M(1,0,0) = 6
# (...all reflections above have the same d(hkl) and so we have overlap of 6 reflections
# (note: here we hard-code multiplicity for cubic system

def Mhkl(h,k,l):
    o = sorted([h,k,l])
    # ...three zero indexes => M=1
    if o[0] == o[1] == o[2] == 0: return(1)
    # ...two zero indexes => M=6
    if o[0] == o[1] == 0: return(6)
    # ...one zero index => M=12/24 (the remaining two the same/different)
    if o[0] == 0:
        if o[1] == o[2] : return(12)
        else            : return(24)
    # ...no zero index, three are the same => M=8
    if o[0] == o[1] == o[2]        : return(8)
    # ...no zero index, two are the same => M=24
    if o[0] == o[1] or o[1] == o[2]: return(24)
    # ...no zero index, all different  => M=48
    return(48)

Check: Mhkl - reflections with more different (h,k,l) have higher multiplicity $\Rightarrow$ intensity

print('Multiplicity of (00l),(0ll),(0kl):', Mhkl(0,0,1), Mhkl(0,1,1), Mhkl(0,1,2))
print('Multiplicity of (hhh),(hkk),(hkl):', Mhkl(1,1,1), Mhkl(1,1,2), Mhkl(1,2,3))

Multiplicity of (00l),(0ll),(0kl): 6 12 24
Multiplicity of (hhh),(hkk),(hkl): 8 24 48

4.3. Calculation of LP(hkl) = combined Lorentz and Polarization factor

# LP = LP-factor = Lorentz-Polarization factor
# (LP-factor = decrease in intensity of Fhkl due to polarization and experiment geometry
# (LP-factor is specific for each experiment geometry:
# (- classics/Weissenberg   : LP = (1 + cos(2*Th)**2) / (2 * sin(2*Th)
# (- standard/Bragg-Brentano: LP = (1 + cos(2*Th)**2) / (2 * sin(Th)**2 * cos(Th))
# (- further complexities occur for polarized light or different experiment geometries...

def LP(h,k,l):
    # For clarity/efficiency of the following formula, we pre-calculate Th = theta
    Th = theta(h,k,l)
    # LP for standard Bragg-Brentano geometry, non-polarized/non-monochromated radiation
    LP = (1 + cos(2*Th)**2) / (2 * sin(Th)**2 * cos(Th))
    return(LP)

Check: LPhkl should decrease with increasing diffraction angle

print('%3s %8s %8s' % ('hkl', 'Th[deg]','LP(hkl)'))
for h in range(1,10):
    print('%d00 %8.2f %8.3f' % (h, theta(h,0,0)*180/pi, LP(h,0,0)))

hkl  Th[deg]  LP(hkl)
100     4.99  130.599
200    10.02   31.564
300    15.13   13.273
400    20.37    6.931
500    25.79    4.067
600    31.47    2.596
700    37.52    1.813
800    44.11    1.439
900    51.54    1.378

5. Final intensities: I(hkl) = calculated diffraction pattern

• This is final calculation; we are going to calculate two arrays:
  • difraction intensities: $I(hkl) \propto |F(hkl)|^2 \cdot T(hkl)^2 \cdot M(hkl) \cdot \mathrm{LP}(hkl)$
  • diffraction vectors: $q = 4\pi \sin\theta/\lambda$ (which correspond to $I(hkl)$ and depend on $d(hkl)$, of course)
• Final plot = calculated diffraction pattern:
  • diffraction intensity as a function of diffraction vector: $I(hkl) = f(q)$

# Max = maximum value of diffraction indexes h,k,l
# Eps = minimal value, below which absFhkl is regarded as zero
Max = 6
Eps = 1e-5

# Calculate final values of diffraction intensities Ihkl
# ...arrays, into which the results will be saved
diff_vectors = []
diff_intensities = []
# ...calculate intensities, save only non-zero values
for h in range(0,Max+1):
    for k in range(h,Max+1):
        for l in range(k,Max+1):
            # Skip F(0,0,0); due to division by zero
            if not(h==0 and k==0 and l==0):
                # Skip F(h,k,l) = 0; systematic extinctions
                absFhkl_value = absFhkl(h,k,l)
                if absFhkl_value > Eps:
                    diff_vectors.append(4*pi*sinTL(h,k,l))
                    diff_intensities.append(
                        absFhkl_value**2 * Thkl(h,k,l)**2 * Mhkl(h,k,l) * LP(h,k,l))         

# Normalize calculated intensities
diff_intensities = np.array(diff_intensities)
max_intensity    = np.max(diff_intensities)
diff_intensities = diff_intensities/max_intensity

# Plot calculated intensities
plt.vlines(diff_vectors,0,diff_intensities, color='red')
plt.title('Calculated diffraction pattern of Au nanocrystals')
plt.xlabel('Diffraction vector q [1/A]')
plt.ylabel('Normalized intensity I(hkl)')
plt.xlim(2,12)
plt.ylim(0,1.05)
plt.show()

Check: Comparison with the diffraction pattern calculated with PowderCell

• PowderCell is a standard program for calculation of powder diffraction patterns.
• Powder diffractogram of Au calculated by PowderCell is shown in EM lectures (lecture EM2/TEM).
• Note that the positions and (relative) intensities in our calculation and PowderCell calculation are similar.
• Direct comparison of calculation and experiment:
  M:\MIREK.PRE\1_PRFUK\EM-AFM\Z_WWW.IMC\SUPPL\x_saed_au_srovnani.pptx.png

6. Conclusions

• We have calculated powder diffraction pattern for Au nanocrystals from the very beginning.
  • We have used several approximations (kinematic diffraction theory, f(Au) from XRD, simple corrections...)
  • Nevertheless, our output is in very good agreement with the output from standard diffraction software!
• What we have learnt from this:
  • Kinematic diffraction theory may be difficult, but it is not unusable.
  • Python/Jupyter is a powerful software, suitable for both simple and complex calculations.
• More practical conclusions:
  • For small crystals, the difference between XRD and SAED is not critical.
  • TEM/SAED can be used for identification of known structures (the same diffraction = the same structures).