{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# This cell is added by sphinx-gallery\n\n%matplotlib inline\n\nimport mrsimulator\nprint(f'You are using mrsimulator v{mrsimulator.__version__}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Amorphous material, 29Si (I=1/2)\n\n29Si (I=1/2) simulation of amorphous material.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "One of the advantages of the ``mrsimulator`` package is that it is a fast NMR\nspectrum simulation library. We can exploit this feature to simulate bulk spectra and\neventually model amorphous materials. In this section, we illustrate how the\n``mrsimulator`` library may be used in simulating the NMR spectrum of amorphous\nmaterials.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom mrsimulator import Simulator\nfrom mrsimulator.methods import BlochDecaySpectrum\nfrom mrsimulator.utils.collection import single_site_system_generator\nfrom scipy.stats import multivariate_normal\n\n# global plot configuration\nmpl.rcParams[\"figure.figsize\"] = [4.5, 3.0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Generating tensor parameter distribution\n\nWe model the amorphous material by assuming a distribution of interaction tensors.\nFor example, a tri-variate normal distribution of the shielding tensor parameters,\n`i.e.`, the isotropic chemical shift, the anisotropy parameter, $\\zeta$, and\nthe asymmetry parameter, $\\eta$. In the following, we use pure NumPy and SciPy\nmethods to generate the three-dimensional distribution, as follows,\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean = [-100, 50, 0.15]  # given as [isotropic chemical shift in ppm, zeta in ppm, eta].\ncovariance = [[3.25, 0, 0], [0, 26.2, 0], [0, 0, 0.002]]  # same order as the mean.\n\n# range of coordinates along the three dimensions\niso_range = np.arange(100) * 0.3055 - 115  # in ppm\nzeta_range = np.arange(30) * 2.5 + 10  # in ppm\neta_range = np.arange(21) / 20\n\n# The coordinates grid\niso, zeta, eta = np.meshgrid(iso_range, zeta_range, eta_range, indexing=\"ij\")\npos = np.asarray([iso, zeta, eta]).T\n\n# Three-dimensional probability distribution function.\npdf = multivariate_normal(mean=mean, cov=covariance).pdf(pos).T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, ``iso``, ``zeta``, and ``eta`` are the isotropic chemical shift, nuclear\nshielding anisotropy, and nuclear shielding asymmetry coordinates of the 3D-grid\nsystem over which the multivariate normal probability distribution is evaluated. The\nmean of the distribution is given by the variable ``mean`` and holds a value of -100\nppm, 50 ppm, and 0.15 for the isotropic chemical shift, nuclear shielding anisotropy,\nand nuclear shielding asymmetry parameter, respectively. Similarly, the variable\n``covariance`` holds the covariance matrix of the multivariate normal distribution.\nThe two-dimensional projections from this three-dimensional distribution are shown\nbelow.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_, ax = plt.subplots(1, 3, figsize=(9, 3))\n\n# isotropic shift v.s. shielding anisotropy\nax[0].contourf(zeta_range, iso_range, pdf.sum(axis=2))\nax[0].set_xlabel(r\"shielding anisotropy, $\\zeta$ / ppm\")\nax[0].set_ylabel(\"isotropic chemical shift / ppm\")\n\n# isotropic shift v.s. shielding asymmetry\nax[1].contourf(eta_range, iso_range, pdf.sum(axis=1))\nax[1].set_xlabel(r\"shielding asymmetry, $\\eta$\")\nax[1].set_ylabel(\"isotropic chemical shift / ppm\")\n\n# shielding anisotropy v.s. shielding asymmetry\nax[2].contourf(eta_range, zeta_range, pdf.sum(axis=0))\nax[2].set_xlabel(r\"shielding asymmetry, $\\eta$\")\nax[2].set_ylabel(r\"shielding anisotropy, $\\zeta$ / ppm\")\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create the Simulator object\n\n**Spin system:**\nLet's create the sites and single-site spin system objects from these parameters.\nUse the :func:`~mrsimulator.utils.collection.single_site_system_generator` utility\nfunction to generate single-site spin systems.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "spin_systems = single_site_system_generator(\n    isotopes=\"29Si\",\n    isotropic_chemical_shifts=iso,\n    shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n    abundance=pdf,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, ``iso``, ``zeta``, and ``eta`` are the array of tensor parameter coordinates,\nand ``pdf`` is the array of corresponding amplitudes.\n\n\n**Method:**\nLet's also create the Bloch decay spectrum method.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "method = BlochDecaySpectrum(\n    channels=[\"29Si\"],\n    spectral_dimensions=[\n        {\"spectral_width\": 25000, \"reference_offset\": -7000}  # values in Hz\n    ],\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above method simulates a static $^{29}\\text{Si}$ spectrum at 9.4 T field\n(default value).\n\n**Simulator:**\nNow, that we have the spin systems and the method, create the simulator object and\nadd the respective objects.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim = Simulator()\nsim.spin_systems = spin_systems  # add the spin systems\nsim.methods += [method]  # add the method"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static spectrum\nObserve the static $^{29}\\text{Si}$ NMR spectrum simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The plot of the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = plt.subplot(projection=\"csdm\")\nax.plot(sim.methods[0].simulation, color=\"black\", linewidth=1)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>The broad spectrum seen in the above figure is a result of spectral averaging\n    of spectra arising from a distribution of shielding tensors. There is no\n    line-broadening filter applied to the spectrum.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Spinning sideband simulation at $90^\\circ$\nHere is an example of a sideband simulation, spinning at a 90-degree angle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim.methods[0] = BlochDecaySpectrum(\n    channels=[\"29Si\"],\n    rotor_frequency=5000,  # in Hz\n    rotor_angle=1.57079,  # in rads, equivalent to 90 deg.\n    spectral_dimensions=[\n        {\"spectral_width\": 25000, \"reference_offset\": -7000}  # values in Hz\n    ],\n)\nsim.config.number_of_sidebands = 8  # eight sidebands are sufficient for this example\nsim.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The plot of the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = plt.subplot(projection=\"csdm\")\nax.plot(sim.methods[0].simulation, color=\"black\", linewidth=1)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Spinning sideband simulation at the magic angle\nHere is another example of a sideband simulation at the magic angle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim.methods[0] = BlochDecaySpectrum(\n    channels=[\"29Si\"],\n    rotor_frequency=1000,  # in Hz\n    rotor_angle=54.735 * np.pi / 180.0,  # in rads\n    spectral_dimensions=[\n        {\"spectral_width\": 25000, \"reference_offset\": -7000}  # values in Hz\n    ],\n)\nsim.config.number_of_sidebands = 16  # sixteen sidebands are sufficient for this example\nsim.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The plot of the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = plt.subplot(projection=\"csdm\")\nax.plot(sim.methods[0].simulation, color=\"black\", linewidth=1)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()"
      ]
    }
  ],
  "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.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}