diffpy.srfit/docs/examples/simplepdftwophase.py at c11547896b0cfe46d6bbcfb5e0fba70e6b560ddc · cadenmyers13/diffpy.srfit · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#!/usr/bin/env python
########################################################################
#
# diffpy.srfit      by DANSE Diffraction group
#                   Simon J. L. Billinge
#                   (c) 2010 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 simplified PDF refinement of two-phase structure."""

from crystalpdftwophase import plotResults
from gaussianrecipe import scipyOptimize
from pyobjcryst import loadCrystal

from diffpy.srfit.fitbase import FitRecipe, FitResults
from diffpy.srfit.pdf import PDFContribution

######
#  Example Code


def makeRecipe(niciffile, siciffile, datname):
    """Create a fitting recipe for crystalline PDF data."""

    # Load data and add it to the profile
    contribution = PDFContribution("nisi")
    contribution.loadData(datname)
    contribution.set_calculation_range(xmax=20)

    stru = loadCrystal(niciffile)
    contribution.addStructure("ni", stru)

    stru = loadCrystal(siciffile)
    contribution.addStructure("si", stru)

    # Make the FitRecipe and add the FitContribution.
    recipe = FitRecipe()
    recipe.add_contribution(contribution)

    # Configure the fit variables
    # Start by configuring the scale factor and resolution factors.
    # We want the sum of the phase scale factors to be 1.
    recipe.create_new_variable("scale_ni", 0.1)
    recipe.add_constraint(contribution.ni.scale, "scale_ni")
    recipe.add_constraint(contribution.si.scale, "1 - scale_ni")
    # We also want the resolution factor to be the same on each. This is done
    # for free by the PDFContribution. We simply need to add it to the recipe.
    recipe.add_variable(contribution.qdamp, 0.03)

    # Vary the global scale as well.
    recipe.add_variable(contribution.scale, 1)

    # Now we can configure the structural parameters. Since we're using
    # ObjCrystCrystalParSets, the space group constraints are automatically
    # applied to each phase. We must selectively vary the free parameters.
    #
    # First the nickel parameters.
    # Note that ni is the name of the PDFGenerator that was automatically
    # created by the PDFContribution. We selected this name in addStructure
    # above.
    phase_ni = contribution.ni.phase
    for par in phase_ni.sgpars:
        recipe.add_variable(par, name=par.name + "_ni")
    recipe.add_variable(contribution.ni.delta2, name="delta2_ni")
    # Next the silicon parameters
    phase_si = contribution.si.phase
    for par in phase_si.sgpars:
        recipe.add_variable(par, name=par.name + "_si")
    recipe.add_variable(contribution.si.delta2, name="delta2_si")

    # We have prior information from the earlier examples so we'll use it here
    # in the form of restraints.
    #
    # The nickel lattice parameter was measured to be 3.527. The uncertainty
    # values are invalid for that measurement, since the data from which it is
    # derived has no uncertainty. Thus, we will tell the recipe to scale the
    # residual, which means that it will be weighted as much as the average
    # data point during the fit.
    recipe.add_soft_bounds(
        "a_ni", lower_bound=3.527, upper_bound=3.527, scaled=True
    )
    # Now we do the same with the delta2 and Biso parameters (remember that
    # Biso = 8*pi**2*Uiso)
    recipe.add_soft_bounds(
        "delta2_ni", lower_bound=2.22, upper_bound=2.22, scaled=True
    )
    recipe.add_soft_bounds(
        "Biso_0_ni", lower_bound=0.454, upper_bound=0.454, scaled=True
    )
    #
    # We can do the same with the silicon values. We haven't done a thorough
    # job of measuring the uncertainties in the results, so we'll scale these
    # as well.
    recipe.add_soft_bounds(
        "a_si", lower_bound=5.430, upper_bound=5.430, scaled=True
    )
    recipe.add_soft_bounds(
        "delta2_si", lower_bound=3.54, upper_bound=3.54, scaled=True
    )
    recipe.add_soft_bounds(
        "Biso_0_si", lower_bound=0.645, upper_bound=0.645, scaled=True
    )

    # Give the recipe away so it can be used!
    return recipe


if __name__ == "__main__":

    # Make the data and the recipe
    niciffile = "data/ni.cif"
    siciffile = "data/si.cif"
    data = "data/si90ni10-q27r60-xray.gr"

    # Make the recipe
    recipe = makeRecipe(niciffile, siciffile, data)

    # Optimize
    scipyOptimize(recipe)

    # Generate and print the FitResults
    res = FitResults(recipe)
    res.print_results()

    # Plot!
    plotResults(recipe)

# End of file