Simulation de diffraction des rayons X avec Python

It’s ! Bonjour. Ceci est le deuxième article. Cette fois, j'ai essayé de simuler la diffraction des rayons X avec Python. Le code complet (utilisant Jupyter Notebook) est sur GitHub [^ 1]

Diffraction des rayons X (XRD)

La diffraction des rayons X est un phénomène dans lequel les rayons X sont diffractés dans un cristal. Les rayons X sont des ondes électromagnétiques avec une longueur d'onde dans une certaine plage (environ 1 nm à 1 pm), et la diffraction signifie que quelque chose avec des propriétés d'onde s'enroule autour des obstacles et se propage. Lorsque les ondes sont diffractées, elles interfèrent les unes avec les autres et s'annulent, entraînant une distribution striée. Dans l'analyse par diffraction des rayons X, la structure cristalline de la cible est identifiée en analysant le motif généré par la diffraction. À propos, le phénomène de diffraction se produit non seulement avec les rayons X, mais aussi avec les faisceaux d'électrons et les rayons neutroniques, et les méthodes d'analyse de la structure cristalline qui les utilisent sont également largement utilisées.

Principe de l'analyse de la structure cristalline par diffraction des rayons X

Dans le diagramme créé par diffraction des rayons X, les taches qui sont renforcées par les interférences dues à la diffraction (où les conditions de diffraction sont satisfaites) sont obtenues. La relation entre la condition de diffraction et la structure cristalline est appelée la condition de Bragg. $2dsin{\theta} = n\lambda$ Il est représenté par. d est la distance entre les plans (plans cristallins) formés par les atomes dans le cristal, $ \ theta $ est l'angle entre le plan cristallin et le rayon X incident, et $ \ lambda $ est la longueur d'onde du rayon X.

Facteur de diffusion atomique et facteur de structure cristalline

Dans la condition de Bragg expliquée plus haut, le point est l'atome qui diffuse les rayons X, mais ce qui diffuse réellement les rayons X est le nuage d'électrons distribué autour de l'atome et l'amplitude des rayons X diffusés. Est proportionnel à la densité électronique à l'emplacement diffusé (s'il s'agit d'une diffusion élastique). Par conséquent, l'amplitude f des rayons X diffusés par les atomes est $f = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ Ce sera. Cette équation équivaut à la transformation de Fourier de la densité électronique $ \ rho (\ vec {r}) $ avec le vecteur de diffusion $ \ vec {k} $. Un vecteur de diffusion est un vecteur défini par des composantes (h, k, l) dans l'espace inverse. Ce f est appelé le facteur de diffusion atomique. La densité électronique est nécessaire pour obtenir le facteur de diffusion atomique, mais comme la densité électronique autour de l'atome ne peut être obtenue analytiquement que dans des cas simples, elle est approximativement calculée par la méthode de Hartley-Fock ou similaire. Puisque les facteurs de diffusion atomique de chaque atome ont déjà été obtenus et stockés dans une base de données, nous les utiliserons dans ce calcul.

De plus, la même idée peut être utilisée pour tout le cristal. Le facteur de structure cristalline F peut être obtenu en transformant de Fourier la densité électronique dans le cristal. $F(\vec{k}) = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$ La formule est la même que celle du facteur de diffusion atomique. Cependant, dans le cristal, les électrons sont distribués uniquement autour de la position atomique, ce qui peut être réécrit par les coordonnées atomiques (x, y, z) et le facteur de diffusion atomique f. $F(hkl) = \sum_{j}{f_jT_je^{2πi(hx_j + ky_j + lz_j)}}$ Ce sera. h, k, l sont des composantes du vecteur de diffusion k et sont également les coordonnées des points de grille (points de grille inversés) dans l'espace inverse. $ T_j $, que j'expliquerai plus tard, est le facteur d'appareil-waller.

programme

Il s'agit de la méthode de mise en œuvre réelle. Commencez par lister les points de grille inverses hkl. À l'origine, seuls les points de la sphère d'Ewald devraient être répertoriés, mais comme une plage assez large est couverte sur la sphère d'Ewald à la longueur d'onde des rayons X, j'en énumérerai autant que possible. Hkl.csv sur GitHub est un fichier hkl énuméré. Vous aurez également besoin des coordonnées atomiques de la cellule unitaire du cristal. Ici, seuls les atomes de la maille élémentaire sont suffisants. Cette fois, les coordonnées du réseau cubique face à centre (fcc) sont utilisées comme exemple. Dans GitHub, il s'agit d'un fichier nommé pos.csv.

