02-Function-Calculators.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function Calculators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [DiffPy-CMI](http://www.diffpy.org/products/diffpycmi/) framework the [diffpy.srreal](https://github.com/diffpy/diffpy.srreal) package provides the majority of features for evaluation of structure models and for simulation of structure-based physical quantities.  The diffpy.srreal calculators address physical quantities deriving from atom-pair interactions.  The calculators use the same procedure for iteration over atom pairs and a uniform interface to structure data.  As a result, the package can be easily extended with user-defined calculators or with support for additional structure representations.\n",
    "\n",
    "The calculators available in diffpy.srreal include:\n",
    "\n",
    "- BondCalculator - evaluation of bond distances\n",
    "- BVSCalculator - evaluation of bond valence sums\n",
    "- OverlapCalculator - calculation of atom overlaps asumming hard-sphere model\n",
    "- PDFCalculator - calculation of PDF (atomic pair distribution function) using real-space summation\n",
    "- DebyePDFCalculator - simulation of PDF from Debye scattering equation\n",
    "\n",
    "This notebook gives an overview of the use and extensibility of the calculators in diffpy.srreal.\n",
    "\n",
    "**Contents**\n",
    "\n",
    "> [Calculation of bond distances](#Calculation-of-bond-distances)<br>\n",
    "> [Basic and optimized evaluation](#Basic-and-optimized-evaluation)<br>\n",
    "> [Parallel calculation](#Parallel-calculation)<br>\n",
    "> [PDF calculators](#PDF-calculators)<br>\n",
    "> [Calculator components and their customization](#Calculator-components-and-their-customization)<br>\n",
    "> [Extending with more structure types](#Extending-with-more-structure-types)<br>\n",
    "> [User-defined Lennard-Jones calculator](#User-defined-Lennard-Jones-calculator)<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of bond distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# notebook setup\n",
    "%matplotlib notebook\n",
    "from __future__ import print_function\n",
    "from matplotlib.pyplot import subplots, rc\n",
    "rc('figure', figsize=(6, 4.5))\n",
    "\n",
    "# use several structure representations in DiffPy-CMI\n",
    "from diffpy.Structure import loadStructure\n",
    "from pyobjcryst.crystal import CreateCrystalFromCIF\n",
    "# load NaCl CIF as diffpy.Structure and pyobjcryst.Crystal object\n",
    "nacl_stru = loadStructure('NaCl.cif')\n",
    "nacl_objc = CreateCrystalFromCIF(open('NaCl.cif'))\n",
    "\n",
    "# create BondCalculator object\n",
    "from diffpy.srreal.bondcalculator import BondCalculator\n",
    "bc = BondCalculator(rmax=3)\n",
    "\n",
    "# The bc calculator works with both structure representations\n",
    "print(bc(nacl_stru))\n",
    "print(bc(nacl_objc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned distance arrays are of different lengths because *diffpy.Structure* works with an expanded unit cell.  The Python variable `bc` is an instance of the *BondCalculator* class.  There may be multiple instances of *BondCalculator* or any other calculator object, each with its own configuration and structure models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create two BondCalculators for Na-Na and Na-Cl pairs only\n",
    "bcNaNa = BondCalculator(rmax=4.5)\n",
    "bcNaNa.setTypeMask('Na1+', 'Na1+', True, others=False)\n",
    "bcNaCl = BondCalculator(rmax=4.5)\n",
    "bcNaCl.setTypeMask('Na1+', 'Cl1-', True, others=False)\n",
    "\n",
    "print(\"Na-Na:\", bcNaNa(nacl_objc))\n",
    "print(\"Na-Cl:\", bcNaCl(nacl_objc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *BondCalculator* object behaves like a function which takes a structure model as its argument and returns a sorted array of atom distances.  Additional details about the evaluated distances are thereafter available as calculator attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tablerows = zip(bcNaCl.distances, bcNaCl.sites0, bcNaCl.sites1, \n",
    "                bcNaCl.types0, bcNaCl.types1, bcNaCl.directions)\n",
    "for r in tablerows:\n",
    "    print(*r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A brief description of the functions and attributes in the *BondCalculator* class is available in its Python docstring and also in the package documentation at http://diffpy.github.io/diffpy.srreal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "help(BondCalculator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic and optimized evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculators in diffpy.srreal can be used in an *OPTIMIZED* mode, where pair contributions are evaluated only for atoms that have moved or otherwise changed since the last evaluation.  A full sum over all atom pairs is then carried only for the first evaluation of a structure or when the structure model changes significantly.  The *BASIC* mode performs a full sum every time.  The actual evaluation mode that was used is indicated with the *evaluatortypeused* attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.pdfcalculator import PDFCalculator\n",
    "\n",
    "nacl_stru.Bisoequiv = 0.5\n",
    "pdfc = PDFCalculator(rmax=10)\n",
    "pdfc.evaluatortype = 'OPTIMIZED'\n",
    "r, g = pdfc(nacl_stru)\n",
    "print(\"first evaluation:\", pdfc.evaluatortypeused)\n",
    "r, g = pdfc(nacl_stru)\n",
    "print(\"second evaluation:\", pdfc.evaluatortypeused)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When calculator objects are run with no arguments they reevaluate the last structure.  Here we use this for a simple benchmark of the BASIC and OPTIMIZED calculation modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdfc.evaluatortype = 'BASIC'\n",
    "%timeit pdfc()\n",
    "pdfc.evaluatortype = 'OPTIMIZED'\n",
    "%timeit pdfc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parallel calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculations in diffpy.srreal can also be sped up by evaluating the atom pairs in parallel and then combining the partial results.  The parallel calculations are configured by specifying a template calculator, the number of partial tasks, and a parallel map function that distributes function calls among parallel workers.\n",
    "\n",
    "The example below sets up 2 partial tasks, but uses the builtin *map* function which performs them in a sequence.  The calculation speed is thus not improved, but we can confirm the same result is obtained by combining 2 partial calculations.  If this notebook session is run with an active IPython cluster, remove the comment symbols to use the parallel *map* from the cluster and conduct a parallel evaluation for real."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.parallel import createParallelCalculator\n",
    "\n",
    "pdfc0 = PDFCalculator()\n",
    "pdfc2 = createParallelCalculator(PDFCalculator(), ncpu=2, pmap=map)\n",
    "## uncomment to use parallel map from IPython parallel cluster\n",
    "#from ipyparallel import Client\n",
    "#engineview = Client()[:]\n",
    "#pdfc2 = createParallelCalculator(PDFCalculator(), len(engineview), engineview.map)\n",
    "\n",
    "# calculate PDF using serial and parallel calculator\n",
    "r0, g0 = pdfc0(nacl_stru)\n",
    "r2, g2 = pdfc2(nacl_stru)\n",
    "\n",
    "# compare serial and parallel results\n",
    "fig1, ax1 = subplots()\n",
    "ax1.plot(r0, g0, 'x', r2, g2, '-r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PDF calculators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "diffpy.srreal provides 2 methods for calculating PDFs:\n",
    "\n",
    "- `PDFCalculator` - real-space summation of peak profiles suitable for periodic structures,\n",
    "- `DebyePDFCalculator` - reciprocal-space evaluation of Debye sum with Fourier transformation to real space - preferred for finite clusters and atoms with different *Q*-dependencies of their scattering factors.\n",
    "\n",
    "Here we demonstrate both calculators on an equally spaced chain of 10 atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.Structure import Structure, Atom\n",
    "\n",
    "# generate a chain of atoms\n",
    "c10 = Structure([Atom('C', (0, 0, z)) for z in range(10)])\n",
    "c10.Bisoequiv = 0.5\n",
    "\n",
    "# create PDF calculators per each method\n",
    "from diffpy.srreal.pdfcalculator import PDFCalculator\n",
    "from diffpy.srreal.pdfcalculator import DebyePDFCalculator\n",
    "\n",
    "pdfc = PDFCalculator(qmax=18, rmax=10)\n",
    "dbpdfc = DebyePDFCalculator(qmax=18, rmax=10)\n",
    "\n",
    "# Try qmin=1.0 values to simulate low-Q bound in measured PDF.\n",
    "qmin = 0\n",
    "r0, g0 = pdfc(c10, qmin=qmin)\n",
    "r2, g2 = dbpdfc(c10, qmin=qmin)\n",
    "\n",
    "# compare results\n",
    "fig2, ax2 = subplots()\n",
    "ax2.plot(r0, g0, 'x', r2, g2, '-r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculator components and their customization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculator objects in diffpy.srreal are composed of other objects for specialized partial tasks, such as lookup of scattering factors or calculation of profile functions.  These objects can be reconfigured or replaced at runtime to tweak the calculation.  For example, a simulation of a neutron PDF is achieved by replacing the *ScatteringFactorTable* object attached to a PDFCalculator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.scatteringfactortable import SFTXray, SFTNeutron\n",
    "\n",
    "rutile = loadStructure('TiO2_rutile.cif')\n",
    "pcx = PDFCalculator()\n",
    "pcx.scatteringfactortable = SFTXray()\n",
    "pcn = PDFCalculator()\n",
    "pcn.scatteringfactortable = SFTNeutron()\n",
    "\n",
    "fig3, ax3 = subplots()\n",
    "ax3.plot(*pcx(rutile), marker='x', label='x-ray')\n",
    "ax3.plot(*pcn(rutile), color='red', label='neutron')\n",
    "ax3.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to using the predefined objects, users can define their own components to further customize the calculation procedure.  As diffpy.srreal is a Python interface to the C++ library libdiffpy, the new components can be defined either in C++ or Python.  The example below defines a new rectangular peak profile for a real-space PDF calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# New profiles must derive from a generic PeakProfile class.\n",
    "from diffpy.srreal.peakprofile import PeakProfile\n",
    "\n",
    "class RectangularProfile(PeakProfile):\n",
    "    \"Rectangular profile with a unit area.\"\n",
    "    \n",
    "    # override functions from the base PeakProfile class\n",
    "    \n",
    "    def __call__(self, x, fwhm):\n",
    "        \"Evaluate zero-centered profile at the given position and width.\"\n",
    "        inside = (-fwhm/2.0 < x < +fwhm/2.0)\n",
    "        y = (1.0 / fwhm) if inside else 0.0\n",
    "        return y\n",
    "\n",
    "    def clone(self):\n",
    "        \"Return a copy of this profile object.\"\n",
    "        import copy\n",
    "        return copy.copy(self)\n",
    "\n",
    "    def create(self):\n",
    "        \"Return a new instance of this profile class.\"\n",
    "        return RectangularProfile()\n",
    "\n",
    "    def type(self):\n",
    "        \"Return a unique string name for this profile function.\"\n",
    "        return \"rectangularprofile\"\n",
    "\n",
    "    def xboundhi(self, fwhm):\n",
    "        \"Upper bound where profile becomes negligibly small.\"\n",
    "        return +0.5 * fwhm\n",
    "\n",
    "    def xboundlo(self, fwhm):\n",
    "        \"Lower bound where profile becomes negligibly small.\"\n",
    "        return -0.5 * fwhm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new `RectangularProfile` class can be added to a global registry of available profiles, so it can be specified by its string name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# deregister to allow repeated evaluations\n",
    "PeakProfile._deregisterType('rectangularprofile')\n",
    "RectangularProfile()._registerThisType()\n",
    "print(PeakProfile.getRegisteredTypes())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new peak profile can now be set in the *PDFCalculator* function using either the *RectangularProfile* object or the register name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdfc = PDFCalculator()\n",
    "\n",
    "pdfc.peakprofile = RectangularProfile()\n",
    "print(\"Profile assigned by instance:\", pdfc.peakprofile)\n",
    "\n",
    "pdfc.peakprofile = 'rectangularprofile'\n",
    "print(\"Profile assigned by name:\", pdfc.peakprofile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, here is a plot of the PDF obtained with the new rectangular profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r, g = pdfc(rutile)\n",
    "fig4, ax4 = subplots()\n",
    "ax4.plot(r, g);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extending with more structure types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculators in diffpy.srreal can use either the diffpy.Structure or pyobjcryst structure representations.  In detail, both of these types are converted to a native StructureAdapter class, which is understood by the calculation internals.  Structure objects from other crystallographic packages can thus be readily supported by providing appropriate conversion functions.\n",
    "\n",
    "Here we define a very simple structure class and then make it work with all calculators in diffpy.srreal.  The example structure type *CarbonChain* stores only the *z* coordinate of carbon atoms and their common isotropic displacement parameter Uiso."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CarbonChain(object):\n",
    "    def __init__(self):\n",
    "        self.Uiso = 0.004\n",
    "        self.z = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define a converter function, which builds a compatible adapter from a CarbonChain object.  The converter function is declared using the *@RegisterStructureAdapter* decorator and a dotted Python name of the converted class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.structureconverters import RegisterStructureAdapter\n",
    "\n",
    "@RegisterStructureAdapter('__main__.CarbonChain')\n",
    "def convertCarbonChain(thechain):\n",
    "    # AtomicStructureAdapter is a simple collection of atom positions\n",
    "    from diffpy.srreal.structureadapter import AtomicStructureAdapter\n",
    "    from diffpy.srreal.structureadapter import Atom as AdptAtom\n",
    "    adpt = AtomicStructureAdapter()\n",
    "    a = AdptAtom()\n",
    "    a.atomtype = 'C'\n",
    "    a.uc11 = a.uc22 = a.uc33 = thechain.Uiso\n",
    "    for z in thechain.z:\n",
    "        a.zc = z\n",
    "        adpt.append(a)\n",
    "    return adpt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now test it with the *createStructureAdapter* function to verify that the conversion works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.structureadapter import createStructureAdapter\n",
    "\n",
    "c5 = CarbonChain()\n",
    "c5.z = range(5)\n",
    "createStructureAdapter(c5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there were no error messages, the CarbonChain class is convertible to an adapter and can be directly used with any of the diffpy.srreal calculators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate bond distances in the c5 CarbonChain\n",
    "bc = BondCalculator()\n",
    "bc(c5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate PDF for the c5 chain\n",
    "pdfc = PDFCalculator(rmax=5)\n",
    "c5.Uiso = 0.006\n",
    "fig5, ax5 = subplots()\n",
    "ax5.plot(*pdfc(c5))\n",
    "c5.Uiso = 0.003\n",
    "ax5.plot(*pdfc(c5));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User-defined Lennard-Jones calculator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All calculator objects in diffpy.srreal are derived from the base *PairQuantity* class and use the same procedures for iteration over atom pairs.  Users can define their own calculators for pair-interaction related quantities by deriving from the `PairQuantity` class and modifying just the function which processes pair contributions.  \n",
    "\n",
    "As an example we define a class that calculates the Lennard-Jones potential, where the energy contribution from each atom pair is\n",
    "$$\n",
    "V_{\\mathit{LJ}} = 4\\epsilon \\left[ \\left(\\frac{\\sigma}{r}\\right)^{12} -\n",
    "\\left(\\frac{\\sigma}{r}\\right)^{6} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffpy.srreal.pairquantity import PairQuantity\n",
    "\n",
    "class LennardJonesCalculator(PairQuantity):\n",
    "\n",
    "    # parameters in the LJ formula\n",
    "    epsilon = 1.0\n",
    "    sigma = 1.0\n",
    "\n",
    "    # override functions from the PairQuantity class\n",
    "    \n",
    "    def __init__(self):\n",
    "        PairQuantity.__init__(self)\n",
    "        self.rmax = 10\n",
    "        return\n",
    "    \n",
    "    def _resetValue(self):\n",
    "        'Clear old results at the beginning of evaluation.'\n",
    "        self._vljtotal = 0.0\n",
    "        return\n",
    "\n",
    "    def _addPairContribution(self, bond, sumscale):\n",
    "        '''Process contribution from one atom pair.\n",
    "\n",
    "        bond     -- BaseBondGenerator object that describes the\n",
    "                    current pair of atom sites.\n",
    "        sumscale -- multiplicity of the current atom pair.\n",
    "        '''\n",
    "        nmrij = bond.distance() / self.sigma\n",
    "        vij = 4 * self.epsilon * (nmrij ** -12  -  nmrij ** -6)\n",
    "        self._vljtotal += 0.5 * vij * sumscale\n",
    "        return\n",
    "\n",
    "    def __call__(self, structure=None):\n",
    "        'Return LJ potential for the specified structure.'\n",
    "        self.eval(structure)\n",
    "        return self._vljtotal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new *LennardJonesCalculator* is aware of all structure conversions just as the pre-defined calculators.  Thus, we can try it out on the 5-atom *CarbonChain* object `c5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ljcalc = LennardJonesCalculator()\n",
    "print(\"ljcalc(c5) =\", ljcalc(c5))\n",
    "ljcalc.epsilon = 2\n",
    "print(\"Recalculate with LJ epsilon =\", ljcalc.epsilon)\n",
    "print(\"ljcalc() =\", ljcalc())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "For more information about the diffpy.srreal package see http://www.diffpy.org/diffpy.srreal."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}