Röntgenbeugungssimulation mit Python

It’s ! Hallo. Dies ist der zweite Beitrag. Dieses Mal habe ich versucht, die Röntgenbeugung mit Python zu simulieren. Der vollständige Code (mit Jupyter Notebook) befindet sich auf GitHub [^ 1]

Röntgenbeugung (XRD)

Röntgenbeugung ist ein Phänomen, bei dem Röntgenstrahlen in einem Kristall gebeugt werden. Röntgenstrahlen sind elektromagnetische Wellen mit einer Wellenlänge innerhalb eines bestimmten Bereichs (etwa 1 nm bis 1 pm), und Beugung bedeutet, dass sich etwas mit Welleneigenschaften um Hindernisse wickelt und sich ausbreitet. Wenn die Wellen gebeugt werden, interferieren sie miteinander und heben sich gegenseitig auf, was zu einer Streifenverteilung führt. Bei der Analyse durch Röntgenbeugung wird die Kristallstruktur des Ziels durch Analyse des durch die Beugung erzeugten Musters identifiziert. Übrigens tritt das Beugungsphänomen nicht nur bei Röntgenstrahlen, sondern auch bei Elektronenstrahlen und Neutronenstrahlen auf, und Kristallstrukturanalyseverfahren, die diese verwenden, sind ebenfalls weit verbreitet.

Prinzip der Kristallstrukturanalyse durch Röntgenbeugung

In dem durch Röntgenbeugung erzeugten Muster werden die Punkte erhalten, wo sie durch Interferenz aufgrund von Beugung verstärkt werden (wo die Beugungsbedingungen erfüllt sind). Die Beziehung zwischen der Beugungsbedingung und der Kristallstruktur wird als Bragg-Bedingung bezeichnet. $2dsin{\theta} = n\lambda$ Es wird vertreten durch. d ist der Abstand zwischen den Ebenen (Kristallebenen), die von den Atomen im Kristall gebildet werden, $ \ theta $ ist der Winkel zwischen der Kristallebene und dem einfallenden Röntgenstrahl und $ \ lambda $ ist die Wellenlänge des Röntgenstrahls.

Atomstreufaktor und Kristallstrukturfaktor

In der zuvor erläuterten Bragg-Bedingung ist der Punkt das Atom, das Röntgenstrahlen streut, aber was Röntgenstrahlen tatsächlich streut, ist die um das Atom verteilte Elektronenwolke und die Amplitude der gestreuten Röntgenstrahlen. Ist proportional zur Elektronendichte am gestreuten Ort (wenn es sich um elastische Streuung handelt). Daher beträgt die Amplitude f des vom Atom gestreuten Röntgenstrahls $f = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ Es wird sein. Diese Gleichung entspricht der Fourier-Transformation der Elektronendichte $ \ rho (\ vec {r}) $ mit dem Streuungsvektor $ \ vec {k} $. Ein Streuvektor ist ein Vektor, der durch Komponenten (h, k, l) im inversen Raum definiert ist. Dieses f wird als atomarer Streufaktor bezeichnet. Die Elektronendichte ist erforderlich, um den Atomstreufaktor zu erhalten, aber da die Elektronendichte um das Atom nur in einfachen Fällen analytisch erhalten werden kann, wird sie grob nach der Hartley-Fock-Methode oder dergleichen berechnet. Da die atomaren Streufaktoren jedes Atoms bereits erhalten und in einer Datenbank gespeichert wurden, werden wir sie für diese Berechnung verwenden.

Darüber hinaus kann dieselbe Idee für den gesamten Kristall verwendet werden. Der Kristallstrukturfaktor F kann durch Fourier-Transformation der Elektronendichte im Kristall erhalten werden. $F(\vec{k}) = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ Die Formel entspricht dem atomaren Streufaktor. Im Kristall sind Elektronen jedoch nur um die Atomposition verteilt, so dass dies durch die Atomkoordinaten (x, y, z) und den Atomstreufaktor f umgeschrieben werden kann. $F(hkl) = \sum_{j}{f_jT_je^{2πi(hx_j + ky_j + lz_j)}}$ Es wird sein. h, k, l sind Komponenten des Streuungsvektors k und auch die Koordinaten der Gitterpunkte (umgekehrte Gitterpunkte) im umgekehrten Raum. $ T_j $, das ich später erläutern werde, ist der Device-Waller-Faktor.