Ensuite, à partir des constantes de réseau a, b, c, $ \ alpha $, $ \ beta $, $ \ gamma $ (que nous avons à l'avance comme paramètres), nous assemblons le tenseur de réseau $ G $ dans l'espace réel. $ 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) $ prochain, $G^{\ast} = G^{-1}$ Par conséquent, en prenant la matrice inverse de ce tenseur de réseau, le tenseur de réseau $ G ^ {\ ast} $ dans l'espace inverse peut être obtenu.

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)

En utilisant ce tenseur de réseau inverse $ G ^ {\ ast} , la longueur du vecteur de réseau inverse par rapport à la réflexion K au point de réseau inverse hkl $|\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}|$ Demander. $ |\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|A partir de la condition de Bragg, l'espacement de surface dk à la réflexion K et l'angle de Bragg θk peuvent être obtenus. $d_k = \frac{1}{|\vec{r}|^2}$ $\theta_k = \arcsin(\frac{λ}{2d_k})$

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

Les valeurs de la base de données sont utilisées pour le facteur de diffusion atomique, mais comme le facteur de diffusion atomique est une fonction qui dépend de sinθ / λ, seul le coefficient est répertorié dans la base de données. Pour assembler à partir de coefficients $f(\frac{\sin\theta}{\lambda}) = \sum_{j=1}^{4}{a_j e^{-b_j(\frac{\sin\theta}{λ})^2}} + c$ ça ira. $ a_j $, $ b_j $, c sont les coefficients. Dans le programme, nous avons supposé Na atome et utilisé les données.

Maintenant, nous devons expliquer le facteur de waller. Le facteur d'appareil-waller est un facteur qui compense les effets des vibrations thermiques. Le facteur debai-waller pour les vibrations thermiques isotropes $T = e^{-B(\frac{\sin\theta}{\lambda})^2}$ Il est représenté par. B est le paramètre de déplacement atomique isotrope. Si les vibrations thermiques ne sont pas isotropes (dites anisotropes), elles sont plus complexes, en utilisant la matrice symétrique $ 3 \ times3 $ $ \ beta . $T = e^{-h^{\mathrm{T}} β h}$$ Ce sera. Cette fois, compte tenu des vibrations thermiques isotropes, B = 2,0. J'ai recherché la valeur du paramètre de déplacement atomique, mais je ne l'ai pas trouvée, c'est donc une valeur appropriée. Excusez-moi.

Le facteur de structure cristalline est calculé sur la base du facteur de diffusion atomique et des coordonnées atomiques de la cellule unitaire du cristal. $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)}}$ Comme vous pouvez le voir, le facteur de structure cristalline est un nombre complexe.

Si vous tracez le carré de la valeur absolue du facteur de structure cristalline obtenu (puisque c'est le carré de la valeur absolue de l'amplitude de diffusion, ce sera l'intensité de diffusion) sur l'axe vertical et $ 2 \ theta $ sur l'axe horizontal, un graphique qui peut être obtenu par diffraction des rayons X Il peut être obtenu. Puisque l'angle du rayon X incident est $ \ theta $ et l'angle du rayon X diffusé est $ \ theta $, il est habituel de prendre l'axe horizontal à $ 2 \ theta $.

résultat

L'intrigue ressemble à ceci:

L'angle de Bragg est arrondi à 2 chiffres. Dans ce graphique, l'intensité de diffusion apparaît dans une fonction delta uniquement à l'angle de réflexion qui satisfait la condition de Bragg, mais en diffraction réelle des rayons X, le pic de l'intensité de diffusion a une largeur. Reproduisons l'étalement du pic.

Lorsque vous voulez donner un spread à la fonction delta, la première chose qui vous vient à l'esprit est la convolution par la fonction gaussienne (gaussienne). C'est la soi-disant convolution gaussienne. Ce qui suit est une convolution de la fonction gaussienne dans le graphique précédent. La dispersion a une valeur qui lui ressemble. Vous pouvez voir qu'il y a un peu de largeur à la base. Comme le sait tous ceux qui ont vu les données de mesure de la diffraction des rayons X, ce graphique semble avoir une base légèrement plus étroite que la base réelle. En fait, la fonction gaussienne est adaptée pour exprimer un comportement tel que la vibration thermique de particules, mais elle ne convient pas pour exprimer le pic d'un tel spectre. Par conséquent, convoluons avec la fonction de Lorentz utilisée pour exprimer le pic du spectre. La fonction Lorentz est la suivante. $L = \frac{1}{1 + (\frac{x}{γ})^2}$ $ \ Gamma $ est appelé la largeur à moitié prix, qui est la largeur du pic à 1/2 de la hauteur du pic. Ce paramètre modifie la façon dont la fonction se propage.

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

Celui qui a été plié est le suivant. Encore une fois, la valeur de la fourchette de prix à moitié est exactement comme ça. Elle a une base plus large que la fonction gaussienne précédente. Cela ressemble au pic de diffraction des rayons X. Cela dit, la fonction de Lorentz seule n'est pas précise. Parce que les atomes du cristal vibrent avec la température (il n'y a aucune vibration même au zéro absolu), les vibrations thermiques avec un comportement fonctionnel gaussien affectent également l'étalement des pics. Par conséquent, la fonction de Lorentz et la fonction de Gauss peuvent être pliées avec un poids arbitraire pour exprimer le pic plus précisément.

