Calculating Theoretical Powder Patterns from CIF Files

The end goal of this part of the plugin is to give presenters an easy way to generate theoretical data that can be visually compared to experimental data. In addition to simulating experimental peak broading, it is also useful to be able to compare relative intensities with identified phase fractions to verify that the experimental peaks are fully accounted for by the identified phases.

Unfortunately, the previously established method (in our research group) for generating theoretical patterns through VESTA results in an arbitrarily normalized pattern with intensities that cannot be compared with other calculated patterns.

There is, however, an already-established algorithm for scaling intensities by phase fractions which is used in Rietveld refinement. With this in mind, the CIF import module in the plugin is constructed to approximate a Rietveld-simulated powder pattern. It is not meant to be a simulation with enough accuracy to refine structures, but it is enough of an approximation to visually compare with experimental data, including relative phase fractions.

Source Code

For reference, the main calculation module is in the calculate_pattern() module in PXRD_cifImp.py:

Click to expand code 
def calculate_pattern(
    cif_path,
    fe_wavelengths,
    fe_weights,
    two_theta_range,
    step,
    U=0.0,
    V=0.0,
    W=0.012,
    X=0.0,
    Y=0.0,
    axial_S=0.015
):
    # Read CIF file using pymatgen's Structure module
    structure = Structure.from_file(cif_path)

    # Build list of B values for each atom site
    atom_B = []
    pi2 = np.pi**2
    for site in structure.sites:
        props = site.properties
        if "B_iso" in props:
            B = props["B_iso"]
        elif "Uiso" in props:
            B = 8 * pi2 * props["Uiso"]
        elif "Ueq" in props:
            B = 8 * pi2 * props["Ueq"]
        else:
            B = 0.0
        atom_B.append(B)

    # Generate 2theta column and empty intensity column of same length.
    tmin, tmax = two_theta_range
    two_theta = np.arange(tmin, tmax + step, step)
    intensity = np.zeros_like(two_theta)
    fe_weights = np.array(fe_weights)/np.sum(fe_weights)
    
    # Generate intensity as the sum of all fe intensities by weights.
    for wl, wt in zip(fe_wavelengths, fe_weights):
        # Basic reflections calculated using pymatgen's XRDCalculator.get_pattern()
        xrd = XRDCalculator(wavelength=wl)
        pattern = xrd.get_pattern(structure, two_theta_range=two_theta_range, scaled=False)
        Z = compute_Z(structure)
        cell_volume = structure.lattice.volume
        mu = absorption_proxy(structure) # This is very fudgy
        # Modify reflections with scattering factors
        for idx, (t0, I0) in enumerate(zip(pattern.x, pattern.y)):
            # Useful constants
            theta = np.radians(t0 / 2)
            sin_th = np.sin(theta)
            cos_th = np.cos(theta)
                
            # Useful constants
            s = sin_th / wl
            pi2 = np.pi**2

            # Extra dampening defined at head of this file.
            B_extra = _B_EXTRA

            # Debye-Waller damping by B-factors
            DW_atoms = np.mean([np.exp(-2 * pi2 * B * s**2) for B in atom_B])
            # Fudge factor to match experimental/VESTA heights. (not currently in use)
            DW_extra = np.exp(-2 * pi2 * B_extra * s**2)

            # Modify base intensity with damping
            I0 *= DW_atoms #* DW_extra
            I0 /= Z
            I0 /= cell_volume

            # Caglioti broadening: Gaussian (H_G) and Lorentzian (HL) hybrid
            H_G = np.sqrt(U*np.tan(theta)**2 + V*np.tan(theta) + W)
            H_L = X*np.tan(theta) + Y/np.cos(theta)

            # Each peak is calculated indepedently as a function of all 2theta values.
            pv = tch_pseudo_voigt(two_theta, t0, H_G, H_L)
            asym = fcj_asymmetry(two_theta, t0, H_G, S=axial_S)

            # Peak is a normalized Voight-shape, with tailing added, multiplied by intensity and wavelength weight.
            # Each individual peak is constructed as a function of all 2theta values. They are all overlaid onto the intensity series.
            intensity += wt * I0 * pv * asym

    # Returns series, not indiviual values
    return two_theta, intensity

Approximation of peak shapes and intensities with the Rietveld model

The classical Rietveld method of calculating diffraction intensity for a given sample is:

