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 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.
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.
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
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.
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 = [[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}
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.
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
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
The crystal structure factor is calculated based on the atomic scattering factor and the atomic coordinates of the unit cell of the crystal.
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 $.
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.
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.
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.
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]
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
Recommended Posts