Calculate Profile Data Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how to utilise the calculate_profile_data method for x-ray diffraction simulations of different crystal structures. The structes demonstrated in this notebook are:\n",
    "    -Fm-3d Silicon\n",
    "    -Bcc Iron\n",
    "    -Fcc Nickel\n",
    "    -Graphite \n",
    "    -Arbitray simple cubic (with Ni atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This functionaility has been checked to run in pyxem-0.11.0 (May 2020). Bugs are always possible, do not trust the code blindly, and if you experience any issues please report them here: https://github.com/pyxem/pyxem-demos/issues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. <a href='#tem'> Importing generators</a>\n",
    "2. <a href='#vec'> Define Structures</a>\n",
    "3. <a href='#vec'> Define the diffraction generator class</a>\n",
    "4. <a href='#vec'> Plotting\n",
    "5. <a href='#vec'> Index of methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='tem'></a> 1. Import Generators\n",
    "\n",
    "Import generators required for simulation and indexation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import diffpy\n",
    "import numpy as np\n",
    "from math import pi\n",
    "from transforms3d.euler import euler2mat\n",
    "import collections\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from diffsims.sims.diffraction_simulation import DiffractionSimulation\n",
    "\n",
    "from diffsims.utils.atomic_scattering_params import ATOMIC_SCATTERING_PARAMS\n",
    "from diffsims.utils.lobato_scattering_params import ATOMIC_SCATTERING_PARAMS_LOBATO\n",
    "from diffsims.generators.diffraction_generator import DiffractionGenerator\n",
    "from diffsims.sims.diffraction_simulation import  DiffractionSimulation\n",
    "\n",
    "from diffsims.utils.atomic_diffraction_generator_support.fourier_transform import (\n",
    "    from_recip,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='tem'></a> 2. Define structures\n",
    "\n",
    "Define or import the structures to be used in this simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_si_structure(lattice_parameter=None):\n",
    "    \"\"\"\n",
    "    We construct an Fd-3m silicon (with lattice parameter 5.431 as a default)\n",
    "    \"\"\"\n",
    "    if lattice_parameter is not None:\n",
    "        a = lattice_parameter\n",
    "    else:\n",
    "        a = 5.431\n",
    "    latt = diffpy.structure.lattice.Lattice(a, a, a, 90, 90, 90)\n",
    "    atom_list = []\n",
    "    for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:\n",
    "        x, y, z = coords[0], coords[1], coords[2]\n",
    "        atom_list.append(\n",
    "            diffpy.structure.atom.Atom(atype=\"Si\", xyz=[x, y, z], lattice=latt)\n",
    "        )  # Motif part A\n",
    "        atom_list.append(\n",
    "            diffpy.structure.atom.Atom(\n",
    "                atype=\"Si\", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt\n",
    "            )\n",
    "        )  # Motif part B\n",
    "    return diffpy.structure.Structure(atoms=atom_list, lattice=latt)\n",
    "\n",
    "def make_fe_structure():\n",
    "    a = 2.856\n",
    "    latt = diffpy.structure.lattice.Lattice(a, a, a, 90, 90, 90)\n",
    "    atom_list = []\n",
    "    for coords in [[0, 0, 0], [0.5, 0.5, 0.5]]:\n",
    "        x, y, z = coords[0], coords[1], coords[2]\n",
    "        atom_list.append(\n",
    "            diffpy.structure.atom.Atom(atype=\"Fe\", xyz=[x, y, z], lattice=latt)\n",
    "        )\n",
    "    return diffpy.structure.Structure(atoms=atom_list, lattice=latt)\n",
    "\n",
    "def make_ni_structure():\n",
    "    a = 3.499\n",
    "    latt = diffpy.structure.lattice.Lattice(a, a, a, 90, 90, 90)\n",
    "    atom_list = []\n",
    "    for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:\n",
    "        x, y, z = coords[0], coords[1], coords[2]\n",
    "        atom_list.append(\n",
    "            diffpy.structure.atom.Atom(atype=\"Ni\", xyz=[x, y, z], lattice=latt)\n",
    "        )\n",
    "    return diffpy.structure.Structure(atoms=atom_list, lattice=latt)\n",
    "\n",
    "def make_graphite_structure():\n",
    "    a = 2.461\n",
    "    c = 6.708\n",
    "    latt = diffpy.structure.lattice.Lattice(a, a, c, 90, 90, 120)\n",
    "    atom_list = []\n",
    "    coords = [0, 0, 0]\n",
    "    x, y, z = coords[0], coords[1], coords[2]\n",
    "    atom_list.append(\n",
    "        diffpy.structure.atom.Atom(atype=\"C\", xyz=[x, y, z], lattice=latt)\n",
    "    )\n",
    "    return diffpy.structure.Structure(atoms=atom_list, lattice=latt)\n",
    "\n",
    "def make_cubic_structure():\n",
    "    a = 3\n",
    "    latt = diffpy.structure.lattice.Lattice(a, a, a, 90, 90, 90)\n",
    "    atom_list = []\n",
    "    coords = [0, 0, 0]\n",
    "    x, y, z = coords[0], coords[1], coords[2]\n",
    "    atom_list.append(\n",
    "        diffpy.structure.atom.Atom(atype=\"Ni\", xyz=[x, y, z], lattice=latt)\n",
    "    )\n",
    "    return diffpy.structure.Structure(atoms=atom_list, lattice=latt)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Si_struct = make_si_structure()\n",
    "Fe_struct = make_fe_structure()\n",
    "Ni_struct = make_ni_structure()\n",
    "graphite_struct = make_graphite_structure()\n",
    "cubic_struct = make_cubic_structure()\n",
    "wz_struct = diffpy.structure.loadStructure('/Users/Wood/Documents/GitHub/cif_files/GaAs_mp-8883_conventional_standard.cif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='veca'></a> 3. Generate diffraction data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the calcualte_profile_data method from the DiffractionGenerator class \n",
    "\n",
    "Once defined can use the caluclate_profile_data method to produce the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_gen = DiffractionGenerator(accelerating_voltage=300,\n",
    "                                max_excitation_error=1/10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "si_prof = diff_gen.calculate_profile_data(\n",
    "                        Si_struct,\n",
    "                        reciprocal_radius=1.0,\n",
    "                        magnitude_tolerance=1e-5,\n",
    "                        minimum_intensity=1e-3,\n",
    "                        )\n",
    "\n",
    "fe_prof = diff_gen.calculate_profile_data(\n",
    "                        Fe_struct,\n",
    "                        reciprocal_radius=1.0,\n",
    "                        magnitude_tolerance=1e-5,\n",
    "                        minimum_intensity=1e-3,\n",
    "                        )\n",
    "\n",
    "ni_prof = diff_gen.calculate_profile_data(\n",
    "                        Ni_struct,\n",
    "                        reciprocal_radius=1.0,\n",
    "                        magnitude_tolerance=1e-5,\n",
    "                        minimum_intensity=1e-3,\n",
    "                        )\n",
    "\n",
    "graphite_prof = diff_gen.calculate_profile_data(\n",
    "                        graphite_struct,\n",
    "                        reciprocal_radius=1.0,\n",
    "                        magnitude_tolerance=1e-5,\n",
    "                        minimum_intensity=1e-3,\n",
    "                        )\n",
    "\n",
    "cubic_prof = diff_gen.calculate_profile_data(\n",
    "                        cubic_struct,\n",
    "                        reciprocal_radius=1.0,\n",
    "                        magnitude_tolerance=1e-5,\n",
    "                        minimum_intensity=1e-3,\n",
    "                        )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='veca'></a> 4. Plotting\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the get_plot method from the DiffractionSimulation class to visualise the data produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "si_plot = si_prof.get_plot(5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fe_plot = fe_prof.get_plot(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ni_plot = ni_prof.get_plot(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grpahite_plot = graphite_prof.get_plot(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cubic_plot = cubic_prof.get_plot(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='veca'></a> 5. Index of methods\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Libary of the importnant methods used in the above simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_intesnities_params_dict(reciprocal_lattice, reciprocal_radius):\n",
    "\n",
    "    spot_indices, _, spot_distances = get_points_in_sphere(reciprocal_lattice, reciprocal_radius)\n",
    "    \n",
    "    dict_i_to_d = {}\n",
    "    for i,d in zip(spot_indices, spot_distances):\n",
    "        dict_i_to_d[tuple(i)] = d\n",
    "    \n",
    "    list_hkls = spot_indices.tolist()\n",
    "    \n",
    "    unique_hkls_dict = get_unique_families(list_hkls)\n",
    "    \n",
    "    multiplicites = np.fromiter(unique_hkls_dict.values(), dtype=float)\n",
    "    unique_hkls1 = np.array(list(unique_hkls_dict))\n",
    "\n",
    "    g_hkls = []\n",
    "    for unique_hkl in unique_hkls1:\n",
    "        g_hkls.append(dict_i_to_d[tuple(unique_hkl)])\n",
    "        \n",
    "    return unique_hkls1, multiplicites, g_hkls\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kinematical_intensities(\n",
    "    structure,\n",
    "    g_indices,\n",
    "    g_hkls_array,\n",
    "    debye_waller_factors,\n",
    "    multiplicites,\n",
    "    scattering_params,\n",
    "    shape_factor,\n",
    "):\n",
    "\n",
    "    \"\"\"Calculates peak intensities.\n",
    "\n",
    "    The peak intensity is a combination of the structure factor for a given\n",
    "    peak and the position the Ewald sphere intersects the relrod. In this\n",
    "    implementation, the intensity scales linearly with proximity.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    structure : Structure\n",
    "        The structure for which to derive the structure factors.\n",
    "    indices : array-like\n",
    "        The fractional coordinates of the peaks for which to calculate the\n",
    "        structure factor.\n",
    "    proximities : array-like\n",
    "        The distances between the Ewald sphere and the peak centers.\n",
    "    shape_factor_use : bool\n",
    "        True if the shape factor correction needs to be used\n",
    "    Returns\n",
    "    -------\n",
    "    peak_intensities : array-like\n",
    "        The intensities of the peaks.\n",
    "\n",
    "    \"\"\"\n",
    "    (\n",
    "        coeffs,\n",
    "        fcoords,\n",
    "        occus,\n",
    "        dwfactors,\n",
    "    ) = get_vectorized_list_for_atomic_scattering_factors(\n",
    "        structure=structure,\n",
    "        debye_waller_factors=debye_waller_factors,\n",
    "        scattering_params=scattering_params,\n",
    "    )\n",
    "\n",
    "    # Store array of g_hkls^2 values since used multiple times.\n",
    "\n",
    "    ##length of the unique hkls\n",
    "    g_hkls_sq = g_hkls_array ** 2\n",
    "\n",
    "    # Create array containing atomic scattering factors.\n",
    "    fs = get_atomic_scattering_factors(g_hkls_sq, coeffs, scattering_params)\n",
    "\n",
    "    # Change the coordinate system of fcoords to align with that of g_indices\n",
    "    fcoords = np.dot(\n",
    "        fcoords,\n",
    "        np.linalg.inv(np.dot(structure.lattice.stdbase, structure.lattice.recbase)),\n",
    "    )\n",
    "\n",
    "    # Calculate structure factors for all excited g-vectors.\n",
    "    # dosnet like doing the dot on g_indices since is a dict keys not just array--qq usingunique_hkls here test!\n",
    "    f_hkls = np.sum(\n",
    "        fs\n",
    "        * occus\n",
    "        * np.exp(\n",
    "            2j * np.pi * np.dot(g_indices, fcoords.T)\n",
    "            - 0.25 * np.outer(g_hkls_sq, dwfactors)\n",
    "        ),\n",
    "        axis=-1,\n",
    "    )\n",
    "\n",
    "    # Define an intensity scaling that is linear with distance from Ewald sphere\n",
    "    # along the beam direction.\n",
    "\n",
    "    prefactor = shape_factor * multiplicites\n",
    "\n",
    "    # Calculate the peak intensities from the structure factor and excitation\n",
    "    peak_intensities = prefactor * (f_hkls * f_hkls.conjugate()).real\n",
    "    return peak_intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_profile_data(\n",
    "        self,\n",
    "        structure,\n",
    "        reciprocal_radius=1.0,\n",
    "        magnitude_tolerance=1e-5,\n",
    "        minimum_intensity=1e-3,\n",
    "    ):\n",
    "        \"\"\"\n",
    "        Calculates a one dimensional diffraction profile for a structure.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        structure : Structure\n",
    "            The structure for which to calculate the diffraction profile.\n",
    "        reciprocal_radius : float\n",
    "            The maximum radius of the sphere of reciprocal space to sample, in\n",
    "            reciprocal angstroms.\n",
    "        magnitude_tolerance : float\n",
    "            The minimum difference between diffraction magnitudes in reciprocal\n",
    "            angstroms for two peaks to be consdiered different.\n",
    "        minimum_intensity : float\n",
    "            The minimum intensity required for a diffraction peak to be\n",
    "            considered real. Deals with numerical precision issues.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        diffsims.ProfileSimulation\n",
    "            The diffraction profile corresponding to this structure and\n",
    "            experimental conditions.\n",
    "        \"\"\"\n",
    "        max_r = reciprocal_radius\n",
    "        wavelength = self.wavelength\n",
    "        scattering_params = self.scattering_params\n",
    "\n",
    "        latt = structure.lattice\n",
    "        is_hex = is_lattice_hexagonal(latt)\n",
    "\n",
    "        # Obtain crystallographic reciprocal lattice points within range\n",
    "        recip_latt = latt.reciprocal()\n",
    "        spot_indices, _, spot_distances = get_points_in_sphere(\n",
    "            recip_latt, reciprocal_radius\n",
    "        )\n",
    "\n",
    "        ##spot_indicies is a numpy.array of the hkls allowd in the recip radius\n",
    "        unique_hkls, multiplicites, g_hkls = get_intensities_params(\n",
    "            recip_latt, reciprocal_radius\n",
    "        )\n",
    "        g_indices = unique_hkls\n",
    "        debye_waller_factors = self.debye_waller_factors\n",
    "        excitation_error = None\n",
    "        max_excitation_error = None\n",
    "        g_hkls_array = np.asarray(g_hkls)\n",
    "\n",
    "        i_hkl = get_kinematical_intensities(\n",
    "            structure,\n",
    "            g_indices,\n",
    "            g_hkls_array,\n",
    "            debye_waller_factors,\n",
    "            multiplicites,\n",
    "            scattering_params,\n",
    "            shape_factor=1,\n",
    "        )\n",
    "\n",
    "        if is_hex:\n",
    "            # Use Miller-Bravais indices for hexagonal lattices.\n",
    "            g_indices = (g_indices[0], g_indices[1], -g_indices[0] - g_indices[1], g_indices[2])\n",
    "\n",
    "        hkls_labels = [\"\".join([str(int(x)) for x in xs]) for xs in unique_hkls]\n",
    "\n",
    "        peaks = {}\n",
    "        for l, i, g in zip(hkls_labels, i_hkl, g_hkls):\n",
    "            peaks[l] = [i, g]\n",
    "\n",
    "        # Scale intensities so that the max intensity is 100.\n",
    "\n",
    "        max_intensity = max([v[0] for v in peaks.values()])\n",
    "        x = []\n",
    "        y = []\n",
    "        hkls = []\n",
    "        for k in peaks.keys():\n",
    "            v = peaks[k]\n",
    "            if v[0] / max_intensity * 100 > minimum_intensity and (k != \"000\"):\n",
    "                x.append(v[1])\n",
    "                y.append(v[0])\n",
    "                hkls.append(k)\n",
    "\n",
    "        y = y / max(y) * 100\n",
    "\n",
    "        return ProfileSimulation(x, y, hkls)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_plot(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):\n",
    "\n",
    "        \"\"\"Plots the diffraction profile simulation for the\n",
    "           calculate_profile_data method in DiffractionGenerator.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        g_max : float\n",
    "            Maximum g-vector magnitude to plot.\n",
    "        annotate_peaks : boolean\n",
    "            If True, peaks are annotaed with hkl information.\n",
    "        with_labels : boolean\n",
    "            If True, xlabels and ylabels are added to the plot.\n",
    "        fontsize : integer\n",
    "            Fontsize for peak labels.\n",
    "        \"\"\"\n",
    "\n",
    "        ax = plt.gca()\n",
    "        for g, i, hkls in zip(self.magnitudes, self.intensities, self.hkls):\n",
    "            label = hkls\n",
    "            ax.plot([g, g], [0, i], color=\"k\", linewidth=3, label=label)\n",
    "            if annotate_peaks:\n",
    "                ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize)\n",
    "\n",
    "            if with_labels:\n",
    "                ax.set_xlabel(\"A ($^{-1}$)\")\n",
    "                ax.set_ylabel(\"Intensities (scaled)\")\n",
    "\n",
    "        return plt"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}