Question about calculator in diffpy.srfit.pdf.pdfcontribution module

[Ruochen Cao]

Hi, everyone,

I'm a new user of CMI. Accoring to the example code in add2019-diffpy-cmi (06-Fit-CdSe-nanoparticle.ipynb), I can calculate the PDF of nanoparticle now. But I'm wondering the calculation method/equation used in this code. How does the PDFContribution module get PDF data from my input structure file? 

I've looked up the Documentation page of CMI，diffpy.srfit and diffpy.srreal part. In diffpy.srfit , pdfcontribution part,  there is '_calculators – A managed dictionary of Calculators'  in PDFContribution class, but I can't get this  _calculater value by 

print(gtot._calculator)        or        print(inspect.getmembers(gtot._calculator))

And in diffpy.srreal.pdfcalculator  there are DebyePDFCalculator and PDFCalculator. I suppose the first calculator uses debye scattering equation and the second uses usual PDF equation. Maybe in srfit.pdf.pdfcontribution module, the calculator is similar to PDFCalculator?

And the  last question, I'm going to calculate PDF of a nanoparticle system, is is more accurate to use DebyePDFCalculator ? 

Thanks

Ruochen

Here is the code part:

import matplotlib

# make figures smaller

matplotlib.rc('figure', figsize=(5, 3.75))

import numpy.linalg

from scipy.optimize import leastsq

from diffpy.srfit.fitbase import FitRecipe, FitResults

from cdse_functions import differenceplot

def new_cdsecluster_fit():

    from diffpy.srfit.pdf import PDFContribution

    from diffpy.structure import loadStructure

    gtot = PDFContribution('gtot')

    gtot.qdamp = c_qdamp

    gtot.setQmin(c_qmin)

    gtot.loadData(c_pdfdatafile)

    cluster = loadStructure(c_clusterfile)

    gtot.addStructure('cdse', cluster, periodic=False)

    gtot.profile.setCalculationRange(*c_fitrange)

    fit = FitRecipe()

    fit.clearFitHooks()

    fit.addContribution(gtot)

    sc0 = best_scale(gtot.evaluate(), gtot.profile.y)

    fit.addVar(gtot.scale, sc0)

    fit.addVar(gtot.cdse.delta2, value=1)

    fit.newVar('expansion', value=0)

    fit.newVar('BisoCd', 0.5)

    fit.newVar('BisoSe', 0.5)

    ph = gtot.cdse.phase

    fit.constrain(ph.lattice.a, '(1 + expansion)')

    fit.constrain(ph.lattice.b, '(1 + expansion)')

    fit.constrain(ph.lattice.c, '(1 + expansion)')

    for pa in gtot.cdse.phase.atoms:

        if pa.atom.element == "Cd":

            fit.constrain(pa.Biso, 'BisoCd')

        if pa.atom.element == "Se":

            fit.constrain(pa.Biso, 'BisoSe')

    return fit

fcluster = new_cdsecluster_fit()

leastsq(fcluster.residual, fcluster.values)

differenceplot(fcluster, title="fit with Cd52Se35 tetrahedral cluster")

[Ruochen Cao]

And I noticed that there are PDFGenerator and  DebyePDFGenerator class in diffpy.srfit. What's the difference between those PDF generators and PDF calculators in diffpy.srreal package? How can I fit data from PDFCalculator or DebyePDFCalculator(srreal package) to my experimental data?

Ruochen

[Mikkel Juelsholt]

Hi Ruochen

That is some good questions. Personally I cannot figure out what is going on in the scripts from ADD either. But I can answer some of your questions. So we might need Pavol to make sure we are correct. 

The PDFgenerators are used for generating PDFs for fitting, while the PDFcalculators only calculates a PDF. 

So you should use the PDFgenerators to fit with. I hope I have attached a script which is build in a different way than the examples from ADD. Hopefully it makes it more clear how to use the two generators for a fit. The script actually makes a fit to a single dataset using a crystalline model in the PDFgenerator and a xyz-model for the DebyePDFgenerator.

Which one to use if you only need one of them is not an easy question to answer. The Debye method is very flexible as you can calculate scattering from any object, but it gives you a large number of parameters. It can also be very difficult to refine stuff like size and size distribution, though it can be done. Also because you have abandoned the unit cell it is not easy to refine symmetry related parameters and so on. But for some systems it might be the only option. 

If you "just" want to fit a crystalline nanoparticle I would start with a simple fit with a cif-file and the PDFgenerator. It is going to be faster and a more stable refinement with fewer parameters. 

