08-Orient-Naphthalene.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Refine orientation of naphthalene molecule\n",
    "\n",
    "Simulate PDF with separate inter and intramolecular displacement parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## initialization\n",
    "\n",
    "We start with some necessary imports from diffpy packages and define a convenience function for plotting the PDF fit.\n",
    "It is safe to skip to the [next section](#global-constants).\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use `print()` function as in Python 3\n",
    "from __future__ import print_function\n",
    "\n",
    "# notebook plotting setup\n",
    "%matplotlib notebook\n",
    "from matplotlib.pyplot import subplots\n",
    "\n",
    "# numerical math and optimization\n",
    "import numpy as np\n",
    "from scipy.optimize import leastsq\n",
    "\n",
    "# srfit PDF fitting objects\n",
    "from diffpy.srfit.pdf import PDFContribution\n",
    "from diffpy.srfit.fitbase import FitRecipe, FitResults, Calculator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## global constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the initial CIF structure\n",
    "g_cif = \"naphthalene.cif\"\n",
    "# offset of difference curve in PDF plot\n",
    "g_diff_baseline = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## define several helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def myLoadCrystal(filename):\n",
    "    \"Load pyobjcryst.Crystal object from a CIF file.\"\n",
    "    from pyobjcryst.crystal import CreateCrystalFromCIF\n",
    "    with open(filename) as fp:\n",
    "        rv = CreateCrystalFromCIF(fp, oneScatteringPowerPerElement=True)\n",
    "    return rv\n",
    "\n",
    "\n",
    "def getCrystalWithMolecule(ciffilename):\n",
    "    \"\"\"\n",
    "    Load the CIF file and group all scatterers as one Molecule.\n",
    "\n",
    "    This assumes the asymmetric unit is a single whole Molecule,\n",
    "    i.e., there are no other molecules or independent atoms in\n",
    "    the crystal and the molecule does not need to be\n",
    "    symmetry-expanded.\n",
    "    \"\"\"\n",
    "    from pyobjcryst.molecule import Molecule\n",
    "    # center of the molecule in fractional coordinates\n",
    "    xyzmol = np.array([0.0, 0.0, 0.0])\n",
    "    crst = myLoadCrystal(ciffilename)\n",
    "    mol = Molecule(crst, \"mol\")\n",
    "    for i in reversed(range(crst.GetNbScatterer())):\n",
    "        a = crst.GetScatterer(i)\n",
    "        xyzm = np.array([a.X, a.Y, a.Z]) - xyzmol\n",
    "        xc, yc, zc = crst.FractionalToOrthonormalCoords(*xyzm)\n",
    "        mol.AddAtom(xc, yc, zc, a.GetScatteringPower(), a.GetName())\n",
    "        crst.RemoveScatterer(a)\n",
    "    crst.AddScatterer(mol)\n",
    "    mol.X, mol.Y, mol.Z = xyzmol\n",
    "    return crst\n",
    "\n",
    "\n",
    "class quatzero(Calculator):\n",
    "    \"\"\"Return the zero-th component of normalized quaternion.\n",
    "    \"\"\"\n",
    "\n",
    "    def __call__(self, q1, q2, q3):\n",
    "        ssq = q1**2 + q2**2 + q3**2\n",
    "        q0 = np.sqrt(1.0 - ssq) if ssq < 1 else 0.0\n",
    "        return q0\n",
    "\n",
    "# end of quatzero\n",
    "\n",
    "\n",
    "def plotrecipe(axis, rec, what='obs, calc, diff'):\n",
    "    choice = set(what.replace(',', ' ').split())\n",
    "    r = rec.pcnt.r.value\n",
    "    gobs = rec.pcnt.y.value\n",
    "    gcalc = rec.pcnt.evaluate()\n",
    "    gdiff = gobs - gcalc\n",
    "    yall = [('obs', gobs),\n",
    "            ('calc', gcalc),\n",
    "            ('diff', gdiff + g_diff_baseline)]\n",
    "    pargs = sum([(r, y) for n, y in yall if n in choice], ())\n",
    "    rv = axis.plot(*pargs)\n",
    "    if 'calc' in choice and 'diff' in choice:\n",
    "        rv[-1].set_color(rv[-2].get_color())\n",
    "    return rv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## define PDF model composed of inter and intra-molecular components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcnt = PDFContribution('pcnt')\n",
    "pcnt.loadData('naphthalene.gr')\n",
    "pcnt.qdamp = 0.06\n",
    "pcnt.setCalculationRange(1.1, 25)\n",
    "# intra-molecular contribution from one molecule:\n",
    "crst = getCrystalWithMolecule(g_cif)\n",
    "mol = crst.GetScatterer(\"mol\")\n",
    "pcnt.addStructure('mol', mol, periodic=False)\n",
    "# inter-molecular contributions (wide)\n",
    "wcrst = getCrystalWithMolecule(g_cif)\n",
    "wmol = wcrst.GetScatterer('mol')\n",
    "pcnt.addStructure('wcrst', wcrst, periodic=True)\n",
    "pcnt.addStructure('wmol', wmol, periodic=False)\n",
    "pcnt.setEquation('scale * (mol + wcrst - wmol)')\n",
    "\n",
    "# speed up simulation by using parallel jobs.\n",
    "pcnt.wcrst.parallel(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## define FitRecipe for this PDF model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define FitRecipe for this organic crystal\n",
    "ocfit = FitRecipe()\n",
    "ocfit.clearFitHooks()\n",
    "ocfit.addContribution(pcnt)\n",
    "\n",
    "# expose the overall data scale\n",
    "ocfit.addVar(ocfit.pcnt.scale, name='scale', value=0.1)\n",
    "ocfit.addVar(ocfit.pcnt.qdamp, name='qdamp')\n",
    "\n",
    "# expose unit cell parameters\n",
    "pwcrst = ocfit.pcnt.wcrst.phase\n",
    "ocfit.addVar(pwcrst.a)\n",
    "ocfit.addVar(pwcrst.b)\n",
    "ocfit.addVar(pwcrst.c)\n",
    "# objcryst uses radians for cell angle beta, but let us refine\n",
    "# in degrees.\n",
    "ocfit.newVar('beta', value=np.degrees(pwcrst.beta.value))\n",
    "ocfit.constrain(pwcrst.beta, 'radians(beta)')\n",
    "\n",
    "# objcryst links all atoms of the same type to one scattering power\n",
    "# with common Biso.  It is thus sufficient to constrain just one of\n",
    "# those.\n",
    "\n",
    "# first let's ensure there are only 2 independent scattering powers\n",
    "# for \"C\" and \"H\" in the CIF file.\n",
    "assert 2 == crst.GetScatteringPowerRegistry().GetNb(), \\\n",
    "    \"unexpected number of atom species\"\n",
    "\n",
    "# expose isotropic displacement parameters for carbon and hydrogen\n",
    "# in the molecule.  CIF file has sites labeled \"C1\", \"H1\" so we can\n",
    "# refer to them by those names.\n",
    "pmol = ocfit.pcnt.mol.phase\n",
    "ocfit.addVar(pmol.C1.Biso, name='bisoC', value=1.0)\n",
    "ocfit.restrain('bisoC', lb=0.02, sig=1e-4)\n",
    "# hydrogen gives negligible contribution, so we keep them constant\n",
    "ocfit.addVar(pmol.H1.Biso, name='bisoH', value=1, fixed=True)\n",
    "\n",
    "# for intra-molecular component we use the same binter value for all atoms:\n",
    "ocfit.addVar(pwcrst.mol.C1.Biso, name='binter', value=10)\n",
    "ocfit.constrain(pwcrst.mol.H1.Biso, 'binter')\n",
    "\n",
    "# expose rotation quaternions for the molecule\n",
    "pwmol = ocfit.pcnt.wcrst.phase.mol\n",
    "ocfit.addVar(pwmol.q1, name='wq1')\n",
    "ocfit.addVar(pwmol.q2, name='wq2')\n",
    "ocfit.addVar(pwmol.q3, name='wq3')\n",
    "\n",
    "# constrain q0 so we have normalized quaternion.\n",
    "# first define a function \"fixq0\" that calculates q0 from [q1, q2, q3]\n",
    "ocfit.registerCalculator(quatzero('fixq0'), argnames=[])\n",
    "ocfit.constrain(pwmol.q0, 'fixq0(wq1, wq2, wq3)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fix all variables and show the initial residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocfit.fix('all')\n",
    "ocfit.free('scale')\n",
    "ocfit.scalarResidual()\n",
    "print(\"== INITIAL ==\\n\")\n",
    "print(FitResults(ocfit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## refine cell parameters and displacement factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocfit.free('a', 'b', 'c')\n",
    "ocfit.free('bisoC', 'binter')\n",
    "leastsq(ocfit.residual, ocfit.values)\n",
    "print(\"== INTERMEDIATE ==\\n\")\n",
    "print(FitResults(ocfit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## refine the orientation quaternion for the molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocfit.free('wq1', 'wq2', 'wq3')\n",
    "leastsq(ocfit.residual, ocfit.values)\n",
    "print(\"== REFINED ==\\n\")\n",
    "print(FitResults(ocfit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## rotate molecule away from the minimum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rotate molecule away from the minimum\n",
    "ocfit.wq1 << 0.1\n",
    "ocfit.fix('a', 'b', 'c', 'binter', 'wq1', 'wq2', 'wq3')\n",
    "leastsq(ocfit.residual, ocfit.values)\n",
    "print(\"== MISORIENTED ==\\n\")\n",
    "print(FitResults(ocfit))\n",
    "fig1, ax1 = subplots()\n",
    "plotrecipe(ax1, ocfit)\n",
    "ax1.set_title('misoriented molecule');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## refine molecule orientation to return to the correct structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocfit.free('wq1', 'wq2', 'wq3')\n",
    "leastsq(ocfit.residual, ocfit.values)\n",
    "print(\"== ORIENTED BACK ==\\n\")\n",
    "print(FitResults(ocfit))\n",
    "fig2, ax2 = subplots()\n",
    "plotrecipe(ax2, ocfit);\n",
    "ax2.set_title('molecule oriented back');"
   ]
  }
 ],
 "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": 2
}