From f4149a2d46731b5e05bf1bfb3196f8b2dea4aedc Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 16 Sep 2026 09:57:03 -0400 Subject: [PATCH 1/3] feat: migrate structure subpackage from diffpy.srfit Port all 8 modules of diffpy.srfit.structure into diffpy.cmipdf.structure so cmipdf owns the structure adapters it depends on, and repoint basepdfgenerator at the local copy instead of diffpy.srfit.structure. Per the deprecation plan, cmipdf has never been released and so carries no deprecation shims: constrainAsSpaceGroup and the three setConst forwarders are dropped along with the diffpy.utils._deprecator imports, and the camelCase API is ported as snake_case only. getValue and addParameter stay camelCase because they override core srfit methods that have not been renamed yet; struToParameterSet keeps its name because srfit is not renaming it. Also add extend_path to src/diffpy/__init__.py, without which diffpy.cmipdf shadows the rest of the diffpy namespace. Tests ported from srfit for diffpyparset, objcrystparset and sgconstraints, minus the constrainAsSpaceGroup deprecation test, plus the LaMnO3.stru fixture and a pyobjcryst_available gate. Suite goes from 12 to 29 passing, matching srfit's 18 for the same modules less the one deprecation test. Co-Authored-By: Claude Opus 5 --- news/add-structure-subpackage.rst | 29 + src/diffpy/__init__.py | 3 + src/diffpy/cmipdf/basepdfgenerator.py | 6 +- src/diffpy/cmipdf/debyepdfgenerator.py | 2 +- src/diffpy/cmipdf/structure/__init__.py | 66 + .../cmipdf/structure/basestructureparset.py | 65 + src/diffpy/cmipdf/structure/bvsrestraint.py | 115 + src/diffpy/cmipdf/structure/cctbxparset.py | 328 +++ src/diffpy/cmipdf/structure/diffpyparset.py | 345 +++ src/diffpy/cmipdf/structure/objcrystparset.py | 1883 +++++++++++++++++ src/diffpy/cmipdf/structure/sgconstraints.py | 827 ++++++++ src/diffpy/cmipdf/structure/srrealparset.py | 99 + tests/conftest.py | 18 + tests/test_diffpyparset.py | 134 ++ tests/test_objcrystparset.py | 678 ++++++ tests/test_sgconstraints.py | 252 +++ tests/testdata/LaMnO3.stru | 129 ++ 17 files changed, 4975 insertions(+), 4 deletions(-) create mode 100644 news/add-structure-subpackage.rst create mode 100644 src/diffpy/cmipdf/structure/__init__.py create mode 100644 src/diffpy/cmipdf/structure/basestructureparset.py create mode 100644 src/diffpy/cmipdf/structure/bvsrestraint.py create mode 100644 src/diffpy/cmipdf/structure/cctbxparset.py create mode 100644 src/diffpy/cmipdf/structure/diffpyparset.py create mode 100644 src/diffpy/cmipdf/structure/objcrystparset.py create mode 100644 src/diffpy/cmipdf/structure/sgconstraints.py create mode 100644 src/diffpy/cmipdf/structure/srrealparset.py create mode 100644 tests/test_diffpyparset.py create mode 100644 tests/test_objcrystparset.py create mode 100644 tests/test_sgconstraints.py create mode 100644 tests/testdata/LaMnO3.stru diff --git a/news/add-structure-subpackage.rst b/news/add-structure-subpackage.rst new file mode 100644 index 0000000..9507765 --- /dev/null +++ b/news/add-structure-subpackage.rst @@ -0,0 +1,29 @@ +**Added:** + +* ``diffpy.cmipdf.structure`` subpackage, migrated from ``diffpy.srfit.structure``. + It adapts ``diffpy.structure``, ``pyobjcryst`` and ``cctbx`` structures to the + ``ParameterSet`` interface and generates space group constraints. The migrated + API is snake_case only; ``getValue`` and ``addParameter`` stay camelCase until + core srfit renames them. + +**Changed:** + +* ``BasePDFGenerator`` now takes ``struToParameterSet`` from + ``diffpy.cmipdf.structure`` instead of ``diffpy.srfit.structure``. + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* ``src/diffpy/__init__.py`` now extends the ``diffpy`` namespace path, so + ``diffpy.cmipdf`` no longer shadows the other ``diffpy`` packages. + +**Security:** + +* diff --git a/src/diffpy/__init__.py b/src/diffpy/__init__.py index 43d76dc..57fc75c 100644 --- a/src/diffpy/__init__.py +++ b/src/diffpy/__init__.py @@ -12,3 +12,6 @@ # See LICENSE.rst for license information. # ############################################################################## +from pkgutil import extend_path + +__path__ = extend_path(__path__, __name__) diff --git a/src/diffpy/cmipdf/basepdfgenerator.py b/src/diffpy/cmipdf/basepdfgenerator.py index 14ea083..13716e7 100644 --- a/src/diffpy/cmipdf/basepdfgenerator.py +++ b/src/diffpy/cmipdf/basepdfgenerator.py @@ -23,10 +23,10 @@ import numpy +from diffpy.cmipdf.structure import struToParameterSet from diffpy.srfit.exceptions import SrFitError from diffpy.srfit.fitbase import ProfileGenerator from diffpy.srfit.fitbase.parameter import ParameterAdapter -from diffpy.srfit.structure import struToParameterSet # FIXME - Parameter creation will have to be smarter once deeper calculator # configuration is enabled. @@ -257,7 +257,7 @@ def set_structure(self, structure, name="phase", periodic=True): This creates a DiffpyStructureParSet, ObjCrystCrystalParSet or ObjCrystMoleculeParSet that adapts structure to a ParameterSet interface. - See those classes (located in diffpy.srfit.structure) for how they are + See those classes (located in diffpy.cmipdf.structure) for how they are used. The resulting ParameterSet will be managed by this generator. Parameters @@ -310,7 +310,7 @@ def set_structure_from_parset(self, parset, periodic=True): self.add_parameter_set(parset) # Set periodicity - self._phase.useSymmetry(periodic) + self._phase.use_symmetry(periodic) return def _prepare(self, r): diff --git a/src/diffpy/cmipdf/debyepdfgenerator.py b/src/diffpy/cmipdf/debyepdfgenerator.py index 1825b5d..9df4493 100644 --- a/src/diffpy/cmipdf/debyepdfgenerator.py +++ b/src/diffpy/cmipdf/debyepdfgenerator.py @@ -93,7 +93,7 @@ def set_structure(self, structure, name="phase", periodic=False): This creates a DiffpyStructureParSet, ObjCrystCrystalParSet or ObjCrystMoleculeParSet that adapts structure to a ParameterSet interface. - See those classes (located in diffpy.srfit.structure) for how they are + See those classes (located in diffpy.cmipdf.structure) for how they are used. The resulting ParameterSet will be managed by this generator. Parameters diff --git a/src/diffpy/cmipdf/structure/__init__.py b/src/diffpy/cmipdf/structure/__init__.py new file mode 100644 index 0000000..77531b5 --- /dev/null +++ b/src/diffpy/cmipdf/structure/__init__.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Modules and classes that adapt structure representations to the +ParameterSet interface and automatic structure constraint generation +from space group information.""" + +from diffpy.cmipdf.structure.sgconstraints import constrain_as_space_group + +__all__ = ["constrain_as_space_group", "struToParameterSet"] + + +def struToParameterSet(name, stru): + """Creates a ParameterSet from an structure. + + This returns a ParameterSet adapted for the structure depending on its + type. + + Parameters + ---------- + stru + a structure object known by this module + name + A name to give the structure. + + Raises TypeError if stru cannot be adapted + """ + from diffpy.cmipdf.structure.diffpyparset import DiffpyStructureParSet + + if DiffpyStructureParSet.can_adapt(stru): + return DiffpyStructureParSet(name, stru) + + from diffpy.cmipdf.structure.objcrystparset import ObjCrystCrystalParSet + + if ObjCrystCrystalParSet.can_adapt(stru): + return ObjCrystCrystalParSet(name, stru) + + from diffpy.cmipdf.structure.objcrystparset import ObjCrystMoleculeParSet + + if ObjCrystMoleculeParSet.can_adapt(stru): + return ObjCrystMoleculeParSet(name, stru) + + from diffpy.cmipdf.structure.cctbxparset import CCTBXCrystalParSet + + if CCTBXCrystalParSet.can_adapt(stru): + return CCTBXCrystalParSet(name, stru) + + raise TypeError("Unadaptable structure format") + + +# silence pyflakes checker +assert constrain_as_space_group + +# End of file diff --git a/src/diffpy/cmipdf/structure/basestructureparset.py b/src/diffpy/cmipdf/structure/basestructureparset.py new file mode 100644 index 0000000..48196ed --- /dev/null +++ b/src/diffpy/cmipdf/structure/basestructureparset.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Base class for adapting structures to a ParameterSet interface. + +The BaseStructureParSet is a ParameterSet with functionality required by +all structure adapters. +""" + +__all__ = ["BaseStructureParSet"] + +from diffpy.srfit.fitbase.parameterset import ParameterSet + + +class BaseStructureParSet(ParameterSet): + """Base class for structure adapters. + + BaseStructureParSet derives from ParameterSet and provides methods that + help interface the ParameterSet with the space group constraint methods in + the sgconstraints module and to ProfileGenerators. + + Attributes + ---------- + stru + The adapted object + """ + + @classmethod + def can_adapt(self, stru): + """Return whether the structure can be adapted by this class.""" + return False + + def get_lattice(self): + """Get a ParameterSet containing the lattice Parameters. + + The returned ParameterSet may contain other Parameters than the + lattice Parameters. It is assumed that the lattice parameters + are named "a", "b", "c", "alpha", "beta", "gamma". + + Lattice must also have the "angunits" attribute, which is either + "deg" or "rad", to signify degrees or radians. + """ + raise NotImplementedError("The must be overloaded") + + def get_scatterers(self): + """Get a list of ParameterSets that represents the scatterers. + + The site positions must be accessible from the list entries via + the names "x", "y", and "z". The ADPs must be accessible as + well, but the name and nature of the ADPs (U-factors, B-factors, + isotropic, anisotropic) depends on the adapted structure. + """ + raise NotImplementedError("The must be overloaded") diff --git a/src/diffpy/cmipdf/structure/bvsrestraint.py b/src/diffpy/cmipdf/structure/bvsrestraint.py new file mode 100644 index 0000000..c62fc23 --- /dev/null +++ b/src/diffpy/cmipdf/structure/bvsrestraint.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Bond-valence sum calculator from SrReal wrapped as a Restraint. + +This can be used as an addition to a cost function during a structure +refinement to keep the bond-valence sum within tolerable limits. +""" + +__all__ = ["BVSRestraint"] + +from diffpy.srfit.exceptions import SrFitError +from diffpy.srfit.fitbase.restraint import Restraint + + +class BVSRestraint(Restraint): + """Wrapping of BVSCalculator.bvmsdiff as a Restraint. + + The restraint penalty is the root-mean-square deviation of the theoretical + and calculated bond-valence sum of a structure. + + Attributes + ---------- + _calc + The SrReal BVSCalculator instance. + _parset + The SrRealParSet that created this BVSRestraint. + sig + The uncertainty on the BVS (default 1). + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + + def __init__(self, parset, sig=1, scaled=False): + """Initialize the Restraint. + + Parameters + ---------- + parset + SrRealParSet that creates this BVSRestraint. + sig + The uncertainty on the BVS (default 1). + scaled + A flag indicating if the restraint is scaled + (multiplied) by the unrestrained point-average chi^2 + (chi^2/numpoints) (bool, default False). + """ + from diffpy.srreal.bvscalculator import BVSCalculator + + self._calc = BVSCalculator() + self._parset = parset + self.sig = float(sig) + self.scaled = bool(scaled) + return + + def penalty(self, w=1.0): + """Calculate the penalty of the restraint. + + Parameters + ---------- + w + The point-average chi^2 which is optionally used to scale the + penalty (float, default 1.0). + """ + # Get the bvms from the BVSCalculator + stru = self._parset._get_srreal_structure() + self._calc.eval(stru) + penalty = self._calc.bvmsdiff + + # Scale by the prefactor + penalty /= self.sig**2 + + # Optionally scale by w + if self.scaled: + penalty *= w + + return penalty + + def _validate(self): + """This evaluates the calculator. + + Raises SrFitError if validation fails. + """ + from numpy import nan + + p = self.penalty() + if p is None or p is nan: + raise SrFitError("Cannot evaluate penalty") + v = self._calc.value + if len(v) > 1 and not v.any(): + emsg = ( + "Bond valence sums are all zero. Check atom symbols in " + "the structure or define custom bond-valence parameters." + ) + raise SrFitError(emsg) + return + + # End of class BVSRestraint + + +# End of file diff --git a/src/diffpy/cmipdf/structure/cctbxparset.py b/src/diffpy/cmipdf/structure/cctbxparset.py new file mode 100644 index 0000000..c799e0e --- /dev/null +++ b/src/diffpy/cmipdf/structure/cctbxparset.py @@ -0,0 +1,328 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Wrappers for interfacing cctbx crystal with SrFit. + +This wraps a cctbx.crystal as a ParameterSet with a similar hierarchy, which +can then be used within a FitRecipe. Note that all manipulations to the +cctbx.crystal should be done before wrapping. Changes made to the cctbx.crystal +object after wrapping may not be reflected within the wrapper, which can have +unpredictable results during a structure refinement. + +The following classes are adapted: + +- `CCTBXCrystalParSet`: wrapper for `cctbx.crystal`. +- `CCTBXUnitCellParSet`: wrapper for the unit cell of `cctbx.crystal`. +- `CCTBXScattererParSet`: wrapper for `cctbx.xray.scatterer`. +""" + +from diffpy.cmipdf.structure.basestructureparset import BaseStructureParSet +from diffpy.srfit.fitbase.parameter import ParameterAdapter +from diffpy.srfit.fitbase.parameterset import ParameterSet + +__all__ = ["CCTBXScattererParSet", "CCTBXUnitCellParSet", "CCTBXCrystalParSet"] + + +class CCTBXScattererParSet(ParameterSet): + """A wrapper for cctbx.xray.scatterer. + + This class derives from ParameterSet. + + Attributes + ---------- + name + Name of the scatterer. The name is always of the form + "%s%i" % (element, number), where the number is the running + index of that element type (starting at 0). + x (y, z) -- Atom position in crystal coordinates (ParameterAdapter) + occupancy + Occupancy of the atom on its crystal location + (ParameterAdapter) + Uiso + Isotropic scattering factor (ParameterAdapter). + """ + + def __init__(self, name, strups, idx): + """Initialize. + + Parameters + ---------- + name + The name of this scatterer. + strups + The CCTBXCrystalParSet that contains the cctbx structure + idx + The index of the scatterer in the structure. + """ + ParameterSet.__init__(self, name) + self.strups = strups + self.idx = idx + + # x, y, z, occupancy + self.addParameter( + ParameterAdapter("x", None, self._xyzgetter(0), self._xyzsetter(0)) + ) + self.addParameter( + ParameterAdapter("y", None, self._xyzgetter(1), self._xyzsetter(1)) + ) + self.addParameter( + ParameterAdapter("z", None, self._xyzgetter(2), self._xyzsetter(2)) + ) + self.addParameter( + ParameterAdapter("occupancy", None, self._getocc, self._setocc) + ) + self.addParameter( + ParameterAdapter("Uiso", None, self._getuiso, self._setuiso) + ) + return + + # Getters and setters + + def _xyzgetter(self, i): + + def f(dummy): + return self.strups.stru.scatterers()[self.idx].site[i] + + return f + + def _xyzsetter(self, i): + + def f(dummy, value): + xyz = list(self.strups.stru.scatterers()[self.idx].site) + xyz[i] = value + self.strups.stru.scatterers()[self.idx].site = tuple(xyz) + return + + return f + + def _getocc(self, dummy): + return self.strups.stru.scatterers()[self.idx].occupancy + + def _setocc(self, dummy, value): + self.strups.stru.scatterers()[self.idx].occupancy = value + return + + def _getuiso(self, dummy): + return self.strups.stru.scatterers()[self.idx].u_iso + + def _setuiso(self, dummy, value): + self.strups.stru.scatterers()[self.idx].u_iso = value + return + + def _getelem(self): + return self.stru.element_symbol() + + element = property(_getelem) + + +# End class CCTBXScattererParSet + + +class CCTBXUnitCellParSet(ParameterSet): + """A wrapper for cctbx unit_cell object. + + Attributes + ---------- + name + Always "unitcell". + a + Unit cell parameters (ParameterAdapter). + b + Unit cell parameters (ParameterAdapter). + c + Unit cell parameters (ParameterAdapter). + alpha + Unit cell parameters (ParameterAdapter). + beta + Unit cell parameters (ParameterAdapter). + gamma + Unit cell parameters (ParameterAdapter). + """ + + def __init__(self, strups): + """Initialize. + + Parameters + ---------- + strups + The CCTBXCrystalParSet that contains the cctbx structure + and the unit cell we're wrapper. + """ + ParameterSet.__init__(self, "unitcell") + self.strups = strups + self._latpars = list(self.strups.stru.unit_cell().parameters()) + + self.addParameter( + ParameterAdapter("a", None, self._latgetter(0), self._latsetter(0)) + ) + self.addParameter( + ParameterAdapter("b", None, self._latgetter(1), self._latsetter(1)) + ) + self.addParameter( + ParameterAdapter("c", None, self._latgetter(2), self._latsetter(2)) + ) + self.addParameter( + ParameterAdapter( + "alpha", None, self._latgetter(3), self._latsetter(3) + ) + ) + self.addParameter( + ParameterAdapter( + "beta", None, self._latgetter(4), self._latsetter(4) + ) + ) + self.addParameter( + ParameterAdapter( + "gamma", None, self._latgetter(5), self._latsetter(5) + ) + ) + + return + + def _latgetter(self, i): + + def f(dummy): + return self._latpars[i] + + return f + + def _latsetter(self, i): + + def f(dummy, value): + self._latpars[i] = value + self.strups._update = True + return + + return f + + +# End class CCTBXUnitCellParSet + +# FIXME - Special positions should be constant. + + +class CCTBXCrystalParSet(BaseStructureParSet): + """A wrapper for CCTBX structure. + + Attributes + ---------- + stru + The adapted cctbx structure object. + scatterers + The list of ScattererParSets. + unitcell + The CCTBXUnitCellParSet for the structure. + """ + + def __init__(self, name, stru): + """Initialize. + + Parameters + ---------- + name + A name for this + stru + A CCTBX structure instance. + """ + ParameterSet.__init__(self, name) + self.stru = stru + self.add_parameter_set(CCTBXUnitCellParSet(self)) + self.scatterers = [] + + self._update = False + + cdict = {} + for s in stru.scatterers(): + el = s.element_symbol() + i = cdict.get(el, 0) + sname = "%s%i" % (el, i) + cdict[el] = i + 1 + scatterer = CCTBXScattererParSet(sname, self, i) + self.add_parameter_set(scatterer) + self.scatterers.append(scatterer) + + # Constrain the lattice + from diffpy.cmipdf.structure.sgconstraints import ( + _constrain_space_group, + ) + + symbol = self.get_space_group() + _constrain_space_group(self, symbol) + + return + + def update(self): + """Update the unit_cell to a change in lattice parameters. + + This remakes the unit cell according to a change in the lattice + parameters. Call this function before using the + CCTBXCrystalParSet. The unit_cell will only be remade if + necessary. + """ + if not self._update: + return + + self._update = False + stru = self.stru + sgn = stru.space_group().match_tabulated_settings().number() + + # Create the symmetry object + from cctbx.crystal import symmetry + + symm = symmetry( + unit_cell=self.unitcell._latpars, space_group_symbol=sgn + ) + + # Now the new structure + newstru = stru.__class__( + crystal_symmetry=symm, scatterers=stru.scatterers() + ) + + self.unitcell._latpars = list(newstru.unit_cell().parameters()) + + self.stru = newstru + return + + @classmethod + def can_adapt(self, stru): + """Return whether the structure can be adapted by this class.""" + try: + from cctbx.crystal import special_position_settings + except ImportError: + return False + return isinstance(stru, special_position_settings) + + def get_lattice(self): + """Get the ParameterSet containing the lattice Parameters.""" + return self.unitcell + + def get_scatterers(self): + """Get a list of ParameterSets that represents the scatterers. + + The site positions must be accessible from the list entries via + the names "x", "y", and "z". The ADPs must be accessible as + well, but the name and nature of the ADPs (U-factors, B-factors, + isotropic, anisotropic) depends on the adapted structure. + """ + return self.scatterers + + def get_space_group(self): + """Get the HM space group symbol for the structure.""" + sg = self.stru.space_group() + t = sg.type() + return t.lookup_symbol() + + +# End class CCTBXCrystalParSet diff --git a/src/diffpy/cmipdf/structure/diffpyparset.py b/src/diffpy/cmipdf/structure/diffpyparset.py new file mode 100644 index 0000000..54f74d9 --- /dev/null +++ b/src/diffpy/cmipdf/structure/diffpyparset.py @@ -0,0 +1,345 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Adapters for interfacing a diffpy.structure.Structure with SrFit. + +A diffpy.structure.Structure object is meant to be passed to a +DiffpyStructureParSet object from this module, which can then be used as a +ParameterSet. (It has other methods for interfacing with SrReal calculator +adapters.) Any change to the lattice or existing atoms will be registered with +the Structure. Changes in the number of atoms will not be recognized. Thus, +the diffpy.structure.Structure object should be fully configured before passing +it to DiffpyStructureParSet. + +The following classes are adapted: + +- `DiffpyStructureParSet`: adapter for `diffpy.structure.Structure`. +- `DiffpyLatticeParSet`: adapter for `diffpy.structure.Lattice`. +- `DiffpyAtomParSet`: adapter for `diffpy.structure.Atom`. +""" + +__all__ = ["DiffpyStructureParSet"] + +from diffpy.cmipdf.structure.srrealparset import SrRealParSet +from diffpy.srfit.fitbase.parameter import ParameterAdapter, ParameterProxy +from diffpy.srfit.fitbase.parameterset import ParameterSet +from diffpy.srfit.util.argbinders import bind2nd + + +# Accessor for xyz of atoms +class _xyzgetter(object): + + def __init__(self, i): + self.i = i + + def __call__(self, atom): + return atom.xyz[self.i] + + +class _xyzsetter(object): + + def __init__(self, i): + self.i = i + + def __call__(self, atom, value): + atom.xyz[self.i] = value + + +class DiffpyAtomParSet(ParameterSet): + """A wrapper for diffpy.structure.Atom. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. See + this class for base attributes. + + Attributes + ---------- + atom + The diffpy.structure.Atom this is adapting + element + The element name (property). + + Managed Parameters + ------------------ + occupancy + Occupancy of the atom on its crystal location + (ParameterAdapter) + occ + Proxy for occupancy (ParameterProxy). + U11, U22, U33, U12, U21, U23, U32, U13, U31 + -- Anisotropic displacement factor for atom (ParameterAdapter + or ParameterProxy). Note that the Uij and Uji parameters + are the same. + Uiso + Isotropic ADP (ParameterAdapter). + B11, B22, B33, B12, B21, B23, B32, B13, B31 + -- Anisotropic displacement factor for atom (ParameterAdapter + or ParameterProxy). Note that the Bij and Bji parameters + are the same. (Bij = 8*pi**2*Uij) + Biso + Isotropic ADP (ParameterAdapter). + """ + + def __init__(self, name, atom): + """Initialize. + + Parameters + ---------- + atom + A diffpy.structure.Atom instance + """ + ParameterSet.__init__(self, name) + self.atom = atom + a = atom + # x, y, z, occupancy + self.addParameter( + ParameterAdapter("x", a, _xyzgetter(0), _xyzsetter(0)) + ) + self.addParameter( + ParameterAdapter("y", a, _xyzgetter(1), _xyzsetter(1)) + ) + self.addParameter( + ParameterAdapter("z", a, _xyzgetter(2), _xyzsetter(2)) + ) + occupancy = ParameterAdapter("occupancy", a, attr="occupancy") + self.addParameter(occupancy) + self.addParameter(ParameterProxy("occ", occupancy)) + # U + self.addParameter(ParameterAdapter("U11", a, attr="U11")) + self.addParameter(ParameterAdapter("U22", a, attr="U22")) + self.addParameter(ParameterAdapter("U33", a, attr="U33")) + U12 = ParameterAdapter("U12", a, attr="U12") + U21 = ParameterProxy("U21", U12) + U13 = ParameterAdapter("U13", a, attr="U13") + U31 = ParameterProxy("U31", U13) + U23 = ParameterAdapter("U23", a, attr="U23") + U32 = ParameterProxy("U32", U23) + self.addParameter(U12) + self.addParameter(U21) + self.addParameter(U13) + self.addParameter(U31) + self.addParameter(U23) + self.addParameter(U32) + self.addParameter(ParameterAdapter("Uiso", a, attr="Uisoequiv")) + # B + self.addParameter(ParameterAdapter("B11", a, attr="B11")) + self.addParameter(ParameterAdapter("B22", a, attr="B22")) + self.addParameter(ParameterAdapter("B33", a, attr="B33")) + B12 = ParameterAdapter("B12", a, attr="B12") + B21 = ParameterProxy("B21", B12) + B13 = ParameterAdapter("B13", a, attr="B13") + B31 = ParameterProxy("B31", B13) + B23 = ParameterAdapter("B23", a, attr="B23") + B32 = ParameterProxy("B32", B23) + self.addParameter(B12) + self.addParameter(B21) + self.addParameter(B13) + self.addParameter(B31) + self.addParameter(B23) + self.addParameter(B32) + self.addParameter(ParameterAdapter("Biso", a, attr="Bisoequiv")) + return + + def __repr__(self): + return repr(self.atom) + + def _getelem(self): + return self.atom.element + + def _setelem(self, el): + self.atom.element = el + + element = property(_getelem, _setelem, "type of atom") + + +# End class DiffpyAtomParSet + + +def _latgetter(par): + return bind2nd(getattr, par) + + +def _latsetter(par): + return bind2nd(setattr, par) + + +class DiffpyLatticeParSet(ParameterSet): + """A wrapper for diffpy.structure.Lattice. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. + See this class for base attributes. + + Attributes + ---------- + lattice + The diffpy.structure.Lattice this is adapting + name + Always "lattice" + angunits + "deg", the units of angle + + Parameters + ---------- + a + Unit cell parameters (ParameterAdapter). + b + Unit cell parameters (ParameterAdapter). + c + Unit cell parameters (ParameterAdapter). + alpha + Unit cell parameters (ParameterAdapter). + beta + Unit cell parameters (ParameterAdapter). + gamma + Unit cell parameters (ParameterAdapter). + """ + + def __init__(self, lattice): + """Initialize. + + Parameters + ---------- + lattice + A diffpy.structure.Lattice instance + """ + ParameterSet.__init__(self, "lattice") + self.angunits = "deg" + self.lattice = lattice + lat = lattice + self.addParameter( + ParameterAdapter("a", lat, _latgetter("a"), _latsetter("a")) + ) + self.addParameter( + ParameterAdapter("b", lat, _latgetter("b"), _latsetter("b")) + ) + self.addParameter( + ParameterAdapter("c", lat, _latgetter("c"), _latsetter("c")) + ) + self.addParameter( + ParameterAdapter( + "alpha", lat, _latgetter("alpha"), _latsetter("alpha") + ) + ) + self.addParameter( + ParameterAdapter( + "beta", lat, _latgetter("beta"), _latsetter("beta") + ) + ) + self.addParameter( + ParameterAdapter( + "gamma", lat, _latgetter("gamma"), _latsetter("gamma") + ) + ) + return + + def __repr__(self): + return repr(self.lattice) + + +# End class DiffpyLatticeParSet + + +class DiffpyStructureParSet(SrRealParSet): + """A wrapper for diffpy.structure.Structure. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. See + this class for base attributes. + + Attributes + ---------- + atoms + The list of DiffpyAtomParSets, provided for convenience. + stru + The diffpy.structure.Structure this is adapting + + Managed ParameterSets + --------------------- + lattice + The managed DiffpyLatticeParSet + + A managed DiffpyAtomParSets. is the atomic element and + is the index of that element in the structure, + starting from zero. Thus, for nickel in P1 symmetry, the + managed DiffpyAtomParSets will be named "Ni0", "Ni1", "Ni2" + and "Ni3". + """ + + def __init__(self, name, stru): + """Initialize. + + Parameters + ---------- + name + A name for the structure + stru + A diffpy.structure.Structure instance + """ + SrRealParSet.__init__(self, name) + self.stru = stru + self.add_parameter_set(DiffpyLatticeParSet(stru.lattice)) + self.atoms = [] + + cdict = {} + for a in stru: + el = a.element.title() + # Try to sanitize the name. + el = el.replace("+", "p") + el = el.replace("-", "m") + i = cdict.get(el, 0) + aname = "%s%i" % (el, i) + cdict[el] = i + 1 + atom = DiffpyAtomParSet(aname, a) + self.add_parameter_set(atom) + self.atoms.append(atom) + + return + + def __repr__(self): + return repr(self.stru) + + def get_lattice(self): + """Get the ParameterSet containing the lattice Parameters.""" + return self.lattice + + @classmethod + def can_adapt(self, stru): + """Return whether the structure can be adapted by this class.""" + from diffpy.structure import Structure + + return isinstance(stru, Structure) + + def get_scatterers(self): + """Get a list of ParameterSets that represents the scatterers. + + The site positions must be accessible from the list entries via + the names "x", "y", and "z". The ADPs must be accessible as + well, but the name and nature of the ADPs (U-factors, B-factors, + isotropic, anisotropic) depends on the adapted structure. + """ + return self.atoms + + def _get_srreal_structure(self): + """Get the structure object for use with SrReal calculators. + + If this is periodic, then return the structure, otherwise, pass + it inside of a nosymmetry wrapper. This takes the extra step of + wrapping the structure in a nometa wrapper. + """ + from diffpy.srreal.structureadapter import nometa + + stru = SrRealParSet._get_srreal_structure(self) + return nometa(stru) + + +# End class DiffpyStructureParSet diff --git a/src/diffpy/cmipdf/structure/objcrystparset.py b/src/diffpy/cmipdf/structure/objcrystparset.py new file mode 100644 index 0000000..5924e86 --- /dev/null +++ b/src/diffpy/cmipdf/structure/objcrystparset.py @@ -0,0 +1,1883 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Wrappers for adapting pyobjcryst.crystal.Crystal to a srfit +ParameterSet. + +This will adapt a Crystal or Molecule object from pyobjcryst into the +ParameterSet interface. The following classes are adapted: + +- `ObjCrystCrystalParSet`: adapter for `pyobjcryst.crystal.Crystal`. +- `ObjCrystAtomParSet`: adapter for `pyobjcryst.atom.Atom`. +- `ObjCrystMoleculeParSet`: adapter for `pyobjcryst.molecule.Molecule`. +- `ObjCrystMolAtomParSet`: adapter for `pyobjcryst.molecule.MolAtom`. + +Related to the adaptation of Molecule and MolAtom, there are adaptors +for specifying molecule restraints: + +- `ObjCrystBondLengthRestraint` +- `ObjCrystBondAngleRestraint` +- `ObjCrystDihedralAngleRestraint` + +There are also Parameters for encapsulating and modifying atoms via +their relative positions. These Parameters can also act like +constraints, and can modify the positions of multiple MolAtoms: + +- `ObjCrystBondLengthParameter` +- `ObjCrystBondAngleParameter` +- `ObjCrystDihedralAngleParameter` +""" + +__all__ = ["ObjCrystMoleculeParSet", "ObjCrystCrystalParSet"] + +import numpy +from pyobjcryst.molecule import ( + GetBondAngle, + GetBondLength, + GetDihedralAngle, + StretchModeBondAngle, + StretchModeBondLength, + StretchModeTorsion, +) + +from diffpy.cmipdf.structure.srrealparset import SrRealParSet +from diffpy.srfit.fitbase.parameter import ( + Parameter, + ParameterAdapter, + ParameterProxy, +) +from diffpy.srfit.fitbase.parameterset import ParameterSet + + +class ObjCrystScattererParSet(ParameterSet): + """A base adaptor for an Objcryst Scatterer. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet and + adapts pyobjcryst.scatterer.Scatterer derivatives (Molecule, Atom) and + objects with a similar interface (MolAtom). See the ParameterSet class for + base attributes. + + Attributes + ---------- + scat + The adapted pyobjcryst object. + parent + The ParameterSet this belongs to + + Managed Parameters + ------------------ + occ + Occupancy of the scatterer on its crystal site + (ParameterWraper) + """ + + def __init__(self, name, scat, parent): + """Initialize. + + Parameters + ---------- + name + The name of the scatterer + scat + The pyobjcryst.Scatterer instance + parent + The ParameterSet this belongs to + """ + ParameterSet.__init__(self, name) + self.scat = scat + self.parent = parent + + # x, y, z, occ + self.addParameter(ParameterAdapter("x", self.scat, attr="X")) + self.addParameter(ParameterAdapter("y", self.scat, attr="Y")) + self.addParameter(ParameterAdapter("z", self.scat, attr="Z")) + self.addParameter(ParameterAdapter("occ", self.scat, attr="Occupancy")) + return + + def is_dummy(self): + """Indicate whether this scatterer is a dummy atom.""" + return False + + def has_scatterers(self): + """Indicate if this scatterer has its own scatterers.""" + return hasattr(self, "get_scatterers") + + +# End class ObjCrystScattererParSet + + +class ObjCrystAtomParSet(ObjCrystScattererParSet): + """A adaptor for a pyobjcryst.Atom. + + This class derives from ObjCrystScattererParSet. + + Attributes + ---------- + scat + The adapted pyobjcryst.atom.Atom. + element + Non-refinable name of the element (property). + parent + The ObjCrystCrystalParSet this belongs to. + + Managed Parameters + ------------------ + occ + Occupancy of the atom on its crystal location + (ParameterAdapter) + Biso + Isotropic scattering factor (ParameterAdapter). + B11, B22, B33, B12, B21, B23, B32, B13, B31 + -- Anisotropic displacement factor for scatterer + (ParameterAdapter or ParameterProxy). Note that the Bij and Bji + parameters are the same. + """ + + def __init__(self, name, atom, parent): + """Initialize. + + Parameters + ---------- + name + The name of the scatterer + scat + The Scatterer instance + parent + The ObjCrystCrystalParSet this belongs to + """ + ObjCrystScattererParSet.__init__(self, name, atom, parent) + sp = atom.GetScatteringPower() + + # The B-parameters + self.addParameter(ParameterAdapter("Biso", sp, attr="Biso")) + self.addParameter(ParameterAdapter("B11", sp, attr="B11")) + self.addParameter(ParameterAdapter("B22", sp, attr="B22")) + self.addParameter(ParameterAdapter("B33", sp, attr="B33")) + B12 = ParameterAdapter("B12", sp, attr="B12") + B21 = ParameterProxy("B21", B12) + B13 = ParameterAdapter("B13", sp, attr="B13") + B31 = ParameterProxy("B31", B13) + B23 = ParameterAdapter("B23", sp, attr="B23") + B32 = ParameterProxy("B32", B23) + self.addParameter(B12) + self.addParameter(B21) + self.addParameter(B13) + self.addParameter(B31) + self.addParameter(B23) + self.addParameter(B32) + + # Give a value to Biso if it doesn't have one, and this is isotropic + if sp.IsIsotropic() and self.Biso.value == 0: + self.Biso.value = 0.5 + return + + def _getelem(self): + """Getter for the element type.""" + return self.scat.GetScatteringPower().GetSymbol() + + element = property(_getelem) + + +# End class ObjCrystAtomParSet + + +class ObjCrystMoleculeParSet(ObjCrystScattererParSet): + """A adaptor for a pyobjcryst.Molecule. + + This class derives from ObjCrystScattererParSet. + + Attributes + ---------- + scat + The adapted pyobjcryst.molecule.Molecule. + stru + The adapted pyobjcryst.molecule.Molecule. + parent + The ObjCrystCrystalParSet this belongs to. + ObjCrystMoleculeParSets can be used on their own, in which + case this is None. + + Managed Parameters + ------------------ + occ + Occupancy of the molecule on its crystal location + (ParameterAdapter) + q0, q1, q2, q3 -- Orientational quaternion (ParameterAdapter) + + + Other attributes are inherited from + diffpy.srfit.fitbase.parameterset.ParameterSet + """ + + def __init__(self, name, molecule, parent=None): + """Initialize. + + Parameters + ---------- + name + The name of the scatterer + molecule + The pyobjcryst.Molecule instance + parent + The ObjCrystCrystalParSet this belongs to (default None). + """ + ObjCrystScattererParSet.__init__(self, name, molecule, parent) + self.stru = molecule + + # Add orientation quaternion + self.addParameter(ParameterAdapter("q0", self.scat, attr="Q0")) + self.addParameter(ParameterAdapter("q1", self.scat, attr="Q1")) + self.addParameter(ParameterAdapter("q2", self.scat, attr="Q2")) + self.addParameter(ParameterAdapter("q3", self.scat, attr="Q3")) + + # Wrap the MolAtoms within the molecule + self.atoms = [] + anames = [] + + for a in molecule: + + name = a.GetName() + if not name: + raise AttributeError("Each MolAtom must have a name") + if name in anames: + raise AttributeError("MolAtom name '%s' is duplicated" % name) + + atom = ObjCrystMolAtomParSet(name, a, self) + atom.molecule = self + self.add_parameter_set(atom) + self.atoms.append(atom) + anames.append(name) + + return + + @classmethod + def can_adapt(self, stru): + """Return whether the structure can be adapted by this class.""" + from pyobjcryst.molecule import Molecule + + return isinstance(stru, Molecule) + + # Part of SrRealParSet interface + def use_symmetry(self, use=True): + """Set this structure to use symmetry. + + This structure object does not support symmetry. + """ + return + + # Part of SrRealParSet interface + def using_symmetry(self): + """Check if symmetry is being used. + + This structure object does not support symmetry. + """ + return False + + # Part of SrRealParSet interface + def _get_srreal_structure(self): + """Get the structure object for use with SrReal calculators. + + Molecule objects are never periodic. Return the object and let + the SrReal adapters do the proper thing. + """ + return self.stru + + def get_lattice(self): + """Get the ParameterSet containing the lattice Parameters.""" + lattice = ParameterSet("lattice") + lattice.newPar("a", 1.0) + lattice.newPar("b", 1.0) + lattice.newPar("c", 1.0) + lattice.newPar("alpha", 90) + lattice.newPar("beta", 90) + lattice.newPar("gamma", 90) + lattice.angunits = "deg" + return lattice + + def get_scatterers(self): + """Get a list of ParameterSets that represents the scatterers. + + The site positions must be accessible from the list entries via + the names "x", "y", and "z". The ADPs must be accessible as + well, but the name and nature of the ADPs (U-factors, B-factors, + isotropic, anisotropic) depends on the adapted structure. + """ + return self.atoms + + def wrap_restraints(self): + """Wrap the restraints implicit to the molecule. + + This will wrap MolBonds, MolBondAngles and MolDihedralAngles of + the Molecule as ObjCrystMoleculeRestraint objects. + """ + # Wrap restraints. Restraints wrapped in this way cannot be modified + # from within this class. + for b in self.scat.GetBondList(): + res = ObjCrystMoleculeRestraint(b) + self._restraints.add(res) + + for ba in self.scat.GetBondAngleList(): + res = ObjCrystMoleculeRestraint(ba) + self._restraints.add(res) + + for da in self.scat.GetDihedralAngleList(): + res = ObjCrystMoleculeRestraint(da) + self._restraints.add(res) + + return + + def wrap_stretch_mode_parameters(self): + """Wrap the stretch modes implicit to the Molecule as + Parameters. + + This will wrap StretchModeBondLengths and StretchModeBondAngles of the + Molecule as Parameters. Note that this requires that the MolBondAtoms + in the Molecule came in with unique names. Torsion angles are not + wrapped, as there is not enough information to determine each MolAtom + in the angle. + + The Parameters will be given the concatenated name of its constituents. + bond lengths: "bl_aname1_aname2" + bond angles: "ba_aname1_aname2_aname3" + """ + for mode in self.scat.GetStretchModeBondLengthList(): + name1 = mode.mpAtom0.GetName() + name2 = mode.mpAtom1.GetName() + + name = "bl_" + "_".join((name1, name2)) + + atom1 = getattr(self, name1) + atom2 = getattr(self, name2) + + par = ObjCrystBondLengthParameter(name, atom1, atom2, mode=mode) + + atoms = [] + for a in mode.GetAtoms(): + name = a.GetName() + atoms.append(getattr(self, name)) + + par.AddAtoms(atoms) + + self.addParameter(par) + + for mode in self.scat.GetStretchModeBondAngleList(): + name1 = mode.mpAtom0.GetName() + name2 = mode.mpAtom1.GetName() + name3 = mode.mpAtom2.GetName() + + name = "ba_" + "_".join((name1, name2, name3)) + + atom1 = getattr(self, name1) + atom2 = getattr(self, name2) + atom3 = getattr(self, name3) + + par = ObjCrystBondAngleParameter( + name, atom1, atom2, atom3, mode=mode + ) + + atoms = [] + for a in mode.GetAtoms(): + name = a.GetName() + atoms.append(getattr(self, name)) + par.AddAtoms(atoms) + + self.addParameter(par) + + return + + def restrain_bond_length( + self, atom1, atom2, length, sigma, delta, scaled=False + ): + """Add a bond length restraint. + + This creates an instance of ObjCrystBondLengthRestraint and adds it to + the ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the bond + atom2 + Second atom (ObjCrystMolAtomParSet) in the bond + length + The length of the bond (Angstroms) + sigma + The uncertainty of the bond length (Angstroms) + delta + The width of the bond (Angstroms) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False) + + Returns + ------- + res + The ObjCrystBondLengthRestraint object for use with the + 'unrestrain' method. + """ + res = ObjCrystBondLengthRestraint( + atom1, atom2, length, sigma, delta, scaled + ) + self._restraints.add(res) + + return res + + def restrain_bond_length_parameter( + self, par, length, sigma, delta, scaled=False + ): + """Add a bond length restraint. + + This creates an instance of ObjCrystBondLengthRestraint and adds it to + the ObjCrystMoleculeParSet. + + Parameters + ---------- + par + A ObjCrystBondLengthParameter (see add_bond_length_parameter) + length + The length of the bond (Angstroms) + sigma + The uncertainty of the bond length (Angstroms) + delta + The width of the bond (Angstroms) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False) + + Returns + ------- + ObjCrystBondLengthRestraint object + Returns the ObjCrystBondLengthRestraint object for use with the + 'unrestrain' method. + """ + return self.restrain_bond_length( + par.atom1, par.atom2, length, sigma, delta, scaled + ) + + def restrain_bond_angle( + self, atom1, atom2, atom3, angle, sigma, delta, scaled=False + ): + """Add a bond angle restraint. + + This creates an instance of ObjCrystBondAngleRestraint and adds it to + the ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the bond angle + atom2 + Second (central) atom (ObjCrystMolAtomParSet) in the bond + angle + atom3 + Third atom (ObjCrystMolAtomParSet) in the bond angle + angle + The bond angle (radians) + sigma + The uncertainty of the bond angle (radians) + delta + The width of the bond angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondAngleRestraint object + Returns the ObjCrystBondAngleRestraint object for use with the + 'unrestrain' method. + """ + res = ObjCrystBondAngleRestraint( + atom1, atom2, atom3, angle, sigma, delta, scaled + ) + self._restraints.add(res) + + return res + + def restrain_bond_angle_parameter( + self, par, angle, sigma, delta, scaled=False + ): + """Add a bond angle restraint. + + This creates an instance of ObjCrystBondAngleRestraint and adds it to + the ObjCrystMoleculeParSet. + + Parameters + ---------- + par + A ObjCrystBondAngleParameter (see add_bond_angle_parameter) + angle + The bond angle (radians) + sigma + The uncertainty of the bond angle (radians) + delta + The width of the bond angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondAngleRestraint object + Returns the ObjCrystBondAngleRestraint object for use with the + 'unrestrain' method. + """ + return self.restrain_bond_angle( + par.atom1, par.atom2, par.atom3, angle, sigma, delta, scaled + ) + + def restrain_dihedral_angle( + self, atom1, atom2, atom3, atom4, angle, sigma, delta, scaled=False + ): + """Add a dihedral angle restraint. + + This creates an instance of ObjCrystDihedralAngleRestraint and adds it + to the ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the angle + atom2 + Second (central) atom (ObjCrystMolAtomParSet) in the angle + atom3 + Third (central) atom (ObjCrystMolAtomParSet) in the angle + atom4 + Fourth atom in the angle (ObjCrystMolAtomParSet) + angle + The dihedral angle (radians) + sigma + The uncertainty of the dihedral angle (radians) + delta + The width of the dihedral angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystDihedralAngleRestraint object + Returns the ObjCrystDihedralAngleRestraint object for use with the + 'unrestrain' method. + """ + res = ObjCrystDihedralAngleRestraint( + atom1, atom2, atom3, atom4, angle, sigma, delta, scaled + ) + self._restraints.add(res) + + return res + + def restrain_dihedral_angle_parameter( + self, par, angle, sigma, delta, scaled=False + ): + """Add a dihedral angle restraint. + + This creates an instance of ObjCrystDihedralAngleRestraint and adds it + to the ObjCrystMoleculeParSet. + + Parameters + ---------- + par + A ObjCrystDihedralAngleParameter (see + add_dihedral_angle_parameter) + angle + The dihedral angle (radians) + sigma + The uncertainty of the dihedral angle (radians) + delta + The width of the dihedral angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystDihedralAngleRestraint object + Returns the ObjCrystDihedralAngleRestraint object for use with the + 'unrestrain' method. + """ + return self.restrain_dihedral_angle( + par.atom1, + par.atom2, + par.atom3, + par.atom4, + angle, + sigma, + delta, + scaled, + ) + + def add_bond_length_parameter( + self, name, atom1, atom2, value=None, const=False + ): + """Add a bond length to the Molecule. + + This creates a ObjCrystBondLengthParameter to the + ObjCrystMoleculeParSet that can be adjusted during the fit. + + Parameters + ---------- + name + The name of the ObjCrystBondLengthParameter + atom1 + The first atom (ObjCrystMolAtomParSet) in the bond + atom2 + The second (mutated) atom (ObjCrystMolAtomParSet) in the + bond + value + An initial value for the bond length. If this is None + (default), then the current distance between the atoms will + be used. + const + A flag indicating whether the Parameter is constant + (default False) + + Returns + ------- + ObjCrystBondLengthParameter object + Returns the new ObjCrystBondLengthParameter. + """ + par = ObjCrystBondLengthParameter(name, atom1, atom2, value, const) + self.addParameter(par) + + return par + + def add_bond_angle_parameter( + self, name, atom1, atom2, atom3, value=None, const=False + ): + """Add a bond angle to the Molecule. + + This creates a ObjCrystBondAngleParameter to the ObjCrystMoleculeParSet + that can be adjusted during the fit. + + Parameters + ---------- + name + The name of the ObjCrystBondAngleParameter + atom1 + The first atom (ObjCrystMolAtomParSet) in the bond angle + atom2 + The second (central) atom (ObjCrystMolAtomParSet) in the + bond angle + atom3 + The third (mutated) atom (ObjCrystMolAtomParSet) in the + bond angle + value + An initial value for the bond angle. If this is None + (default), then the current bond angle between the atoms + will be used. + const + A flag indicating whether the Parameter is constant + (default False). + + Returns + ------- + ObjCrystBondAngleParameter object + Returns the new ObjCrystBondAngleParameter. + """ + par = ObjCrystBondAngleParameter( + name, atom1, atom2, atom3, value, const + ) + self.addParameter(par) + + return par + + def add_dihedral_angle_parameter( + self, name, atom1, atom2, atom3, atom4, value=None, const=False + ): + """Add a dihedral angle to the Molecule. + + This creates a ObjCrystDihedralAngleParameter to the + ObjCrystMoleculeParSet that can be adjusted during the fit. + + Parameters + ---------- + name + The name of the ObjCrystDihedralAngleParameter. + atom1 + The first atom (ObjCrystMolAtomParSet) in the dihderal + angle. + atom2 + The second (central) atom (ObjCrystMolAtomParSet) in the + dihderal angle + atom3 + The third (central) atom (ObjCrystMolAtomParSet) in the + dihderal angle + atom4 + The fourth (mutated) atom (ObjCrystMolAtomParSet) in the + dihderal angle + value + An initial value for the dihedral angle. If this is None + (default), then the current dihedral angle between atoms + will be used. + const + A flag indicating whether the Parameter is constant + (default False). + + Returns + ------- + ObjCrystDihedralAngleParameter object + Returns the new ObjCrystDihedralAngleParameter. + """ + par = ObjCrystDihedralAngleParameter( + name, atom1, atom2, atom3, atom4, value, const + ) + self.addParameter(par) + + return par + + +# End class ObjCrystMoleculeParSet + + +class ObjCrystMolAtomParSet(ObjCrystScattererParSet): + """A adaptor for an pyobjcryst.molecule.MolAtom. + + This class derives from srfit.fitbase.parameterset.ParameterSet. Note that + MolAtom does not derive from Scatterer, but the relevant interface is the + same within pyobjcryst. See the ParameterSet class for base attributes. + + Attributes + ---------- + scat + The adapted pyobjcryst.molecule.MolAtom. + parent + The ObjCrystCrystalParSet this belongs to + element + Non-refinable name of the element (property). + + Managed Parameters + ------------------ + occ + Occupancy of the atom on its crystal location + (ParameterAdapter) + Biso + Isotropic scattering factor (ParameterAdapter). This does + not exist for dummy atoms. See the 'is_dummy' method. + B11, B22, B33, B12, B21, B23, B32, B13, B31 + -- Anisotropic displacement factor for scatterer + (ParameterAdapter or ParameterProxy). Note that the Bij and Bji + parameters are the same. + """ + + def __init__(self, name, scat, parent): + """Initialize. + + Parameters + ---------- + name + The name of the scatterer + scat + The Scatterer instance + parent + The ObjCrystCrystalParSet this belongs to + """ + ObjCrystScattererParSet.__init__(self, name, scat, parent) + sp = scat.GetScatteringPower() + + # Only wrap this if there is a scattering power + if sp is not None: + self.addParameter(ParameterAdapter("Biso", sp, attr="Biso")) + self.addParameter(ParameterAdapter("B11", sp, attr="B11")) + self.addParameter(ParameterAdapter("B22", sp, attr="B22")) + self.addParameter(ParameterAdapter("B33", sp, attr="B33")) + B12 = ParameterAdapter("B12", sp, attr="B12") + B21 = ParameterProxy("B21", B12) + B13 = ParameterAdapter("B13", sp, attr="B13") + B31 = ParameterProxy("B31", B13) + B23 = ParameterAdapter("B23", sp, attr="B23") + B32 = ParameterProxy("B32", B23) + self.addParameter(B12) + self.addParameter(B21) + self.addParameter(B13) + self.addParameter(B31) + self.addParameter(B23) + self.addParameter(B32) + + return + + def _getelem(self): + """Getter for the element type.""" + sp = self.scat.GetScatteringPower() + if sp: + return sp.GetSymbol() + else: + return "dummy" + + element = property(_getelem) + + def is_dummy(self): + """Indicate whether this atom is a dummy atom.""" + return self.scat.IsDummy() + + +# End class ObjCrystMolAtomParSet + + +class ObjCrystMoleculeRestraint(object): + """Base class for adapting pyobjcryst Molecule restraints to srfit. + + The 'penalty' method calls 'GetLogLikelihood' of the pyobjcryst restraint. + This implements the 'penalty' method from + diffpy.srfit.fitbase.restraint.Restraint. The 'restrain' method is not + needed or implemented. + + Attributes + ---------- + res + The pyobjcryst Molecule restraint. + scaled + A flag indicating if the restraint is scaled (multiplied) by + the unrestrained point-average chi^2 (chi^2/numpoints) (default + False). + """ + + def __init__(self, res, scaled=False): + """Create a Restraint-like from a pyobjcryst Molecule restraint. + + Parameters + ---------- + res + The pyobjcryst Molecule restraint. + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + self.res = res + self.scaled = scaled + return + + def penalty(self, w=1.0): + """Calculate the penalty of the restraint. + + Parameters + ---------- + w + The point-average chi^2 which is optionally used to scale the + penalty (default 1.0). + """ + penalty = self.res.GetLogLikelihood() + if self.scaled: + penalty *= w + return penalty + + +# End class ObjCrystMoleculeRestraint + + +class ObjCrystBondLengthRestraint(ObjCrystMoleculeRestraint): + """Restrain the distance between two atoms. + + Attributes + ---------- + atom1 + The first atom in the bond (ObjCrystMolAtomParSet) + atom2 + The second atom in the bond (ObjCrystMolAtomParSet) + length + The length of the bond (Angstroms) + sigma + The uncertainty of the bond length (Angstroms) + delta + The width of the bond (Angstroms) + res + The pyobjcryst BondLength restraint + scaled + A flag indicating if the restraint is scaled (multiplied) by + the unrestrained point-average chi^2 (chi^2/numpoints) (default + False) + """ + + def __init__(self, atom1, atom2, length, sigma, delta, scaled=False): + """Create a bond length restraint. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the bond + atom2 + Second atom (ObjCrystMolAtomParSet) in the bond + length + The length of the bond (Angstroms) + sigma + The uncertainty of the bond length (Angstroms) + delta + The width of the bond (Angstroms) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False) + """ + self.atom1 = atom1 + self.atom2 = atom2 + + m = self.atom1.scat.GetMolecule() + res = m.AddBond(atom1.scat, atom2.scat, length, sigma, delta) + + ObjCrystMoleculeRestraint.__init__(self, res, scaled) + return + + # Give access to the parameters of the restraint + length = property( + lambda self: self.res.GetLength0(), + lambda self, val: self.res.SetLength0(val), + ) + sigma = property( + lambda self: self.res.GetLengthSigma(), + lambda self, val: self.res.SetLengthSigma(val), + ) + delta = property( + lambda self: self.res.GetLengthDelta(), + lambda self, val: self.res.SetLengthDelta(val), + ) + + +# End class ObjCrystBondLengthRestraint + + +class ObjCrystBondAngleRestraint(ObjCrystMoleculeRestraint): + """Restrain the angle defined by three atoms. + + Attributes + ---------- + atom1 + The first atom in the angle (ObjCrystMolAtomParSet) + atom2 + The second atom in the angle (ObjCrystMolAtomParSet) + atom3 + The third atom in the angle (ObjCrystMolAtomParSet) + angle + The bond angle (radians) + sigma + The uncertainty of the bond angle (radians) + delta + The width of the bond angle (radians) + res + The pyobjcryst BondAngle restraint + scaled + A flag indicating if the restraint is scaled (multiplied) by + the unrestrained point-average chi^2 (chi^2/numpoints) (default + False) + """ + + def __init__(self, atom1, atom2, atom3, angle, sigma, delta, scaled=False): + """Create a bond angle restraint. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the bond angle + atom2 + Second (central) atom (ObjCrystMolAtomParSet) in the bond + angle + atom3 + Third atom (ObjCrystMolAtomParSet) in the bond angle + angle + The bond angle (radians) + sigma + The uncertainty of the bond angle (radians) + delta + The width of the bond angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + self.atom1 = atom1 + self.atom2 = atom2 + self.atom3 = atom3 + + m = self.atom1.scat.GetMolecule() + res = m.AddBondAngle( + atom1.scat, atom2.scat, atom3.scat, angle, sigma, delta + ) + + ObjCrystMoleculeRestraint.__init__(self, res, scaled) + return + + # Give access to the parameters of the restraint + angle = property( + lambda self: self.res.GetAngle0(), + lambda self, val: self.res.SetAngle0(val), + ) + sigma = property( + lambda self: self.res.GetAngleSigma(), + lambda self, val: self.res.SetAngleSigma(val), + ) + delta = property( + lambda self: self.res.GetAngleDelta(), + lambda self, val: self.res.SetAngleDelta(val), + ) + + +# End class ObjCrystBondAngleRestraint + + +class ObjCrystDihedralAngleRestraint(ObjCrystMoleculeRestraint): + """Restrain the dihedral (torsion) angle defined by four atoms. + + Attributes + ---------- + atom1 + The first atom in the angle (ObjCrystMolAtomParSet) + atom2 + The second (central) atom in the angle (ObjCrystMolAtomParSet) + atom3 + The third (central) atom in the angle (ObjCrystMolAtomParSet) + atom4 + The fourth atom in the angle (ObjCrystMolAtomParSet) + angle + The dihedral angle (radians) + sigma + The uncertainty of the dihedral angle (radians) + delta + The width of the dihedral angle (radians) + res + The pyobjcryst DihedralAngle restraint + scaled + A flag indicating if the restraint is scaled (multiplied) by + the unrestrained point-average chi^2 (chi^2/numpoints) (default + False) + """ + + def __init__( + self, atom1, atom2, atom3, atom4, angle, sigma, delta, scaled=False + ): + """Create a dihedral angle restraint. + + Parameters + ---------- + atom1 + First atom (ObjCrystMolAtomParSet) in the angle + atom2 + Second (central) atom (ObjCrystMolAtomParSet) in the angle + atom3 + Third (central) atom (ObjCrystMolAtomParSet) in the angle + atom4 + Fourth atom in the angle (ObjCrystMolAtomParSet) + angle + The dihedral angle (radians) + sigma + The uncertainty of the dihedral angle (radians) + delta + The width of the dihedral angle (radians) + scaled + A flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + self.atom1 = atom1 + self.atom2 = atom2 + self.atom3 = atom3 + self.atom4 = atom4 + + m = self.atom1.scat.GetMolecule() + res = m.AddDihedralAngle( + atom1.scat, atom2.scat, atom3.scat, atom4.scat, angle, sigma, delta + ) + + ObjCrystMoleculeRestraint.__init__(self, res, scaled) + return + + # Give access to the parameters of the restraint + angle = property( + lambda self: self.res.GetAngle0(), + lambda self, val: self.res.SetAngle0(val), + ) + sigma = property( + lambda self: self.res.GetAngleSigma(), + lambda self, val: self.res.SetAngleSigma(val), + ) + delta = property( + lambda self: self.res.GetAngleDelta(), + lambda self, val: self.res.SetAngleDelta(val), + ) + + +# End class ObjCrystDihedralAngleRestraint + + +class StretchModeParameter(Parameter): + """Partial Parameter class encapsulating pyobjcryst stretch modes. + + This class relies upon attributes that do not belong to it. Do not + instantiate this class. + + Required attributes + ------------------- + matoms + The set of all mutated AtomParSets + molecule + The ObjCrystMoleculeParSet the atoms belong to + mode + The pyobjcryst.molecule.StretchMode used to change atomic + positions. + keepcenter + Flag indicating whether to keep the center of mass of the + molecule stationary within the crystal when changing the + value of the parameter (bool, default True). + """ + + def __init__(self, name, value=None, const=False): + """Initialization. + + Parameters + ---------- + name + The name of this Parameter (must be a valid attribute + identifier) + value + The initial value of this Parameter (default 0). + const + A flag inticating whether the Parameter is a constant (like + pi). + + Raises ValueError if the name is not a valid attribute identifier + """ + Parameter.__init__(self, name, value, const) + self.keepcenter = True + + def set_value(self, val): + """Change the value of the Parameter.""" + curval = self.getValue() + val = float(val) + + if val == curval: + return self + + # The StretchMode expects the change in mutated value. + delta = val - curval + self.mode.Stretch(delta, self.keepcenter) + + # Let Parameter take care of the general details + Parameter.set_value(self, val) + + return self + + def add_atoms(self, atomlist): + """Associate ObjCrystMolAtomParSets with the Parameter. + + This will associate additional ObjCrystMolAtomParSets with the + Parameter. These will be mutated in the exact same way as the + primary mutated ObjCrystMolAtomParSet. This is useful when a + group of atoms should move rigidly in response to a change in a + bond property. + """ + if not hasattr(atomlist, "__iter__"): + atomlist = [atomlist] + # Record the added atoms in the Parameter + self.matoms.update(atomlist) + # Make sure we're observing these atoms + for a in atomlist: + a.x.addObserver(self._flush) + a.y.addObserver(self._flush) + a.z.addObserver(self._flush) + + # Record the added atoms in the StretchMode + scatlist = [a.scat for a in atomlist] + self.mode.AddAtoms(scatlist) + return self + + def notify(self, other=()): + """Notify all mutated Parameters and observers. + + Some of the mutated parameters will be observing us. At the same + time we need to observe them. Observable won't let us do both, + so we notify the Parameters that we mutate directly. + """ + noneother = () + # Notify the atoms that have moved + for a in self.matoms: + a.x._flush(noneother) + a.y._flush(noneother) + a.z._flush(noneother) + # Notify the molecule position + self.molecule.x._flush(noneother) + self.molecule.y._flush(noneother) + self.molecule.z._flush(noneother) + + # Notify observers + Parameter.notify(self, other) + return + + +# End class StretchModeParameter + + +class ObjCrystBondLengthParameter(StretchModeParameter): + """Class for abstracting a bond length in a Molecule to a Parameter. + + This wraps up a pyobjcryst.molecule.StretchModeBondLength object so that + the distance between two MolAtoms in a Molecule can be used as an + adjustable Parameter. When a bond length is adjusted, the second MolAtom is + moved, and the absolute position of the Molecule is altered to preserve the + location of the center of mass within the Crystal. Thus, the x, y and z + Parameters of the MolAtom and its parent Molecule are altered. This can be + changed by setting the 'keepcenter' attribute of the parameter to False. + + This Parameter makes it possible to mutate a MolAtom multiple times in a + single refinement step. If these mutations are not orthogonal, then this + could lead to nonconvergence of a fit, depending on the optimizer. Consider + mutating atom2 of a bond directly, and via a ObjCrystBondLengthParameter. + The two mutations of atom2 may be determined independently by the + optimizer, in which case the composed mutation will have an unexpected + effect on the residual. It is best practice to either modify MolAtom + positions directly, or thorough BondLengthParameters, BondAngleParameters + and DihedralAngleParameters (which are mutually orthogonal). + + Note that by making a ObjCrystBondLengthParameter constant it also makes + the underlying ObjCrystMolAtomParSets constant. When setting it as + nonconstant, each ObjCrystMolAtomParSet is set nonconstant. Changing the + bond length changes the position of the second MolAtom and Molecule, even + if either is set as constant. + + Attributes + ---------- + atom1 + The first ObjCrystMolAtomParSet in the bond + atom2 + The second (mutated) ObjCrystMolAtomParSet in the bond + matoms + The set of all mutated ObjCrystMolAtomParSets + molecule + The ObjCrystMoleculeParSet the ObjCrystMolAtomParSets + belong to + mode + The pyobjcryst.molecule.StretchModeBondLength for the bond + + Inherited Attributes + -------------------- + name + A name for this Parameter. + const + A flag indicating whether this is considered a constant. + _value + The value of the Parameter. Modified with 'set_value'. + value + Property for 'getValue' and 'set_value'. + constraint + A callable that calculates the value of this Parameter. If + this is None (None), the the Parameter is responsible for its + own value. The callable takes no arguments. + bounds + A 2-list defining the bounds on the Parameter. This can be + used by some optimizers when the Parameter is varied. + """ + + def __init__(self, name, atom1, atom2, value=None, const=False, mode=None): + """Create a ObjCrystBondLengthParameter. + + Parameters + ---------- + name + The name of the ObjCrystBondLengthParameter + atom1 + The first atom (ObjCrystMolAtomParSet) in the bond + atom2 + The second (mutated) atom (ObjCrystMolAtomParSet) in the + bond + value + An initial value for the bond length. If this is None + (default), then the current distance between the atoms will + be used. + const + A flag indicating whether the Parameter is constant + (default False) + mode + An extant pyobjcryst.molecule.StretchModeBondLength to use. + If this is None (default), then a new StretchModeBondLength + will be built. + """ + # Create the mode + self.mode = mode + if mode is None: + self.mode = StretchModeBondLength(atom1.scat, atom2.scat, None) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom2.scat) + self.matoms = set([atom2]) + + # Observe the atom positions + for a in [atom1, atom2]: + a.x.addObserver(self._flush) + a.y.addObserver(self._flush) + a.z.addObserver(self._flush) + + self.atom1 = atom1 + self.atom2 = atom2 + self.molecule = atom1.parent + + # We do this last so the atoms are defined before we set any values. + if value is None: + value = GetBondLength(atom1.scat, atom2.scat) + StretchModeParameter.__init__(self, name, value, const) + self.set_constant(const) + + return + + def set_constant(self, is_constant=True, value=None): + """Toggle the Parameter as constant. + + This sets the underlying ObjCrystMolAtomParSet positions + constant as well. + + Parameters + ---------- + is_constant + Flag indicating if the Parameter is constant (default + True). + value + An optional value for the Parameter (default None). If this + is not None, then the Parameter will get a new value, + constant or otherwise. + + Return + ------ + self + Returns self so that mutators can be chained. + """ + StretchModeParameter.set_constant(self, is_constant, value) + + for a in [self.atom1, self.atom2]: + a.x.set_constant(is_constant) + a.y.set_constant(is_constant) + a.z.set_constant(is_constant) + return self + + def getValue(self): + """This calculates the value if it might have been changed. + + There is no guarantee that the ObjCrystMolAtomParSets underlying + the bond won't change, so the bond length is calculated if + necessary each time this is called. + """ + if self._value is None: + val = GetBondLength(self.atom1.scat, self.atom2.scat) + Parameter.set_value(self, val) + + return self._value + + +# End class ObjCrystBondLengthParameter + + +class ObjCrystBondAngleParameter(StretchModeParameter): + """Class for abstracting a bond angle in a Molecule to a Parameter. + + This wraps up a pyobjcryst.molecule.StretchModeBondAngle object so that the + angle defined by three MolAtoms in a Molecule can be used as an adjustable + Parameter. When a bond angle is adjusted, the third MolAtom is moved, and + the absolute position of the Molecule is altered to preserve the location + of the center of mass within the crystal. This can be changed by setting + the 'keepcenter' attribute of the parameter to False. + + See precautions in the ObjCrystBondLengthParameter class. + + Attributes + ---------- + atom1 + The first ObjCrystAtomParSet in the bond angle + atom2 + The second (central) ObjCrystMolAtomParSet in the bond angle + atom3 + The third (mutated) ObjCrystMolAtomParSet in the bond angle + matoms + The set of all mutated ObjCrystMolAtomParSets + molecule + The ObjCrystMoleculeParSet the ObjCrystMolAtomParSets + belong to + mode + The pyobjcryst.molecule.StretchModeBondAngle for the bond angle + + Inherited Attributes + -------------------- + name + A name for this Parameter. + const + A flag indicating whether this is considered a constant. + _value + The value of the Parameter. Modified with 'set_value'. + value + Property for 'getValue' and 'set_value'. + constraint + A callable that calculates the value of this Parameter. If + this is None (None), the the Parameter is responsible for its + own value. The callable takes no arguments. + bounds + A 2-list defining the bounds on the Parameter. This can be + used by some optimizers when the Parameter is varied. + """ + + def __init__( + self, name, atom1, atom2, atom3, value=None, const=False, mode=None + ): + """Create a ObjCrystBondAngleParameter. + + Parameters + ---------- + name + The name of the ObjCrystBondAngleParameter. + atom1 + The first atom (ObjCrystMolAtomParSet) in the bond angle + atom2 + The second (central) atom (ObjCrystMolAtomParSet) in the + bond angle + atom3 + The third (mutated) atom (ObjCrystMolAtomParSet) in the + bond angle + value + An initial value for the bond length. If this is None + (default), then the current bond angle between the atoms + will be used. + const + A flag indicating whether the Parameter is constant + (default False). + mode + A pre-built mode to place in this Parameter. If this is + None (default), then a StretchMode will be built. + """ + # Create the stretch mode + self.mode = mode + if mode is None: + self.mode = StretchModeBondAngle( + atom1.scat, atom2.scat, atom3.scat, None + ) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom3.scat) + self.matoms = set([atom3]) + + # Observe the atom positions + for a in [atom1, atom2, atom3]: + a.x.addObserver(self._flush) + a.y.addObserver(self._flush) + a.z.addObserver(self._flush) + + self.atom1 = atom1 + self.atom2 = atom2 + self.atom3 = atom3 + self.molecule = atom1.parent + + # We do this last so the atoms are defined before we set any values. + if value is None: + value = GetBondAngle(atom1.scat, atom2.scat, atom3.scat) + StretchModeParameter.__init__(self, name, value, const) + self.set_constant(const) + + return + + def set_constant(self, is_constant=True, value=None): + """Toggle the Parameter as constant. + + This sets the underlying ObjCrystMolAtomParSet positions + constant as well. + + Parameters + ---------- + is_constant + Flag indicating if the Parameter is constant (default + True). + value + An optional value for the Parameter (default None). If this + is not None, then the Parameter will get a new value, + constant or otherwise. + + Return + ------ + self + Returns self so that mutators can be chained. + """ + StretchModeParameter.set_constant(self, is_constant, value) + for a in [self.atom1, self.atom2, self.atom3]: + a.x.set_constant(is_constant) + a.y.set_constant(is_constant) + a.z.set_constant(is_constant) + return self + + def getValue(self): + """This calculates the value if it might have been changed. + + There is no guarantee that the MolAtoms underlying the bond + angle won't change, so the bond angle is calculated if necessary + each time this is called. + """ + if self._value is None: + val = GetBondAngle( + self.atom1.scat, self.atom2.scat, self.atom3.scat + ) + Parameter.set_value(self, val) + + return self._value + + +# End class ObjCrystBondAngleParameter + + +class ObjCrystDihedralAngleParameter(StretchModeParameter): + """Class for abstracting a dihedral angle in a Molecule to a + Parameter. + + This wraps up a pyobjcryst.molecule.StretchModeTorsion object so that the + angle defined by four MolAtoms ([a1-a2].[a3-a4]) in a Molecule can be used + as an adjustable parameter. When a dihedral angle is adjusted, the fourth + MolAtom is moved, and the absolute position of the Molecule is altered to + preserve the location of the center of mass within the crystal. This can + be changed by setting the 'keepcenter' attribute of the parameter to False. + + See precautions in the ObjCrystBondLengthParameter class. + + Attributes + ---------- + atom1 + The first ObjCrystMolAtomParSet in the dihedral angle + atom2 + The second (central) ObjCrystMolAtomParSet in the dihedral + angle + atom3 + The third (central) ObjCrystMolAtomParSet in the dihedral angle + atom4 + The fourth (mutated) ObjCrystMolAtomParSet in the dihedral + angle + matoms + The set of all mutated ObjCrystMolAtomParSets + molecule + The ObjCrystMoleculeParSet the atoms belong to + mode + The pyobjcryst.molecule.StretchModeTorsion for the dihedral + angle + + Inherited Attributes + -------------------- + name + A name for this Parameter. + const + A flag indicating whether this is considered a constant. + _value + The value of the Parameter. Modified with 'set_value'. + value + Property for 'getValue' and 'set_value'. + constraint + A callable that calculates the value of this Parameter. If + this is None (None), the the Parameter is responsible for its + own value. The callable takes no arguments. + bounds + A 2-list defining the bounds on the Parameter. This can be + used by some optimizers when the Parameter is varied. + """ + + def __init__( + self, + name, + atom1, + atom2, + atom3, + atom4, + value=None, + const=False, + mode=None, + ): + """Create a ObjCrystDihedralAngleParameter. + + Parameters + ---------- + name + The name of the ObjCrystDihedralAngleParameter + atom1 + The first atom (ObjCrystMolAtomParSet) in the dihderal + angle + atom2 + The second (central) atom (ObjCrystMolAtomParSet) in the + dihderal angle + atom3 + The third (central) atom (ObjCrystMolAtomParSet) in the + dihderal angle + atom4 + The fourth (mutated) atom (ObjCrystMolAtomParSet) in the + dihderal angle + value + An initial value for the bond length. If this is None + (default), then the current dihedral angle between atoms + will be used. + const + A flag indicating whether the Parameter is constant + (default False). + mode + A pre-built mode to place in this Parameter. If this is + None (default), then a StretchMode will be built. + """ + # Create the stretch mode + self.mode = mode + if mode is None: + self.mode = StretchModeTorsion(atom2.scat, atom3.scat, None) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom4.scat) + self.matoms = set([atom4]) + + # Observe the atom positions + for a in [atom1, atom2, atom3, atom4]: + a.x.addObserver(self._flush) + a.y.addObserver(self._flush) + a.z.addObserver(self._flush) + + self.atom1 = atom1 + self.atom2 = atom2 + self.atom3 = atom3 + self.atom4 = atom4 + self.molecule = atom1.parent + + # We do this last so the atoms are defined before we set any values. + if value is None: + value = GetDihedralAngle( + atom1.scat, atom2.scat, atom3.scat, atom4.scat + ) + StretchModeParameter.__init__(self, name, value, const) + self.set_constant(const) + + return + + def set_constant(self, is_constant=True, value=None): + """Toggle the Parameter as constant. + + This sets the underlying ObjCrystMolAtomParSet positions const as well. + + Parameters + ---------- + is_constant + Flag indicating if the Parameter is constant (default + True). + value + An optional value for the Parameter (default None). If this + is not None, then the Parameter will get a new value, + constant or otherwise. + + Return + ------ + self + Returns self so that mutators can be chained. + """ + StretchModeParameter.set_constant(self, is_constant, value) + for a in [self.atom1, self.atom2, self.atom3, self.atom4]: + a.x.set_constant(is_constant) + a.y.set_constant(is_constant) + a.z.set_constant(is_constant) + return self + + def getValue(self): + """This calculates the value if it might have been changed. + + There is no guarantee that the ObjCrystMolAtomParSets underlying + the dihedral angle won't change from some other Parameter, so + the value is recalculated each time. + """ + if self._value is None: + val = GetDihedralAngle( + self.atom1.scat, + self.atom2.scat, + self.atom3.scat, + self.atom4.scat, + ) + Parameter.set_value(self, val) + + return self._value + + +# End class ObjCrystDihedralAngleParameter + + +class ObjCrystCrystalParSet(SrRealParSet): + """A adaptor for pyobjcryst.crystal.Crystal instance. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. + See this class for base attributes. + + Attributes + ---------- + stru + The adapted pyobjcryst.Crystal. + scatterers + The list of aggregated ScattererParSets (either + ObjCrystAtomParSet or ObjCrystMoleculeParSet), provided for + convenience. + _sgpars + A BaseSpaceGroupParameters object containing free structure + Parameters. See the diffpy.cmipdf.structure.sgconstraints + module. + sgpars + property that creates _sgpars when it is needed. + angunits + "rad", the units of angle + + Parameters + ---------- + x + Scatterer position in crystal coordinates (ParameterWraper) + y + Scatterer position in crystal coordinates (ParameterWraper) + z + Scatterer position in crystal coordinates (ParameterWraper) + occ + Occupancy of the scatterer on its crystal site (ParameterWraper) + """ + + def __init__(self, name, cryst): + """Initialize. + + Parameters + ---------- + name + A name for this ParameterSet + cryst + An pyobjcryst.Crystal instance. + """ + SrRealParSet.__init__(self, name) + self.angunits = "rad" + self.stru = cryst + self._sgpars = None + + self.addParameter(ParameterAdapter("a", self.stru, attr="a")) + self.addParameter(ParameterAdapter("b", self.stru, attr="b")) + self.addParameter(ParameterAdapter("c", self.stru, attr="c")) + self.addParameter(ParameterAdapter("alpha", self.stru, attr="alpha")) + self.addParameter(ParameterAdapter("beta", self.stru, attr="beta")) + self.addParameter(ParameterAdapter("gamma", self.stru, attr="gamma")) + + # Now we must loop over the scatterers and create parameter sets from + # them. + self.scatterers = [] + snames = [] + + for j in range(self.stru.GetNbScatterer()): + s = self.stru.GetScatt(j) + name = s.GetName() + if not name: + raise ValueError("Each Scatterer must have a name") + if name in snames: + raise ValueError("Scatterer name '%s' is duplicated" % name) + + # Now create the proper object + cname = s.GetClassName() + if cname == "Atom": + parset = ObjCrystAtomParSet(name, s, self) + elif cname == "Molecule": + parset = ObjCrystMoleculeParSet(name, s, self) + else: + raise TypeError("Unrecognized scatterer '%s'" % cname) + + self.add_parameter_set(parset) + self.scatterers.append(parset) + snames.append(name) + + return + + def _constrain_space_group(self): + """Constrain the space group.""" + if self._sgpars is not None: + return self._sgpars + sg = self._create_space_group(self.stru.GetSpaceGroup()) + from diffpy.cmipdf.structure.sgconstraints import ( + _constrain_as_space_group, + ) + + adpsymbols = ["B11", "B22", "B33", "B12", "B13", "B23"] + isosymbol = "Biso" + sgoffset = [0, 0, 0] + self._sgpars = _constrain_as_space_group( + self, + sg, + self.scatterers, + sgoffset, + adpsymbols=adpsymbols, + isosymbol=isosymbol, + ) + return self._sgpars + + sgpars = property(_constrain_space_group) + + @staticmethod + def _create_space_group(sgobjcryst): + """Create a diffpy.structure SpaceGroup object from pyobjcryst. + + Parameters + ---------- + sgobjcryst + A pyobjcryst.spacegroup.SpaceGroup instance. + + This uses the actual space group operations from the + pyobjcryst.spacegroup.SpaceGroup instance so there is no ambiguity + about the actual space group. + """ + import copy + + from diffpy.structure.spacegroups import GetSpaceGroup, SymOp + + name = sgobjcryst.GetName() + extnstr = ":%s" % sgobjcryst.GetExtension() + if name.endswith(extnstr): + name = name[: -len(extnstr)] + + # Get whatever spacegroup we can get by name. This will set the proper + # crystal system. Creating a copy of the singleton from GetSpaceGroup, + # as this function messes with sg.symop_list. + sg = copy.copy(GetSpaceGroup(name)) + + # Replace the symmetry operations to guarantee that we get it right. + symops = sgobjcryst.GetSymmetryOperations() + tranops = sgobjcryst.GetTranslationVectors() + sg.symop_list = [] + + for trans in tranops: + for shift, rot in symops: + tv = trans + shift + tv -= numpy.floor(tv) + sg.symop_list.append(SymOp(rot, tv)) + + if sgobjcryst.IsCentrosymmetric(): + center = sgobjcryst.GetInversionCenter() + for trans in tranops: + for shift, rot in symops: + tv = center - trans - shift + tv -= numpy.floor(tv) + sg.symop_list.append(SymOp(-rot, tv)) + + return sg + + @classmethod + def can_adapt(self, stru): + """Return whether the structure can be adapted by this class.""" + from pyobjcryst.crystal import Crystal + + return isinstance(stru, Crystal) + + def get_lattice(self): + """Get the ParameterSet containing the lattice Parameters.""" + return self + + def get_scatterers(self): + """Get a list of ParameterSets that represents the scatterers. + + The site positions must be accessible from the list entries via + the names "x", "y", and "z". The ADPs must be accessible as + well, but the name and nature of the ADPs (U-factors, B-factors, + isotropic, anisotropic) depends on the adapted structure. + """ + return self.scatterers + + +# End class ObjCrystCrystalParSet diff --git a/src/diffpy/cmipdf/structure/sgconstraints.py b/src/diffpy/cmipdf/structure/sgconstraints.py new file mode 100644 index 0000000..5e0647c --- /dev/null +++ b/src/diffpy/cmipdf/structure/sgconstraints.py @@ -0,0 +1,827 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Code to set space group constraints for a crystal structure.""" + +import re + +import numpy + +from diffpy.srfit.fitbase.parameter import ParameterProxy +from diffpy.srfit.fitbase.recipeorganizer import RecipeContainer + +__all__ = ["constrain_as_space_group"] + + +def constrain_as_space_group( + phase, + spacegroup, + scatterers=None, + sgoffset=[0, 0, 0], + constrainlat=True, + constrainadps=True, + adpsymbols=None, + isosymbol="Uiso", +): + """Constrain the structure to the space group. + + This applies space group constraints to a StructureParSet with P1 + symmetry. Passed scatterers are explicitly constrained to the + specified space group. The ADPs and lattice may be constrained as well. + + Parameters + ---------- + phase + A BaseStructure object. + spacegroup + The space group number, symbol or an instance of + SpaceGroup class from diffpy.structure package. + sgoffset + Optional offset for sg origin (default [0, 0, 0]). + scatterers + The scatterer ParameterSets to constrain. If scatterers + is None (default), then all scatterers accessible from + phase.get_scatterers will be constrained. + constrainlat + Flag indicating whether to constrain the lattice + (default True). + constrainadps + Flag indicating whether to constrain the ADPs + (default True). + adpsymbols + A list of the ADP names. By default this is equal to + diffpy.structure.symmetryutilities.stdUsymbols (U11, + U22, etc.). The names must be given in the same order + as stdUsymbols. + isosymbol + Symbol for isotropic ADP (default "Uiso"). If None, + isotropic ADPs will be constrained via the anisotropic ADPs. + + + New Parameters that are used in constraints are created within a + SpaceGroupParameters object, which is returned from this function. + Constraints are created in ParameterSet that contains the constrained + Parameter. This will erase any constraints or constant flags on the + scatterers, lattice or ADPs if they are to be constrained. + + The lattice constraints are applied as following. + + Crystal System + Triclinic + No constraints. + Monoclinic + alpha and beta are fixed to 90 unless alpha != beta and + alpha == gamma, in which case alpha and gamma are fixed + to 90. + Orthorhombic + alpha, beta and gamma are fixed to 90. + Tetragonal + b is constrained to a and alpha, beta and gamma are + fixed to 90. + Trigonal + If gamma == 120, then b is constrained to a, alpha + and beta are fixed to 90 and gamma is fixed to 120. + Otherwise, b and c are constrained to a, beta and gamma + are fixed to alpha. + Hexagonal + b is constrained to a, alpha and beta are fixed to 90 + and gamma is fixed to 120. + Cubic + b and c are constrained to a, and alpha, beta and + gamma are fixed to 90. + """ + from diffpy.structure.spacegroups import GetSpaceGroup, SpaceGroup + + sg = spacegroup + if not isinstance(spacegroup, SpaceGroup): + sg = GetSpaceGroup(spacegroup) + sgp = _constrain_as_space_group( + phase, + sg, + scatterers, + sgoffset, + constrainlat, + constrainadps, + adpsymbols, + isosymbol, + ) + + return sgp + + +def _constrain_as_space_group( + phase, + sg, + scatterers=None, + sgoffset=[0, 0, 0], + constrainlat=True, + constrainadps=True, + adpsymbols=None, + isosymbol="Uiso", +): + """Restricted interface to constrain_as_space_group. + + Arguments: As constrain_as_space_group, except + ----------------------------------------------- + sg + diffpy.structure.spacegroups.SpaceGroup instance + """ + from diffpy.structure.symmetryutilities import stdUsymbols + + if scatterers is None: + scatterers = phase.get_scatterers() + if adpsymbols is None: + adpsymbols = stdUsymbols + + sgp = SpaceGroupParameters( + phase, + sg, + scatterers, + sgoffset, + constrainlat, + constrainadps, + adpsymbols, + isosymbol, + ) + + return sgp + + +# End constrain_as_space_group + + +class BaseSpaceGroupParameters(RecipeContainer): + """Base class for holding space group Parameters. + + This class is used to store the variable Parameters of a structure, leaving + out those that constrained or fixed due to space group. This class has the + same Parameter attribute access of a ParameterSet. The purpose of this + class is to make it easy to access the free variables of a structure for + scripting purposes. + + Attributes + ---------- + name + "sgpars" + """ + + def __init__(self, name="sgpars"): + """Create the BaseSpaceGroupParameters object. + + This initializes the attributes. + """ + RecipeContainer.__init__(self, name) + return + + def addParameter(self, par, check=True): + """Store a Parameter. + + Parameters + ---------- + par + The Parameter to be stored. + check + If True (default), a ValueError is raised a Parameter of + the specified name has already been inserted. + + Raises ValueError if the Parameter has no name. + """ + # Store the Parameter + RecipeContainer._add_object(self, par, self._parameters, check) + return + + +# End class BaseSpaceGroupParameters + + +class SpaceGroupParameters(BaseSpaceGroupParameters): + """Class for holding and creating space group Parameters. + + This class is used to store the variable Parameters of a structure, leaving + out those that constrained or fixed due to space group. This does the work + of the constrain_as_space_group method. This class has the same Parameter + attribute access of a ParameterSet. + + Attributes + ---------- + name + "sgpars" + phase + The constrained BaseStructure object. + sg + The diffpy.structure.spacegroups.SpaceGroup object + corresponding to the space group. + sgoffset + Optional offset for the space group origin. + scatterers + The constrained scatterer ParameterSets. + constrainlat + Flag indicating whether the lattice is constrained. + constrainadps + Flag indicating whether the ADPs are constrained. + adpsymbols + A list of the ADP names. + _xyzpars + BaseSpaceGroupParameters of free xyz Parameters that are + constrained to. + xyzpars + Property that populates _xyzpars. + _latpars + BaseSpaceGroupParameters of free lattice Parameters that + are constrained to. + latpars + Property that populates _latpars. + _adppars + BaseSpaceGroupParameters of free ADPs that are constrained + to. + adppars + Property that populates _adppars. + """ + + def __init__( + self, + phase, + sg, + scatterers, + sgoffset, + constrainlat, + constrainadps, + adpsymbols, + isosymbol, + ): + """Create the SpaceGroupParameters object. + + Parameters + ---------- + phase + A BaseStructure object to be constrained. + sg + The space group number or symbol (compatible with + diffpy.structure.spacegroups.GetSpaceGroup. + sgoffset + Optional offset for sg origin. + scatterers + The scatterer ParameterSets to constrain. If scatterers + is None, then all scatterers accessible from + phase.get_scatterers will be constrained. + constrainlat + Flag indicating whether to constrain the lattice. + constrainadps + Flag indicating whether to constrain the ADPs. + adpsymbols + A list of the ADP names. The names must be given in the + same order as + diffpy.structure.symmetryutilities.stdUsymbols. + isosymbol + Symbol for isotropic ADP (default "Uiso"). If None, + isotropic ADPs will be constrained via the anisotropic + ADPs. + """ + BaseSpaceGroupParameters.__init__(self) + self._latpars = None + self._xyzpars = None + self._adppars = None + + self._parsets = {} + self._manage(self._parsets) + + self.phase = phase + self.sg = sg + self.sgoffset = sgoffset + self.scatterers = scatterers + self.constrainlat = constrainlat + self.constrainadps = constrainadps + self.adpsymbols = adpsymbols + self.isosymbol = isosymbol + + return + + def __iter__(self): + """Iterate over top-level parameters.""" + if ( + self._latpars is None + or self._xyzpars is None + or self._adppars is None + ): + self._make_constraints() + return RecipeContainer.__iter__(self) + + latpars = property(lambda self: self._get_lat_pars()) + + def _get_lat_pars(self): + """Accessor for _latpars.""" + if self._latpars is None: + self._constrain_lattice() + return self._latpars + + xyzpars = property(lambda self: self._get_xyz_pars()) + + def _get_xyz_pars(self): + """Accessor for _xyzpars.""" + positions = [] + for scatterer in self.scatterers: + xyz = [scatterer.x, scatterer.y, scatterer.z] + positions.append([p.value for p in xyz]) + if self._xyzpars is None: + self._constrain_xyzs(positions) + return self._xyzpars + + adppars = property(lambda self: self._get_adp_pars()) + + def _get_adp_pars(self): + """Accessor for _adppars.""" + positions = [] + for scatterer in self.scatterers: + xyz = [scatterer.x, scatterer.y, scatterer.z] + positions.append([p.value for p in xyz]) + if self._adppars is None: + self._constrain_adps(positions) + return self._adppars + + def _make_constraints(self): + """Constrain the structure to the space group. + + This works as described by the constrain_as_space_group method. + """ + # Start by clearing the constraints + self._clear_constraints() + + scatterers = self.scatterers + + # Prepare positions + positions = [] + for scatterer in scatterers: + xyz = [scatterer.x, scatterer.y, scatterer.z] + positions.append([p.value for p in xyz]) + + self._constrain_lattice() + self._constrain_xyzs(positions) + self._constrain_adps(positions) + + return + + def _clear_constraints(self): + """Clear old constraints. + + This only clears constraints where new ones are going to be + applied. + """ + phase = self.phase + scatterers = self.scatterers + isosymbol = self.isosymbol + adpsymbols = self.adpsymbols + + # Clear xyz + for scatterer in scatterers: + + for par in [scatterer.x, scatterer.y, scatterer.z]: + if scatterer.is_constrained(par): + scatterer.remove_constraint(par) + par.set_constant(False) + + # Clear the lattice + if self.constrainlat: + + lattice = phase.get_lattice() + latpars = [ + lattice.a, + lattice.b, + lattice.c, + lattice.alpha, + lattice.beta, + lattice.gamma, + ] + for par in latpars: + if lattice.is_constrained(par): + lattice.remove_constraint(par) + par.set_constant(False) + + # Clear ADPs + if self.constrainadps: + for scatterer in scatterers: + if isosymbol: + par = scatterer.get(isosymbol) + if par is not None: + if scatterer.is_constrained(par): + scatterer.remove_constraint(par) + par.set_constant(False) + + for pname in adpsymbols: + par = scatterer.get(pname) + if par is not None: + if scatterer.is_constrained(par): + scatterer.remove_constraint(par) + par.set_constant(False) + + return + + def _constrain_lattice(self): + """Constrain the lattice parameters.""" + if not self.constrainlat: + return + + phase = self.phase + sg = self.sg + + lattice = phase.get_lattice() + system = sg.crystal_system + if not system: + system = "Triclinic" + system = system.title() + # This makes the constraints + f = _constraintMap[system] + f(lattice) + + # Now get the unconstrained, non-constant lattice pars and store them. + self._latpars = BaseSpaceGroupParameters("latpars") + latpars = [ + lattice.a, + lattice.b, + lattice.c, + lattice.alpha, + lattice.beta, + lattice.gamma, + ] + pars = [p for p in latpars if not p.const and not p.constrained] + for par in pars: + # FIXME - the original parameter will still appear as + # constrained. + newpar = self.__add_par(par.name, par) + self._latpars.addParameter(newpar) + + return + + def _constrain_xyzs(self, positions): + """Constrain the positions. + + Parameters + ---------- + positions + The coordinates of the scatterers. + """ + from diffpy.structure.symmetryutilities import SymmetryConstraints + + sg = self.sg + sgoffset = self.sgoffset + + # We do this without ADPs here so we can skip much complication. See + # the _constrain_adps method for details. + g = SymmetryConstraints(sg, positions, sgoffset=sgoffset) + + scatterers = self.scatterers + self._xyzpars = BaseSpaceGroupParameters("xyzpars") + + # Make proxies to the free xyz parameters + xyznames = [name[:1] + "_" + name[1:] for name, val in g.pospars] + for pname in xyznames: + name, idx = pname.rsplit("_", 1) + idx = int(idx) + par = scatterers[idx].get(name) + newpar = self.__add_par(pname, par) + self._xyzpars.addParameter(newpar) + + # Constrain non-free xyz parameters + fpos = g.positionFormulas(xyznames) + for idx, tmp in enumerate(zip(scatterers, fpos)): + scatterer, fp = tmp + + # Extract the constraint equation from the formula + for parname, formula in fp.items(): + _makeconstraint( + parname, formula, scatterer, idx, self._parameters + ) + + return + + def _constrain_adps(self, positions): + """Constrain the ADPs. + + Parameters + ---------- + positions + The coordinates of the scatterers. + """ + from diffpy.structure.symmetryutilities import ( + SymmetryConstraints, + stdUsymbols, + ) + + if not self.constrainadps: + return + + sg = self.sg + sgoffset = self.sgoffset + scatterers = self.scatterers + isosymbol = self.isosymbol + adpsymbols = self.adpsymbols + adpmap = dict(zip(stdUsymbols, adpsymbols)) + self._adppars = BaseSpaceGroupParameters("adppars") + + # Prepare ADPs. Note that not all scatterers have constrainable ADPs. + # For example, MoleculeParSet from objcryststructure does not. We + # discard those. + nonadps = [] + Uijs = [] + for sidx, scatterer in enumerate(scatterers): + + pars = [scatterer.get(symb) for symb in adpsymbols] + + if None in pars: + nonadps.append(sidx) + continue + + Uij = numpy.zeros((3, 3), dtype=float) + for idx, par in enumerate(pars): + i, j = _idxtoij[idx] + Uij[i, j] = Uij[j, i] = par.getValue() + + Uijs.append(Uij) + + # Discard any positions for the nonadps + positions = list(positions) + nonadps.reverse() + [positions.pop(idx) for idx in nonadps] + + # Now we can create symmetry constraints without having to worry about + # the nonadps + g = SymmetryConstraints(sg, positions, Uijs, sgoffset=sgoffset) + + adpnames = [adpmap[name[:3]] + "_" + name[3:] for name, val in g.Upars] + + # Make proxies to the free adp parameters. We start by filtering out + # the isotropic ones so we can use the isotropic parameter. + isoidx = [] + isonames = [] + for pname in adpnames: + name, idx = pname.rsplit("_", 1) + idx = int(idx) + # Check for isotropic ADPs + scatterer = scatterers[idx] + if isosymbol and g.Uisotropy[idx] and idx not in isoidx: + isoidx.append(idx) + par = scatterer.get(isosymbol) + if par is not None: + parname = "%s_%i" % (isosymbol, idx) + newpar = self.__add_par(parname, par) + self._adppars.addParameter(newpar) + isonames.append(newpar.name) + else: + par = scatterer.get(name) + if par is not None: + newpar = self.__add_par(pname, par) + self._adppars.addParameter(newpar) + + # Constrain dependent isotropics + for idx, isoname in zip(isoidx[:], isonames): + for j in g.coremap[idx]: + if j == idx: + continue + isoidx.append(j) + scatterer = scatterers[j] + scatterer.add_constraint( + isosymbol, isoname, params=self._parameters + ) + + fadp = g.UFormulas(adpnames) + + # Constrain dependent anisotropics. We use the fact that an + # anisotropic cannot be dependent on an isotropic. + for idx, tmp in enumerate(zip(scatterers, fadp)): + if idx in isoidx: + continue + scatterer, fa = tmp + # Extract the constraint equation from the formula + for stdparname, formula in fa.items(): + pname = adpmap[stdparname] + _makeconstraint( + pname, formula, scatterer, idx, self._parameters + ) + + def __add_par(self, parname, par): + """Constrain a parameter via proxy with a specified name. + + Parameters + ---------- + par + Parameter to constrain + idx + Index to identify scatterer from which par comes + """ + newpar = ParameterProxy(parname, par) + self.addParameter(newpar) + return newpar + + +# End class SpaceGroupParameters + +# crystal system rules +# ref: Benjamin, W. A., Introduction to crystallography, +# New York (1969), p.60 + + +def _constrain_triclinic(lattice): + """Make constraints for Triclinic systems.""" + return + + +def _constrain_monoclinic(lattice): + """Make constraints for Monoclinic systems. + + alpha and beta are fixed to 90 unless alpha != beta and alpha == + gamma, in which case alpha and gamma are constrained to 90. + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + lattice.alpha.set_constant(True, ang90) + beta = lattice.beta.getValue() + gamma = lattice.gamma.getValue() + + if ang90 != beta and ang90 == gamma: + lattice.gamma.set_constant(True, ang90) + else: + lattice.beta.set_constant(True, ang90) + return + + +def _constrain_orthorhombic(lattice): + """Make constraints for Orthorhombic systems. + + alpha, beta and gamma are constrained to 90 + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + lattice.alpha.set_constant(True, ang90) + lattice.beta.set_constant(True, ang90) + lattice.gamma.set_constant(True, ang90) + return + + +def _constrain_tetragonal(lattice): + """Make constraints for Tetragonal systems. + + b is constrained to a and alpha, beta and gamma are constrained to + 90. + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + lattice.alpha.set_constant(True, ang90) + lattice.beta.set_constant(True, ang90) + lattice.gamma.set_constant(True, ang90) + lattice.add_constraint(lattice.b, lattice.a) + return + + +def _constrain_trigonal(lattice): + """Make constraints for Trigonal systems. + + If gamma == 120, then b is constrained to a, alpha and beta are + constrained to 90 and gamma is constrained to 120. Otherwise, b and + c are constrained to a, beta and gamma are constrained to alpha. + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + ang120 = 120.0 * afactor + if lattice.gamma.getValue() == ang120: + lattice.add_constraint(lattice.b, lattice.a) + lattice.alpha.set_constant(True, ang90) + lattice.beta.set_constant(True, ang90) + lattice.gamma.set_constant(True, ang120) + else: + lattice.add_constraint(lattice.b, lattice.a) + lattice.add_constraint(lattice.c, lattice.a) + lattice.add_constraint(lattice.beta, lattice.alpha) + lattice.add_constraint(lattice.gamma, lattice.alpha) + return + + +def _constrain_hexagonal(lattice): + """Make constraints for Hexagonal systems. + + b is constrained to a, alpha and beta are constrained to 90 and + gamma is constrained to 120. + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + ang120 = 120.0 * afactor + lattice.add_constraint(lattice.b, lattice.a) + lattice.alpha.set_constant(True, ang90) + lattice.beta.set_constant(True, ang90) + lattice.gamma.set_constant(True, ang120) + return + + +def _constrain_cubic(lattice): + """Make constraints for Cubic systems. + + b and c are constrained to a, alpha, beta and gamma are constrained + to 90. + """ + afactor = 1 + if lattice.angunits == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + lattice.add_constraint(lattice.b, lattice.a) + lattice.add_constraint(lattice.c, lattice.a) + lattice.alpha.set_constant(True, ang90) + lattice.beta.set_constant(True, ang90) + lattice.gamma.set_constant(True, ang90) + return + + +# This is used to map the correct crystal system to the proper constraint +# function. +_constraintMap = { + "Triclinic": _constrain_triclinic, + "Monoclinic": _constrain_monoclinic, + "Orthorhombic": _constrain_orthorhombic, + "Tetragonal": _constrain_tetragonal, + "Trigonal": _constrain_trigonal, + "Hexagonal": _constrain_hexagonal, + "Cubic": _constrain_cubic, +} + + +def _makeconstraint(parname, formula, scatterer, idx, ns={}): + """Constrain a parameter according to a formula. + + Parameters + ---------- + parname + Name of parameter + formula + Constraint formula + scatterer + scatterer containing par of parname + idx + Index to identify scatterer from which par comes + ns + namespace to draw extra names from (default {}) + + Returns + ------- + par + Returns the parameter if it is free. + """ + par = scatterer.get(parname) + + if par is None: + return + + compname = "%s_%i" % (parname, idx) + + # Check to see if this parameter is free + pat = r"%s *([+-] *\d+)?$" % compname + if re.match(pat, formula): + return par + + # Check to see if it is a constant + fval = _get_float(formula) + if fval is not None: + par.set_constant() + return + + # If we got here, then we have a constraint equation + # Fix any division issues + formula = formula.replace("/", "*1.0/") + scatterer.add_constraint(par, formula, params=ns) + return + + +def _get_float(formula): + """Get a float from a formula string, or None if this is not + possible.""" + try: + return eval(formula) + except NameError: + return None + + +# Constants needed above +_idxtoij = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)] +deg2rad = numpy.pi / 180 +rad2deg = 1.0 / deg2rad + + +# End of file diff --git a/src/diffpy/cmipdf/structure/srrealparset.py b/src/diffpy/cmipdf/structure/srrealparset.py new file mode 100644 index 0000000..8eb8403 --- /dev/null +++ b/src/diffpy/cmipdf/structure/srrealparset.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Structure wrapper class for structures compatible with SrReal.""" + +__all__ = ["SrRealParSet"] + +from diffpy.cmipdf.structure.basestructureparset import BaseStructureParSet +from diffpy.cmipdf.structure.bvsrestraint import BVSRestraint + + +class SrRealParSet(BaseStructureParSet): + """Base class for SrReal-compatible structure adapters. + + This derives from BaseStructureParSet and provides some extended + functionality provided by SrReal. + + Attributes + ---------- + stru + The adapted object + _usesymmetry + A flag indicating if SrReal calculators that operate on + this object should use symmetry. By default this is + True. + """ + + def __init__(self, *args, **kw): + BaseStructureParSet.__init__(self, *args, **kw) + self._usesymmetry = True + self.stru = None + return + + def restrain_bvs(self, sig=1, scaled=False): + """Restrain the bond-valence sum to zero. + + This adds a penalty to the cost function equal to + bvmsdiff / sig**2 + where bvmsdiff is the mean-squared difference between the calculated + and expected bond valence sums for the structure. If scaled is True, + this is also scaled by the current point-averaged chi^2 value so the + restraint is roughly equally weighted in the fit. + + Parameters + ---------- + sig + The uncertainty on the BVS (default 1). + scaled + A flag indicating if the restraint is scaled + (multiplied) by the unrestrained point-average chi^2 + (chi^2/numpoints) (default False). + + Returns the BVSRestraint object for use with the 'unrestrain' method. + """ + # Create the Restraint object + res = BVSRestraint(self, sig, scaled) + # Add it to the _restraints set + self._restraints.add(res) + # Our configuration changed. Notify observers. + self._update_configuration() + # Return the Restraint object + return res + + def use_symmetry(self, use=True): + """Set this structure to use symmetry. + + This determines how the structure is treated by SrReal + calculators. + """ + self._usesymmetry = bool(use) + return + + def using_symmetry(self): + """Check if symmetry is being used.""" + return self._usesymmetry + + def _get_srreal_structure(self): + """Get the structure object for use with SrReal calculators. + + If this is periodic, then return the structure, otherwise, pass + it inside of a nosymmetry wrapper. + """ + from diffpy.srreal.structureadapter import nosymmetry + + if self._usesymmetry: + return self.stru + return nosymmetry(self.stru) diff --git a/tests/conftest.py b/tests/conftest.py index 9245c31..8863646 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -38,6 +38,19 @@ def has_diffpy_structure(): return False +# pyobjcryst +@lru_cache() +def has_pyobjcryst(): + try: + import pyobjcryst as m + + del m + return True + except ImportError: + logger.warning("Cannot import pyobjcryst, pyobjcryst tests skipped.") + return False + + # diffpy.srreal @@ -63,6 +76,11 @@ def diffpy_srreal_available(): return has_diffpy_srreal() +@pytest.fixture(scope="session") +def pyobjcryst_available(): + return has_pyobjcryst() + + @pytest.fixture(scope="session") def datafile(): """Fixture to load a test data file from the testdata package diff --git a/tests/test_diffpyparset.py b/tests/test_diffpyparset.py new file mode 100644 index 0000000..f4b69c5 --- /dev/null +++ b/tests/test_diffpyparset.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Tests for diffpy.cmipdf.structure package.""" + +import pickle +import unittest + +import numpy as np + +from diffpy.cmipdf.structure.diffpyparset import DiffpyStructureParSet + + +def testDiffpyStructureParSet(): + """Test the structure conversion.""" + from diffpy.structure import Atom, Lattice, Structure + + a1 = Atom("Cu", xyz=np.array([0.0, 0.1, 0.2]), Uisoequiv=0.003) + a2 = Atom("Ag", xyz=np.array([0.3, 0.4, 0.5]), Uisoequiv=0.002) + lattice = Lattice(2.5, 2.5, 2.5, 90, 90, 90) + + dsstru = Structure([a1, a2], lattice) + # Structure makes copies + a1 = dsstru[0] + a2 = dsstru[1] + + s = DiffpyStructureParSet("CuAg", dsstru) + + assert s.name == "CuAg" + + def _testAtoms(): + # Check the atoms thoroughly + assert a1.element == s.Cu0.element + assert a2.element == s.Ag0.element + assert a1.Uisoequiv == s.Cu0.Uiso.getValue() + assert a2.Uisoequiv == s.Ag0.Uiso.getValue() + assert a1.Bisoequiv == s.Cu0.Biso.getValue() + assert a2.Bisoequiv == s.Ag0.Biso.getValue() + for i in range(1, 4): + for j in range(i, 4): + uijstru = getattr(a1, "U%i%i" % (i, j)) + uij = getattr(s.Cu0, "U%i%i" % (i, j)).getValue() + uji = getattr(s.Cu0, "U%i%i" % (j, i)).getValue() + assert uijstru == uij + assert uijstru == uji + bijstru = getattr(a1, "B%i%i" % (i, j)) + bij = getattr(s.Cu0, "B%i%i" % (i, j)).getValue() + bji = getattr(s.Cu0, "B%i%i" % (j, i)).getValue() + assert bijstru == bij + assert bijstru == bji + + assert a1.xyz[0] == s.Cu0.x.getValue() + assert a1.xyz[1] == s.Cu0.y.getValue() + assert a1.xyz[2] == s.Cu0.z.getValue() + return + + def _testLattice(): + + # Test the lattice + assert dsstru.lattice.a == s.lattice.a.getValue() + assert dsstru.lattice.b == s.lattice.b.getValue() + assert dsstru.lattice.c == s.lattice.c.getValue() + assert dsstru.lattice.alpha == s.lattice.alpha.getValue() + assert dsstru.lattice.beta == s.lattice.beta.getValue() + assert dsstru.lattice.gamma == s.lattice.gamma.getValue() + + _testAtoms() + _testLattice() + + # Now change some values from the diffpy Structure + a1.xyz[1] = 0.123 + a1.U11 = 0.321 + a1.B32 = 0.111 + dsstru.lattice.setLatPar(a=3.0, gamma=121) + _testAtoms() + _testLattice() + + # Now change values from the DiffpyStructureParSet + s.Cu0.x.set_value(0.456) + s.Cu0.U22.set_value(0.441) + s.Cu0.B13.set_value(0.550) + d = dsstru.lattice.dist(a1.xyz, a2.xyz) + s.lattice.b.set_value(4.6) + s.lattice.alpha.set_value(91.3) + _testAtoms() + _testLattice() + # Make sure the distance changed + assert d != dsstru.lattice.dist(a1.xyz, a2.xyz) + return + + +def test___repr__(): + """Test representation of DiffpyStructureParSet objects.""" + from diffpy.structure import Atom, Lattice, Structure + + lat = Lattice(3, 3, 2, 90, 90, 90) + atom = Atom("C", [0, 0.2, 0.5]) + stru = Structure([atom], lattice=lat) + dsps = DiffpyStructureParSet("dsps", stru) + assert repr(stru) == repr(dsps) + assert repr(lat) == repr(dsps.lattice) + assert repr(atom) == repr(dsps.atoms[0]) + return + + +def test_pickling(): + """Test pickling of DiffpyStructureParSet.""" + from diffpy.structure import Atom, Structure + + stru = Structure([Atom("C", [0, 0.2, 0.5])]) + dsps = DiffpyStructureParSet("dsps", stru) + data = pickle.dumps(dsps) + dsps2 = pickle.loads(data) + assert 1 == len(dsps2.atoms) + assert 0.2 == dsps2.atoms[0].y.value + return + + +# End of class TestParameterAdapter + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_objcrystparset.py b/tests/test_objcrystparset.py new file mode 100644 index 0000000..afa74ab --- /dev/null +++ b/tests/test_objcrystparset.py @@ -0,0 +1,678 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Tests for diffpy.cmipdf.structure package.""" + +import unittest + +import numpy +import pytest + +# Global variables to be assigned in setUp +ObjCrystCrystalParSet = spacegroups = None +Crystal = Atom = Molecule = ScatteringPowerAtom = None + + +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 a crystal containing the C60 molecule using pyobjcryst.""" + pi = numpy.pi + c = Crystal(100, 100, 100, "P1") + c.SetName("c60frame") + m = Molecule(c, "c60") + + c.AddScatterer(m) + + sp = ScatteringPowerAtom("C", "C") + sp.SetBiso(8 * pi * pi * 0.003) + # c.AddScatteringPower(sp) + + for i, l in enumerate(c60xyz.strip().splitlines()): + x, y, z = map(float, l.split()) + m.AddAtom(x, y, z, sp, "C%i" % i) + + return c + + +# ---------------------------------------------------------------------------- + + +class TestParameterAdapter: + @pytest.fixture(autouse=True) + def setup(self, pyobjcryst_available): + # shared setup + if not pyobjcryst_available: + pytest.skip("pyobjcryst package not available") + + global ObjCrystCrystalParSet, Crystal, Atom, Molecule + global ScatteringPowerAtom + from pyobjcryst.atom import Atom + from pyobjcryst.crystal import Crystal + from pyobjcryst.molecule import Molecule + from pyobjcryst.scatteringpower import ScatteringPowerAtom + + from diffpy.cmipdf.structure.objcrystparset import ( + ObjCrystCrystalParSet, + ) + + self.occryst = makeC60() + self.ocmol = self.occryst.GetScatterer("c60") + return + + def tearDown(self): + del self.occryst + del self.ocmol + return + + def testImplicitBondAngleRestraints(self): + """Test the structure with implicit bond angles.""" + occryst = self.occryst + ocmol = self.ocmol + + # Add some bond angles to the molecule + ocmol.AddBondAngle(ocmol[0], ocmol[5], ocmol[8], 1.1, 0.1, 0.1) + ocmol.AddBondAngle(ocmol[0], ocmol[7], ocmol[44], 1.3, 0.1, 0.1) + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + m.wrap_restraints() + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + res0, res1 = m._restraints + p0 = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + p1 = set([angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()]) + assert p0 == p1 + + return + + def testObjCrystParSet(self): + """Test the structure conversion.""" + occryst = self.occryst + ocmol = self.ocmol + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + assert cryst.name == "bucky" + + def _testCrystal(): + # Test the lattice + assert occryst.a == pytest.approx(cryst.a.value) + assert occryst.b == pytest.approx(cryst.b.getValue()) + assert occryst.c == pytest.approx(cryst.c.getValue()) + assert occryst.alpha == pytest.approx(cryst.alpha.getValue()) + assert occryst.beta == pytest.approx(cryst.beta.getValue()) + assert occryst.gamma == pytest.approx(cryst.gamma.getValue()) + return + + def _testMolecule(): + + # Test position / occupancy + assert ocmol.X == pytest.approx(m.x.getValue()) + assert ocmol.Y == pytest.approx(m.y.getValue()) + assert ocmol.Z == pytest.approx(m.z.getValue()) + assert ocmol.Occupancy == pytest.approx(m.occ.getValue()) + + # Test orientation + assert ocmol.Q0 == pytest.approx(m.q0.getValue()) + assert ocmol.Q1 == pytest.approx(m.q1.getValue()) + assert ocmol.Q2 == pytest.approx(m.q2.getValue()) + assert ocmol.Q3 == pytest.approx(m.q3.getValue()) + + # Check the atoms thoroughly + for i in range(len(ocmol)): + oca = ocmol[i] + ocsp = oca.GetScatteringPower() + a = m.atoms[i] + assert ocsp.GetSymbol() == a.element + assert oca.X == pytest.approx(a.x.getValue()) + assert oca.Y == pytest.approx(a.y.getValue()) + assert oca.Z == pytest.approx(a.z.getValue()) + assert oca.Occupancy == pytest.approx(a.occ.getValue()) + assert ocsp.Biso == pytest.approx(a.Biso.getValue()) + return + + _testCrystal() + _testMolecule() + + # Now change some values from ObjCryst + ocmol[0].X *= 1.1 + ocmol[0].Occupancy *= 1.1 + ocmol[0].GetScatteringPower().Biso *= 1.1 + ocmol.Q0 *= 1.1 + occryst.a *= 1.1 + + _testCrystal() + _testMolecule() + + # Now change values from the StructureParSet + cryst.c60.C44.x.set_value(1.1) + cryst.c60.C44.occ.set_value(1.1) + cryst.c60.C44.Biso.set_value(1.1) + cryst.c60.q3.set_value(1.1) + cryst.a.set_value(1.1) + + _testCrystal() + _testMolecule() + return + + def testImplicitBondLengthRestraints(self): + """Test the structure with implicit bond lengths.""" + occryst = self.occryst + ocmol = self.ocmol + + # Add some bonds to the molecule + ocmol.AddBond(ocmol[0], ocmol[5], 3.3, 0.1, 0.1) + ocmol.AddBond(ocmol[0], ocmol[7], 3.3, 0.1, 0.1) + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + m.wrap_restraints() + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + res0, res1 = m._restraints + p0 = set([res0.penalty(), res1.penalty()]) + bonds = ocmol.GetBondList() + p1 = set([bonds[0].GetLogLikelihood(), bonds[1].GetLogLikelihood()]) + assert p0 == p1 + + return + + def testImplicitDihedralAngleRestraints(self): + """Test the structure with implicit dihedral angles.""" + occryst = self.occryst + ocmol = self.ocmol + + # Add some bond angles to the molecule + ocmol.AddDihedralAngle( + ocmol[0], ocmol[5], ocmol[8], ocmol[41], 1.1, 0.1, 0.1 + ) + ocmol.AddDihedralAngle( + ocmol[0], ocmol[7], ocmol[44], ocmol[2], 1.3, 0.1, 0.1 + ) + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + m.wrap_restraints() + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + res0, res1 = m._restraints + p0 = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + p1 = set([angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()]) + assert p0 == p1 + + return + + def testImplicitStretchModes(self): + """Test the molecule with implicit stretch modes.""" + # Not sure how to make this happen. + pass + + def testExplicitBondLengthRestraints(self): + """Test the structure with explicit bond lengths.""" + occryst = self.occryst + ocmol = self.ocmol + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + # make some bond angle restraints + res0 = m.restrain_bond_length(m.atoms[0], m.atoms[5], 3.3, 0.1, 0.1) + res1 = m.restrain_bond_length(m.atoms[0], m.atoms[7], 3.3, 0.1, 0.1) + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + p0 = [res0.penalty(), res1.penalty()] + bonds = ocmol.GetBondList() + assert 2 == len(bonds) + p1 = [b.GetLogLikelihood() for b in bonds] + assert p0 == p1 + + return + + def testExplicitBondAngleRestraints(self): + """Test the structure with explicit bond angles. + + Note that this cannot work with co-linear points as the + direction of rotation cannot be defined in this case. + """ + occryst = self.occryst + ocmol = self.ocmol + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + # restrain some bond angles + res0 = m.restrain_bond_angle( + m.atoms[0], m.atoms[5], m.atoms[8], 3.3, 0.1, 0.1 + ) + res1 = m.restrain_bond_angle( + m.atoms[0], m.atoms[7], m.atoms[44], 3.3, 0.1, 0.1 + ) + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + p0 = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + p1 = set([angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()]) + assert p0 == p1 + + return + + def testExplicitDihedralAngleRestraints(self): + """Test the structure with explicit dihedral angles.""" + occryst = self.occryst + ocmol = self.ocmol + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + # Restrain some dihedral angles. + res0 = m.restrain_dihedral_angle( + m.atoms[0], m.atoms[5], m.atoms[8], m.atoms[41], 1.1, 0.1, 0.1 + ) + res1 = m.restrain_dihedral_angle( + m.atoms[0], m.atoms[7], m.atoms[44], m.atoms[2], 1.1, 0.1, 0.1 + ) + + # make sure that we have some restraints in the molecule + assert 2 == len(m._restraints) + + # make sure these evaluate to whatver we get from objcryst + p0 = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + p1 = set([angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()]) + assert p0 == p1 + + return + + def testExplicitBondLengthParameter(self): + """Test adding bond length parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + + # Add a parameter + p1 = m.add_bond_length_parameter("C07", a0, a7) + # Have another atom tag along for the ride + p1.add_atoms([a20]) + + xyz0 = numpy.array([a0.x.getValue(), a0.y.getValue(), a0.z.getValue()]) + xyz7 = numpy.array([a7.x.getValue(), a7.y.getValue(), a7.z.getValue()]) + xyz20 = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + + dd = xyz0 - xyz7 + d0 = numpy.dot(dd, dd) ** 0.5 + assert d0 == pytest.approx(p1.getValue(), abs=1e-6) + + # Record the unit direction of change for later + u = dd / d0 + + # Change the value + scale = 1.05 + p1.set_value(scale * d0) + + # Verify that it has changed. + assert scale * d0 == pytest.approx(p1.getValue(), abs=1e-6) + + xyz0a = numpy.array( + [a0.x.getValue(), a0.y.getValue(), a0.z.getValue()] + ) + xyz7a = numpy.array( + [a7.x.getValue(), a7.y.getValue(), a7.z.getValue()] + ) + xyz20a = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + + dda = xyz0a - xyz7a + d1 = numpy.dot(dda, dda) ** 0.5 + + assert scale * d0 == pytest.approx(d1, abs=1e-6) + + # Verify that only the second and third atoms have moved. + + assert numpy.array_equal(xyz0, xyz0a) + + xyz7calc = xyz7 + (1 - scale) * d0 * u + for i in range(3): + assert xyz7a[i] == pytest.approx(xyz7calc[i], abs=1e-5) + + xyz20calc = xyz20 + (1 - scale) * d0 * u + for i in range(3): + assert xyz20a[i] == pytest.approx(xyz20calc[i], abs=1e-6) + + return + + def testExplicitBondAngleParameter(self): + """Test adding bond angle parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + a25 = m.atoms[25] + + xyz0 = numpy.array([a0.x.getValue(), a0.y.getValue(), a0.z.getValue()]) + xyz7 = numpy.array([a7.x.getValue(), a7.y.getValue(), a7.z.getValue()]) + xyz20 = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + xyz25 = numpy.array( + [a25.x.getValue(), a25.y.getValue(), a25.z.getValue()] + ) + + v1 = xyz7 - xyz0 + d1 = numpy.dot(v1, v1) ** 0.5 + v2 = xyz7 - xyz20 + d2 = numpy.dot(v2, v2) ** 0.5 + + angle0 = numpy.arccos(numpy.dot(v1, v2) / (d1 * d2)) + + # Add a parameter + p1 = m.add_bond_angle_parameter("C0720", a0, a7, a20) + # Have another atom tag along for the ride + p1.add_atoms([a25]) + + assert angle0 == pytest.approx(p1.getValue(), abs=1e-6) + + # Change the value + scale = 1.05 + p1.set_value(scale * angle0) + + # Verify that it has changed. + assert scale * angle0 == pytest.approx(p1.getValue(), abs=1e-6) + + xyz0a = numpy.array( + [a0.x.getValue(), a0.y.getValue(), a0.z.getValue()] + ) + xyz7a = numpy.array( + [a7.x.getValue(), a7.y.getValue(), a7.z.getValue()] + ) + xyz20a = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + xyz25a = numpy.array( + [a25.x.getValue(), a25.y.getValue(), a25.z.getValue()] + ) + + v1a = xyz7a - xyz0a + d1a = numpy.dot(v1a, v1a) ** 0.5 + v2a = xyz7a - xyz20a + d2a = numpy.dot(v2a, v2a) ** 0.5 + + angle1 = numpy.arccos(numpy.dot(v1a, v2a) / (d1a * d2a)) + + assert scale * angle0 == pytest.approx(angle1, abs=1e-6) + + # Verify that only the last two atoms have moved. + + assert numpy.array_equal(xyz0, xyz0a) + assert numpy.array_equal(xyz7, xyz7a) + assert not numpy.array_equal(xyz20, xyz20a) + assert not numpy.array_equal(xyz25, xyz25a) + + return + + def testExplicitDihedralAngleParameter(self): + """Test adding dihedral angle parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + cryst = ObjCrystCrystalParSet("bucky", occryst) + m = cryst.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + a25 = m.atoms[25] + a33 = m.atoms[33] + + xyz0 = numpy.array([a0.x.getValue(), a0.y.getValue(), a0.z.getValue()]) + xyz7 = numpy.array([a7.x.getValue(), a7.y.getValue(), a7.z.getValue()]) + xyz20 = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + xyz25 = numpy.array( + [a25.x.getValue(), a25.y.getValue(), a25.z.getValue()] + ) + xyz33 = numpy.array( + [a33.x.getValue(), a33.y.getValue(), a33.z.getValue()] + ) + + v12 = xyz0 - xyz7 + v23 = xyz7 - xyz20 + v34 = xyz20 - xyz25 + v123 = numpy.cross(v12, v23) + v234 = numpy.cross(v23, v34) + + d123 = numpy.dot(v123, v123) ** 0.5 + d234 = numpy.dot(v234, v234) ** 0.5 + angle0 = -numpy.arccos(numpy.dot(v123, v234) / (d123 * d234)) + + # Add a parameter + p1 = m.add_dihedral_angle_parameter("C072025", a0, a7, a20, a25) + # Have another atom tag along for the ride + p1.add_atoms([a33]) + + assert angle0 == pytest.approx(p1.getValue(), abs=1e-6) + + # Change the value + scale = 1.05 + p1.set_value(scale * angle0) + + # Verify that it has changed. + assert scale * angle0 == pytest.approx(p1.getValue(), abs=1e-6) + + xyz0a = numpy.array( + [a0.x.getValue(), a0.y.getValue(), a0.z.getValue()] + ) + xyz7a = numpy.array( + [a7.x.getValue(), a7.y.getValue(), a7.z.getValue()] + ) + xyz20a = numpy.array( + [a20.x.getValue(), a20.y.getValue(), a20.z.getValue()] + ) + xyz25a = numpy.array( + [a25.x.getValue(), a25.y.getValue(), a25.z.getValue()] + ) + xyz33a = numpy.array( + [a33.x.getValue(), a33.y.getValue(), a33.z.getValue()] + ) + + v12a = xyz0a - xyz7a + v23a = xyz7a - xyz20a + v34a = xyz20a - xyz25a + v123a = numpy.cross(v12a, v23a) + v234a = numpy.cross(v23a, v34a) + + d123a = numpy.dot(v123a, v123a) ** 0.5 + d234a = numpy.dot(v234a, v234a) ** 0.5 + angle1 = -numpy.arccos(numpy.dot(v123a, v234a) / (d123a * d234a)) + assert scale * angle0 == pytest.approx(angle1, abs=1e-6) + + # Verify that only the last two atoms have moved. + + assert numpy.array_equal(xyz0, xyz0a) + assert numpy.array_equal(xyz7, xyz7a) + assert numpy.array_equal(xyz20, xyz20a) + assert not numpy.array_equal(xyz25, xyz25a) + assert not numpy.array_equal(xyz33, xyz33a) + + return + + +class TestCreateSpaceGroup: + """Test space group creation from pyobjcryst structures. + + This makes sure that the space groups created by the structure + parameter set are correct. + """ + + @pytest.fixture(autouse=True) + def setup(self, diffpy_structure_available, pyobjcryst_available): + # shared setup + if not diffpy_structure_available: + pytest.skip("diffpy.structure package not available") + if not pyobjcryst_available: + pytest.skip("pyobjcryst package not available") + + global ObjCrystCrystalParSet, spacegroups + from diffpy.cmipdf.structure.objcrystparset import ( + ObjCrystCrystalParSet, + ) + from diffpy.structure import spacegroups + + @staticmethod + def getObjCrystParSetSpaceGroup(sg): + """Make an ObjCrystCrystalParSet with the proper space group.""" + from pyobjcryst.spacegroup import SpaceGroup + + sgobjcryst = SpaceGroup(sg.short_name) + sgnew = ObjCrystCrystalParSet._create_space_group(sgobjcryst) + return sgnew + + @staticmethod + def hashDiffPySpaceGroup(sg): + lines = [str(sg.number % 1000)] + sorted(map(str, sg.iter_symops())) + s = "\n".join(lines) + return s + + def sgsEquivalent(self, sg1, sg2): + """Check to see if two space group objects are the same.""" + hash1 = self.hashDiffPySpaceGroup(sg1) + hash2 = self.hashDiffPySpaceGroup(sg2) + return hash1 == hash2 + + # FIXME: only about 50% of the spacegroups pass the assertion + # test disabled even if cctbx is installed + def xtestCreateSpaceGroup(self): + """Check all sgtbx space groups for proper conversion to + SpaceGroup.""" + try: + from cctbx import sgtbx + except ImportError: + return + + for smbls in sgtbx.space_group_symbol_iterator(): + shn = smbls.hermann_mauguin() + short_name = shn.replace(" ", "") + if spacegroups.IsSpaceGroupIdentifier(short_name): + sg = spacegroups.GetSpaceGroup(shn) + sgnew = self.getObjCrystParSetSpaceGroup(sg) + # print("dbsg: " + repr(self.sgsEquivalent(sg, sgnew))) + assert self.sgsEquivalent(sg, sgnew) + return + + +# End of class TestCreateSpaceGroup + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_sgconstraints.py b/tests/test_sgconstraints.py new file mode 100644 index 0000000..1d3cc6f --- /dev/null +++ b/tests/test_sgconstraints.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2025 Simon Billinge. +# All rights reserved. +# +# File coded by: Caden Myers, Simon Billinge, and members of the Billinge +# group. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmipdf/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Tests space group constraints.""" + +import unittest + +import numpy +import pytest + +# ---------------------------------------------------------------------------- + + +def test_ObjCryst_constrain_space_group(pyobjcryst_available): + """Make sure that all Parameters are constrained properly. + + This tests constrainSpaceGroup from + diffpy.cmipdf.structure.sgconstraints, which is performed + automatically when an ObjCrystCrystalParSet is created. + """ + if not pyobjcryst_available: + pytest.skip("pyobjcrysta package not available") + + from diffpy.cmipdf.structure.objcrystparset import ObjCrystCrystalParSet + + pi = numpy.pi + + occryst = makeLaMnO3() + stru = ObjCrystCrystalParSet(occryst.GetName(), occryst) + # Make sure we actually create the constraints + stru._constrain_space_group() + # Make the space group parameters individually + stru.sgpars.latpars + stru.sgpars.xyzpars + stru.sgpars.adppars + + # Check the orthorhombic lattice + lattice = stru.get_lattice() + assert lattice.alpha.const + assert lattice.beta.const + assert lattice.gamma.const + assert pi / 2 == lattice.alpha.getValue() + assert pi / 2 == lattice.beta.getValue() + assert pi / 2 == lattice.gamma.getValue() + + assert not lattice.a.const + assert not lattice.b.const + assert not lattice.c.const + assert 0 == len(lattice._constraints) + + # Now make sure the scatterers are constrained properly + scatterers = stru.get_scatterers() + la = scatterers[0] + assert not la.x.const + assert not la.y.const + assert la.z.const + assert 0 == len(la._constraints) + + mn = scatterers[1] + assert mn.x.const + assert mn.y.const + assert mn.z.const + assert 0 == len(mn._constraints) + + o1 = scatterers[2] + assert not o1.x.const + assert not o1.y.const + assert o1.z.const + assert 0 == len(o1._constraints) + + o2 = scatterers[3] + assert not o2.x.const + assert not o2.y.const + assert not o2.z.const + assert 0 == len(o2._constraints) + + # Make sure we can't constrain these + with pytest.raises(ValueError): + mn.add_constraint(mn.x, "y") + + with pytest.raises(ValueError): + mn.add_constraint(mn.y, "z") + + with pytest.raises(ValueError): + mn.add_constraint(mn.z, "x") + + # Nor can we make them into variables + from diffpy.srfit.fitbase.fitrecipe import FitRecipe + + f = FitRecipe() + with pytest.raises(ValueError): + f.add_variable(mn.x) + + return + + +def test_DiffPy_constrain_as_space_group(datafile, pyobjcryst_available): + """Test the constrain_as_space_group function.""" + if not pyobjcryst_available: + pytest.skip("pyobjcrysta package not available") + + from diffpy.cmipdf.structure.diffpyparset import DiffpyStructureParSet + from diffpy.cmipdf.structure.sgconstraints import constrain_as_space_group + + stru = makeLaMnO3_P1(datafile) + parset = DiffpyStructureParSet("LaMnO3", stru) + + sgpars = constrain_as_space_group( + parset, + "P b n m", + scatterers=parset.get_scatterers()[::2], + constrainadps=True, + ) + + # Make sure that the new parameters were created + for par in sgpars: + assert par is not None + assert par.getValue() is not None + + # Test the unconstrained atoms + for scatterer in parset.get_scatterers()[1::2]: + assert not scatterer.x.const + assert not scatterer.y.const + assert not scatterer.z.const + assert not scatterer.U11.const + assert not scatterer.U22.const + assert not scatterer.U33.const + assert not scatterer.U12.const + assert not scatterer.U13.const + assert not scatterer.U23.const + assert 0 == len(scatterer._constraints) + + proxied = [p.par for p in sgpars] + + def _consttest(par): + return par.const + + def _constrainedtest(par): + return par.constrained + + def _proxytest(par): + return par in proxied + + def _alltests(par): + return _consttest(par) or _constrainedtest(par) or _proxytest(par) + + for idx, scatterer in enumerate(parset.get_scatterers()[::2]): + # Under this scheme, atom 6 is free to vary + test = False + for par in [scatterer.x, scatterer.y, scatterer.z]: + test |= _alltests(par) + assert test + + test = False + for par in [ + scatterer.U11, + scatterer.U22, + scatterer.U33, + scatterer.U12, + scatterer.U13, + scatterer.U23, + ]: + test |= _alltests(par) + + assert test + + return + + +def test_constrain_as_space_group_args(pyobjcryst_available, datafile): + """Test the arguments processing of constrain_as_space_group + function.""" + if not pyobjcryst_available: + pytest.skip("pyobjcrysta package not available") + + from diffpy.cmipdf.structure.diffpyparset import DiffpyStructureParSet + from diffpy.cmipdf.structure.sgconstraints import constrain_as_space_group + from diffpy.structure.spacegroups import GetSpaceGroup + + stru = makeLaMnO3_P1(datafile) + parset = DiffpyStructureParSet("LaMnO3", stru) + sgpars = constrain_as_space_group(parset, "P b n m") + sg = GetSpaceGroup("P b n m") + parset2 = DiffpyStructureParSet("LMO", makeLaMnO3_P1(datafile)) + sgpars2 = constrain_as_space_group(parset2, sg) + list(sgpars) + list(sgpars2) + assert sgpars.names == sgpars2.names + return + + +def makeLaMnO3_P1(datafile): + from diffpy.structure import Structure + + stru = Structure() + stru.read(datafile("LaMnO3.stru")) + return stru + + +def makeLaMnO3(): + from pyobjcryst.atom import Atom + from pyobjcryst.crystal import Crystal + from pyobjcryst.scatteringpower import ScatteringPowerAtom + + pi = numpy.pi + # It appears that ObjCryst only supports standard symbols + crystal = Crystal(5.486341, 5.619215, 7.628206, "P b n m") + crystal.SetName("LaMnO3") + # La1 + sp = ScatteringPowerAtom("La1", "La") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.996096, 0.0321494, 0.25, "La1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # Mn1 + sp = ScatteringPowerAtom("Mn1", "Mn") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0, 0.5, 0, "Mn1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O1 + sp = ScatteringPowerAtom("O1", "O") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.0595746, 0.496164, 0.25, "O1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O2 + sp = ScatteringPowerAtom("O2", "O") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.720052, 0.289387, 0.0311126, "O2", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + + return crystal + + +# ---------------------------------------------------------------------------- + +if __name__ == "__main__": + unittest.main() diff --git a/tests/testdata/LaMnO3.stru b/tests/testdata/LaMnO3.stru new file mode 100644 index 0000000..044869a --- /dev/null +++ b/tests/testdata/LaMnO3.stru @@ -0,0 +1,129 @@ +title Cell structure file of LaMnO3.0 +format pdffit +scale 1.000000 +sharp 0.000000, 0.000000, 1.000000, 3.500000 +spcgr Pbnm +cell 5.486341, 5.619215, 7.628206, 90.000000, 90.000000, 90.000000 +dcell 0.000118, 0.000156, 0.000118, 0.000000, 0.000000, 0.000000 +ncell 1, 1, 1, 20 +atoms +LA 0.99609631 0.03214940 0.25000000 1.0000 + 0.00003041 0.00000852 0.00000000 0.0000 + 0.00253993 0.00253993 0.00253993 + 0.00000214 0.00000214 0.00000214 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +LA 0.49609631 0.46785060 0.75000000 1.0000 + 0.00003041 0.00000852 0.00000000 0.0000 + 0.00253993 0.00253993 0.00253993 + 0.00000214 0.00000214 0.00000214 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +LA 0.00390369 0.96785063 0.75000000 1.0000 + 0.00003041 0.00000852 0.00000000 0.0000 + 0.00253993 0.00253993 0.00253993 + 0.00000214 0.00000214 0.00000214 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +LA 0.50390369 0.53214937 0.25000000 1.0000 + 0.00003041 0.00000852 0.00000000 0.0000 + 0.00253993 0.00253993 0.00253993 + 0.00000214 0.00000214 0.00000214 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +MN 0.00000000 0.50000000 0.00000000 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00065337 0.00065337 0.00065337 + 0.00000165 0.00000165 0.00000165 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +MN 0.50000000 0.00000000 0.00000000 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00065337 0.00065337 0.00065337 + 0.00000165 0.00000165 0.00000165 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +MN 0.00000000 0.50000000 0.50000000 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00065337 0.00065337 0.00065337 + 0.00000165 0.00000165 0.00000165 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +MN 0.50000000 0.00000000 0.50000000 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00065337 0.00065337 0.00065337 + 0.00000165 0.00000165 0.00000165 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.05957463 0.49616399 0.25000000 1.0000 + 0.00001546 0.00001610 0.00000000 0.0000 + 0.00082010 0.00082010 0.00082010 + 0.00000137 0.00000137 0.00000137 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.55957460 0.00383601 0.75000000 1.0000 + 0.00001546 0.00001610 0.00000000 0.0000 + 0.00082010 0.00082010 0.00082010 + 0.00000137 0.00000137 0.00000137 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.94042540 0.50383604 0.75000000 1.0000 + 0.00001546 0.00001610 0.00000000 0.0000 + 0.00082010 0.00082010 0.00082010 + 0.00000137 0.00000137 0.00000137 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.44042537 0.99616396 0.25000000 1.0000 + 0.00001546 0.00001610 0.00000000 0.0000 + 0.00082010 0.00082010 0.00082010 + 0.00000137 0.00000137 0.00000137 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.72005206 0.28938726 0.03111255 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.22005206 0.21061274 0.96888745 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.27994794 0.71061277 0.53111255 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.77994794 0.78938723 0.46888745 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.27994794 0.71061277 0.96888745 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.77994794 0.78938723 0.03111255 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.72005206 0.28938726 0.46888745 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +O 0.22005206 0.21061274 0.53111255 1.0000 + 0.00001528 0.00001560 0.00002506 0.0000 + 0.00512371 0.00512371 0.00512371 + 0.00000153 0.00000153 0.00000153 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 From 4f41b5facb9cde5fd6f3c6f9d74980b19f2e1a2a Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Wed, 16 Sep 2026 10:37:45 -0400 Subject: [PATCH 2/3] add pyobjcryst to test.txt file --- requirements/tests.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements/tests.txt b/requirements/tests.txt index a727786..d888d99 100644 --- a/requirements/tests.txt +++ b/requirements/tests.txt @@ -4,3 +4,4 @@ codecov coverage pytest-cov pytest-env +pyobjcryst From 8748233e4ef1ff6b8830fe73bc2d4f9e196af29b Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Thu, 17 Sep 2026 14:36:35 -0400 Subject: [PATCH 3/3] news --- news/add-structure-subpackage.rst | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/news/add-structure-subpackage.rst b/news/add-structure-subpackage.rst index 9507765..413fc2d 100644 --- a/news/add-structure-subpackage.rst +++ b/news/add-structure-subpackage.rst @@ -1,15 +1,10 @@ **Added:** -* ``diffpy.cmipdf.structure`` subpackage, migrated from ``diffpy.srfit.structure``. - It adapts ``diffpy.structure``, ``pyobjcryst`` and ``cctbx`` structures to the - ``ParameterSet`` interface and generates space group constraints. The migrated - API is snake_case only; ``getValue`` and ``addParameter`` stay camelCase until - core srfit renames them. +* No news needed: add the ``structure`` subpackage. **Changed:** -* ``BasePDFGenerator`` now takes ``struToParameterSet`` from - ``diffpy.cmipdf.structure`` instead of ``diffpy.srfit.structure``. +* **Deprecated:** @@ -21,8 +16,7 @@ **Fixed:** -* ``src/diffpy/__init__.py`` now extends the ``diffpy`` namespace path, so - ``diffpy.cmipdf`` no longer shadows the other ``diffpy`` packages. +* **Security:**