\[I_{calc}(2\theta) = S_F \sum_{k}^{phases}\left[\frac{\phi_k}{V_k^2}\sum_{j}^{peaks}\left(m_{k,j}L_k\left|F_{k,j}\right|^2S_j\left(2\theta-2\theta_{k,j}\right)P_{k,j}A_{j}\right)\right] + bkg\]

Reflection Intensities

First, we focus on the peak-scaled intensities outlined in the expression:

\[\sum_{j}^{peaks}\left(m_{k,j}L_k\left|F_{k,j}\right|^2 S_j\left(2\theta-2\theta_{k,j}\right)P_{k,j}A_{j}\right)\]

Where:

  • $m_{k,j}$ is the peak multiplicity
  • $L_k$ is the Lorentz-Polarization factor
  • $\left\vert F_{k,j}\right\vert^2$ is the structure factor
  • $S_j\left(2\theta-2\theta_{k,j}\right)$ is the shape function for a peak centered at $2\theta_{k,j}$.
  • $P_{k,j}$ is the March-Dollase factor for preferred orientation
  • $A_{j}$ is an absorption factor.

In the python code, the base peak intensities are generated by making use of pymatgen’s XRDCalculator.get_pattern() method:

xrd = XRDCalculator(wavelength=wl)
pattern = xrd.get_pattern(structure, two_theta_range=two_theta_range, scaled=False)
# In the get_pattern() method, setting scaled=False prevents normalization to [0,100]

This generates a set of sharp peak intensities (pattern = $I_{peak}$) that account for multiplicity, Lorentz-Polarization and structure factor only:

\[I_{peak} = m_{k,j}L_k\left|F_{k,j}\right|^2\]

Since the patterns we import are intended for visualization, not refinement, we also choose to neglect preferred orientation and absorption effects:

\[\sum_{j}^{peaks}\left(m_{k,j}L_k\left|F_{k,j}\right|^2 S_j\left(2\theta-2\theta_{k,j}\right)P_{k,j}A_{j}\right) \approx \sum_{j}^{peaks}\left(m_{k,j}L_k\left|F_{k,j}\right|^2 S_j\left(2\theta-2\theta_{k,j}\right)\right) = \sum_{j}^{peaks}\left(I_{peak}S_j\left(2\theta-2\theta_{k,j}\right)\right)\]

Peak Shapes and Broadening

Now we must apply the shape function, $S_j\left(2\theta-2\theta_{k,j}\right)$, which converts theoretical peak intensities into realistically-broadened peak shapes. In this plugin, I opted to construct the peak function as a combination of:

  • Phase-averaged Debye-Waller Peak Dampening (sometimes included in $\left\vert F_{k,j}\right\vert^2$, but not in the case of pymatgen’s get_pattern())
  • Caglioti/Pseudo-Voight Peak Broadening (Hybrid of Gaussian and Lorentzian peak shape)
  • Finger-Cox-Jephcoat axial divergence asymmetry for peak tailing

In the python code, the output data starts as an empty data set:

two_theta = np.arange(tmin, tmax + step, step)
intensity = np.zeros_like(two_theta)

Then, it loops through every peak intensity in pattern:

for idx, (t0, I0) in enumerate(zip(pattern.x, pattern.y)):

For each peak intensity, it first multiplies intensity by an approximated Debye-Waller dampening factor averaged from all atoms in the structure (strictly speaking, dampening is unique to each atom, but for visualizing patterns, the average approximation is sufficient):

sin_th = np.sin(theta)
s = sin_th / wl # wl = wavelength
pi2 = np.pi**2
DW_atoms = np.mean([np.exp(-2 * pi2 * B * s**2) for B in atom_B])
I0 *= DW_atoms

Next, it generates ($2\theta$, Intensity) peak shape data sets for each individual peak:

H_G = np.sqrt(U*np.tan(theta)**2 + V*np.tan(theta) + W)
H_L = X*np.tan(theta) + Y/np.cos(theta)

pv = tch_pseudo_voigt(two_theta, t0, H_G, H_L)
asym = fcj_asymmetry(two_theta, t0, H_G, S=axial_S)
# Note that the variable two_theta was previously defined as a series of all two_theta values, so each peak is generated as an individual function for the entire x-axis.

