It’s ! Hello. This is the second post. This time, I tried to simulate X-ray diffraction with Python. The full code (using Jupyter Notebook) is on GitHub [^ 1]

X-Ray Diffraction (XRD)

X-ray diffraction is a phenomenon in which X-rays are diffracted in a crystal. X-rays are electromagnetic waves with wavelengths within a certain range (about 1 nm to 1 pm), and diffraction means that something with wave properties wraps around obstacles and propagates. When the waves diffract, they interfere with each other and cancel each other out, resulting in a striped distribution. In the analysis by X-ray diffraction, the crystal structure of the object is identified by analyzing the pattern generated by the diffraction. By the way, the diffraction phenomenon occurs not only with X-rays but also with electron beams and neutron rays, and crystal structure analysis methods using them are also widely used.

Principle of crystal structure analysis by X-ray diffraction

In the pattern created by X-ray diffraction, the spots that are strengthened by the interference caused by diffraction (where the diffraction conditions are satisfied) are obtained. The relationship between the diffraction condition and the crystal structure is called the Bragg condition. $2dsin{\theta} = n\lambda$ It is represented by. d is the distance between the planes (crystal planes) formed by the atoms in the crystal, $ \ theta $ is the angle between the crystal plane and the incident X-ray, and $ \ lambda $ is the wavelength of the X-ray.

Atomic scattering factor and crystal structure factor

In the Bragg condition explained above, the atom that scatters X-rays is used as a point, but it is the electron cloud distributed around the atom that actually scatters X-rays, and the amplitude of the scattered X-rays. Is proportional to the electron density at the scattered location (if it is elastic scattering). Therefore, the amplitude f of the X-rays scattered by the atoms is $f = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ It will be. This equation is equivalent to Fourier transforming the electron density $ \ rho (\ vec {r}) $ with the scattering vector $ \ vec {k} $. A scattering vector is a vector defined by components (h, k, l) in reciprocal space. This f is called an atomic scattering factor. The electron density is required to obtain the atomic scattering factor, but since the electron density around the atom cannot be calculated analytically except in simple cases, it is approximated by the Hartree-Fock method. Since the atomic scattering factors of each atom have already been obtained and stored in a database, we will use them in this calculation.

Moreover, the same idea can be used for the entire crystal. Crystal structure factor F can be obtained by Fourier transforming the electron density in the crystal. $F(\vec{k}) = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ The formula is the same as the atomic scattering factor. However, in a crystal, electrons are distributed only around the atomic position, so this can be rewritten with the atomic coordinates (x, y, z) and the atomic scattering factor f. $F(hkl) = \sum_{j}{f_jT_je^{2πi(hx_j + ky_j + lz_j)}}$ It will be. h, k, l are the components of the scattering vector k and are also the coordinates of the lattice points (reciprocal lattice points) in the reciprocal space. $ T_j $, which I will explain later, is the Debye-Waller factor.

program

This is the actual implementation method. First, list the reciprocal lattice points hkl. Originally, only the points on the Ewald sphere should be listed, but since a fairly wide range is covered on the Ewald sphere at the wavelength of X-rays, I will list as many as possible. Hkl.csv on GitHub is an enumerated hkl file. You will also need the atomic coordinates of the unit cell of the crystal. For this, only the atoms in the unit cell are sufficient. This time, the coordinates of the face-centered cubic lattice (fcc) are used as an example. On GitHub, it's a file named pos.csv.

Next, from the lattice constants a, b, c, $ \ alpha $, $ \ beta $, $ \ gamma $ (which we have in advance as parameters), we assemble the lattice tensor $ G $ in real space. $ G = \left( \begin{matrix} a^{2} & ab\cos\gamma & ac\cos\beta\\\ ba\cos\gamma & b^{2} & ba\cos\alpha\\\ ca\cos\beta & cb\cos\alpha & c^{2} \end{matrix} \right) $ next, $G^{\ast} = G^{-1}$ Therefore, if we take the inverse matrix of this lattice tensor, we can find the lattice tensor $ G ^ {\ ast} $ in the reciprocal space.

