#!/usr/bin/env python
########################################################################
#
# diffpy.srfit by DANSE Diffraction group
# Simon J. L. Billinge
# (c) 2009 The Trustees of Columbia University
# in the City of New York. All rights reserved.
#
# File coded by: Chris Farrow
#
# See AUTHORS.txt for a list of people who contributed.
# See LICENSE_DANSE.txt for license information.
#
########################################################################
"""Example of a PDF refinement of two-phase structure.
This example uses PDFGenerator to refine a the two phase nickel-silicon
structure to all the available data.
"""
import numpy
from pyobjcryst import loadCrystal
from diffpy.srfit.pdf import PDFGenerator, PDFParser
from diffpy.srfit.fitbase import Profile
from diffpy.srfit.fitbase import FitContribution, FitRecipe
from diffpy.srfit.fitbase import FitResults
from gaussianrecipe import scipyOptimize
####### Example Code
def makeProfile(datafile):
"""Make an place data within a Profile."""
profile = Profile()
parser = PDFParser()
parser.parseFile(datafile)
profile.loadParsedData(parser)
profile.setCalculationRange(xmax = 20)
return profile
def makeContribution(name, generator, profile):
"""Make a FitContribution and add a generator and profile."""
contribution = FitContribution(name)
contribution.addProfileGenerator(generator)
contribution.setProfile(profile, xname = "r")
return contribution
def makeRecipe(ciffile_ni, ciffile_si, xdata_ni, ndata_ni, xdata_si,
xdata_sini):
"""Create a fitting recipe for crystalline PDF data."""
## The Profiles
# We need a profile for each data set.
xprofile_ni = makeProfile(xdata_ni)
xprofile_si = makeProfile(xdata_si)
nprofile_ni = makeProfile(ndata_ni)
xprofile_sini = makeProfile(xdata_sini)
## The ProfileGenerators
# We create one for each phase and share the phases.
xgenerator_ni = PDFGenerator("xG_ni")
stru = loadCrystal(ciffile_ni)
xgenerator_ni.setStructure(stru)
phase_ni = xgenerator_ni.phase
xgenerator_si = PDFGenerator("xG_si")
stru = loadCrystal(ciffile_si)
xgenerator_si.setStructure(stru)
phase_si = xgenerator_si.phase
ngenerator_ni = PDFGenerator("nG_ni")
ngenerator_ni.setPhase(phase_ni)
xgenerator_sini_ni = PDFGenerator("xG_sini_ni")
xgenerator_sini_ni.setPhase(phase_ni)
xgenerator_sini_si = PDFGenerator("xG_sini_si")
xgenerator_sini_si.setPhase(phase_si)
## The FitContributions
# We one of these for each data set.
xcontribution_ni = makeContribution("xnickel", xgenerator_ni, xprofile_ni)
xcontribution_si = makeContribution("xsilicon", xgenerator_si, xprofile_si)
ncontribution_ni = makeContribution("nnickel", ngenerator_ni, nprofile_ni)
xcontribution_sini = makeContribution("xsini", xgenerator_sini_ni,
xprofile_sini)
xcontribution_sini.addProfileGenerator(xgenerator_sini_si)
xcontribution_sini.setEquation("scale * (xG_sini_ni + xG_sini_si)")
# As explained in another example, we want to minimize using Rw^2.
xcontribution_ni.setResidualEquation("resv")
xcontribution_si.setResidualEquation("resv")
ncontribution_ni.setResidualEquation("resv")
xcontribution_sini.setResidualEquation("resv")
# Make the FitRecipe and add the FitContributions.
recipe = FitRecipe()
recipe.addContribution(xcontribution_ni)
recipe.addContribution(xcontribution_si)
recipe.addContribution(ncontribution_ni)
recipe.addContribution(xcontribution_sini)
# Now we vary and constrain Parameters as before.
for par in phase_ni.sgpars:
recipe.addVar(par, name = par.name + "_ni")
delta2_ni = recipe.newVar("delta2_ni", 2.5)
recipe.constrain(xgenerator_ni.delta2, delta2_ni)
recipe.constrain(ngenerator_ni.delta2, delta2_ni)
recipe.constrain(xgenerator_sini_ni.delta2, delta2_ni)
for par in phase_si.sgpars:
recipe.addVar(par, name = par.name + "_si")
delta2_si = recipe.newVar("delta2_si", 2.5)
recipe.constrain(xgenerator_si.delta2, delta2_si)
recipe.constrain(xgenerator_sini_si.delta2, delta2_si)
# Now the experimental parameters
recipe.addVar(xgenerator_ni.scale, name = "xscale_ni")
recipe.addVar(xgenerator_si.scale, name = "xscale_si")
recipe.addVar(ngenerator_ni.scale, name = "nscale_ni")
recipe.addVar(xcontribution_sini.scale, 1.0, "xscale_sini")
recipe.newVar("pscale_sini_ni", 0.8)
recipe.constrain(xgenerator_sini_ni.scale, "pscale_sini_ni")
recipe.constrain(xgenerator_sini_si.scale, "1 - pscale_sini_ni")
# The qdamp parameters are too correlated to vary so we fix them based on
# previous measurments.
xgenerator_ni.qdamp.value = 0.055
xgenerator_si.qdamp.value = 0.051
ngenerator_ni.qdamp.value = 0.030
xgenerator_sini_ni.qdamp.value = 0.052
xgenerator_sini_si.qdamp.value = 0.052
# Give the recipe away so it can be used!
return recipe
def plotResults(recipe):
"""Plot the results contained within a refined FitRecipe."""
# All this should be pretty familiar by now.
xnickel = recipe.xnickel
xr_ni = xnickel.profile.x
xg_ni = xnickel.profile.y
xgcalc_ni = xnickel.profile.ycalc
xdiffzero_ni = -0.8 * max(xg_ni) * numpy.ones_like(xg_ni)
xdiff_ni = xg_ni - xgcalc_ni + xdiffzero_ni
xsilicon = recipe.xsilicon
xr_si = xsilicon.profile.x
xg_si = xsilicon.profile.y
xgcalc_si = xsilicon.profile.ycalc
xdiffzero_si = -0.8 * max(xg_si) * numpy.ones_like(xg_si)
xdiff_si = xg_si - xgcalc_si + xdiffzero_si
nnickel = recipe.nnickel
nr_ni = nnickel.profile.x
ng_ni = nnickel.profile.y
ngcalc_ni = nnickel.profile.ycalc
ndiffzero_ni = -0.8 * max(ng_ni) * numpy.ones_like(ng_ni)
ndiff_ni = ng_ni - ngcalc_ni + ndiffzero_ni
xsini = recipe.xsini
xr_sini = xsini.profile.x
xg_sini = xsini.profile.y
xgcalc_sini = xsini.profile.ycalc
xdiffzero_sini = -0.8 * max(xg_sini) * numpy.ones_like(xg_sini)
xdiff_sini = xg_sini - xgcalc_sini + xdiffzero_sini
import pylab
pylab.subplot(2, 2, 1)
pylab.plot(xr_ni,xg_ni,'bo',label="G(r) x-ray nickel Data")
pylab.plot(xr_ni,xgcalc_ni,'r-',label="G(r) x-ray nickel Fit")
pylab.plot(xr_ni,xdiff_ni,'g-',label="G(r) x-ray nickel diff")
pylab.plot(xr_ni,xdiffzero_ni,'k-')
pylab.xlabel(r"$r (\AA)$")
pylab.ylabel(r"$G (\AA^{-2})$")
pylab.legend(loc=1)
pylab.subplot(2, 2, 2)
pylab.plot(xr_si,xg_si,'bo',label="G(r) x-ray silicon Data")
pylab.plot(xr_si,xgcalc_si,'r-',label="G(r) x-ray silicon Fit")
pylab.plot(xr_si,xdiff_si,'g-',label="G(r) x-ray silicon diff")
pylab.plot(xr_si,xdiffzero_si,'k-')
pylab.legend(loc=1)
pylab.subplot(2, 2, 3)
pylab.plot(nr_ni,ng_ni,'bo',label="G(r) neutron nickel Data")
pylab.plot(nr_ni,ngcalc_ni,'r-',label="G(r) neutron nickel Fit")
pylab.plot(nr_ni,ndiff_ni,'g-',label="G(r) neutron nickel diff")
pylab.plot(nr_ni,ndiffzero_ni,'k-')
pylab.legend(loc=1)
pylab.subplot(2, 2, 4)
pylab.plot(xr_sini,xg_sini,'bo',label="G(r) x-ray sini Data")
pylab.plot(xr_sini,xgcalc_sini,'r-',label="G(r) x-ray sini Fit")
pylab.plot(xr_sini,xdiff_sini,'g-',label="G(r) x-ray sini diff")
pylab.plot(xr_sini,xdiffzero_sini,'k-')
pylab.legend(loc=1)
pylab.show()
return
if __name__ == "__main__":
# Make the data and the recipe
ciffile_ni = "data/ni.cif"
ciffile_si = "data/si.cif"
xdata_ni = "data/ni-q27r60-xray.gr"
ndata_ni = "data/ni-q27r100-neutron.gr"
xdata_si = "data/si-q27r60-xray.gr"
xdata_sini = "data/si90ni10-q27r60-xray.gr"
# Make the recipe
recipe = makeRecipe(ciffile_ni, ciffile_si, xdata_ni, ndata_ni, xdata_si,
xdata_sini)
# Optimize
scipyOptimize(recipe)
# Generate and print the FitResults
res = FitResults(recipe)
res.printResults()
# Plot!
plotResults(recipe)
# End of file