Before finally multiplying the peak shapes by the dampened intensity and adding each peak to the intially empty data set:

intensity += wt * I0 * pv * asym
#wt refers to the relative weight of the current wavelength (for doublet splitting)

This generates a full diffraction pattern that approximates a single phase per unit cell without relative scaling (two_theta, intensity = $I_{k}\left(2\theta\right)$ ).

\[I_{k}\left(2\theta\right) = \sum_{j}^{peaks}\left(I_{peak}S_j\left(2\theta-2\theta_{k,j}\right)\right)\]

Scaling by Molar Phase Fraction

Now, we only have to figure out what factor to multiply each phase pattern by in order for it to be imported on a molar basis. Substituting our calculations and approximations so far back into our Rietveld model, we get:

\[I_{calc}(2\theta) = S_F \sum_{k}^{phases}\left[\frac{\phi_k}{V_k^2}I_{k}\left(2\theta\right)\right] + bkg\]

For theoretical patterns, we can neglect background. Also, since we will be normalizing our phases after relative scaling, we can also set $S_F = 1$:

\[I_{calc}(2\theta) =\sum_{k}^{phases}\left[\frac{\phi_k}{V_k^2}I_{k}\left(2\theta\right)\right]\]

Now, the scale factor $\frac{\phi_k}{V_k^2}$ is what we must account for, where:

  • $\phi_k$ is the phase volume fraction.
  • $V_k$ is the phase unit cell volume.

If we express volume fraction as a function of a molar phase fraction by formula unit, $x_k$:

\[\phi_k \propto x_k \frac{V_k}{Z_k}\]

Substituting back in, we get:

\[I_{calc}(2\theta)\propto\sum_{k}^{phases}\left[\frac{x_k}{Z_k V_k}I_{k}\left(2\theta\right)\right]\]

Since we will be normalizing to match experimental data, maintaining proportionality is sufficient. We will also be multiplying by $x_k$ after import to allow for dynamically changing phase fractions, so we should isolate that from our imported pattern:

\(I_{calc}(2\theta)\propto\sum_{k}^{phases}\left[x_k I_{k,imported}(2\theta)\right]\) \(I_{k,imported}(2\theta) = \frac{I_{k}\left(2\theta\right)}{Z_k V_k}\)

To modify our calculated $I_k$ patterns into $I_{k,imported}$, we only need to divide our peak intensities by $Z_k$ and $V_k$ before adding them to our data set.

I0 /= Z
I0 /= cell_volume

'''
Peak shape functions pv and asym defined here
'''

intensity += wt * I0 * pv * asym

Finally, this results in a simulated pattern per formula unit for each imported phase. When importing CIF patterns using the “Phase Fraction Analysis” normalization mode, there is an additional parameter row for each phase where the user can specify molar fractions. This automatically scales the resulting patterns before normalizing the whole data set.

Choosing Your Phase Fraction

To simplify writing the dynamic formulas to Origin, CIF patterns will always be imported for molar scaling. However, the user can now specify which type of phase fraction they want to adjust during the import menu. Whichever option is selected will be the row that is visible in Origin. The molar row is then auto-calculated before being used for pattern scaling. Since patterns are normalized after phase-scaling, the conversions to molar basis only need to be proportional, not normalized to 100%.

Mole Fraction

This mode is useful if you already know the relative molar fractions of each phase. The patterns are already imported on a molar basis, so scaling is as simple as multiplying by the user’s input mole fraction (done automatically by column formula)

Weight Fraction

Some Reitveld applications (like Match! software’s embedded FullProf) output phase fractions directly as weight fractions ($w_k$). A user input row for weight fraction is added, and then the hidden mole fraction row is calculated using molecular weight:

\[x_k\propto\frac{w_k}{MW_k}\]

Cell Fraction (GSAS “Phase Fraction” Parameter)

In GSAS-II, the refined parameter for phase fraction is actually cell fraction, $c_k$; that is, the fraction of unit cells contributed by each phase in a given space:

\[c_k = \frac{\phi_k}{V_k}\]

This is done because the cell fraction is useful in other calculations during refinement (whereas volume fraction is not), so it simplifies the number of calculated variables when refining.

When using this option, the hidden molar fraction row is calculated using Z:

\[x_k\propto Z c_k\]
This site uses Just the Docs, a documentation theme for Jekyll.