|
| 1 | +#!/usr/bin/env python |
| 2 | +# |
| 3 | +# Copyright (c) 2025 Authors and contributors |
| 4 | +# (see the AUTHORS.rst file for the full list of names) |
| 5 | +# |
| 6 | +# Released under the GNU Public Licence, v3 or any higher version |
| 7 | +# SPDX-License-Identifier: GPL-3.0-or-later |
| 8 | +""".. _howto-saxs: |
| 9 | +
|
| 10 | +Small-angle X-ray scattering |
| 11 | +============================ |
| 12 | +
|
| 13 | +Small-angle X-ray scattering (SAXS) can be extracted using Scatterkit. To follow this |
| 14 | +how-to guide, you should download the :download:`topology <water.tpr>` and the |
| 15 | +:download:`trajectory <water.trr>` files of the water system. |
| 16 | +
|
| 17 | +For more details on the theory see :ref:`saxs-explanations`. |
| 18 | +
|
| 19 | +First, we import Matplotlib, MDAnalysis, NumPy and Scatterkit: |
| 20 | +""" # noqa: D415 |
| 21 | +# %% |
| 22 | + |
| 23 | +import matplotlib.pyplot as plt |
| 24 | +import MDAnalysis as mda |
| 25 | +from MDAnalysis.analysis.rdf import InterRDF |
| 26 | + |
| 27 | +import scatterkit |
| 28 | +from scatterkit.lib.math import atomic_form_factor, rdf_structure_factor |
| 29 | + |
| 30 | +# %% |
| 31 | +# The `water` system consists of 510 water molecules in the liquid state. The |
| 32 | +# molecules are placed in a periodic cubic cell with an extent of :math:`25 \times 25 |
| 33 | +# \times 25\,\textrm{Å}^3`. |
| 34 | +# |
| 35 | +# Load Simulation Data |
| 36 | +# -------------------- |
| 37 | +# |
| 38 | +# Create a :class:`MDAnalysis.core.universe.Universe` and define a group containing only |
| 39 | +# the oxygen atoms and a group containing only the hydrogen atoms: |
| 40 | + |
| 41 | +u = mda.Universe("water.tpr", "water.trr") |
| 42 | + |
| 43 | +group_O = u.select_atoms("type O*") |
| 44 | +group_H = u.select_atoms("type H*") |
| 45 | + |
| 46 | +# %% |
| 47 | +# Extract small angle x-ray scattering (SAXS) intensities |
| 48 | +# ------------------------------------------------------- |
| 49 | +# |
| 50 | +# Let us use the :class:`scatterkit.Saxs` class of Scatterkit and apply it to all atoms |
| 51 | +# in the system: |
| 52 | + |
| 53 | +saxs = scatterkit.Saxs(u.atoms).run(stop=30) |
| 54 | + |
| 55 | +# %% |
| 56 | +# .. Note:: |
| 57 | +# SAXS computations are extensive calculations. Here, to get an overview of the |
| 58 | +# scattering intensities, we reduce the number of frames to be analyzed from ``101`` |
| 59 | +# to ``30``, by adding the ``stop = 30`` parameter to the ``run`` method. Due to the |
| 60 | +# small number of analyzed frames, the scattering intensities shown in this tutorial |
| 61 | +# should not be used to draw any conclusions from the data. |
| 62 | +# |
| 63 | +# Extract the scattering vectors and the averaged structure factor and SAXS scattering |
| 64 | +# intensities from the ``results`` attribute: |
| 65 | + |
| 66 | +scattering_values = saxs.results.scattering_vectors |
| 67 | +structure_factors = saxs.results.structure_factors |
| 68 | +scattering_intensities = saxs.results.scattering_intensities |
| 69 | + |
| 70 | +# %% |
| 71 | +# The scattering intensities (and structure factors) are given as a 1D array, let us |
| 72 | +# look at the 10 first lines: |
| 73 | + |
| 74 | +print(scattering_intensities[:10]) |
| 75 | + |
| 76 | +# %% |
| 77 | +# By default, the binwidth in the recipocal :math:`(q)` space is :math:`0.1 Å^{-1}`. |
| 78 | +# |
| 79 | +# We now plot the structure factor as well a the scattering intensities together. |
| 80 | + |
| 81 | +fig1, ax1 = plt.subplots(nrows=2, sharex=True, layout="constrained") |
| 82 | + |
| 83 | +ax1[0].plot(scattering_values, structure_factors) |
| 84 | +ax1[1].plot(scattering_values, scattering_intensities) |
| 85 | + |
| 86 | +ax1[-1].set_xlabel(r"q (1/Å)") |
| 87 | + |
| 88 | +ax1[0].set_ylabel(r"structure factor $S(q)$") |
| 89 | +ax1[1].set_ylabel(r"scattering intensities $I(q)$") |
| 90 | +fig1.align_labels() |
| 91 | + |
| 92 | +fig1.show() |
| 93 | + |
| 94 | + |
| 95 | +# %% |
| 96 | +# The structure factor :math:`S(q)` and the scattering intensities :math:`I(q)` are |
| 97 | +# related via |
| 98 | +# |
| 99 | +# .. math:: |
| 100 | +# I(q) = [f(q)]^2 S(q) |
| 101 | +# |
| 102 | +# where :math:`f(q)` are the atomic form factors. We will investigate the relation below |
| 103 | +# in more details. |
| 104 | +# |
| 105 | +# Computing oxygen and hydrogen contributions |
| 106 | +# ------------------------------------------- |
| 107 | +# |
| 108 | +# An advantage of full atomistic simulations is their ability to investigate atomic |
| 109 | +# contributions individually. Let us calculate both oxygen and hydrogen contributions, |
| 110 | +# respectively: |
| 111 | + |
| 112 | +saxs_O = scatterkit.Saxs(group_O).run(stop=30) |
| 113 | +saxs_H = scatterkit.Saxs(group_H).run(stop=30) |
| 114 | + |
| 115 | +# %% |
| 116 | +# Let us plot the results for the structure factor, the squared atomic form factor as |
| 117 | +# well scattering intensities together. For computing the atomic form factor we will use |
| 118 | +# :func:`scatterkit.lib.math.atomic_form_factor`. Note that for the |
| 119 | +# ``structure_factors`` and the ``scattering_intensities` we access the results directly |
| 120 | +# from the ``results`` attribute without storing them in individual variables as before: |
| 121 | + |
| 122 | +fig2, ax2 = plt.subplots(nrows=3, sharex=True, layout="constrained") |
| 123 | + |
| 124 | +# structure factors |
| 125 | +ax2[0].plot( |
| 126 | + saxs_O.results.scattering_vectors, |
| 127 | + saxs_O.results.structure_factors, |
| 128 | + label="Oxygen", |
| 129 | +) |
| 130 | +ax2[0].plot( |
| 131 | + saxs_H.results.scattering_vectors, |
| 132 | + saxs_H.results.structure_factors, |
| 133 | + label="Hydrogen", |
| 134 | +) |
| 135 | + |
| 136 | +# atomic form factors |
| 137 | +ax2[1].plot( |
| 138 | + saxs_O.results.scattering_vectors, |
| 139 | + atomic_form_factor(saxs_O.results.scattering_vectors, "O") ** 2, |
| 140 | +) |
| 141 | +ax2[1].plot( |
| 142 | + saxs_H.results.scattering_vectors, |
| 143 | + atomic_form_factor(saxs_H.results.scattering_vectors, "H") ** 2, |
| 144 | +) |
| 145 | + |
| 146 | +# scattering intensities |
| 147 | +ax2[2].plot(saxs_O.results.scattering_vectors, saxs_O.results.scattering_intensities) |
| 148 | +ax2[2].plot(saxs_H.results.scattering_vectors, saxs_H.results.scattering_intensities) |
| 149 | + |
| 150 | +ax2[-1].set_xlabel(r"q (1/Å)") |
| 151 | +ax2[0].set_ylabel(r"$S(q)$") |
| 152 | +ax2[1].set_ylabel(r"$f(q)^2$") |
| 153 | +ax2[2].set_ylabel(r"$I(q)$") |
| 154 | + |
| 155 | +ax2[0].legend() |
| 156 | +fig2.align_labels() |
| 157 | + |
| 158 | +fig2.show() |
| 159 | + |
| 160 | +# %% |
| 161 | +# The figure above nicely shows that multiplying the structure factor :math:`S(q)` and |
| 162 | +# the squared atomic form factor :math:`f(q)^2` results in the scattering intensity |
| 163 | +# :math:`I(q)`. |
| 164 | +# |
| 165 | +# The atomic form factors are monotonically decaying as a function of :math:`q` and |
| 166 | +# their value for :math:`q=0` is the same number of electrons for the element. Also, it |
| 167 | +# is worth to notice that due to small atomic form factor of hydrogen there is basically |
| 168 | +# no contribution of the hydrogen atoms to the total scattering intensity of water. |
| 169 | +# |
| 170 | +# Connection of the structure factor to the radial distribution function |
| 171 | +# ---------------------------------------------------------------------- |
| 172 | +# |
| 173 | +# As in details explained in :ref:`saxs-explanations`, the structure factor can be |
| 174 | +# related to the radial distribution function (RDF). We denote this structure factor by |
| 175 | +# :math:`S^\mathrm{FT}(q)` since it is based on Fourier transforming the RDF. The |
| 176 | +# structure factor which can be directly obtained from the trajectory is denoted by |
| 177 | +# :math:`S^\mathrm{D}(q)`. |
| 178 | +# |
| 179 | +# To relate these two we first calculate the oxygen-oxygen RDF up to half the box length |
| 180 | +# using :class:`MDAnalysis.analysis.rdf.InterRDF` and save the result in |
| 181 | +# variables for an easier access. |
| 182 | + |
| 183 | +box_lengh = u.dimensions[0] |
| 184 | + |
| 185 | +oo_inter_rdf = InterRDF( |
| 186 | + g1=group_O, g2=group_O, range=(0, box_lengh / 2), exclude_same="residue" |
| 187 | +).run() |
| 188 | + |
| 189 | +r_oo = oo_inter_rdf.results.bins |
| 190 | +rdf_oo = oo_inter_rdf.results.rdf |
| 191 | + |
| 192 | +# %% |
| 193 | +# We use ``exclude_same="residue"`` to exclude atomic self contributions resulting in a |
| 194 | +# large peak at 0. Next, we convert the RDF into a structure factor using |
| 195 | +# :func:`scatterkit.lib.math.rdf_structure_factor` and the number density of the |
| 196 | +# oxygens. |
| 197 | + |
| 198 | +density = group_O.n_atoms / u.trajectory.ts.volume |
| 199 | + |
| 200 | +q_rdf, struct_factor_rdf = rdf_structure_factor(rdf=rdf_oo, r=r_oo, density=density) |
| 201 | + |
| 202 | +# %% |
| 203 | +# Now we can plot everything together and find that the direct evaluation from above and |
| 204 | +# the transformed RDF give the same structure factor. |
| 205 | + |
| 206 | +fig3, ax3 = plt.subplots(2, layout="constrained") |
| 207 | + |
| 208 | +ax3[0].axhline(1, c="gray", ls="dashed") |
| 209 | +ax3[0].plot(r_oo, rdf_oo, label="Oxygen-Oxygen") |
| 210 | +ax3[0].set_xlabel("r (Å)") |
| 211 | +ax3[0].set_ylabel("radial distribution function") |
| 212 | +ax3[0].set_xlim(0, 10) |
| 213 | + |
| 214 | +ax3[1].plot(q_rdf, struct_factor_rdf, label=r"$S^\mathrm{FT}$") |
| 215 | +ax3[1].plot( |
| 216 | + saxs_O.results.scattering_vectors, |
| 217 | + saxs_O.results.structure_factors, |
| 218 | + label=r"$S^\mathrm{D}$", |
| 219 | + ls="dashed", |
| 220 | +) |
| 221 | + |
| 222 | +ax3[1].set_xlabel("q (1/Å)") |
| 223 | +ax3[1].set_ylabel("structure factor $S(q)$") |
| 224 | +ax3[1].set_xlim(0, 7) |
| 225 | + |
| 226 | +ax3[1].legend() |
| 227 | +ax3[0].legend() |
| 228 | +fig3.align_labels() |
| 229 | + |
| 230 | +fig3.show() |
0 commit comments