Is the Debye-Waller correction in XRDCalculator correct?

[Higomon]

Describe the bug

I am not an expert in XRD.

An XRD expert on Twitter(https://bit.ly/3jhhoDi) pointed out that rough X-ray diffraction intensity calculations ignoring Debye-Waller factors are not valid.

Therefore, I tried simulating by changing the Debye-Waller factors as shown in the attached figure. Even when I raised the Debye-Waller factor of iron, the intensity seems too high on the high-angle side compared to the VESTA results.

Is the target formula (line 243) of xrd.py correct?

By referring to the expert's book on XRD, I modified it as follows and got reasonable results.

pi2 = np.pi ** 2
243: dw_correction = np.exp(-2 * pi2 * ndwfactors * s2)

I would like the expert to verify this.

pymatgen.analysis.diffraction.xrd module

xrd.py

243:  dw_correction = np.exp(-dwfactors * s2)

Environment

• OS: Windows 10 Pro 22H2, build 19045.2486
• Python: 3.9.13
• pymatgen: 2022.0.17

CIF file

gamma-Fe2O3_fiz15840.zip

---

Code

from os.path import dirname, basename
import pandas as pd
import numpy as np
from scipy.special import wofz

from pymatgen.io.cif import CifParser
from pymatgen.analysis.diffraction.xrd import XRDCalculator

# setting --------------------------------
x_str = 10  # /deg (start of X range)
x_end = 120  # /deg (end of X range)

# Debye-Waller factor
DWF = {"Fe": 0.15, "O": 0.6}
# ----------------------------------------


def voigt(x, amp, pos, fwhm, shape):
    """
    x: x data
    amp: amplitude(peak height)
    pos: peak center
    fwhm: Full Width at Half Maximum
    shape:
    """
    tmp = 1 / wofz(np.zeros((len(x))) + 1j * np.sqrt(np.log(2.0)) * shape).real
    tmp = (
        tmp
        * amp
        * wofz(
            2 * np.sqrt(np.log(2.0)) * (x - pos) / fwhm
            + 1j * np.sqrt(np.log(2.0)) * shape
        ).real
    )
    return tmp


def sum_voigt_peaks(pat, x_str, x_end):
    x = np.arange(x_str, x_end, 0.01)
    calc_y = np.zeros(len(x))

    N = len(pat.x)
    fwhm = 0.06
    shape = 0.5
    for i in range(N):
        pos = pat.x[i]
        amp = pat.y[i]
        calc_y += voigt(x, amp, pos, fwhm, shape)

    return x, calc_y


def output_hklTable(pat, f):
    d_space = pat.d_hkls
    _hkl = pat.hkls

    hkl_index = [i[0]["hkl"] for i in _hkl]
    m = [i[0]["multiplicity"] for i in _hkl]
    dict = {
        "hkl": hkl_index,
        "multiplicity": m,
        "d": d_space,
        "2theta": pat.x,
        "Int": pat.y,
    }
    df = pd.DataFrame(dict)
    df.index = df.index + 1
    d_name = dirname(f)
    f_name = basename(f).replace(".cif", ".csv")
    df.to_csv(
        d_name + "\\hkl_" + f_name,
        encoding="utf-8",
        index=True,
    )


def output_XRDpattern(pat, f):
    x, y = sum_voigt_peaks(pat, x_str, x_end)
    dict = {"2theta": x, "Int": y}
    df = pd.DataFrame(dict)
    d_name = dirname(f)
    f_name = basename(f).replace(".cif", ".csv")
    df.to_csv(
        d_name + "\\prof_" + f_name,
        encoding="utf-8",
        index=False,
    )


def main():
    f = r"D:\gamma-Fe2O3_fiz15840.cif"
    parser = CifParser(f)
    structure = parser.get_structures()[0]
    lat = XRDCalculator(
        wavelength=1.54056,
        debye_waller_factors=DWF,
    )
    pat = lat.get_pattern(structure, scaled=True, two_theta_range=(x_str, x_end))

    output_hklTable(pat, f)

    output_XRDpattern(pat, f)


if __name__ == "__main__":
    main()

---

[Andrew-S-Rosen]

Could you share a reference for the expression in the exponential? I'm not an expert in XRD but can help out if I have a resource to look at.

[Higomon]

Thank you for your reply.

I only have the text written in Japanese by Fujio Izumi.

Instead, I asked ChatGPT4, and the answer is below. I'm no expert in crystallography, but the following seems to be a reasonable answer.

---

Isotropic Debye-Waller factor

In X-ray diffraction, the Debye-Waller factor (also known as the temperature factor) accounts for the reduction in diffracted intensity due to atomic vibrations. It's typically denoted as B, and it's related to the mean square displacement of atoms, <u²>, via the equation:

B = 8π²<u²>

This factor is applied to the scattering intensity of X-rays using the formula:

I = I0 exp(-2Bsin²θ/λ²)

where:

• I is the scattered intensity
• I0 is the initial intensity (at 0°, or when the temperature effects are ignored)
• B is the Debye-Waller factor
• θ is the scattering angle
• λ is the wavelength of the X-ray

Here, B assumes isotropic vibrations, i.e., atoms vibrate equally in all directions, which may not be the case in many real-world applications.

Anisotropic Debye-Waller factor

The Debye-Waller factor (B) assumes isotropic vibrations, i.e., that atoms vibrate equally in all directions. However, in many real-world cases, atomic vibrations can be anisotropic, meaning they vary depending on the direction.

In cases of anisotropic vibrations, the Debye-Waller factor is replaced by an anisotropic displacement parameter (ADP), also known as anisotropic temperature factor. This anisotropic factor is usually expressed as a symmetric 3x3 tensor (matrix):

U =
| U11 U12 U13 |
| U12 U22 U23 |
| U13 U23 U33 |

Each element of this matrix represents mean-square displacement along specific crystallographic directions.

When this anisotropic displacement parameter is incorporated, the intensity equation becomes:

I = I0 exp[-2π²(h²a²U11 + k²b²U22 + l²c²U33 + 2hkabU12 + 2hlccU13 + 2klbc*U23)]

where:

• h, k, l are Miller indices
• a*, b*, c* are reciprocal lattice vectors

This allows for more accurate modeling of atomic vibrations in X-ray diffraction studies, particularly in real-world cases where these vibrations are likely to be anisotropic.

Please note that the above form is a general form and may take different forms depending on the specific crystal structure or experimental conditions.