#!/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 nanoparticle PDF refinement using DebyePDFGenerator.
This example is similar to crystalpdfobjcryst.py, except that it uses
the DebyePDFGenerator from SrReal to refine a pyobjcryst Molecule.
"""
from pathlib import Path
from diffpy.srfit.fitbase import (
FitContribution,
FitRecipe,
FitResults,
Profile,
)
from diffpy.srfit.pdf import DebyePDFGenerator
######
# Example Code
def makeRecipe(molecule, datname):
"""Create a recipe that uses the DebyePDFGenerator."""
# The Profile
profile = Profile()
# Load data and add it to the profile
profile.loadtxt(datname)
profile.set_calculation_range(xmin=1.2, xmax=8)
# The ProfileGenerator
# Create a DebyePDFGenerator named "G".
generator = DebyePDFGenerator("G")
generator.setStructure(molecule)
# These are metadata needed by the generator
generator.setQmin(0.68)
generator.setQmax(22)
# The FitContribution
contribution = FitContribution("bucky")
contribution.add_profile_generator(generator)
contribution.set_profile(profile, xname="r")
# Make a FitRecipe.
recipe = FitRecipe()
recipe.add_contribution(contribution)
# Specify which parameters we want to refine. We'll be using the
# MoleculeParSet within the generator, so let's get a handle to it. See the
# diffpy.srfit.structure.objcryststructure module for more information
# about the MoleculeParSet hierarchy.
c60 = generator.phase
# First, the isotropic thermal displacement factor.
Biso = recipe.create_new_variable("Biso")
for atom in c60.getScatterers():
# We have defined a 'center' atom that is a dummy, which means that it
# has no scattering power. It is only used as a reference point for
# our bond length. We don't want to constrain it.
if not atom.isDummy():
recipe.add_constraint(atom.Biso, Biso)
# We need to let the molecule expand. If we were modeling it as a crystal,
# we could let the unit cell expand. For instruction purposes, we use a
# Molecule to model C60, and molecules have different modeling options than
# crystals. To make the molecule expand from a central point, we will
# constrain the distance from each atom to a dummy center atom that was
# created with the molecule, and allow that distance to vary. (We could
# also let the nearest-neighbor bond lengths vary, but that would be much
# more difficult to set up.)
center = c60.center
# Create a new Parameter that represents the radius of the molecule. Note
# that we don't give it an initial value. Since the variable is being
# directly constrained to further below, its initial value will be inferred
# from the constraint.
radius = recipe.create_new_variable("radius")
for i, atom in enumerate(c60.getScatterers()):
if atom.isDummy():
continue
# This creates a Parameter that moves the second atom according to the
# bond length. Note that each Parameter needs a unique name.
par = c60.addBondLengthParameter("rad%i" % i, center, atom)
recipe.add_constraint(par, radius)
# Add the correlation term, scale. The scale is too short to effectively
# determine qdamp.
recipe.add_variable(generator.delta2, 2)
recipe.add_variable(generator.scale, 1.3e4)
# Give the recipe away so it can be used!
return recipe
plot_styles = {
"xlabel": r"$r (\AA)$",
"ylabel": r"$G (\AA^{-2})$",
}
def main():
molecule = makeC60()
# Make the data and the recipe
recipe = makeRecipe(molecule, Path(__file__).parent / "data/C60.gr")
# Tell the fithook that we want very verbose output.
recipe.fithooks[0].verbose = 3
# Optimize
from scipy.optimize import leastsq
leastsq(recipe.residual, recipe.get_values())
# Print results
res = FitResults(recipe)
res.print_results()
# Plot results
recipe.plot_recipe(**plot_styles)
return
c60xyz = """
3.451266498 0.685000000 0.000000000
3.451266498 -0.685000000 0.000000000
-3.451266498 0.685000000 0.000000000
-3.451266498 -0.685000000 0.000000000
0.685000000 0.000000000 3.451266498
-0.685000000 0.000000000 3.451266498
0.685000000 0.000000000 -3.451266498
-0.685000000 0.000000000 -3.451266498
0.000000000 3.451266498 0.685000000
0.000000000 3.451266498 -0.685000000
0.000000000 -3.451266498 0.685000000
0.000000000 -3.451266498 -0.685000000
3.003809890 1.409000000 1.171456608
3.003809890 1.409000000 -1.171456608
3.003809890 -1.409000000 1.171456608
3.003809890 -1.409000000 -1.171456608
-3.003809890 1.409000000 1.171456608
-3.003809890 1.409000000 -1.171456608
-3.003809890 -1.409000000 1.171456608
-3.003809890 -1.409000000 -1.171456608
1.409000000 1.171456608 3.003809890
1.409000000 -1.171456608 3.003809890
-1.409000000 1.171456608 3.003809890
-1.409000000 -1.171456608 3.003809890
1.409000000 1.171456608 -3.003809890
1.409000000 -1.171456608 -3.003809890
-1.409000000 1.171456608 -3.003809890
-1.409000000 -1.171456608 -3.003809890
1.171456608 3.003809890 1.409000000
-1.171456608 3.003809890 1.409000000
1.171456608 3.003809890 -1.409000000
-1.171456608 3.003809890 -1.409000000
1.171456608 -3.003809890 1.409000000
-1.171456608 -3.003809890 1.409000000
1.171456608 -3.003809890 -1.409000000
-1.171456608 -3.003809890 -1.409000000
2.580456608 0.724000000 2.279809890
2.580456608 0.724000000 -2.279809890
2.580456608 -0.724000000 2.279809890
2.580456608 -0.724000000 -2.279809890
-2.580456608 0.724000000 2.279809890
-2.580456608 0.724000000 -2.279809890
-2.580456608 -0.724000000 2.279809890
-2.580456608 -0.724000000 -2.279809890
0.724000000 2.279809890 2.580456608
0.724000000 -2.279809890 2.580456608
-0.724000000 2.279809890 2.580456608
-0.724000000 -2.279809890 2.580456608
0.724000000 2.279809890 -2.580456608
0.724000000 -2.279809890 -2.580456608
-0.724000000 2.279809890 -2.580456608
-0.724000000 -2.279809890 -2.580456608
2.279809890 2.580456608 0.724000000
-2.279809890 2.580456608 0.724000000
2.279809890 2.580456608 -0.724000000
-2.279809890 2.580456608 -0.724000000
2.279809890 -2.580456608 0.724000000
-2.279809890 -2.580456608 0.724000000
2.279809890 -2.580456608 -0.724000000
-2.279809890 -2.580456608 -0.724000000
"""
def makeC60():
"""Make the C60 molecule using pyobjcryst."""
from pyobjcryst.crystal import Crystal
from pyobjcryst.molecule import Molecule
from pyobjcryst.scatteringpower import ScatteringPowerAtom
# make a crystal box to put the molecule in
c = Crystal(1, 1, 1, "P1")
c.SetName("c60frame")
# put a molecule inside the box
m = Molecule(c, "c60")
c.AddScatterer(m)
# Create a dummy atom at the center.
m.AddAtom(0, 0, 0, None, "center")
# Create the scattering power object for the carbon atoms
sp = ScatteringPowerAtom("C", "C")
c.AddScatteringPower(sp)
sp.SetBiso(0.25)
# Add the other atoms. They will be named C1, C2, ..., C60.
for i, l in enumerate(c60xyz.strip().splitlines()):
x, y, z = map(float, l.split())
m.AddAtom(x, y, z, sp, "C%i" % (i + 1))
return m
if __name__ == "__main__":
main()
# End of file