Hope this help, even though it does not answer the underlying questions about Diffpy-CMI classes and packages

Mikkel

[Mikkel Juelsholt]

The script looks like this. I am not allowed ot upload the notebook itself.

# Import the packages we'll need

%matplotlib notebook

from matplotlib.pyplot import *

import string

from collections import OrderedDict

import numpy as np

from pyobjcryst import loadCrystal

from diffpy.srfit.pdf import DebyePDFGenerator, PDFGenerator, PDFParser

from diffpy.srfit.fitbase import Profile

from diffpy.srfit.fitbase import FitContribution, FitRecipe

from diffpy.srfit.fitbase import FitResults, initializeRecipe

from diffpy.Structure import Structure

from diffpy.Structure import loadStructure

# Define functions for optimization and plotting that we will use later

def scipyOptimize(recipe):

    from scipy.optimize.minpack import leastsq

    print("Fit using scipy's LM optimizer")

    leastsq(recipe.residual, recipe.getValues())

    return

def plotRecipe(recipe, ax=None):

    """Plot recipe in the specified axes `ax`.

    If `ax` is None, create a new figure.

    Return `ax`

    """

    if ax is None:

        fig, ax = subplots()

    r = recipe.pdf.profile.x

    g = recipe.pdf.profile.y

    gcalc = recipe.pdf.evaluate()

    diffzero = -0.8 * max(g) * np.ones_like(g)

    diff = g - gcalc + diffzero

    ax.plot(r,g,'o',label="G(r) Data", mfc='none', mec='blue')

    ax.plot(r, gcalc,'r-',label="G(r) Fit")

    ax.plot(r,diff,'g-',label="G(r) diff")

    ax.plot(r, diffzero,'k-')

    ax.set_xlabel("$r (\AA)$")

    ax.set_ylabel("$G (\AA^{-2})$")

    ax.legend(loc=1)

    return ax

# Import the data and make it a PDFprofile. 

grdata =  'ZrO2_JDR_85_20171001-225809_6bd5cf_001q23q22.5.gr'

pdfprofile = Profile()

pdfparser = PDFParser()

pdfparser.parseFile(grdata)

pdfprofile.loadParsedData(pdfparser)

# Setup the PDFgenerator 1 that calculates the PDF from the CIF-file

pdfgenerator_crystal = PDFGenerator("G1")

pdfgenerator_crystal.setQmax(22.0)

pdfgenerator_crystal.setQmin(0.6)

# Load structure from the CIF file.

ciffile = "tet_ZrO2.cif"

crystal_structure = loadCrystal(ciffile)

pdfgenerator_crystal.setStructure(crystal_structure)

# Setup the PDFgenerator 2 that calculates the PDF from the xyz-file.

pdfgenerator_cluster = DebyePDFGenerator("G2")

pdfgenerator_cluster.setQmax(22.0)

pdfgenerator_cluster.setQmin(0.6)

# Load structure from the xyz-file

xyzfile = "ZrPO.xyz"

xyz_structure = loadStructure(xyzfile)

pdfgenerator_cluster.setStructure(xyz_structure, periodic = False)

# Add the profile and both generators to the PDFcontribution

pdfcontribution = FitContribution("pdf")

pdfcontribution.setProfile(pdfprofile, xname="r") 

pdfcontribution.addProfileGenerator(pdfgenerator_crystal)

pdfcontribution.addProfileGenerator(pdfgenerator_cluster)

# Setup correction scaling due to finite spherical shape. SP diameter from PDFgui

from diffpy.srfit.pdf.characteristicfunctions import sphericalCF

pdfcontribution.registerFunction(sphericalCF, name="f1", argnames=['r', 'psize1'])

pdfcontribution.psize1 << 70;

# Use scaling factors proportional to molar content

pdfcontribution.setEquation('mc1*G1*f1 + mc2*G2')

# Define the recipe to do the fit and add it to the PDFcontribution

recipe = FitRecipe()

recipe.addContribution(pdfcontribution)

# Avoid too much output during fitting 

recipe.clearFitHooks()

# Add the two scale factors.

recipe.addVar(pdfcontribution.mc1, 1.0, tag = "scale")

recipe.addVar(pdfcontribution.mc2, 1.0, tag = 'scale')

# set qdamp, qbroad for the instrument

qdamp = 0.04

qbroad = 0.01

# Add the instrumental parameters to the two generators

pdfgenerator_crystal.qdamp.value = qdamp

pdfgenerator_crystal.qbroad.value = qbroad

