Simulating powder x-ray diffraction in Python with Pymatgen and Xrayutilities.


Student: “Professor! I made a material!”

Professor: “That’s great! What did you make?”

Student: “I don’t know!”

Professor: “…”

X-ray diffraction (XRD) is probably the most widely-used characterization technique in materials science.  It can be used to probe the atomic structure of unknown compounds and to confirm that you have synthesized a particular phase of interest. Powder XRD, like the name suggests, can be used to characterize the crystal structures of powdered and polycrystalline samples, making it an extremely popular variant of the technique.

It essentially works by shooting x-rays with a known wavelength at a sample and watching for x-rays that come off of it at various angles. These angles correspond to distances between planes of atoms in the crystal structure(s) through a physical principle called Bragg’s Law and the intensity of the peaks at each angle give us information about the actual atoms that the sample is made of.


Figure 1: Illustration of Bragg’s Law. Image source: Hydrargyrum – Own work, CC BY-SA 3.0, 

Now, actually determining a structure from an XRD pattern is incredibly difficult in practice and only the most expert crystallographers have any confidence in their ability to do so. However solving the opposite problem of turning a structure into a simulated XRD pattern isn’t too hard, especially if you have a computer do it for you!

Currently, the Materials Project database includes tens of thousands of simulated XRD patterns. Unfortunately, these calculated patterns (stick plots) don’t do a great job of replicating what we see in the real-world (peaks broadening, shifting, etc.) As you can see below, these calculated patterns aren’t a dead-ringer for the real thing.

Screen Shot 2017-04-19 at 1.13.00 AM
Screen Shot 2017-04-19 at 1.21.29 AM

Figure 2: Real XRD pattern vs calculated pattern for \beta-CsCl (mp-22865) from

One of the projects that I’m working on involves doing more advanced simulations of XRD patterns that actually simulates the physics of a diffractometer. This is called the Fundamental Parameters Approach (FPA), an implementation of which is described in a great paper out of NIST. Using FPA to simulate diffraction patterns is useful because it produces extremely accurate patterns that replicate many of the features that we see in real diffraction patterns. The NIST team was also nice enough to release their implementation out into the wild for free, which means starving young researchers like me can play around with it!

While I’m doing shout-outs, I should also thank Eugen Wintersberger and Dominik Kriegner for writing a great open-source library for python called xrayutilities, which includes both an FPA code and a ton of other useful tools for working XRD data. Their library is by far the easiest way to get your hands on a FPA code.

Alright, now to start simulating some XRD patterns!

Let’s use the same high-temperature phase β-CsCl (Fm\overline{3}m) that I showed above as our example. We’ll first make the structure using Pymatgen and then simulate it’s diffraction pattern using xrayutilities.

In [1]:

# Set up some imports that we will need
from pymatgen import Lattice, Structure
import numpy
import xrayutilities as xru
from xrayutilities.materials.cif import CIFFile
from xrayutilities.materials.material import Crystal
from IPython.display import Image, display
from tempfile import NamedTemporaryFile
import matplotlib.pyplot as plt 
import matplotlib.ticker as ticker
%matplotlib inline
matplotlib.rcParams.update({'font.size': 18})
fig_size = [15, 12]
plt.rcParams["figure.figsize"] = fig_size
In [2]:
# Create beta-CsCl structure
a = 6.923 #Angstrom
latt = Lattice.cubic(a)
structure = Structure(latt, ["Cs", "Cs", "Cs", "Cs", "Cl", "Cl", 
"Cl", "Cl"], [[0, 0, 0], [0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0.5], [0, 0, 0.5], [0, 0.5, 0], [0.5, 0, 0]])

temp_cif = NamedTemporaryFile(delete=False)"cif",
xu_cif = CIFFile(
xu_crystal = Crystal(name="b-CsCl", lat=xu_cif.SGLattice())

two_theta = numpy.arange(0, 80, 0.01)

powder = xru.simpack.smaterials.Powder(xu_crystal, 1)
pm = xru.simpack.PowderModel(powder, I0=100)
intensities = pm.simulate(two_theta)
ax = plt.axes()
plt.title("XRD Pattern for " +
plt.xlabel("2 theta (degrees)")
Our simulated diffraction pattern is shown below:
Screen Shot 2017-04-19 at 1.12.53 AMFigure 3: Simulated XRD pattern of β-CsCl (mp-22865) using xrayutilities.
When we compare it with the experimental XRD pattern for β-CsCl below, we can see that our simulation pretty much nailed it!

Screen Shot 2017-04-19 at 1.13.00 AMFigure 4: Real XRD pattern of β-CsCl (mp-22865).