La convolution de la fonction de Lorenz et de la fonction de Gauss est appelée fonction de Vogt. Les fonctions Vogt sont plus rapides à créer à l'aide de fonctions d'erreur complexes qu'à simplement convoluer des fonctions de Lorentz et de Gauss. Je ne suis pas sûr en raison du manque d'étude, mais il semble que la partie réelle de la fonction d'erreur complexe ressemble à une fonction Voigt. [^ 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

Ce qui suit est une fonction de Vogt qui exprime l'étalement du pic. C'est bien de l'essayer, mais il semble que le spread ne s'exprime que par la fonction Lorenz. La fonction Voigt nécessite une fourchette à moitié prix pour chacune des parties fonction Gauss et fonction Lorentz en tant que paramètres, mais cette fois, elle est définie de manière appropriée. En approchant théoriquement la moitié de la fourchette de prix, il peut être plus proche de l'écart réel du pic.

Image de diffraction

En diffraction des rayons X, en plus du tracé ci-dessus, une image appelée image de diffraction, qui est une copie de la diffraction des rayons X telle quelle, peut être obtenue. Puisque le réseau inverse peut être vu dans l'image de diffraction, je le ferai en prime cette fois.

En tant qu'implémentation, numpy de toute taille.zeros()Uniquement lorsque la condition de réflexion est satisfaite en générant un tableau de|F|^Entrez une valeur de 2 et pliez la fonction gaussienne pour la faire ressembler à un point.

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

Dans l'image de diffraction, seuls les points inverses du réseau perpendiculaires à la direction du rayon X incident apparaissent sous forme de points. Par conséquent, en plus des conditions existantes, la direction incidente des rayons X doit également être prise en compte. Avec la direction incidente des rayons X comme [001], seules les réflexions qui prennent le produit interne et deviennent 0 sont extraites. Ce processus est facile s'il est [001], mais si la direction d'incident est [110] ou [111], un traitement compliqué sera nécessaire. En Python, vous pouvez utiliser matplotlib.pyplot.imshow () pour créer une image telle qu'elle est, alors faites-en une image avec un contraste noir et blanc.

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)

C'est pratique. L'image dans fcc est comme ci-dessous.

Ça ressemble à ça. Je n'ai dit ce genre de chose que depuis un moment. Je l'ai également essayé avec un réseau cubique centré sur le corps (cci). C'est plus amusant que fcc.

Résumé

J'ai essayé de simuler la diffraction des rayons X. Bien que la théorie de la diffraction soit compliquée, le traitement lui-même par un ordinateur est facile, il peut donc être un bon matériel pédagogique pour comprendre la diffraction des rayons X. A vrai dire, pour réaliser une simulation précise, il faut considérer plus de facteurs que ceux mis en œuvre cette fois (Lorentz / facteur de polarisation, fonction d'orientation sélective, etc.). De plus, les paramètres définis de manière appropriée cette fois (paramètre de déplacement atomique, largeur de demi-valeur de la fonction Vogt) doivent également utiliser des valeurs précises. Alors, ces paramètres peuvent-ils être déterminés théoriquement? Ces paramètres sont constitués d'une combinaison complexe de divers facteurs, donc même si vous pouvez déduire des facteurs dominants, il sera difficile de les déterminer analytiquement. Afin d'obtenir des valeurs précises pour les paramètres ci-dessus, nous utilisons une méthode de mise en correspondance des valeurs expérimentales avec des simulations pour réduire les erreurs. Donc, l'article précédent [^ 3]

Les références

Cet article est basé sur les articles et livres suivants. Je vais le présenter sans permission.

--Pratique de l'analyse aux rayons X sur poudre Introduction à la méthode Reetbelt (édité par la conférence de la table ronde de recherche sur l'analyse aux rayons X de la Japan Society for Analytical Chemistry, édité par Izumi Nakai, édité par Fujio Izumi)

• Comment tracer la fonction Voigt en python

code

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