G = [[a**2, a*b*np.cos(gamma), a*c*np.cos(beta)], [b*a*np.cos(gamma), b**2, b*a*np.cos(alpha)], [c*a*np.cos(beta), c*b*np.cos(alpha), c**2]]
invG = np.linalg.inv(G)

Using this reciprocal lattice tensor $ G ^ {\ ast} , the length of the reciprocal lattice vector with respect to the reflection K at the reciprocal lattice point hkl $|\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}|$ Ask for. $ |\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}| = h^{\mathrm{T}}G^{\ast}h$$

hkl = [h[i], k[i], l[i]]
thkl = np.transpose(hkl)
dk = 1/((np.matmul(np.matmul(invG, thkl), hkl))**(1/2))

|r|From the Bragg condition, the surface spacing dk at the reflection K and the Bragg angle θk can be obtained. $d_k = \frac{1}{|\vec{r}|^2}$ $\theta_k = \arcsin(\frac{λ}{2d_k})$

sinthetak = lamb/(2*dk)
thetak = np.rad2deg(np.arcsin(sinthetak))

The values in the database are used for the atomic scattering factor, but since the atomic scattering factor is a function that depends on sinθ / λ, only the coefficient is listed in the database. To assemble from coefficients $f(\frac{\sin\theta}{\lambda}) = \sum_{j=1}^{4}{a_j e^{-b_j(\frac{\sin\theta}{λ})^2}} + c$ will do. $ a_j $, $ b_j $, c are the coefficients. In the program, we assumed Na atom and used the data.

Now we have to explain the Debye-Waller factor. The Debye-Waller factor is a factor that compensates for the effects of thermal vibration. The Debye-Waller factor for isotropic thermal vibrations $T = e^{-B(\frac{\sin\theta}{\lambda})^2}$ It is represented by. B is the isotropic atomic displacement parameter. If the thermal vibration is not isotropic (called anisotropy), it is more complex, using the $ 3 \ times3 $ symmetric matrix $ \ beta . $T = e^{-h^{\mathrm{T}} β h}$$ It will be. This time, considering isotropic thermal vibration, B = 2.0. I searched for the value of the atomic displacement parameter, but I could not find it, so it is an appropriate value. Excuse me.

The crystal structure factor is calculated based on the atomic scattering factor and the atomic coordinates of the unit cell of the crystal. $F_i(x, y, z, h_i, k_i, l_i, f_i) = f_i\sum_{j}{e^{2πi(hx_j + ky_j + lz_j)}}$ As you can see, the crystal structure factor is a complex number.

If you plot the square of the absolute value of the obtained crystal structure factor (the square of the absolute value of the scattering amplitude, so it becomes the scattering intensity) on the vertical axis and $ 2 \ theta $ on the horizontal axis, a plot that can be obtained by X-ray diffraction It can be obtained. Since the angle of incident X-rays is $ \ theta $ and the angle of scattered X-rays is $ \ theta $, it is customary to take the horizontal axis at $ 2 \ theta $.

result

The plot looks like this:

Bragg angle is rounded to two decimal places. In this plot, the scattering intensity appears delta-functionally only at the angle of reflection that satisfies the Bragg condition, but in actual X-ray diffraction, the peak of the scattering intensity has a width. Let's reproduce the spread of the peak.

When you want to give the delta function a spread, the first thing that comes to mind is the convolution by the Gaussian function (Gaussian). It's the so-called Gaussian Convolution. The following is a convolution of the Gaussian function in the previous plot. The variance has a value that looks like it. You can see that there is a little width at the base. As anyone who has seen the X-ray diffraction measurement data knows, this plot seems to have a slightly narrower base than the actual one. In fact, the Gaussian function is suitable for expressing the behavior of particles such as thermal vibration, but it is not suitable for expressing the peak of such a spectrum. Therefore, let's convolve with the Lorentz function used to express the peak of the spectrum. The Lorentz function is as follows. $L = \frac{1}{1 + (\frac{x}{γ})^2}$ $ \ Gamma $ is called the half width, which is the peak width at half the peak height. This parameter changes how the function spreads.