Programm

Dies ist die eigentliche Implementierungsmethode. Listen Sie zunächst die inversen Gitterpunkte hkl auf. Ursprünglich sollten nur die Punkte auf der Ewald-Kugel aufgelistet werden, aber da bei der Röntgenwellenlänge ein ziemlich großer Bereich auf der Ewald-Kugel abgedeckt ist, werde ich so viele wie möglich auflisten. Hkl.csv auf GitHub ist eine aufgezählte hkl-Datei. Sie benötigen auch die Atomkoordinaten der Elementarzelle des Kristalls. Hier reichen nur die Atome in der Elementarzelle aus. Dieses Mal werden die Koordinaten des kubischen Gitters von Fläche zu Mitte (fcc) als Beispiel verwendet. In GitHub ist es eine Datei mit dem Namen pos.csv.

Als nächstes setzen wir aus den Gitterkonstanten a, b, c, $ \ alpha $, $ \ beta $, $ \ gamma $ (die wir im Voraus als Parameter haben) den Gittertensor $ G $ im realen Raum zusammen. $ 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) $ Nächster, $G^{\ast} = G^{-1}$ Daher kann durch Nehmen der inversen Matrix dieses Gittertensors der Gittertensor $ G ^ {\ ast} $ im inversen Raum erhalten werden.

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)

Unter Verwendung dieses inversen Gittertensors $ G ^ {\ ast} $ wird die Länge des inversen Gittervektors in Bezug auf die Reflexion K am inversen Gitterpunkt hkl $|\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}|$ Nachfragen. $ |\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|Aus der Bragg-Bedingung kann der Oberflächenabstand dk bei der Reflexion K und der Bragg-Winkel & thgr; k erhalten werden. $d_k = \frac{1}{|\vec{r}|^2}$ $\theta_k = \arcsin(\frac{λ}{2d_k})$

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

Die Werte in der Datenbank werden für den Atomstreufaktor verwendet. Da der Atomstreufaktor jedoch eine Funktion ist, die von sinθ / λ abhängt, wird nur der Koeffizient in der Datenbank aufgeführt. Aus Koeffizienten zusammensetzen $f(\frac{\sin\theta}{\lambda}) = \sum_{j=1}^{4}{a_j e^{-b_j(\frac{\sin\theta}{λ})^2}} + c$ Wird besorgt. $ a_j $, $ b_j $, c sind die Koeffizienten. Im Programm haben wir das Na-Atom angenommen und die Daten verwendet.

Jetzt müssen wir den Debay-Waller-Faktor erklären. Der Device-Waller-Faktor ist ein Faktor, der die Auswirkungen thermischer Vibrationen kompensiert. Der Debai-Waller-Faktor für isotrope thermische Schwingungen $T = e^{-B(\frac{\sin\theta}{\lambda})^2}$ Es wird vertreten durch. B ist der isotrope Atomverschiebungsparameter. Wenn die thermischen Schwingungen nicht isotrop sind (Anisotropie genannt), sind sie unter Verwendung der symmetrischen Matrix $ 3 \ times3 $ $ \ beta $ komplexer. $T = e^{-h^{\mathrm{T}} β h}$ Es wird sein. Diesmal ist unter Berücksichtigung der isotropen thermischen Schwingung B = 2,0. Ich habe nach dem Wert des Atomverschiebungsparameters gesucht, konnte ihn jedoch nicht finden, daher ist er ein geeigneter Wert. Entschuldigung.

