Chapter 2: Diffraction
Plotting Diffraction Pattern¶
part of
MSE672: Introduction to Transmission Electron Microscopy
Spring 2026
by Gerd Duscher
Microscopy Facilities
Institute of Advanced Materials & Manufacturing
Materials Science & Engineering
The University of Tennessee, Knoxville
Background and methods to analysis and quantification of data acquired with transmission electron microscopes.
import sys
import importlib.metadata
def test_package(package_name):
"""Test if package exists and returns version or -1"""
try:
version = importlib.metadata.version(package_name)
except importlib.metadata.PackageNotFoundError:
version = '-1'
return version
if test_package('pyTEMlib') < '0.2026.1.1':
print('installing pyTEMlib')
!{sys.executable} -m pip install pyTEMlib -q --upgrade
print('done')installing pyTEMlib
done
Load the plotting and figure packages¶
Import the python packages that we will use:
Beside the basic numerical (numpy) and plotting (pylab of matplotlib) libraries,
three dimensional plotting
and some libraries from the pyTEMlib
diffraction_tools
crystal_tools
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
import sys
if 'google.colab' in sys.modules:
from google.colab import output
output.enable_custom_widget_manager()
# additional package
import itertools
import scipy
# Import libraries from the pyTEMlib
import pyTEMlib
__notebook_version__ = '2026.01.08'
print('pyTEM version: ', pyTEMlib.__version__)
print('notebook version: ', __notebook_version__)pyTEM version: 0.2026.1.0
notebook version: 2026.01.08
Define Silicon crystal¶
the crystal structure is defined in an ase.Atoms format
atoms = pyTEMlib.crystal_tools.structure_by_name('Silicon')
print(atoms)
#Reciprocal Lattice
# We use the linear algebra package of numpy to invert the unit_cell "matrix"
reciprocal_unit_cell = atoms.cell.reciprocal() # np.linalg.inv(atoms.cell).T # transposed of inverted unit_cell
Lattice(symbols='Si8', pbc=True, cell=[5.43088, 5.43088, 5.43088])
All Possible and Allowed Reflections¶
see
Kinematic Scattering Geometry for details
Because the nuclei of the atoms in a material are positive, there is a small acceleration of the electron within the material.
We add this to the contribution to the magnitude of the incident wave-vector. The acceleration potential is called inner potential and is calculated from the scattering factors according to (Kirkland Eq 6.10:
which reduces for to:
which forms an uniform potential inside the material
# --------------- INPUT ------------------------
zone_hkl = np.array([1, 1, 0])
hkl_max = 35 # maximum allowed Miller index
Sg_max = 0.03 # 1/Ang maximum allowed excitation error
acceleration_voltage_V = 200.0 * 1000.0 #V
# -------------------------------------------
wave_length_A = pyTEMlib.diffraction_tools.get_wavelength(acceleration_voltage_V, unit='A')
U_0 = 0 # in (Ang)
# atom form factor of zero reflection angle is the inner potential in 1/A
for atom in atoms:
U_0 += pyTEMlib.diffraction_tools.form_factor(atom.symbol, 0)
scattering_factor_to_volts = (scipy.constants.h*1e10)**2 / (2 * np.pi * scipy.constants.m_e * scipy.constants.e) / atoms.cell.volume
# Conversion of inner potential to Volts
U_0 = U_0 * scattering_factor_to_volts
print('The inner potential is {0:.2f}V'.format(U_0))
incident_wave_vector_vacuum = 1/wave_length_A
K0_magnitude = incident_wave_vector = np.sqrt(1/wave_length_A**2 + U_0 )#1/Ang
cent = np.dot(zone_hkl,reciprocal_unit_cell)
cent = cent /np.linalg.norm(cent)* incident_wave_vector
# zone axis in global coordinate system
zone_vector = np.dot(zone_hkl,reciprocal_unit_cell)
h = np.linspace(-hkl_max,hkl_max, 2*hkl_max+1) # all evaluated single Miller Indices
hkl = np.array(list(itertools.product(h, h, h) )) # all evaluated Miller indices
g = np.dot(hkl, reciprocal_unit_cell) # all evaluated reciprocal lattice points
# Calculate excitation errors for all reciprocal lattice points
## Zuo and Spence, 'Adv TEM', 2017 -- Eq 3:14
S = (K0_magnitude**2-np.linalg.norm(g - cent, axis =1)**2)/(2*K0_magnitude)
print(cent)
# Determine reciprocal lattice points with excitation error less than the maximum allowed one: Sg_max
reflections = abs(S)< Sg_max
Sg = S[reflections]
g_hkl = g[reflections]
hkl = hkl[reflections]
print ('Of the {0} tested reciprocal lattice points, {1} have an excitation error less than {2:.2f} 1/nm'.format( len(g) , len(g_hkl), Sg_max*10))
# Calculate Structure Factors
structure_factors = []
base = atoms.positions
for j in range(len(g_hkl)):
F = 0
for atom in atoms:
f = pyTEMlib.diffraction_tools.form_factor(atom.symbol, np.linalg.norm(g_hkl[j]))
F += f * np.exp(-2*np.pi*1j*(g_hkl[j]*atom.position).sum())
structure_factors.append(F)
F = structure_factors = np.array(structure_factors)
allowed = np.absolute(structure_factors) > 0.000001
print('Of the {0} possible reflection {1} are allowed.'.format(hkl.shape[0], allowed.sum() ) )
# information of allowed reflections
Sg_allowed = Sg[allowed]
hkl_allowed = hkl[allowed][:]
g_allowed = g_hkl[allowed, :]
F_allowed = F[allowed]
intensities = np.absolute(F_allowed)**2
# forbidden reflections:
forbidden = np.logical_not(allowed)
g_forbidden = g_hkl[forbidden, :]
# Determine Laue Zone for reflections
zone_axis = np.dot(zone_hkl, atoms.cell.array)
Laue_Zone = np.floor(abs(np.dot(hkl_allowed, zone_hkl))) # Weiss zone law: works only for all crystal systems
ZOLZ = Laue_Zone ==0
HOLZ = Laue_Zone > 1
print ('Of those, there are {0} in ZOLZ and {1} in HOLZ'.format(ZOLZ.sum(),HOLZ.sum()))
intensities = np.absolute(F_allowed)**2
K0_magnitude
The inner potential is 13.90V
[28.31779972 28.31779972 0. ]
Of the 357911 tested reciprocal lattice points, 2109 have an excitation error less than 0.30 1/nm
Of the 2109 possible reflection 405 are allowed.
Of those, there are 57 in ZOLZ and 348 in HOLZ
np.float64(40.047416413731916)Plotting of Ewald Sphere¶
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(cent[0], cent[1], cent[2], c='green')
ax.plot((cent[0], 0), (cent[1], 0), (cent[2], 0))
# draw sphere
u, v = np.mgrid[0:2*np.pi:80j, 0:np.pi:40j]
x = np.cos(u)*np.sin(v)*K0_magnitude+cent[0]
y = np.sin(u)*np.sin(v)*K0_magnitude+cent[1]
z = np.cos(v)*K0_magnitude+cent[2]
ax.scatter(g_allowed[ZOLZ,0], g_allowed[ZOLZ,1], g_allowed[ZOLZ,2], c='red')
ax.scatter(g_allowed[HOLZ,0], g_allowed[HOLZ,1], g_allowed[HOLZ,2], c='blue')
ax.plot_surface( x, y, z, rstride=1, cstride=1, color='c', alpha=0.3, linewidth=0)
ax.set_xlim([-10, 10])
ax.set_ylim([-10, 10])
ax.set_zlim([-10, 10])
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
ax.set_aspect("equal");
plt.tight_layout();without sphere
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(g_allowed[ZOLZ,0], g_allowed[ZOLZ,1], g_allowed[ZOLZ,2], c='red')
ax.scatter(g_allowed[HOLZ,0], g_allowed[HOLZ,1], g_allowed[HOLZ,2], c='blue')
ax.scatter(cent[0], cent[1], cent[2], c='green')
ax.plot((cent[0], 0), (cent[1], 0), (cent[2],0))
33
ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([-5, 5])
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
ax.set_zlabel(r'scattering angle in (1/$\AA$)')
ax.set_aspect("equal");
plt.tight_layout(); Projection¶
We need to project the active reciprocal points onto the plane perpendicular to the zone axis.
For that we tilt the zone axis onto the z-axis
We use the spherical coordinates to determine the rotation matrix
def get_rotation_matrix(zone):
#spherical coordinates of zone
r = np.linalg.norm(zone)
theta = np.arccos(zone[2]/r)
if zone[0] < 0:
theta = -theta
if zone[0] == 0:
phi= np.pi/2
else:
phi = (np.arctan(zone[1]/zone[0]))
print('Rotation theta ',np.degrees(theta),' phi ',np.degrees(phi))
#first we rotate phi about z-axis
c, s = np.cos(phi), np.sin(phi)
rotz = np.array([[c, -s , 0],[s,c,0],[0,0,1]])
# second we rotate theta about y-axis
c, s = np.cos(theta), np.sin(theta)
roty = np.array([[c, 0 ,s],[0,1,0],[-s,0,c]])
# the rotation now makes z-axis coincide with plane normal
return np.dot(rotz,roty), np.degrees(theta), np.degrees(phi)
# zone axis in global coordinate system
zone_vector = np.dot(zone_hkl,reciprocal_unit_cell)
rotation_matrix, theta, phi = get_rotation_matrix(zone_vector)
print(np.round(rotation_matrix, 3))
print('\n Zone axis can now be rotated parallel to z axis')
print(np.round(np.dot(zone_vector,rotation_matrix), 4))Rotation theta 90.0 phi 45.0
[[ 0. -0.707 0.707]
[ 0. 0.707 0.707]
[-1. 0. 0. ]]
Zone axis can now be rotated parallel to z axis
[0. 0. 0.2604]
Rotation¶
We use the rotation matrix to rotate all lattice vectors
K0_vector_rotated = np.dot(cent, rotation_matrix)
cent_rotated = K0_vector_rotated
print(np.round(cent_rotated, 3))
g_hkl_rotated = np.dot(g_allowed, rotation_matrix)[ 0. 0. 40.047]
Now we can plot these diffraction spots in 2D by just setting the z-coordinate to zero. That is our projection procedure.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(g_hkl_rotated[ZOLZ,0], g_hkl_rotated[ZOLZ,1], g_hkl_rotated[ZOLZ,2]+10, c='red')
ax.scatter(g_hkl_rotated[HOLZ,0], g_hkl_rotated[HOLZ,1], g_hkl_rotated[HOLZ,2]+10, c='blue')
ax.scatter(g_hkl_rotated[ZOLZ,0], g_hkl_rotated[ZOLZ,1], 0, c='red')
ax.scatter(g_hkl_rotated[HOLZ,0], g_hkl_rotated[HOLZ,1], 0, c='blue')
ax.set_xlim([-10,10])
ax.set_ylim([-10,10])
ax.set_zlim([0,12])
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
ax.set_zlabel(r'scattering angle in (1/$\AA$)')
#ax.set_aspect("equal");
plt.tight_layout(); plt.show()In the above graph the z-axis is much smaller than the x-, y-axis.
In perspective to the incident beam wave vector the ewald sphere appears almost flat.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(g_hkl_rotated[ZOLZ,0], g_hkl_rotated[ZOLZ,1], g_hkl_rotated[ZOLZ,2]+10, c='red')
ax.scatter(g_hkl_rotated[HOLZ,0], g_hkl_rotated[HOLZ,1], g_hkl_rotated[HOLZ,2]+10, c='blue')
ax.scatter(cent_rotated[0], cent_rotated[1], cent_rotated[2], c='green')
ax.plot((cent_rotated[0], 0), (cent_rotated[1], 0), (cent_rotated[2],0))
#ax.set_aspect("equal");
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
ax.set_zlabel(r'scattering angle in (1/$\AA$)')
plt.tight_layout();2D projection¶
To compare this simulation with an experimental diffraction pattern, we only need the 2D projection.
Please note that all calculations were done without rotation. Only the plotting requires that. So rotation and setting the z-coordinate to zero is our projection procedure.
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(g_hkl_rotated[ZOLZ,0], g_hkl_rotated[ZOLZ,1], c='red')
ax.scatter(g_hkl_rotated[HOLZ,0], g_hkl_rotated[HOLZ,1], c='blue')
## Add indices
#for i in range (ZOLZ_plot.sum()):
# ax.text(g_hk_rotated[ZOLZ_plot][i,0]-2.5, g_hk_rotated[ZOLZ_plot][i,1]+1.2,hkl[ZOLZ_plot][i], fontsize=10)
#for i in range (HOLZ_plot.sum()):
# ax.text(g_hk_rotated[HOLZ_plot][i,0]-2.5, g_hk_rotated[HOLZ_plot][i,1]+1.2,hkl[HOLZ_plot][i], fontsize=10)
ax.set_aspect("equal")
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
ax.set_xlim([-6, 6])
ax.set_ylim([-6, 6])
plt.tight_layout();Error in 2D projection¶
For this projection, there is only a small error for high acceleration volatges because the Ewald sphere is very flat. But really we need to calculate where the rays from the center of the Ewald sphere cross the detector plane.
A easy way to accomplish this more accurate and universal projection is using spherical coordinates. The distance from center is related to
the polar angle and
the azimutal angle is directly related to the in plane distribution.
For the calculation, we need to put the origin in the center of the Ewald sphere, which is an addition to the z-variable for the rotated reciprocal lattice points.
And then with :
The angles with respect to the different quandrants are take care with the numpy function np.arctan2, which gives the angle to the positive x axis.
r= np.linalg.norm(g_hkl_rotated, axis=1)
g_spherical = g_hkl_rotated[r > 0]
r_spherical = r[r > 0]
theta = np.arctan( r_spherical/ K0_magnitude)
phi = np.arctan2(g_spherical[:, 0], g_spherical[:, 1])So now we can plot the diffraction pattern in mrad
x = theta*np.cos(phi)*1000 # in mrad
y = theta*np.sin(phi)*1000
plt.figure()
plt.scatter(x[ZOLZ[r > 0]], y[ZOLZ[r > 0]])
plt.scatter(x[HOLZ[r > 0]], y[HOLZ[r > 0]])
plt.axis('equal')
plt.xlabel('angle (mrad)')or in 1/nm with the trigonometeric determination of the radial distance with the magnitude of the incident wavevector
# ---- Input -----
rotation = np.radians(30)
# ----------------
x = r_spherical *np.cos(phi)*10 # 1/A*10 = 1/nm
y = r_spherical * np.sin(phi)*10
plt.figure()
plt.scatter(x[ZOLZ[r>0]], y[ZOLZ[r>0]])
plt.axis('equal')
plt.xlabel('reciprocal distance (1/nm)')As we will see later on, there are many advantages to keep the reciprocal reflections in polar coordinates.
One advantage is that to rotate the diffraction pattern we only need to add a value to the the azimutal angle
# ---- Input -----
rotation = np.radians(30)
# ----------------
x = r_spherical *np.cos(phi+rotation)*10
y = r_spherical * np.sin(phi+rotation)*10
plt.figure()
plt.scatter(x[ZOLZ[r>0]], y[ZOLZ[r>0]])
plt.axis('equal')
plt.xlabel('reciprocal distance (1/nm)')Forbidden Reflections¶
We can now pot the forbidden reflections into this diffraction pattern.
As can be seen the first order Laue Zone (SOLZ) is completely forbidden as are all the other odd Laue Zones.
fig = plt.figure()
ax = fig.add_subplot(111)
forbidden = np.dot(g_forbidden, rotation_matrix)
ax.scatter(forbidden[:,0], forbidden[:,1], c='green', alpha = 0.5, label='forbidden')
ax.scatter(g_hkl_rotated[ZOLZ,0], g_hkl_rotated[ZOLZ,1], c='red', label='ZOLZ')
ax.scatter(g_hkl_rotated[HOLZ,0], g_hkl_rotated[HOLZ,1], c='blue', label='HOLZ')
ax.set_aspect('equal')
ax.set_xlabel(r'scattering angle in (1/$\AA$)')
ax.set_ylabel(r'scattering angle in (1/$\AA$)')
plt.legend()
ax.set_xlim([-6, 6])
ax.set_ylim([-6, 6])
plt.tight_layout()Kinematic Scattering with pyTEMlib¶
the kinematic scattering plots can be calculated with the provided kinematic_scattering module of pyTEMlib. All data are gathered in a python dictionary named atoms.info .
All experimental parameter are stored in atoms.info[‘experimental’] dictionary and serve as an input to the kinematic scattering calculations.
All results from the kinematic scattering library are stored in the diff_dict dictionary
Any additional parameter for the output is stored in the diff_dict[‘output’] dictionary
# copy experimental conditions into the tags directory
atoms.info['experimental'] = {}
atoms.info['experimental']['acceleration_voltage_V'] = acceleration_voltage_V
atoms.info['experimental']['convergence_angle_mrad'] = 0
atoms.info['experimental']['zone_hkl'] = zone_hkl # incident nearest zone axis: defines Laue Zones!!!!
atoms.info['experimental']['mistilt'] = np.array([0, 0, 0]) # mistilt in degrees
atoms.info['experimental']['Sg_max'] = Sg_max # 1/Ang maximum allowed excitation error ; This parameter is related to the thickness
atoms.info['experimental']['hkl_max'] = hkl_max # Highest evaluated Miller indices
# calculate kinematic scattering data
diff_dict = {}
diff_dict = pyTEMlib.diffraction_tools.get_bragg_reflections(atoms, atoms.info['experimental'], verbose=True)
atoms.info['experimental']['plot_FOV'] = 4
# plot diffraction pattern
diff_dict['output'] = pyTEMlib.diffraction_tools.plot_saed_parameter()
diff_dict['output']['color_reflections'] = None
diff_dict['output']['plot_labels'] = False
diff_dict['output']['plot_Kikuchi'] = False
diff_dict['output']['plot_dynamically_allowed'] = False
v = pyTEMlib.diffraction_tools.plot_diffraction_pattern(diff_dict)Of the 1982 possible reflection 1920 are allowed.
Of those, there are 46 in ZOLZ and 1874 in HOLZ
Of the 54 forbidden reflection in ZOLZ 4 can be dynamically activated.
C:\Users\gduscher\AppData\Local\anaconda3\Lib\site-packages\pyTEMlib\diffraction_tools\diffraction_plot.py:793: UserWarning: No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
plt.legend()
Generating Spot Pattern¶
What needs to be done for a spot pattern
Generate the reciprocal lattice
Rotate zone axis to coincide with z-axis; Rotate reciprocal lattice the same way.
Select all reciprocal lattice points that have an excitation error smaller than S.
Determine the allowed reflections of those reciprocal lattice points.
Plot these spots (projection is done by ignoring z-axis value, which is the Laue zone).
Dynamically Activated Reflections¶
Forbidden Reflection that are a linear combination of two scattering events are called dynamically activated.
For example, the combination of reflection gives the in silicon forbidden reflection, visible in any diffraction pattern of silicon.
Here we test all combinations of any pair of allowed reflections to see whether they are a forbidden one, in which case they are dynamically activated or dynamically allowed
hkl_allowed = diff_dict['allowed']['hkl']
hkl_forbidden = diff_dict['forbidden']['hkl'].tolist()
indices = range(len(hkl_allowed))
combinations = [list(x) for x in itertools.permutations(indices, 2)]
dynamically_allowed = np.zeros(len(hkl_forbidden), dtype=bool)
for [i, j] in combinations:
possible = (hkl_allowed[i] + hkl_allowed[j]).tolist()
if possible in hkl_forbidden:
dynamically_allowed[hkl_forbidden.index(possible)] = True
fig = plt.figure()
points = pyTEMlib.diffraction_tools.plotting_coordinates(diff_dict['allowed']['g'])
plt.scatter(points[:,0], points[:,1], label='allowed', color='blue')
#plt.scatter(points[:,0], points[:,1])
points = pyTEMlib.diffraction_tools.plotting_coordinates(diff_dict['forbidden']['g'])
plt.scatter(points[np.logical_not(dynamically_allowed),0], points[np.logical_not(dynamically_allowed), 1],
color='green', alpha=0.2, label='forbidden')
plt.scatter(points[dynamically_allowed, 0], points[dynamically_allowed, 1],color='orange', alpha=0.5, label='dynamically allowed')
plt.legend()
plt.axis('equal');With kinematic_scattering library of pyTEMlib¶
diff_dict = pyTEMlib.diffraction_tools.get_bragg_reflections(atoms, atoms.info['experimental'], verbose=True)
# plot diffraction pattern
diff_dict['output'] = pyTEMlib.diffraction_tools.plot_saed_parameter()
# diff_dict['output']['color_reflections'] = None
diff_dict['output']['plot_labels'] = False
diff_dict['output']['plot_Kikuchi'] = False
diff_dict['output']['plot_dynamically_allowed'] = True
diff_dict['output']['plot_forbidden'] = False
v = pyTEMlib.diffraction_tools.plot_diffraction_pattern(diff_dict)
plt.gca().set_aspect('equal')Of the 1982 possible reflection 1920 are allowed.
Of those, there are 46 in ZOLZ and 1874 in HOLZ
Of the 54 forbidden reflection in ZOLZ 4 can be dynamically activated.
with forbidden reflections, that are dynamicall allowed
diff_dict['output'] = pyTEMlib.diffraction_tools.plot_saed_parameter()
# diff_dict['output']['color_reflections'] = None
diff_dict['output']['plot_labels'] = False
diff_dict['output']['plot_Kikuchi'] = False
diff_dict['output']['plot_dynamically_allowed'] = True
diff_dict['output']['plot_forbidden'] = True
v = pyTEMlib.diffraction_tools.plot_diffraction_pattern(diff_dict)
plt.gca().set_aspect('equal')Conclusion¶
The scattering geometry provides all the tools to determine which reciprocal lattice points are possible and which of them are allowed.
The diffraction pattern is a projection onto the plane perpendicular to the zone axis. For an easy projection we tilt everything so that the x,y plane is our projection plane.
Navigation¶
Chapter 2: Diffraction
List of Content: Front