diff --git a/news/add-structure-subpackage.rst b/news/add-structure-subpackage.rst new file mode 100644 index 0000000..413fc2d --- /dev/null +++ b/news/add-structure-subpackage.rst @@ -0,0 +1,23 @@ +**Added:** + +* No news needed: add the ``structure`` subpackage. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* 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 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