From Bragg's Law to Measured Data Analysis
By studying this chapter, you will acquire the following knowledge and skills:
X-rays are electromagnetic waves with wavelengths ranging from approximately 0.01 to 10 nm (10 Å). In materials science, primarily Cu Kα radiation (λ = 1.5406 Å) and Mo Kα radiation (λ = 0.7107 Å) are used.
Since the wavelength of X-rays is on the same order as interatomic distances (several Å), diffraction phenomena occur due to crystal lattices.
| X-ray Source | Wavelength (Å) | Energy (keV) | Primary Applications |
|---|---|---|---|
| Cu Kα | 1.5406 | 8.05 | Powder XRD, general structure analysis |
| Mo Kα | 0.7107 | 17.48 | Single crystal XRD, when short wavelength needed |
| Co Kα | 1.7902 | 6.93 | Iron-containing samples (avoiding fluorescence) |
| Synchrotron | Variable | Variable | High-brilliance, high-resolution measurements |
When X-rays enter a crystal, the following process occurs:
Each atom in the crystal scatters X-rays, and the scattered waves interfere. When the phases of scattered waves align at specific angles, they constructively interfere, observed as diffraction peaks.
In 1913, William Lawrence Bragg and his father William Henry Bragg proposed a simple yet powerful model that treats X-ray diffraction as specular reflection from crystal planes.
Where:
Consider X-rays incident at angle θ on crystal planes spaced at distance d. If the path difference between waves reflected from the upper and lower planes equals an integer multiple of wavelength λ, constructive interference occurs.
The path difference is geometrically calculated as:
The condition for constructive interference is:
Bragg's law is a necessary condition for diffraction to occur, but not sufficient. In practice, the structure factor Fhkl must also be non-zero (discussed later).
Calculate diffraction angles from major crystal planes of silicon (Si):
import numpy as np
import matplotlib.pyplot as plt
def cubic_d_spacing(a, h, k, l):
"""
Calculate interplanar spacing for cubic crystal
Parameters:
-----------
a : float
Lattice parameter (Å)
h, k, l : int
Miller indices
Returns:
--------
float : Interplanar spacing (Å)
"""
return a / np.sqrt(h**2 + k**2 + l**2)
def bragg_angle(d_hkl, wavelength, n=1):
"""
Calculate diffraction angle 2θ from Bragg's law
Parameters:
-----------
d_hkl : float
Interplanar spacing (Å)
wavelength : float
X-ray wavelength (Å)
n : int
Order of reflection
Returns:
--------
float : Diffraction angle 2θ (degrees), or None if diffraction not possible
"""
sin_theta = n * wavelength / (2 * d_hkl)
if abs(sin_theta) > 1:
return None # Diffraction condition not satisfied
theta = np.arcsin(sin_theta)
return np.degrees(2 * theta) # Return 2θ
# Silicon (Si) parameters
a_Si = 5.4310 # Å
wavelength_CuKa = 1.5406 # Å
print("=== X-ray Diffraction Pattern Prediction for Silicon (Si) ===")
print(f"Lattice parameter: a = {a_Si} Å")
print(f"X-ray wavelength: λ = {wavelength_CuKa} Å (Cu Kα)\n")
# Major crystal planes
planes = [
(1, 1, 1), (2, 2, 0), (3, 1, 1),
(4, 0, 0), (3, 3, 1), (4, 2, 2)
]
print(f"{'(hkl)':<10} {'d (Å)':<12} {'2θ (deg)':<12}")
print("-" * 40)
results = []
for hkl in planes:
h, k, l = hkl
d = cubic_d_spacing(a_Si, h, k, l)
two_theta = bragg_angle(d, wavelength_CuKa)
if two_theta is not None:
print(f"({h}{k}{l}){'':<8} {d:8.4f} {two_theta:8.3f}")
results.append((hkl, two_theta))
# Visualize peak pattern in graph
fig, ax = plt.subplots(figsize=(12, 6))
for (h, k, l), two_theta in results:
ax.axvline(two_theta, color='red', linewidth=2, alpha=0.7)
ax.text(two_theta, 1.05, f'({h}{k}{l})',
rotation=90, va='bottom', ha='right', fontsize=9)
ax.set_xlim(20, 100)
ax.set_ylim(0, 1.2)
ax.set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
ax.set_ylabel('Relative Intensity (arbitrary units)', fontsize=14, fontweight='bold')
ax.set_title('Theoretical XRD Pattern of Silicon (Si)', fontsize=16, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig('si_xrd_pattern.png', dpi=150, bbox_inches='tight')
plt.show()
print("\nXRD pattern saved: si_xrd_pattern.png")
Compare how diffraction angles change for the same crystal plane with Cu Kα and Mo Kα radiation:
def compare_xray_sources():
"""Compare diffraction patterns with different X-ray sources"""
# X-ray source parameters
sources = {
'Cu Kα': 1.5406,
'Mo Kα': 0.7107,
'Co Kα': 1.7902
}
# Aluminum (Al) parameters
a_Al = 4.0495 # Å
planes = [(1, 1, 1), (2, 0, 0), (2, 2, 0), (3, 1, 1)]
print("=== Comparison of Diffraction Angles with Different X-ray Sources (Al) ===\n")
print(f"Lattice parameter: a = {a_Al} Å\n")
# Calculations for each X-ray source
fig, ax = plt.subplots(figsize=(14, 8))
colors = {'Cu Kα': 'red', 'Mo Kα': 'blue', 'Co Kα': 'green'}
for source_name, wavelength in sources.items():
print(f"--- {source_name} (λ = {wavelength} Å) ---")
print(f"{'(hkl)':<10} {'d (Å)':<12} {'2θ (deg)':<12}")
print("-" * 40)
y_offset = list(sources.keys()).index(source_name) * 0.3
for hkl in planes:
h, k, l = hkl
d = cubic_d_spacing(a_Al, h, k, l)
two_theta = bragg_angle(d, wavelength)
if two_theta is not None:
print(f"({h}{k}{l}){'':<8} {d:8.4f} {two_theta:8.3f}")
# Display as bar graph
ax.plot([two_theta, two_theta], [y_offset, y_offset + 0.25],
color=colors[source_name], linewidth=3)
if y_offset == 0: # Display labels only for first source
ax.text(two_theta, y_offset + 0.28, f'({h}{k}{l})',
rotation=90, va='bottom', ha='center', fontsize=8)
print()
# Graph decoration
ax.set_xlim(0, 150)
ax.set_ylim(-0.1, 1.0)
ax.set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
ax.set_yticks([0.125, 0.425, 0.725])
ax.set_yticklabels(['Cu Kα', 'Mo Kα', 'Co Kα'])
ax.set_title('Comparison of Aluminum (Al) Diffraction Patterns with Different X-ray Sources',
fontsize=16, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
# Legend
from matplotlib.lines import Line2D
legend_elements = [Line2D([0], [0], color=color, lw=3, label=name)
for name, color in colors.items()]
ax.legend(handles=legend_elements, loc='upper right', fontsize=12)
plt.tight_layout()
plt.savefig('xray_source_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
print("Comparison graph saved: xray_source_comparison.png")
# Execute
compare_xray_sources()
Even if Bragg's law is satisfied, not all reflections are observed. The structure factor Fhkl determines the actual diffraction intensity.
Where:
Important: When Fhkl = 0, diffraction does not occur even if Bragg's law is satisfied. This is called systematic absence.
import numpy as np
import pandas as pd
def structure_factor(positions, f_atoms, h, k, l):
"""
Calculate structure factor F_hkl
Parameters:
-----------
positions : list of tuples
Fractional coordinates of atoms in unit cell [(x1,y1,z1), (x2,y2,z2), ...]
f_atoms : list of float
Atomic scattering factor for each atom
h, k, l : int
Miller indices
Returns:
--------
complex : Structure factor F_hkl
"""
F = 0 + 0j
for (x, y, z), f in zip(positions, f_atoms):
phase = 2 * np.pi * (h*x + k*y + l*z)
F += f * np.exp(1j * phase)
return F
def analyze_structure_factors():
"""Analyze structure factors for different lattice types"""
# Atomic positions for each lattice type
structures = {
'SC (Simple Cubic)': [
(0, 0, 0)
],
'BCC (Body-Centered Cubic)': [
(0, 0, 0),
(0.5, 0.5, 0.5)
],
'FCC (Face-Centered Cubic)': [
(0, 0, 0),
(0.5, 0.5, 0),
(0.5, 0, 0.5),
(0, 0.5, 0.5)
],
'Diamond Structure': [
(0, 0, 0),
(0.5, 0.5, 0),
(0.5, 0, 0.5),
(0, 0.5, 0.5),
(0.25, 0.25, 0.25),
(0.75, 0.75, 0.25),
(0.75, 0.25, 0.75),
(0.25, 0.75, 0.75)
]
}
# List of Miller indices
planes = [
(1, 0, 0), (1, 1, 0), (1, 1, 1),
(2, 0, 0), (2, 2, 0), (3, 1, 1),
(2, 2, 2), (4, 0, 0), (3, 3, 1)
]
print("=== Structure Factors and Extinction Rules for Different Lattice Types ===\n")
for structure_name, positions in structures.items():
print(f"\n【{structure_name}】")
print(f"Number of atoms in unit cell: {len(positions)}\n")
print(f"{'(hkl)':<10} {'|F_hkl|^2':<15} {'Observable':<10} {'Note'}")
print("-" * 60)
# Assume all atoms are identical with atomic scattering factor f = 1
f_atoms = [1.0] * len(positions)
for hkl in planes:
h, k, l = hkl
F = structure_factor(positions, f_atoms, h, k, l)
F_squared = abs(F)**2
# Determine observability (|F|^2 > 0.01 considered observable)
observable = "○" if F_squared > 0.01 else "×"
# Explanation of extinction conditions
remarks = ""
if structure_name == 'BCC (Body-Centered Cubic)':
if (h + k + l) % 2 != 0:
remarks = "h+k+l is odd → extinct"
elif structure_name == 'FCC (Face-Centered Cubic)':
if not (h % 2 == k % 2 == l % 2):
remarks = "h,k,l mixed → extinct"
print(f"({h}{k}{l}){'':<8} {F_squared:12.4f} {observable:<10} {remarks}")
print("\n" + "="*60)
print("Summary of Extinction Rules:")
print(" SC : All reflections are observed")
print(" BCC: Only h+k+l even are observed")
print(" FCC: Only all even or all odd h,k,l are observed")
print(" Diamond: FCC + additionally only h+k+l=4n (n:integer) show strong reflections")
print("="*60)
# Execute
analyze_structure_factors()
Extinction rules are extremely important for determining crystal structures. For example, if the (100) peak is not observed, you know it is not simple cubic (SC) but BCC or FCC.
The actual diffraction intensity Ihkl depends on many factors beyond the structure factor:
| Factor | Name | Physical Meaning |
|---|---|---|
| |Fhkl|2 | Structure factor | Effect of atomic arrangement in unit cell |
| mhkl | Multiplicity | Number of symmetrically equivalent planes |
| L | Lorentz factor | Geometric effect of crystal |
| P | Polarization factor | Effect of X-ray polarization state |
| A | Absorption factor | X-ray absorption by sample |
| exp(-2M) | Temperature factor (Debye-Waller factor) | Diffuse scattering due to atomic thermal vibration |
from itertools import permutations, product
def multiplicity(h, k, l, crystal_system='cubic'):
"""
Calculate multiplicity (number of equivalent planes)
Parameters:
-----------
h, k, l : int
Miller indices
crystal_system : str
Crystal system
Returns:
--------
int : Multiplicity
"""
planes = set()
if crystal_system == 'cubic':
for perm in permutations([abs(h), abs(k), abs(l)]):
for signs in product([1, -1], repeat=3):
plane = tuple(s * p for s, p in zip(signs, perm))
if plane != (0, 0, 0):
planes.add(plane)
return len(planes)
def lorentz_polarization_factor(two_theta):
"""
Lorentz-polarization factor for powder X-ray diffraction
Parameters:
-----------
two_theta : float
Diffraction angle 2θ (radians)
Returns:
--------
float : LP factor
"""
theta = two_theta / 2
LP = (1 + np.cos(two_theta)**2) / (np.sin(theta)**2 * np.cos(theta))
return LP
def calculate_intensity_with_factors():
"""Calculate diffraction intensity considering various factors"""
# Aluminum (FCC, a = 4.0495 Å)
a_Al = 4.0495
wavelength = 1.5406 # Cu Kα
# Atomic positions for FCC structure
fcc_positions = [
(0, 0, 0),
(0.5, 0.5, 0),
(0.5, 0, 0.5),
(0, 0.5, 0.5)
]
f_Al = 13 # Aluminum atomic number (simplified as atomic scattering factor)
planes = [
(1, 1, 1), (2, 0, 0), (2, 2, 0),
(3, 1, 1), (2, 2, 2), (4, 0, 0)
]
print("=== XRD Intensity Calculation with Various Factors (Al, FCC) ===\n")
print(f"{'(hkl)':<10} {'d (Å)':<10} {'2θ':<10} {'|F|^2':<12} {'m':<6} {'LP':<10} {'I_rel':<10}")
print("-" * 80)
intensities = []
for hkl in planes:
h, k, l = hkl
# Interplanar spacing
d = cubic_d_spacing(a_Al, h, k, l)
# Bragg angle
two_theta_deg = bragg_angle(d, wavelength)
if two_theta_deg is None:
continue
two_theta_rad = np.radians(two_theta_deg)
# Structure factor
F = structure_factor(fcc_positions, [f_Al]*4, h, k, l)
F_squared = abs(F)**2
# Check extinction rule (mixed indices extinct for FCC)
if F_squared < 0.01:
continue
# Multiplicity
m = multiplicity(h, k, l, 'cubic')
# Lorentz-polarization factor
LP = lorentz_polarization_factor(two_theta_rad)
# Relative intensity (temperature factor and absorption omitted)
I_rel = F_squared * m * LP
intensities.append((hkl, I_rel))
print(f"({h}{k}{l}){'':<8} {d:8.4f} {two_theta_deg:8.2f} {F_squared:10.2f} {m:<6} {LP:8.4f} {I_rel:8.2f}")
# Normalize by maximum intensity
max_intensity = max(I for _, I in intensities)
print("\n--- Normalized Relative Intensity (maximum = 100) ---")
print(f"{'(hkl)':<10} {'Relative Intensity':<15}")
print("-" * 30)
for hkl, I in intensities:
I_normalized = 100 * I / max_intensity
print(f"({hkl[0]}{hkl[1]}{hkl[2]}){'':<8} {I_normalized:8.1f}")
# Execute
calculate_intensity_with_factors()
Powder X-ray Diffraction (PXRD) is a measurement method for samples consisting of fine crystalline grains oriented in random directions. It is one of the most frequently used structural characterization techniques in materials science.
Generate realistic XRD patterns including peak shapes (Gaussian functions) and background:
def gaussian_peak(two_theta, center, intensity, fwhm):
"""
Gaussian peak function
Parameters:
-----------
two_theta : array
Array of 2θ values
center : float
Peak center position
intensity : float
Peak intensity
fwhm : float
Full Width at Half Maximum
Returns:
--------
array : Intensity of Gaussian peak
"""
sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
return intensity * np.exp(-((two_theta - center)**2) / (2 * sigma**2))
def simulate_xrd_pattern(material_name, a, c=None, crystal_system='cubic',
positions=None, f_atoms=None,
wavelength=1.5406, two_theta_range=(20, 100),
fwhm=0.2, background=50):
"""
Simulate complete XRD pattern
Parameters:
-----------
material_name : str
Material name
a, c : float
Lattice parameters
crystal_system : str
Crystal system
positions : list
Fractional coordinates of atoms in unit cell
f_atoms : list
Atomic scattering factors
wavelength : float
X-ray wavelength
two_theta_range : tuple
Measurement range of 2θ
fwhm : float
Full width at half maximum of peaks
background : float
Background intensity
Returns:
--------
two_theta, intensity : arrays
"""
# Array of 2θ values
two_theta = np.linspace(two_theta_range[0], two_theta_range[1], 4000)
intensity = np.ones_like(two_theta) * background # Background
# Calculate peaks
max_hkl = 5
peaks_info = []
for h in range(max_hkl + 1):
for k in range(h, max_hkl + 1):
for l in range(k, max_hkl + 1):
if h == 0 and k == 0 and l == 0:
continue
# Interplanar spacing
if crystal_system == 'cubic':
d = cubic_d_spacing(a, h, k, l)
# Bragg angle
two_theta_peak = bragg_angle(d, wavelength)
if two_theta_peak is None or two_theta_peak > two_theta_range[1]:
continue
# Structure factor
if positions is not None and f_atoms is not None:
F = structure_factor(positions, f_atoms, h, k, l)
F_squared = abs(F)**2
if F_squared < 0.01:
continue
else:
F_squared = 1.0
# Multiplicity
m = multiplicity(h, k, l, crystal_system)
# Lorentz-polarization factor
LP = lorentz_polarization_factor(np.radians(two_theta_peak))
# Peak intensity
I_peak = F_squared * m * LP
# Add Gaussian peak
intensity += gaussian_peak(two_theta, two_theta_peak, I_peak, fwhm)
peaks_info.append(((h, k, l), two_theta_peak, I_peak))
# Normalize intensity
intensity = (intensity / intensity.max()) * 1000
return two_theta, intensity, peaks_info
# Simulate XRD pattern for silicon
a_Si = 5.4310
diamond_positions = [
(0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5),
(0.25, 0.25, 0.25), (0.75, 0.75, 0.25),
(0.75, 0.25, 0.75), (0.25, 0.75, 0.75)
]
f_Si = [14] * 8 # Silicon atomic number
two_theta, intensity, peaks = simulate_xrd_pattern(
'Silicon (Si)',
a=a_Si,
crystal_system='cubic',
positions=diamond_positions,
f_atoms=f_Si,
fwhm=0.15
)
# Display graph
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(two_theta, intensity, 'b-', linewidth=1.5, label='Simulation')
ax.fill_between(two_theta, 0, intensity, alpha=0.2, color='blue')
# Label major peaks
peaks_sorted = sorted(peaks, key=lambda x: x[2], reverse=True)[:6]
for (h, k, l), pos, I in peaks_sorted:
ax.annotate(f'({h}{k}{l})',
xy=(pos, I * 1000 / intensity.max()),
xytext=(pos, I * 1000 / intensity.max() + 80),
ha='center', fontsize=10,
arrowprops=dict(arrowstyle='->', color='red', lw=1))
ax.set_xlim(20, 100)
ax.set_ylim(0, 1100)
ax.set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
ax.set_ylabel('Intensity (arbitrary units)', fontsize=14, fontweight='bold')
ax.set_title('Powder XRD Pattern of Silicon (Si) (Simulation)',
fontsize=16, fontweight='bold')
ax.grid(axis='both', alpha=0.3)
ax.legend(fontsize=12)
plt.tight_layout()
plt.savefig('si_powder_xrd_simulation.png', dpi=150, bbox_inches='tight')
plt.show()
print("Powder XRD simulation saved: si_powder_xrd_simulation.png")
Read actual XRD measurement data (text file) and automatically detect peaks:
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
def read_xrd_data(filename):
"""
Read XRD data file
Common format:
2theta Intensity
20.0 150.2
20.1 152.3
...
Parameters:
-----------
filename : str
Path to data file
Returns:
--------
two_theta, intensity : arrays
"""
try:
data = np.loadtxt(filename, skiprows=1) # Skip header row
two_theta = data[:, 0]
intensity = data[:, 1]
return two_theta, intensity
except FileNotFoundError:
print(f"File {filename} not found.")
print("Generating sample data.")
# Generate sample data
return simulate_xrd_pattern('Sample', a=5.0, fwhm=0.3)[:2]
def detect_peaks_in_xrd(two_theta, intensity, prominence=50, distance=10):
"""
Detect peaks in XRD pattern
Parameters:
-----------
two_theta : array
2θ data
intensity : array
Intensity data
prominence : float
Threshold for peak detection (prominence)
distance : int
Minimum distance between peaks (number of data points)
Returns:
--------
peak_positions, peak_intensities : arrays
"""
peaks_idx, properties = find_peaks(intensity,
prominence=prominence,
distance=distance)
peak_positions = two_theta[peaks_idx]
peak_intensities = intensity[peaks_idx]
peak_prominences = properties['prominences']
return peak_positions, peak_intensities, peak_prominences
def analyze_xrd_data():
"""Demonstration of measured XRD data analysis"""
# Read data (simulate if actual file doesn't exist)
two_theta, intensity = read_xrd_data('sample_xrd.txt')
# Peak detection
peak_pos, peak_int, peak_prom = detect_peaks_in_xrd(
two_theta, intensity,
prominence=100,
distance=20
)
print("=== Detected Peaks ===\n")
print(f"{'Peak Number':<12} {'2θ (deg)':<12} {'Intensity':<15} {'d-spacing (Å)':<12}")
print("-" * 60)
wavelength = 1.5406 # Cu Kα
for i, (pos, intensity_val) in enumerate(zip(peak_pos, peak_int), 1):
# Calculate d-spacing from Bragg's law
theta = np.radians(pos / 2)
d_spacing = wavelength / (2 * np.sin(theta))
print(f"{i:<12} {pos:10.2f} {intensity_val:12.1f} {d_spacing:10.4f}")
# Display graph
fig, ax = plt.subplots(figsize=(14, 7))
# Plot XRD pattern
ax.plot(two_theta, intensity, 'b-', linewidth=1.5, label='Measured Data')
# Mark detected peaks
ax.plot(peak_pos, peak_int, 'ro', markersize=8, label='Detected Peaks')
# Display peak numbers
for i, (pos, int_val) in enumerate(zip(peak_pos, peak_int), 1):
ax.annotate(f'{i}',
xy=(pos, int_val),
xytext=(pos, int_val + 80),
ha='center', fontsize=10, fontweight='bold',
bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.7))
ax.set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
ax.set_ylabel('Intensity (arbitrary units)', fontsize=14, fontweight='bold')
ax.set_title('XRD Pattern and Peak Detection', fontsize=16, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(axis='both', alpha=0.3)
plt.tight_layout()
plt.savefig('xrd_peak_detection.png', dpi=150, bbox_inches='tight')
plt.show()
print("\nPeak detection results saved: xrd_peak_detection.png")
# Execute
analyze_xrd_data()
Mathematically fit peak shapes to determine accurate peak positions and widths:
def pseudo_voigt(x, amplitude, center, fwhm, eta):
"""
Pseudo-Voigt function (mixture of Gaussian and Lorentzian components)
Parameters:
-----------
x : array
Data points
amplitude : float
Peak amplitude
center : float
Peak center
fwhm : float
Full width at half maximum
eta : float
Fraction of Lorentzian component (0: Gaussian, 1: Lorentzian)
Returns:
--------
array : Peak shape
"""
# Gaussian component
sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
gaussian = np.exp(-((x - center)**2) / (2 * sigma**2))
# Lorentzian component
gamma = fwhm / 2
lorentzian = gamma**2 / ((x - center)**2 + gamma**2)
# Mixture
return amplitude * (eta * lorentzian + (1 - eta) * gaussian)
def fit_single_peak(two_theta, intensity, peak_center, window=2.0):
"""
Fit single peak
Parameters:
-----------
two_theta : array
2θ data
intensity : array
Intensity data
peak_center : float
Approximate peak center position
window : float
Fitting range (±window degrees)
Returns:
--------
popt : array
Optimized parameters [amplitude, center, fwhm, eta]
pcov : array
Covariance matrix
"""
# Extract fitting range
mask = (two_theta >= peak_center - window) & (two_theta <= peak_center + window)
x_data = two_theta[mask]
y_data = intensity[mask]
# Initial estimates
amplitude_init = np.max(y_data) - np.min(y_data)
center_init = peak_center
fwhm_init = 0.2
eta_init = 0.5
p0 = [amplitude_init, center_init, fwhm_init, eta_init]
# Boundary conditions
bounds = ([0, peak_center - 1, 0.05, 0],
[amplitude_init * 2, peak_center + 1, 1.0, 1])
try:
popt, pcov = curve_fit(pseudo_voigt, x_data, y_data, p0=p0, bounds=bounds)
return popt, pcov, x_data, y_data
except RuntimeError:
print(f"Fitting failed for peak at {peak_center:.2f}°.")
return None, None, x_data, y_data
def demo_peak_fitting():
"""Demonstration of peak fitting"""
# Generate sample data
two_theta, intensity = simulate_xrd_pattern('Sample', a=5.4, fwhm=0.2)[:2]
# Peak detection
peak_pos, _, _ = detect_peaks_in_xrd(two_theta, intensity, prominence=100)
# Fit three strongest peaks
strongest_peaks = sorted(zip(peak_pos, intensity[np.isin(two_theta, peak_pos)]),
key=lambda x: x[1], reverse=True)[:3]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
for ax, (peak_center, _) in zip(axes, strongest_peaks):
result = fit_single_peak(two_theta, intensity, peak_center, window=2.5)
if result[0] is not None:
popt, pcov, x_fit, y_fit = result
amplitude, center, fwhm, eta = popt
# Fitting curve
x_fine = np.linspace(x_fit.min(), x_fit.max(), 500)
y_fine = pseudo_voigt(x_fine, *popt)
# Plot
ax.plot(x_fit, y_fit, 'bo', markersize=4, label='Measured Data')
ax.plot(x_fine, y_fine, 'r-', linewidth=2, label='Fitting')
# Display results
textstr = f'Center: {center:.3f}°\nFWHM: {fwhm:.3f}°\nη: {eta:.2f}'
ax.text(0.05, 0.95, textstr, transform=ax.transAxes,
fontsize=10, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
ax.set_xlabel('2θ (degrees)', fontsize=12, fontweight='bold')
ax.set_ylabel('Intensity', fontsize=12, fontweight='bold')
ax.set_title(f'Peak @ {center:.1f}°', fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('xrd_peak_fitting.png', dpi=150, bbox_inches='tight')
plt.show()
print("Peak fitting results saved: xrd_peak_fitting.png")
# Execute
demo_peak_fitting()
pymatgen is a powerful Python library for materials science that can automatically generate XRD patterns from crystal structures.
try:
from pymatgen.core import Structure, Lattice
from pymatgen.analysis.diffraction.xrd import XRDCalculator
PYMATGEN_AVAILABLE = True
except ImportError:
print("pymatgen is not installed.")
print("Install with: pip install pymatgen")
PYMATGEN_AVAILABLE = False
def generate_xrd_with_pymatgen():
"""Generate XRD pattern using pymatgen"""
if not PYMATGEN_AVAILABLE:
print("This example is skipped because pymatgen is not installed.")
return
# Define silicon structure
lattice = Lattice.cubic(5.4310)
species = ['Si', 'Si', 'Si', 'Si', 'Si', 'Si', 'Si', 'Si']
coords = [
[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5],
[0.25, 0.25, 0.25], [0.75, 0.75, 0.25],
[0.75, 0.25, 0.75], [0.25, 0.75, 0.75]
]
si_structure = Structure(lattice, species, coords)
print("=== XRD Pattern Generation with pymatgen ===\n")
print(f"Crystal structure: {si_structure.composition}")
print(f"Space group: {si_structure.get_space_group_info()}\n")
# Initialize XRD calculator
calculator = XRDCalculator(wavelength='CuKa') # Cu Kα radiation
# Calculate XRD pattern
pattern = calculator.get_pattern(si_structure, two_theta_range=(20, 100))
print(f"{'2θ (deg)':<12} {'d-spacing (Å)':<15} {'(hkl)':<15} {'Relative Int.':<12}")
print("-" * 60)
for i in range(len(pattern.x)):
two_theta = pattern.x[i]
intensity = pattern.y[i]
hkl = pattern.hkls[i][0]['hkl'] # Get first hkl
d_spacing = pattern.d_hkls[i]
print(f"{two_theta:10.2f} {d_spacing:12.4f} {str(hkl):<15} {intensity:10.1f}")
# Display graph
fig, ax = plt.subplots(figsize=(14, 7))
# Plot as bar graph
ax.vlines(pattern.x, 0, pattern.y, colors='blue', linewidth=2, label='pymatgen')
# Label peaks with hkl
for i, (two_theta, intensity, hkls_data) in enumerate(zip(pattern.x, pattern.y, pattern.hkls)):
if intensity > 20: # Display label only for peaks with intensity > 20
hkl = hkls_data[0]['hkl']
ax.text(two_theta, intensity + 5, f'({hkl[0]}{hkl[1]}{hkl[2]})',
rotation=90, va='bottom', ha='center', fontsize=9)
ax.set_xlim(20, 100)
ax.set_ylim(0, max(pattern.y) * 1.15)
ax.set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
ax.set_ylabel('Relative Intensity', fontsize=14, fontweight='bold')
ax.set_title('XRD Pattern of Silicon (Si) - pymatgen Generated',
fontsize=16, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(axis='both', alpha=0.3)
plt.tight_layout()
plt.savefig('si_xrd_pymatgen.png', dpi=150, bbox_inches='tight')
plt.show()
print("\npymatgen XRD pattern saved: si_xrd_pymatgen.png")
# Compare multiple materials
compare_materials_xrd()
def compare_materials_xrd():
"""Compare XRD patterns of multiple materials"""
if not PYMATGEN_AVAILABLE:
return
# Define materials
materials = {
'Si (Diamond)': Structure(
Lattice.cubic(5.4310),
['Si']*8,
[[0,0,0], [0.5,0.5,0], [0.5,0,0.5], [0,0.5,0.5],
[0.25,0.25,0.25], [0.75,0.75,0.25], [0.75,0.25,0.75], [0.25,0.75,0.75]]
),
'Al (FCC)': Structure(
Lattice.cubic(4.0495),
['Al']*4,
[[0,0,0], [0.5,0.5,0], [0.5,0,0.5], [0,0.5,0.5]]
),
'Fe (BCC)': Structure(
Lattice.cubic(2.8665),
['Fe']*2,
[[0,0,0], [0.5,0.5,0.5]]
)
}
calculator = XRDCalculator(wavelength='CuKa')
fig, axes = plt.subplots(3, 1, figsize=(14, 12))
for ax, (name, structure) in zip(axes, materials.items()):
pattern = calculator.get_pattern(structure, two_theta_range=(20, 100))
# Bar graph
ax.vlines(pattern.x, 0, pattern.y, colors='darkblue', linewidth=2.5)
# Peak labels
for two_theta, intensity, hkls_data in zip(pattern.x, pattern.y, pattern.hkls):
if intensity > 15:
hkl = hkls_data[0]['hkl']
ax.text(two_theta, intensity + 3, f'({hkl[0]}{hkl[1]}{hkl[2]})',
rotation=90, va='bottom', ha='center', fontsize=9)
ax.set_xlim(20, 100)
ax.set_ylim(0, 110)
ax.set_ylabel('Relative Intensity', fontsize=12, fontweight='bold')
ax.set_title(name, fontsize=14, fontweight='bold', loc='left')
ax.grid(axis='x', alpha=0.3)
axes[-1].set_xlabel('2θ (degrees)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('materials_xrd_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
print("Material comparison graph saved: materials_xrd_comparison.png")
# Execute
generate_xrd_with_pymatgen()
Rietveld refinement is a technique that fits the entire powder XRD pattern with a crystal structure model. It was developed by Hugo Rietveld in 1969.
The Rietveld method optimizes the following function by least squares:
Where:
The calculated intensity is expressed as:
The Rietveld method is a technique that requires a structure model. If the initial structure model is significantly incorrect, it will not converge to the correct solution. Typically, initial models are created from known similar structures or single crystal XRD data.
For the face-centered cubic (FCC) structure of copper (Cu) with lattice parameter a = 3.615 Å. Using Cu Kα radiation (λ = 1.5406 Å) for powder XRD measurement, at what angle (2θ) will diffraction peaks from the following (hkl) planes be observed?
Also, explain whether the (100) peak is observed according to the FCC extinction rule.
Calculation of interplanar spacings:
Calculation of Bragg angles (λ = 2d sinθ):
About the (100) plane:
The FCC extinction rule is "only all even or all odd h, k, l are observed." (100) has h=1 (odd), k=0 (even), l=0 (even), so it is mixed, and the structure factor F100 = 0. Therefore it is not observed.
The following peaks were observed in an XRD measurement: (110), (200), (211), (220), (310), (222), (321), (400)
Is this material simple cubic (SC), body-centered cubic (BCC), or face-centered cubic (FCC)? Determine from the extinction rules.
Check the sum h+k+l for each index:
All peaks have even h+k+l. This matches the BCC (body-centered cubic) extinction rule.
If it were FCC, mixed indices like (210) or (221) would be extinct, but there is no such regularity in the observed peaks. If it were SC, all peaks including (100) should be observed.
Answer: BCC (body-centered cubic)
From the XRD pattern of an unknown sample, the following d-spacings (Å) were obtained: 3.35, 2.46, 2.13, 1.91, 1.80
Which of the following materials could this sample be? Determine from lattice parameters.
Calculate representative d-spacings for each material:
A) NaCl (FCC, a=5.64 Å):
B) Si (a=5.43 Å):
C) Graphite (hexagonal, a=2.46, c=6.71 Å):
Answer: C) Graphite
The maximum d-spacing of 3.35 Å matches the characteristic (002) plane spacing of graphite. This corresponds to the important peak for graphite interlayer distance.
Under the following conditions, create a program that simulates the XRD pattern of a fictional material and displays the positions and relative intensities of the major peaks (top 5):
Use a combination of Code Examples 1 and 3. Implement FCC atomic positions and structure factor calculations, Bragg angle calculations using Bragg's law, and intensity calculations considering multiplicity.
The completed program should display peaks like (111), (200), (220), (311), (222) at correct angles and relative intensities.
In this chapter, we learned the principles of X-ray diffraction and practical analysis methods:
In the next chapter, we will learn crystal structure visualization and analysis, and acquire practical skills to retrieve and analyze structures from actual materials databases.