pdfgenerator_cluster.qdamp.value = qdamp

pdfgenerator_cluster.qbroad.value = qbroad

# Add the delta2 parameters for the two phases, and make sure they cannot take unphysical values

recipe.addVar(pdfgenerator_crystal.delta2, 2, name = "delta2_crystal", tag = "delta2")

#Start by have delta2 for the cluster 0

recipe.addVar(pdfgenerator_cluster.delta2, 0, name = "delta2_cluster", tag = "delta2")

# Uncommment below to contrain the delta2s to tke the same value

#recipe.constrain("delta2_crystal","delta2_cluster")

# Add the psize variable for diameter of the spherical nanoparticle.

recipe.addVar(pdfcontribution.psize1);

# Add the structural paramters using space group constaints for the crystal

# Ignore the ADP-s which will be set later.

phase_crystal = pdfgenerator_crystal.phase

sgpars_crystal = phase_crystal.sgpars

for par in sgpars_crystal.latpars:

    recipe.addVar(par, name=par.name+'_s', tag='cell')

for par in sgpars_crystal.xyzpars:

    lclabel = par.par.obj.GetName().lower()

    lcsymbol = lclabel.rstrip(string.digits)

    # name the variables

    name="{}_{}".format(par.par.name, lclabel)

    tags = ['xyz', 'xyz_' + lclabel, 'xyz_' + lcsymbol]

    recipe.addVar(par, name=name, tags=tags)

phase_crystal = pdfgenerator_crystal.phase

sgpars_crystal = phase_crystal.sgpars

# Use the same isotropic displacement parameters for all Zr in the crystal

initial_Biso_Zr = 1.0

recipe.newVar(name='Biso_Zr', value=initial_Biso_Zr, tags=['adp_zr'])

oxypars = [s for s in phase_crystal.scatterers if s.element.startswith('Zr')]

for a in oxypars:

    recipe.constrain(a.Biso, 'Biso_Zr')

# Use the same isotropic displacement parameters for all oxygens

initial_Biso_O = 1.0

recipe.newVar(name='Biso_O', value=initial_Biso_O, tags=['adp_o'])

oxypars = [s for s in phase_crystal.scatterers if s.element.startswith('O')]

for a in oxypars:

    recipe.constrain(a.Biso, 'Biso_O')

# Add ADP and "cell" for the cluster

phase_cluster = pdfgenerator_cluster.phase

atoms = phase_cluster.getScatterers()

#Add "zoomscale" to the cluster. Zoomscale1 is a factor that is multiplied on all x-coordinates and so on. 

# Can change the distances between the atoms in the cluster

lat = phase_cluster.getLattice()

recipe.newVar("zoomscale1", 1.0, tag = "lat")

recipe.newVar("zoomscale2", 1.0, tag = "lat")

recipe.newVar("zoomscale3", 1,0, tag = "lat")

recipe.constrain(lat.a, 'zoomscale1')

recipe.constrain(lat.b, 'zoomscale2')

recipe.constrain(lat.c, 'zoomscale3')

# Set the ADP of the cluster atoms. For now the O and C have the same ADP

Zr_cluster = recipe.newVar("Zr_Biso_cluster", 0.5, tag = 'adp_zr')

O_cluster = recipe.newVar("O_Biso_cluster", 0.1, tag = 'adp_o')

for atom in atoms:

    if atom.element.title() == "Zr":

        recipe.constrain(atom.Biso, Zr_cluster)

    elif atom.element.title() == "O" or "C":

        recipe.constrain(atom.Biso, O_cluster)

# Uncomment the line below to contrain the O in the crystal and O in the cluster to be the same.

#recipe.constrain(O_cluster,'Biso_O_s')

pos_limit1 = 0.04

#Add refinable parameters for the Zr atoms in the xyz-structure

for atom in atoms:

    if atom.element == 'Zr':

        recipe.addVar(atom.x, name=atom.name+'_x1', tag='xyz1')

        recipe.addVar(atom.y, name=atom.name+'_y1', tag='xyz1')

        recipe.addVar(atom.z, name=atom.name+'_z1', tag='xyz1')       

        #restrain atomic positions

        lbx = atom.x.value-pos_limit1

        ubx = atom.x.value+pos_limit1

        recipe.restrain(atom.x, lb=lbx, ub = ubx, sig=0.001)

        lby = atom.y.value-pos_limit1

        uby = atom.y.value+pos_limit1

        recipe.restrain(atom.y, lb=lby, ub = uby, sig=0.001)

        lbz = atom.z.value-pos_limit1

        ubz = atom.z.value+pos_limit1

        recipe.restrain(atom.y, lb=lby, ub = uby, sig=0.001)