def lorentzian(theta, gamma):
    lorentz = 1/(1 + (theta/gamma)**2)
    return lorentz

The one that was actually folded is as follows. Once again, the value of the half width is an appropriate value. It has a wider base than the Gaussian function. It looks like the peak of X-ray diffraction. That said, the Lorentz function alone is not accurate. Because the atoms in the crystal vibrate with temperature (there is zero-point vibration even at absolute zero), thermal vibrations with Gaussian functional behavior also affect peak spread. Therefore, the Lorentz function and Gaussian function can be convoluted with arbitrary weight to express the peak more accurately.

The fold of the Lorentz function and the Gaussian function is called the Voigt function. The Voigt function is faster to create using the complex error function than to simply convolve the Lorentz and Gaussian functions. I'm not sure because of lack of study, but it seems that the real part of the complex error function looks like a Voigt function. [^ 2]

def voigt(theta, Hg, Hl):
    z = (theta + 1j*Hl)/(Hg * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    v = (w.real)/(Hg * np.sqrt(2.0*np.pi))   
    return v

The following is a Voigt function that expresses the spread of the peak. It's good to try it, but it seems that the spread is expressed only by the Lorenz function. The Voigt function requires a half-value width for each of the Gaussian function part and the Lorenz function part as parameters, but this time it is set appropriately. By theoretically approaching the half width, it may approach the actual spread of the peak.

Diffraction image

In X-ray diffraction, in addition to the above plot, an image called a diffraction image, which is a direct copy of the X-ray diffraction, can be obtained. Since the reciprocal lattice can be seen in the diffraction image, I will make it as a bonus this time.

As an implementation, numpy of any size.zeros()Only where the reflection condition is satisfied by generating an array of|F|^Enter a value of 2 and convolve the Gaussian function so that it looks like a spot.

def gaussian(theta, lamb):
    gauss = (1/(2*np.pi*sigma))*np.exp(-(theta**2 + lamb**2)/(2*sigma))
    return gauss

In the diffraction image, only the reciprocal lattice points perpendicular to the direction of the incident X-ray appear as spots. Therefore, in addition to the existing conditions, the incident X-ray direction must also be considered. With the incident X-ray direction as [001], only the reflection that takes the inner product and becomes 0 is extracted. This process is easy if it is [001], but if the incident direction is [110] or [111], complicated processing will be required. In Python, the array can be made into an image as it is with matplotlib.pyplot.imshow (), so make it an image with black and white contrast.

t, l = np.meshgrid(np.arange(-1, 1, 0.1), np.arange(-1, 1, 0.1))
g = gaussian(t, l)
xrdimage = scipy.signal.convolve2d(I, g, boundary='wrap', mode='same')
plt.imshow(xrdimage, interpolation="nearest", cmap=plt.cm.gray)

It's convenient. The image in fcc is as below.

It looks like that. I've only said that kind of thing since a while ago. I also tried it with body-centered cubic lattice (bcc). It's more fun than fcc.

Summary

I tried to simulate X-ray diffraction. Diffraction theory is complicated, but the processing itself by computer is easy, so it may be a good teaching material to understand X-ray diffraction. To tell the truth, in order to perform a precise simulation, we have to consider more factors than those implemented this time (Lorentz / polarization factor, selective orientation function, etc.). In addition, the parameters set appropriately this time (atomic displacement parameter, half width of Voigt function) should also use accurate values. So can these parameters be theoretically determined? Since these parameters are made up of various factors that influence each other in a complex manner, it will be difficult to obtain them analytically, even if they can be inferred from the dominant factors. In order to obtain accurate values for the above parameters, we use a method of matching experimental values with simulations to reduce errors. So, the previous article [^ 3]

References

This article is based on the following articles and books. I will introduce it without permission.

--Practice of powder X-ray analysis Introduction to the Reetbelt method (edited by the Japan Society for Analytical Chemistry X-ray analysis research round-table conference, edited by Izumi Nakai, edited by Fujio Izumi)

-How to plot the Voigt function in python

code

https://github.com/PicricAcid/XRD

X-ray diffraction simulation in Python