search

materialsproject / pymatgen

Python Materials Genomics (pymatgen) is a robust materials analysis code that defines classes for structures and molecules with support for many electronic structure codes. It powers the Materials Project.
https://pymatgen.org
Other
1.49k stars 857 forks source link

Is the Debye-Waller correction in XRDCalculator correct? #2827

Open Higomon opened 1 year ago

Higomon commented 1 year ago

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)

Graph2

Environment

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 commented 1 year ago

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 commented 1 year ago

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:

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:

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.