# Set the r-range to fit

pdfprofile.setCalculationRange(xmin = 1.5, xmax = 60)

# First refine only the scale

recipe.fix('all') # Fix the parameters written in the bracket, in this case all of them

recipe.free('scale', 'cell', 'lat') #Free the parameters with the tag 'scale' and 'xyz'

scipyOptimize(recipe) # Do the least-square refinement

recipe.free('adp_zr') #Free the parameters with the tag 'adp_zr'

scipyOptimize(recipe) # Do the least-square refinement

recipe.free('delta2_crystal') ##Free the parameter delta2_crystal

scipyOptimize(recipe) # Do the least-square refinement again with the extra freed parameter

recipe.free('psize1') ##Free the parameter psize1

scipyOptimize(recipe) # Do the least-square refinement again with the extra freed parameter

recipe.free('Biso_O') ##Free the parameter Biso_O

scipyOptimize(recipe) # Do the least-square refinement again with the extra freed parameter

recipe.free('xyz') ##Free the parameters with the tag 'xyz'

scipyOptimize(recipe) # Do the least-square refinement again with the extra freed parameter

print(FitResults(recipe)) # Print the results from the fit

plotRecipe(recipe); # Plot the fit

# Export the fit for plotting and plot the two contributions to the fit independently

r = recipe.pdf.profile.x

g = recipe.pdf.profile.y

gcalc = recipe.pdf.evaluate()

diffzero = -0.8 * max(g) * np.ones_like(g)

diff = g - gcalc + diffzero

g1 = recipe.pdf.evaluateEquation('mc1 * G1 * f1')

g2 = recipe.pdf.evaluateEquation('mc2 * G2*10')

#Do the actual plotting

fig, ax = subplots()

plotRecipe(recipe, ax);

ax.plot(r, g1+10, label='crystal')

ax.plot(r, g2+15, label='cluster')

ax.legend();

[Ruochen Cao]

Hi Mikkel, 

Thanks for your reply. Althongh question about Diffpy-CMI classes and packages is not clear yet , the script does help a lot. I'll try to use PDFgenerators to calculate and fit with my system.

Ruochen

[He Lin]

Hi: 

     I  have this related question that if the  DebyePDFGenerator  will always calculate the fourier transform in each step of refinement ad compared with G_exp(r)? or actually what is refined is the fourier transformed G_exp(r), that is F(Q)?  Then  maybe we can use S_exp(Q) or F_exp(Q) as input directly?

Is there fast algorithm to do fourier transform of F_cal(Q) in each refinement step? Otherwize it seems to take time?

    Thanks,

                   he

 

[Simon Billinge]

Thanks He.

This is probably a question for Pavol, though we can dig into the code too.  Since we are computing the DSE anyway, it seems like it would make sense to allow it to be refined directly if people want to.   

Other programs like DEBUSSY (http://scripts.iucr.org/cgi-bin/paper?S0021889810041889) will also do this, but they are generally used for larger models than we are interested in in PDF-land.

S

--
You received this message because you are subscribed to the Google Groups "diffpy-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to diffpy-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/diffpy-users/9775132b-cfd4-4fc7-8662-1b07daea2445n%40googlegroups.com.

--

Professor Simon Billinge

Columbia University

[He Lin]

Thanks, Simon,

     We will look into the DEBUSSY software, I believe Pavol has though about this, or maybe it is already implemented in CMI, it would be really interesting to see both Q and r space fitting,  :-) 

                       he

[Pavol Juhas]

Hi He,

There was a similar question about F(Q), S(Q) calculations last year -
https://groups.google.com/d/msg/diffpy-users/kKSHs5D_2wg/NuRxs15RBQAJ -
it might be relevant to your questions.

DebyePDFGenerator uses DebyePDFCalculator to do the computations.
It calculates F(Q) and converts it to G(r) in every refinement step.
The FFT function that does the conversion is exposed to Python and
can be used on its own.
https://www.diffpy.org/diffpy.srreal/api/diffpy.srreal.html#diffpy.srreal.pdfcalculator.fftftog

Hope this helps,

Pavol

PS: My apologies ahead, I pretty much have no bandwidth to follow up here.

[He Lin]

Hi, Pavol:

     Many thanks. Yes, we noticed the last year's  discussion and download the SAScalculator as well before.  We don't quite get why the division by zero issue can be avoided using it. We will revisit the discussion and codes then.

     Best,

             he