Der Kristallstrukturfaktor wird basierend auf dem Atomstreufaktor und den Atomkoordinaten der Elementarzelle des Kristalls berechnet. $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)}}$ Wie Sie sehen können, ist der Kristallstrukturfaktor eine komplexe Zahl.

Wenn Sie das Quadrat des Absolutwerts des erhaltenen Kristallstrukturfaktors (da es das Quadrat des Absolutwerts der Streuamplitude ist, ist es die Streuintensität) auf der vertikalen Achse und $ 2 \ theta $ auf der horizontalen Achse darstellen, ein Diagramm, das durch Röntgenbeugung erhalten werden kann Es kann erhalten werden. Da der Winkel des einfallenden Röntgenstrahls $ \ theta $ und der Winkel des gestreuten Röntgenstrahls $ \ theta $ beträgt, ist es üblich, die horizontale Achse bei $ 2 \ theta $ zu nehmen.

Ergebnis

Die Handlung sieht folgendermaßen aus:

Der Bragg-Winkel wird auf 2 Stellen aufgerundet. In diesem Diagramm erscheint die Streuintensität in einer Delta-Funktion nur unter dem Reflexionswinkel, der die Bragg-Bedingung erfüllt, aber bei der tatsächlichen Röntgenbeugung hat der Peak der Streuintensität eine Breite. Lassen Sie uns die Ausbreitung des Peaks reproduzieren.

Wenn Sie der Delta-Funktion einen Spread geben möchten, fällt Ihnen als Erstes die Faltung durch die Gauß-Funktion (Gauß-Funktion) ein. Es ist die sogenannte Gaußsche Faltung. Das Folgende ist eine Faltung der Gaußschen Funktion in der vorherigen Darstellung. Die Dispersion hat einen Wert, der so aussieht. Sie können sehen, dass es an der Basis eine kleine Breite gibt. Wie jeder weiß, der die Messdaten der Röntgenbeugung gesehen hat, scheint diese Darstellung eine etwas engere Basis zu haben als die tatsächliche. Tatsächlich ist die Gaußsche Funktion geeignet, um Verhalten wie die thermische Schwingung von Partikeln auszudrücken, aber sie ist nicht geeignet, um den Peak eines solchen Spektrums auszudrücken. Falten wir uns daher mit der Lorentz-Funktion zusammen, mit der der Peak des Spektrums ausgedrückt wird. Die Lorentz-Funktion ist wie folgt. $L = \frac{1}{1 + (\frac{x}{γ})^2}$ $ \ Gamma $ wird als halbe Preisbreite bezeichnet, dh die Peakbreite bei 1/2 der Peakhöhe. Dieser Parameter ändert die Verbreitung der Funktion.

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

Die tatsächlich gefaltete ist wie folgt. Auch hier ist der Wert der halben Preisspanne einfach so. Es hat eine breitere Basis als die vorherige Gaußsche Funktion. Es sieht aus wie der Peak der Röntgenbeugung. Die Lorentz-Funktion allein ist jedoch nicht genau. Da die Atome im Kristall mit der Temperatur schwingen (selbst bei absolutem Nullpunkt gibt es keine Schwingung), beeinflussen thermische Schwingungen mit Gaußschem Funktionsverhalten auch die Ausbreitung des Peaks. Daher können die Lorentz-Funktion und die Gauß-Funktion mit einem beliebigen Gewicht gefaltet werden, um den Peak genauer auszudrücken.

