Updated mPDF tutorial materials

cmi_scripts/mpdf/example_corefinement.ipynb

@@ -0,0 +1,307 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Atomic and magnetic PDF co-refinement using SrFit"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook provides an example of a simultaneous atomic + magnetic PDF refinement using SrFit. We will again use MnO as our test case."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Import necessary modules"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "### Import the necessary libraries\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from diffpy.mpdf import *\n",
+    "from diffpy.Structure import loadStructure\n",
+    "from scipy.optimize import least_squares\n",
+    "from diffpy.Structure.Parsers import getParser\n",
+    "from diffpy.srfit.pdf import PDFGenerator, PDFParser\n",
+    "from diffpy.srfit.fitbase import FitRecipe, FitResults\n",
+    "from diffpy.srfit.fitbase import Profile, FitContribution\n",
+    "\n",
+    "### Set all plots to be inline\n",
+    "%matplotlib notebook"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Set up nuclear/magnetic PDF structures and calculators"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Files containing our experimental data and structure file\n",
+    "dataFile = \"npdf_07334.gr\"\n",
+    "structureFile = \"MnO_R-3m.cif\"\n",
+    "\n",
+    "# load structure and space group from the CIF file\n",
+    "pcif = getParser('cif')\n",
+    "mno = pcif.parseFile(structureFile)\n",
+    "\n",
+    "# prepare profile object with experimental data\n",
+    "profile = Profile()\n",
+    "parser = PDFParser()\n",
+    "parser.parseFile(dataFile)\n",
+    "profile.loadParsedData(parser)\n",
+    "\n",
+    "# define range for pdf calculation\n",
+    "rmin = 0.01\n",
+    "rmax = 20\n",
+    "rstep = 0.01\n",
+    "\n",
+    "# setup calculation range for the PDF simulation\n",
+    "profile.setCalculationRange(xmin=rmin, xmax=rmax, dx=rstep)\n",
+    "\n",
+    "# prepare nucpdf function that simulates the nuclear PDF\n",
+    "nucpdf = PDFGenerator(\"nucpdf\")\n",
+    "nucpdf.setStructure(mno)\n",
+    "nucpdf.setProfile(profile)\n",
+    "\n",
+    "# prepare mpdf function that simulates the magnetic PDF\n",
+    "\n",
+    "# Create the Mn2+ magnetic species. Note that we are using a different\n",
+    "# setting for the unit cell, so the propagation vector is actually (0,0,3/2)\n",
+    "# instead of the (1/2,1/2,1/2) as it was for the pseudocubic setting.\n",
+    "mn2p = MagSpecies(struc=mno, label='Mn2+', magIdxs=[0,1,2],\n",
+    "                  basisvecs=np.array([1,0,0]), kvecs=np.array([0,0,1.5]),\n",
+    "                  ffparamkey='Mn2')\n",
+    "\n",
+    "\n",
+    "# Create and prep the magnetic structure\n",
+    "mstr = MagStructure()\n",
+    "mstr.loadSpecies(mn2p)\n",
+    "mstr.makeAll()\n",
+    "\n",
+    "# Set up the mPDF calculator.\n",
+    "\n",
+    "mc=MPDFcalculator(magstruc=mstr,rmin=rmin,rmax=rmax,rstep=rstep)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Create a \"total PDF\" calculator"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "totpdf = FitContribution('totpdf')\n",
+    "totpdf.addProfileGenerator(nucpdf)\n",
+    "totpdf.setProfile(profile)\n",
+    "\n",
+    "# Define an mPDF function that will be added to the total PDF calculator\n",
+    "def mpdf(parascale, ordscale, damp):\n",
+    "    mc.paraScale, mc.ordScale, mc.dampRate = parascale, ordscale, damp\n",
+    "    mc.magstruc.makeAll()\n",
+    "    dr = mc.calc(both=True)[2]\n",
+    "    return dr\n",
+    "\n",
+    "#Add mPDF to the FitContribution\n",
+    "totpdf.registerFunction(mpdf)\n",
+    "totpdf.setEquation(\"nucscale * nucpdf + mpdf(parascale, ordscale, damp)\")\n",
+    "\n",
+    "# Make magnetic PDF depend on any changes to the atomic structure.\n",
+    "# A structure change will now trigger reevaluation of the mPDF.\n",
+    "nucpdf.phase.addObserver(totpdf.ordscale.notify)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Set up the fit recipe and add the parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# The FitRecipe does the work of calculating the PDF with the fit variable\n",
+    "# that we give it.\n",
+    "mnofit = FitRecipe()\n",
+    "\n",
+    "# give the PDFContribution to the FitRecipe\n",
+    "mnofit.addContribution(totpdf)\n",
+    "\n",
+    "# Configure the fit variables and give them to the recipe.  We can use the\n",
+    "# srfit function constrainAsSpaceGroup to constrain the lattice and ADP\n",
+    "# parameters according to the CIF-loaded space group.\n",
+    "from diffpy.srfit.structure import constrainAsSpaceGroup\n",
+    "sgpars = constrainAsSpaceGroup(nucpdf.phase, pcif.spacegroup.short_name)\n",
+    "\n",
+    "# We can now cycle through the parameters and activate them as variables\n",
+    "for par in sgpars.latpars:\n",
+    "    mnofit.addVar(par)\n",
+    "# Set initial value for the ADP parameters, because CIF had no ADP data.\n",
+    "for par in sgpars.adppars:\n",
+    "    mnofit.addVar(par, value=0.003)\n",
+    "\n",
+    "# As usual, we add variables for the overall scale of the PDF and a delta2\n",
+    "# parameter for correlated motion of neighboring atoms.\n",
+    "mnofit.addVar(totpdf.nucscale, 1)\n",
+    "mnofit.addVar(nucpdf.delta2, 1.5)\n",
+    "\n",
+    "# We fix Qdamp based on prior information about our beamline.\n",
+    "mnofit.addVar(nucpdf.qdamp, 0.03, fixed=True)\n",
+    "\n",
+    "# add the mPDF variables\n",
+    "mnofit.addVar(totpdf.parascale, 4)\n",
+    "mnofit.addVar(totpdf.ordscale, 1.5)\n",
+    "mnofit.addVar(totpdf.damp, 0.01)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Perform the refinement"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Turn off printout of iteration number.\n",
+    "mnofit.clearFitHooks()\n",
+    "\n",
+    "# Initial structural fit\n",
+    "print \"Refine PDF using scipy's least-squares optimizer:\"\n",
+    "print \"  variables:\", mnofit.names\n",
+    "print \"  initial values:\", mnofit.values\n",
+    "least_squares(mnofit.residual, mnofit.values)\n",
+    "print \"  final values:\", mnofit.values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plot the result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "r = mnofit.totpdf.profile.x\n",
+    "gobs = mnofit.totpdf.profile.y\n",
+    "\n",
+    "# Get the calculated PDF and compute the difference between the calculated and\n",
+    "# measured PDF\n",
+    "gcalc = mnofit.totpdf.evaluate()\n",
+    "gnuc = mnofit.totpdf.evaluateEquation('nucscale * nucpdf')\n",
+    "gmag = mnofit.totpdf.evaluateEquation('mpdf')\n",
+    "\n",
+    "baseline = 1.1 * gobs.min()\n",
+    "gdiff = gobs - gcalc\n",
+    "baseline2 = 1.1 * (gdiff+baseline).min()\n",
+    "\n",
+    "# Plot the best-fit total PDF\n",
+    "ax=plt.figure().add_subplot(111)\n",
+    "ax.plot(r, gobs, 'bo', label=\"G(r) data\", markerfacecolor='none', markeredgecolor='b')\n",
+    "ax.plot(r, gcalc, 'r-', lw=2, label=\"G(r) fit\")\n",
+    "ax.plot(r, gdiff + baseline,'g-')\n",
+    "ax.plot(r, np.zeros_like(r) + baseline, 'k:')\n",
+    "ax.set_xlabel(r\"r ($\\AA$)\")\n",
+    "ax.set_ylabel(r\"G ($\\AA^{-2}$)\")\n",
+    "ax.set_xlim(xmax=mc.rmax)\n",
+    "plt.legend()\n",
+    "\n",
+    "plt.show()\n",
+    "\n",
+    "\n",
+    "# Plot just the mPDF\n",
+    "ax=plt.figure().add_subplot(111)\n",
+    "ax.plot(r, gobs-gnuc, 'bo', label=\"mPDF data\", markerfacecolor='none', markeredgecolor='b')\n",
+    "ax.plot(r, gmag, 'r-', lw=2, label=\"mPDF fit\")\n",
+    "ax.set_xlabel(r\"r ($\\AA$)\")\n",
+    "ax.set_ylabel(r\"D ($\\AA^{-2}$)\")\n",
+    "ax.set_xlim(xmax=mc.rmax)\n",
+    "plt.legend(loc=1)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "plt.close('all')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "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.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}