05-Fit-NaCl-BVS-restraint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo PDF refinement with BVS restraint\n",
    "\n",
    "In this example we setup a refinement of NaCl structure to experimental X-ray PDF.  The structure model will be also restrained by agreement of bond-valence sums with the expected atom valences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize plotting within notebook.\n",
    "%matplotlib notebook\n",
    "\n",
    "from matplotlib.pyplot import *\n",
    "rc('figure', figsize=(8, 6))\n",
    "\n",
    "# DiffPy-CMI functions for loading data and building a fitting recipe\n",
    "from diffpy.Structure import loadStructure\n",
    "from diffpy.srfit.pdf import PDFContribution\n",
    "from diffpy.srfit.structure import constrainAsSpaceGroup\n",
    "from diffpy.srfit.fitbase import FitRecipe, FitResults\n",
    "\n",
    "# A least squares fitting algorithm from scipy\n",
    "from scipy.optimize import leastsq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Files containing our experimental data and structure file\n",
    "dataFile = \"NaCl.gr\"\n",
    "structureFile = \"NaCl.cif\"\n",
    "spaceGroup = \"F m -3 m\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The first thing to construct is a contribution object which associates\n",
    "# observed data and numerical model.  PDFContribution is a specialized\n",
    "# contribution designed for PDF refinement of structure models.\n",
    "# Here we create a new PDFContribution named \"cpdf\".\n",
    "cpdf = PDFContribution(\"cpdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the PDF data and set the r-range over which we'll fit.\n",
    "cpdf.loadData(dataFile)\n",
    "cpdf.setCalculationRange(xmin=1, xmax=30, dx=0.02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a structure model that will be used for PDF calculation.\n",
    "nacl = loadStructure(structureFile)\n",
    "cpdf.addStructure(\"nacl\", nacl);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now cpdf.nacl now handles parameters for PDF calculation.\n",
    "# cpdf.nacl.phase contains parameters relevant for the structure model.\n",
    "# We can use the srfit function constrainAsSpaceGroup to constrain\n",
    "# the lattice and ADP according to the relevant space group.\n",
    "sgpars = constrainAsSpaceGroup(cpdf.nacl.phase, spaceGroup)\n",
    "print(\"Space group parameters are \" + \", \".join(p.name for p in sgpars) + \".\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cpdf.nacl.phase also provides a restrainBVS function, which defines\n",
    "# a soft restraint for agreement between the expected and calculated valences.\n",
    "# restrainBVS returns the active Restraint object.  We save it so we can\n",
    "# later manipulate its weight in the cost function.\n",
    "rbv = cpdf.nacl.phase.restrainBVS()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The FitRecipe does the work of managing refined variables and calculating\n",
    "# residuals from all contributions and restraints.\n",
    "thefit = FitRecipe()\n",
    "# Turn off printing of iteration number.\n",
    "thefit.clearFitHooks()\n",
    "\n",
    "# We give our PDF model to the fit to be optimized.\n",
    "thefit.addContribution(cpdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We now link various model parameters to the fit variables that\n",
    "# will be refined.  We will start with PDF scale, resolution damping\n",
    "# factor qdamp and the peak sharpening coefficient delta2.\n",
    "thefit.addVar(cpdf.scale, value=1)\n",
    "thefit.addVar(cpdf.qdamp, value=0.03)\n",
    "thefit.addVar(cpdf.nacl.delta2, value=5)\n",
    "\n",
    "# We will also refine independent structure parameters that were found\n",
    "# for our space group and atom coordinates.\n",
    "for par in sgpars.latpars:\n",
    "    thefit.addVar(par)\n",
    "# Here we set the initial value for the anisotropic displacement\n",
    "# parameters, because CIF had no ADP data.\n",
    "for par in sgpars.adppars:\n",
    "    thefit.addVar(par, value=0.005)\n",
    "# Position parameters can be also constrained.  This does nothing\n",
    "# for NaCl, because all atoms are at a special positions.\n",
    "for par in sgpars.xyzpars:\n",
    "    thefit.addVar(par)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now execute the fit using scipy's least square optimizer.\n",
    "# Let's define a few functions so it is easier to rerun the fit later.\n",
    "\n",
    "def namesandvalues():\n",
    "    \"Format names and values of the active fit variables.\"\n",
    "    return ' '.join(\"%s=%g\" % nv for nv in zip(thefit.names, thefit.values))\n",
    "\n",
    "def chattyfit():\n",
    "    print(\"Refine PDF using scipy's least-squares optimizer:\")\n",
    "    print(\"  initial: \" + namesandvalues())\n",
    "    rv = leastsq(thefit.residual, thefit.values)\n",
    "    print(\"  final: \" +  namesandvalues())\n",
    "    print('')\n",
    "    return rv\n",
    "\n",
    "def plotthefit(ax):\n",
    "    # Get the experimental data from the recipe\n",
    "    r = thefit.cpdf.profile.x\n",
    "    gobs = thefit.cpdf.profile.y\n",
    "    gcalc = thefit.cpdf.evaluate()\n",
    "    baseline = 1.1 * gobs.min()\n",
    "    gdiff = gobs - gcalc\n",
    "    ax.plot(r, gobs, 'bo', label=\"G(r) data\",\n",
    "            markerfacecolor='none', markeredgecolor='b')\n",
    "    ax.plot(r, gcalc, 'r-', label=\"G(r) fit\")\n",
    "    ax.plot(r, gdiff + baseline, 'g-', label=\"G(r) diff\")\n",
    "    ax.plot(r, 0.0 * r + baseline, 'k:')\n",
    "    ax.set_xlim(0, 30)\n",
    "    ax.set_xlabel(u\"r (Å)\")\n",
    "    ax.set_ylabel(u\"G (Å$^{-2}$)\")\n",
    "    ax.legend()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform the fit and plot it now.\n",
    "chattyfit();\n",
    "fig1, ax1 = subplots()\n",
    "plotthefit(ax1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Report fit results:\n",
    "results = FitResults(thefit)\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore convergence of the BVS-restrained fit\n",
    "\n",
    "Let's try to run fit with a far-off value of the lattice parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thefit.fix('delta2', 'Uiso_0', 'Uiso_4', 'qdamp')\n",
    "thefit.a = 4\n",
    "chattyfit()\n",
    "fig2, ax2 = subplots()\n",
    "plotthefit(ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit looks terrible, however it has still converged to a local minimum.  Note the second peak in the simulated PDF lines up with the first peak in the data.  To move out of the local minimum, we need to increase the weight of the BVS restraint *rvb*.  This can be done by adjusting its *sig* attribute, which provides a scale for the difference in valences.  The penalty contribution due to BVS is $[(V_{exp} - V_{sim}) / sig]^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase weight of the BVS restraint and refine again\n",
    "rbv.sig = 0.1\n",
    "chattyfit()\n",
    "fig3, ax3 = subplots()\n",
    "plotthefit(ax3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore effect of restraint weight on the lattice parameter\n",
    "\n",
    "We will run the refinement at different values of the restraint parameter sig to explore where does it overweight the lattice parameter refined from PDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "sigvalues = numpy.logspace(-4, 0)\n",
    "avalues = []\n",
    "for rbv.sig in sigvalues:\n",
    "    leastsq(thefit.residual, thefit.values)\n",
    "    avalues.append(thefit.a.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig4, ax4 = subplots()\n",
    "ax4.semilogx(sigvalues, avalues)\n",
    "ax4.set_xlabel('sig parameter of the BVS restraint')\n",
    "ax4.set_ylabel(u'cell parameter a (Å)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more tight restraint at *rvb.sig = 0.1* allowed the refinement to escape a far-off local minimum and recover the original fit.  The restraint is however still at the PDF-dominant region, so there is no significant change of the refined cell parameter."
   ]
  }
 ],
 "metadata": {
  "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
}