Die Faltung der Lorenz-Funktion und der Gauß-Funktion wird als Vogt-Funktion bezeichnet. Vogt-Funktionen lassen sich mit komplexen Fehlerfunktionen schneller erstellen als einfach Lorentz- und Gauß-Funktionen. Ich bin mir wegen mangelnder Studien nicht sicher, aber es scheint, dass der Realteil der komplexen Fehlerfunktion wie eine Voigt-Funktion aussieht. [^ 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

Das Folgende ist eine Vogt-Funktion, die die Ausbreitung des Peaks ausdrückt. Es ist gut, es zu versuchen, aber es scheint, dass die Ausbreitung nur durch die Lorenz-Funktion ausgedrückt wird. Die Voigt-Funktion erfordert einen halben Preisbereich für jeden Gauß-Funktionsteil und Lorentz-Funktionsteil als Parameter, diesmal wird sie jedoch entsprechend eingestellt. Wenn Sie sich theoretisch der halben Preisspanne nähern, kann dies näher an der tatsächlichen Streuung des Peaks liegen.

Beugungsbild

Bei der Röntgenbeugung kann zusätzlich zu dem obigen Diagramm ein Bild erhalten werden, das als Beugungsbild bezeichnet wird und eine Kopie der Röntgenbeugung ist, wie sie ist. Da das inverse Gitter im Beugungsbild zu sehen ist, werde ich es diesmal als Bonus machen.

Als Implementierung numpy jeder Größe.zeros()Nur wenn die Reflexionsbedingung durch Erzeugen eines Arrays von erfüllt ist|F|^Geben Sie den Wert 2 ein und falten Sie die Gaußsche Funktion so, dass sie wie ein Punkt aussieht.

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

Im Beugungsbild erscheinen nur die inversen Gitterpunkte senkrecht zur Richtung des einfallenden Röntgenstrahls als Flecken. Daher muss zusätzlich zu den bestehenden Bedingungen auch die einfallende Röntgenrichtung berücksichtigt werden. Mit der einfallenden Röntgenrichtung als [001] werden nur die Reflexionen extrahiert, die das innere Produkt aufnehmen und 0 werden. Dieser Prozess ist einfach, wenn es [001] ist, aber wenn die Einfallsrichtung [110] oder [111] ist, ist eine komplizierte Verarbeitung erforderlich. In Python kann das Array mit matplotlib.pyplot.imshow () so wie es ist in ein Bild konvertiert werden, sodass es in ein Bild mit Schwarz-Weiß-Kontrast konvertiert wird.

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)

Es ist bequem. Das Bild in fcc ist wie folgt.

Es sieht so aus. Ich habe so etwas erst seit einiger Zeit gesagt. Ich habe es auch mit einem körperzentrierten kubischen Gitter (bcc) versucht. Es macht mehr Spaß als fcc.

Zusammenfassung

Ich habe versucht, Röntgenbeugung zu simulieren. Obwohl die Theorie der Beugung kompliziert ist, ist die Verarbeitung selbst durch einen Computer einfach, so dass sie ein gutes Lehrmaterial zum Verständnis der Röntgenbeugung sein kann. Um die Wahrheit zu sagen, müssen wir zur Durchführung einer präzisen Simulation mehr Faktoren berücksichtigen als die diesmal implementierten (Lorentz / Polarisationsfaktor, selektive Orientierungsfunktion usw.). Darüber hinaus sollten die diesmal entsprechend eingestellten Parameter (Atomverschiebungsparameter, Halbwertsbreite der Vogt-Funktion) auch genaue Werte verwenden. Können diese Parameter also theoretisch bestimmt werden? Diese Parameter bestehen aus einer komplexen Kombination verschiedener Faktoren. Selbst wenn Sie aus den dominanten Faktoren schließen können, ist es schwierig, sie analytisch zu bestimmen. Um genaue Werte für die obigen Parameter zu erhalten, verwenden wir eine Methode zum Abgleichen von experimentellen Werten mit Simulationen, um Fehler zu reduzieren. Also, der vorherige Artikel [^ 3]

Verweise

Dieser Artikel basiert auf den folgenden Artikeln und Büchern. Ich werde es ohne Erlaubnis vorstellen.

• Praxis der Pulverröntgenanalyse Einführung in die Reetbelt-Methode (herausgegeben von der Japan Society for Analytical Chemistry, Forschungsrunde der Röntgenanalyse, herausgegeben von Izumi Nakai, herausgegeben von Fujio Izumi)

• So zeichnen Sie die Voigt-Funktion in Python

Code

• https://github.com/PicricAcid/XRD