diff --git a/.github/workflows/check-news-item.yml b/.github/workflows/check-news-item.yml index 20440f6..aa0fc9d 100644 --- a/.github/workflows/check-news-item.yml +++ b/.github/workflows/check-news-item.yml @@ -3,7 +3,7 @@ name: Check for News on: pull_request_target: branches: - - main # GitHub does not evaluate expressions in trigger filters; edit this value if your base branch is not main + - main # GitHub does not evaluate expressions in trigger filters; edit this value if your base branch is not main jobs: check-news-item: diff --git a/LICENSE_DANSE.rst b/LICENSE_DANSE.rst new file mode 100644 index 0000000..c92ca2d --- /dev/null +++ b/LICENSE_DANSE.rst @@ -0,0 +1,34 @@ +This program is part of the DiffPy and DANSE open-source projects at Columbia +University and is available subject to the conditions and terms laid out below. + +Copyright (c) 2008-2011, The Trustees of Columbia University in +the City of New York. All rights reserved. + +For more information please visit the diffpy web-page at + http://www.diffpy.org +or email Prof. Simon Billinge at sb2896@columbia.edu. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither the names of COLUMBIA UNIVERSITY, MICHIGAN STATE UNIVERSITY nor the + names of their contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/news/migrate-structure.rst b/news/migrate-structure.rst new file mode 100644 index 0000000..de444e5 --- /dev/null +++ b/news/migrate-structure.rst @@ -0,0 +1,23 @@ +**Added:** + +* No news added. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/pyproject.toml b/pyproject.toml index 8fd77f0..72f1f89 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ maintainers = [ description = "diffpy.cmi package for doing refinements with structure objects" keywords = ['diffraction', 'PDF', 'X-ray', 'neutron'] readme = "README.rst" -requires-python = ">=3.12, <3.15" +requires-python = ">=3.11, <3.14" classifiers = [ 'Development Status :: 5 - Production/Stable', 'Environment :: Console', @@ -25,9 +25,9 @@ classifiers = [ 'Operating System :: Microsoft :: Windows', 'Operating System :: POSIX', 'Operating System :: Unix', + 'Programming Language :: Python :: 3.11', 'Programming Language :: Python :: 3.12', 'Programming Language :: Python :: 3.13', - 'Programming Language :: Python :: 3.14', 'Topic :: Scientific/Engineering :: Physics', 'Topic :: Scientific/Engineering :: Chemistry', ] @@ -49,7 +49,7 @@ exclude = [] # exclude packages matching these glob patterns (empty by default) namespaces = false # to disable scanning PEP 420 namespaces (true by default) [project.scripts] -diffpy-cmistructure = "diffpy.cmistructure.app:main" +diffpy-cmistructure = "diffpy.cmistructure.cmistructure_app:main" [tool.setuptools.dynamic] dependencies = {file = ["requirements/pip.txt"]} diff --git a/requirements/conda.txt b/requirements/conda.txt index 24ce15a..6b1daa8 100644 --- a/requirements/conda.txt +++ b/requirements/conda.txt @@ -1 +1,4 @@ numpy +diffpy.srreal +diffpy.srfit +diffpy.structure diff --git a/requirements/pip.txt b/requirements/pip.txt index 24ce15a..6b1daa8 100644 --- a/requirements/pip.txt +++ b/requirements/pip.txt @@ -1 +1,4 @@ numpy +diffpy.srreal +diffpy.srfit +diffpy.structure 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 bfbc26c..e4af1d0 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/cmistructure/__init__.py b/src/diffpy/cmistructure/__init__.py index 77affe7..22451be 100644 --- a/src/diffpy/cmistructure/__init__.py +++ b/src/diffpy/cmistructure/__init__.py @@ -12,11 +12,65 @@ # See LICENSE.rst for license information. # ############################################################################## -"""diffpy.cmi package for doing refinements with structure objects""" +"""Modules and classes that adapt structure representations to the +ParameterSet interface and automatic structure constraint generation +from space group information.""" + +from diffpy.cmistructure.sgconstraints import constrain_as_space_group # package version from diffpy.cmistructure.version import __version__ # noqa +__all__ = ["constrain_as_space_group", "structure_to_parameter_set"] + + +def structure_to_parameter_set(name, structure): + """Create a ParameterSet adapted to a structure object. + + The adapter is chosen from the type of `structure`. Supported types are + diffpy.structure.Structure, pyobjcryst.crystal.Crystal, + pyobjcryst.molecule.Molecule and cctbx.crystal.special_position_settings. + + Parameters + ---------- + name : str + The name to give the structure. + structure : object + The structure object to adapt. + + Returns + ------- + BaseStructureParSet + The ParameterSet adapting `structure`. + + Raises + ------ + TypeError + If `structure` is not one of the supported structure types. + """ + from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet + + if DiffpyStructureParSet.can_adapt(structure): + return DiffpyStructureParSet(name, structure) + + from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet + + if ObjCrystCrystalParSet.can_adapt(structure): + return ObjCrystCrystalParSet(name, structure) + + from diffpy.cmistructure.objcrystparset import ObjCrystMoleculeParSet + + if ObjCrystMoleculeParSet.can_adapt(structure): + return ObjCrystMoleculeParSet(name, structure) + + from diffpy.cmistructure.cctbxparset import CCTBXCrystalParSet + + if CCTBXCrystalParSet.can_adapt(structure): + return CCTBXCrystalParSet(name, structure) + + raise TypeError("Unadaptable structure format") + + # silence the pyflakes syntax checker assert __version__ or True diff --git a/src/diffpy/cmistructure/basestructureparset.py b/src/diffpy/cmistructure/basestructureparset.py new file mode 100644 index 0000000..1dcae91 --- /dev/null +++ b/src/diffpy/cmistructure/basestructureparset.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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 + ---------- + structure : object + The adapted structure object. + """ + + @classmethod + def can_adapt(self, structure): + """Return whether the structure can be adapted by this class. + + Parameters + ---------- + structure : object + The structure object to check. + + Returns + ------- + bool + The flag indicating if `structure` can be adapted. The base class + always returns False. + """ + return False + + def get_lattice(self): + """Return the 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". + + The lattice must also have the "angle_units" attribute, which is + either "deg" or "rad", to signify degrees or radians. + + Returns + ------- + ParameterSet + The ParameterSet holding the lattice Parameters. + + Raises + ------ + NotImplementedError + If the subclass does not override this method. + """ + raise NotImplementedError("The must be overloaded") + + def get_scatterers(self): + """Return the list of ParameterSets that represent 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. + + Returns + ------- + list of ParameterSet + The ParameterSets of the scatterers in the structure. + + Raises + ------ + NotImplementedError + If the subclass does not override this method. + """ + raise NotImplementedError("The must be overloaded") diff --git a/src/diffpy/cmistructure/bvsrestraint.py b/src/diffpy/cmistructure/bvsrestraint.py new file mode 100644 index 0000000..9d3bc85 --- /dev/null +++ b/src/diffpy/cmistructure/bvsrestraint.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2010 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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 : BVSCalculator + The SrReal BVSCalculator instance. + _parameter_set : SrRealParSet + The SrRealParSet that created this BVSRestraint. + sig : float + The uncertainty on the BVS (default 1). + scaled : bool + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + + def __init__(self, parameter_set, sig=1, scaled=False): + """Initialize the Restraint. + + Parameters + ---------- + parameter_set : SrRealParSet + The SrRealParSet that creates this BVSRestraint. + sig : float, optional + The uncertainty on the BVS (default 1). + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + from diffpy.srreal.bvscalculator import BVSCalculator + + self._calc = BVSCalculator() + self._parameter_set = parameter_set + self.sig = float(sig) + self.scaled = bool(scaled) + return + + def penalty(self, w=1.0): + """Calculate the penalty of the restraint. + + Parameters + ---------- + w : float, optional + The point-average chi^2 which is optionally used to scale the + penalty (default 1.0). + + Returns + ------- + float + The bond-valence penalty. + """ + # Get the bvms from the BVSCalculator + structure = self._parameter_set._get_srreal_structure() + self._calc.eval(structure) + 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/cmistructure/cctbxparset.py b/src/diffpy/cmistructure/cctbxparset.py new file mode 100644 index 0000000..34161cd --- /dev/null +++ b/src/diffpy/cmistructure/cctbxparset.py @@ -0,0 +1,364 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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.cmistructure.basestructureparset import BaseStructureParSet +from diffpy.srfit.fitbase.parameter import ParameterAdapter +from diffpy.srfit.fitbase.parameterset import ParameterSet + +__all__ = ["CCTBXScattererParSet", "CCTBXUnitCellParSet", "CCTBXCrystalParSet"] + + +class CCTBXScattererParSet(ParameterSet): + """Adapt a cctbx.xray.scatterer to the ParameterSet interface. + + This class derives from ParameterSet. + + Attributes + ---------- + name : str + The 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 : ParameterAdapter + The atom position in crystal coordinates. + occupancy : ParameterAdapter + The occupancy of the atom on its crystal location. + Uiso : ParameterAdapter + The isotropic displacement factor of the atom. + """ + + def __init__(self, name, structure_parameter_set, index): + """Initialize the scatterer ParameterSet. + + Parameters + ---------- + name : str + The name of this scatterer. + structure_parameter_set : CCTBXCrystalParSet + The CCTBXCrystalParSet that contains the cctbx structure. + index : int + The index of the scatterer in the structure. + """ + ParameterSet.__init__(self, name) + self.structure_parameter_set = structure_parameter_set + self.index = index + + # x, y, z, occupancy + self.add_parameter( + ParameterAdapter("x", None, self._xyzgetter(0), self._xyzsetter(0)) + ) + self.add_parameter( + ParameterAdapter("y", None, self._xyzgetter(1), self._xyzsetter(1)) + ) + self.add_parameter( + ParameterAdapter("z", None, self._xyzgetter(2), self._xyzsetter(2)) + ) + self.add_parameter( + ParameterAdapter("occupancy", None, self._getocc, self._setocc) + ) + self.add_parameter( + ParameterAdapter("Uiso", None, self._getuiso, self._setuiso) + ) + return + + # Getters and setters + + def _xyzgetter(self, i): + + def f(dummy): + return self.structure_parameter_set.structure.scatterers()[ + self.index + ].site[i] + + return f + + def _xyzsetter(self, i): + + def f(dummy, value): + xyz = list( + self.structure_parameter_set.structure.scatterers()[ + self.index + ].site + ) + xyz[i] = value + self.structure_parameter_set.structure.scatterers()[ + self.index + ].site = tuple(xyz) + return + + return f + + def _getocc(self, dummy): + return self.structure_parameter_set.structure.scatterers()[ + self.index + ].occupancy + + def _setocc(self, dummy, value): + self.structure_parameter_set.structure.scatterers()[ + self.index + ].occupancy = value + return + + def _getuiso(self, dummy): + return self.structure_parameter_set.structure.scatterers()[ + self.index + ].u_iso + + def _setuiso(self, dummy, value): + self.structure_parameter_set.structure.scatterers()[ + self.index + ].u_iso = value + return + + def _getelem(self): + return self.structure.element_symbol() + + element = property(_getelem) + + +# End class CCTBXScattererParSet + + +class CCTBXUnitCellParSet(ParameterSet): + """Adapt a cctbx unit_cell to the ParameterSet interface. + + Attributes + ---------- + name : str + The name of this ParameterSet, always "unitcell". + a, b, c, alpha, beta, gamma : ParameterAdapter + The unit cell parameters. + """ + + def __init__(self, structure_parameter_set): + """Initialize the unit cell ParameterSet. + + Parameters + ---------- + structure_parameter_set : CCTBXCrystalParSet + The CCTBXCrystalParSet that contains the cctbx structure + and the unit cell being wrapped. + """ + ParameterSet.__init__(self, "unitcell") + self.structure_parameter_set = structure_parameter_set + self._lattice_parameters = list( + self.structure_parameter_set.structure.unit_cell().parameters() + ) + + self.add_parameter( + ParameterAdapter("a", None, self._latgetter(0), self._latsetter(0)) + ) + self.add_parameter( + ParameterAdapter("b", None, self._latgetter(1), self._latsetter(1)) + ) + self.add_parameter( + ParameterAdapter("c", None, self._latgetter(2), self._latsetter(2)) + ) + self.add_parameter( + ParameterAdapter( + "alpha", None, self._latgetter(3), self._latsetter(3) + ) + ) + self.add_parameter( + ParameterAdapter( + "beta", None, self._latgetter(4), self._latsetter(4) + ) + ) + self.add_parameter( + ParameterAdapter( + "gamma", None, self._latgetter(5), self._latsetter(5) + ) + ) + + return + + def _latgetter(self, i): + + def f(dummy): + return self._lattice_parameters[i] + + return f + + def _latsetter(self, i): + + def f(dummy, value): + self._lattice_parameters[i] = value + self.structure_parameter_set._update = True + return + + return f + + +# End class CCTBXUnitCellParSet + +# FIXME - Special positions should be constant. + + +class CCTBXCrystalParSet(BaseStructureParSet): + """Adapt a cctbx structure to the ParameterSet interface. + + Attributes + ---------- + structure : cctbx.crystal.special_position_settings + The adapted cctbx structure object. + scatterers : list of CCTBXScattererParSet + The scatterer ParameterSets. + unitcell : CCTBXUnitCellParSet + The unit cell ParameterSet for the structure. + """ + + def __init__(self, name, structure): + """Initialize the crystal ParameterSet. + + Parameters + ---------- + name : str + The name of this ParameterSet. + structure : cctbx.crystal.special_position_settings + The cctbx structure to adapt. + """ + ParameterSet.__init__(self, name) + self.structure = structure + self.add_parameter_set(CCTBXUnitCellParSet(self)) + self.scatterers = [] + + self._update = False + + cdict = {} + for s in structure.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.cmistructure.sgconstraints import _constrain_space_group + + symbol = self.get_space_group() + _constrain_space_group(self, symbol) + + return + + def update(self): + """Rebuild the unit cell after a change in lattice parameters. + + Call this function before using the CCTBXCrystalParSet. The unit + cell is only remade if a lattice parameter has changed. + """ + if not self._update: + return + + self._update = False + structure = self.structure + sgn = structure.space_group().match_tabulated_settings().number() + + # Create the symmetry object + from cctbx.crystal import symmetry + + symm = symmetry( + unit_cell=self.unitcell._lattice_parameters, space_group_symbol=sgn + ) + + # Now the new structure + newstru = structure.__class__( + crystal_symmetry=symm, scatterers=structure.scatterers() + ) + + self.unitcell._lattice_parameters = list( + newstru.unit_cell().parameters() + ) + + self.structure = newstru + return + + @classmethod + def can_adapt(self, structure): + """Return whether the structure can be adapted by this class. + + Parameters + ---------- + structure : object + The structure object to check. + + Returns + ------- + bool + The flag indicating if `structure` is a + cctbx.crystal.special_position_settings. False if cctbx is + not installed. + """ + try: + from cctbx.crystal import special_position_settings + except ImportError: + return False + return isinstance(structure, special_position_settings) + + def get_lattice(self): + """Return the ParameterSet containing the lattice Parameters. + + Returns + ------- + CCTBXUnitCellParSet + The unit cell ParameterSet of the structure. + """ + return self.unitcell + + def get_scatterers(self): + """Return the list of ParameterSets that represent the + scatterers. + + Returns + ------- + list of CCTBXScattererParSet + The scatterer ParameterSets of the structure. + """ + return self.scatterers + + def get_space_group(self): + """Return the Hermann-Mauguin space group symbol of the + structure. + + Returns + ------- + str + The Hermann-Mauguin space group symbol. + """ + space_group = self.structure.space_group() + t = space_group.type() + return t.lookup_symbol() + + +# End class CCTBXCrystalParSet diff --git a/src/diffpy/cmistructure/cmistructure_app.py b/src/diffpy/cmistructure/cmistructure_app.py index 3bea673..9d2d518 100644 --- a/src/diffpy/cmistructure/cmistructure_app.py +++ b/src/diffpy/cmistructure/cmistructure_app.py @@ -7,7 +7,8 @@ def main(): parser = argparse.ArgumentParser( prog="diffpy.cmistructure", description=( - "diffpy.cmi package for doing refinements with structure objects\n\n" + "diffpy.cmi package for doing refinements with " + "structure objects\n\n" "For more information, visit: " "https://github.com/diffpy/diffpy.cmistructure/" ), diff --git a/src/diffpy/cmistructure/diffpyparset.py b/src/diffpy/cmistructure/diffpyparset.py new file mode 100644 index 0000000..36a318d --- /dev/null +++ b/src/diffpy/cmistructure/diffpyparset.py @@ -0,0 +1,349 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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.cmistructure.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): + """Adapt a diffpy.structure.Atom to the ParameterSet interface. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. + See that class for base attributes. + + Attributes + ---------- + atom : diffpy.structure.Atom + The atom this is adapting. + element : str + The element name (property). + x, y, z : ParameterAdapter + The fractional coordinates of the atom. + occupancy : ParameterAdapter + The occupancy of the atom on its crystal location. + occ : ParameterProxy + The proxy for `occupancy`. + Uij : ParameterAdapter or ParameterProxy + The anisotropic displacement factors U11, U22, U33, U12, U21, U13, + U31, U23 and U32 of the atom. The Uij and Uji parameters are the + same. + Uiso : ParameterAdapter + The isotropic displacement factor of the atom. + Bij : ParameterAdapter or ParameterProxy + The anisotropic displacement factors B11, B22, B33, B12, B21, B13, + B31, B23 and B32 of the atom, with Bij = 8*pi**2*Uij. The Bij and + Bji parameters are the same. + Biso : ParameterAdapter + The isotropic displacement factor of the atom, as a B-factor. + """ + + def __init__(self, name, atom): + """Initialize the atom ParameterSet. + + Parameters + ---------- + name : str + The name of this ParameterSet. + atom : diffpy.structure.Atom + The atom to adapt. + """ + ParameterSet.__init__(self, name) + self.atom = atom + a = atom + # x, y, z, occupancy + self.add_parameter( + ParameterAdapter("x", a, _xyzgetter(0), _xyzsetter(0)) + ) + self.add_parameter( + ParameterAdapter("y", a, _xyzgetter(1), _xyzsetter(1)) + ) + self.add_parameter( + ParameterAdapter("z", a, _xyzgetter(2), _xyzsetter(2)) + ) + occupancy = ParameterAdapter("occupancy", a, attr="occupancy") + self.add_parameter(occupancy) + self.add_parameter(ParameterProxy("occ", occupancy)) + # U + self.add_parameter(ParameterAdapter("U11", a, attr="U11")) + self.add_parameter(ParameterAdapter("U22", a, attr="U22")) + self.add_parameter(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.add_parameter(U12) + self.add_parameter(U21) + self.add_parameter(U13) + self.add_parameter(U31) + self.add_parameter(U23) + self.add_parameter(U32) + self.add_parameter(ParameterAdapter("Uiso", a, attr="Uisoequiv")) + # B + self.add_parameter(ParameterAdapter("B11", a, attr="B11")) + self.add_parameter(ParameterAdapter("B22", a, attr="B22")) + self.add_parameter(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.add_parameter(B12) + self.add_parameter(B21) + self.add_parameter(B13) + self.add_parameter(B31) + self.add_parameter(B23) + self.add_parameter(B32) + self.add_parameter(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(parameter): + return bind2nd(getattr, parameter) + + +def _latsetter(parameter): + return bind2nd(setattr, parameter) + + +class DiffpyLatticeParSet(ParameterSet): + """Adapt a diffpy.structure.Lattice to the ParameterSet interface. + + This class derives from diffpy.srfit.fitbase.parameterset.ParameterSet. + See that class for base attributes. + + Attributes + ---------- + lattice : diffpy.structure.Lattice + The lattice this is adapting. + name : str + The name of this ParameterSet, always "lattice". + angle_units : str + The units of the lattice angles, always "deg". + a, b, c, alpha, beta, gamma : ParameterAdapter + The unit cell parameters. + """ + + def __init__(self, lattice): + """Initialize the lattice ParameterSet. + + Parameters + ---------- + lattice : diffpy.structure.Lattice + The lattice to adapt. + """ + ParameterSet.__init__(self, "lattice") + self.angle_units = "deg" + self.lattice = lattice + lat = lattice + self.add_parameter( + ParameterAdapter("a", lat, _latgetter("a"), _latsetter("a")) + ) + self.add_parameter( + ParameterAdapter("b", lat, _latgetter("b"), _latsetter("b")) + ) + self.add_parameter( + ParameterAdapter("c", lat, _latgetter("c"), _latsetter("c")) + ) + self.add_parameter( + ParameterAdapter( + "alpha", lat, _latgetter("alpha"), _latsetter("alpha") + ) + ) + self.add_parameter( + ParameterAdapter( + "beta", lat, _latgetter("beta"), _latsetter("beta") + ) + ) + self.add_parameter( + ParameterAdapter( + "gamma", lat, _latgetter("gamma"), _latsetter("gamma") + ) + ) + return + + def __repr__(self): + return repr(self.lattice) + + +# End class DiffpyLatticeParSet + + +class DiffpyStructureParSet(SrRealParSet): + """Adapt a diffpy.structure.Structure to the ParameterSet interface. + + This class derives from SrRealParSet. See that class for base + attributes. + + Attributes + ---------- + atoms : list of DiffpyAtomParSet + The atom ParameterSets, provided for convenience. + structure : diffpy.structure.Structure + The structure this is adapting. + lattice : DiffpyLatticeParSet + The managed lattice ParameterSet. + : DiffpyAtomParSet + The managed atom ParameterSets. is the atomic element and + is the index of that element in the structure, starting + from zero. For nickel in P1 symmetry, the managed + DiffpyAtomParSets are named "Ni0", "Ni1", "Ni2" and "Ni3". + """ + + def __init__(self, name, structure): + """Initialize the structure ParameterSet. + + Parameters + ---------- + name : str + The name of the structure. + structure : diffpy.structure.Structure + The structure to adapt. + """ + SrRealParSet.__init__(self, name) + self.structure = structure + self.add_parameter_set(DiffpyLatticeParSet(structure.lattice)) + self.atoms = [] + + cdict = {} + for a in structure: + 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.structure) + + def get_lattice(self): + """Return the ParameterSet containing the lattice Parameters. + + Returns + ------- + DiffpyLatticeParSet + The lattice ParameterSet of the structure. + """ + return self.lattice + + @classmethod + def can_adapt(self, structure): + """Return whether the structure can be adapted by this class. + + Parameters + ---------- + structure : object + The structure object to check. + + Returns + ------- + bool + The flag indicating if `structure` is a diffpy.structure.Structure. + """ + from diffpy.structure import Structure + + return isinstance(structure, Structure) + + def get_scatterers(self): + """Return the list of ParameterSets that represent the + scatterers. + + Returns + ------- + list of DiffpyAtomParSet + The atom ParameterSets of the 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 + + structure = SrRealParSet._get_srreal_structure(self) + return nometa(structure) + + +# End class DiffpyStructureParSet diff --git a/src/diffpy/cmistructure/functions.py b/src/diffpy/cmistructure/functions.py deleted file mode 100644 index e7e2c8e..0000000 --- a/src/diffpy/cmistructure/functions.py +++ /dev/null @@ -1,31 +0,0 @@ -import numpy as np - - -def dot_product(a, b): - """Compute the dot product of two vectors of any size. - - Ensure that the inputs, a and b, are of the same size. - The supported types are "array_like" objects, which can - be converted to a NumPy array. Examples include lists and tuples. - - Parameters - ---------- - a : array_like - The first input vector. - b : array_like - The second input vector. - - Returns - ------- - float - The dot product of the two vectors. - - Examples - -------- - Compute the dot product of two lists: - >>> a = [1, 2, 3] - >>> b = [4, 5, 6] - >>> dot_product(a, b) - 32.0 - """ - return float(np.dot(a, b)) diff --git a/src/diffpy/cmistructure/objcrystparset.py b/src/diffpy/cmistructure/objcrystparset.py new file mode 100644 index 0000000..ef3ed1b --- /dev/null +++ b/src/diffpy/cmistructure/objcrystparset.py @@ -0,0 +1,1996 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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.cmistructure.srrealparset import SrRealParSet +from diffpy.srfit.fitbase.parameter import ( + Parameter, + ParameterAdapter, + ParameterProxy, +) +from diffpy.srfit.fitbase.parameterset import ParameterSet + + +class ObjCrystScattererParSet(ParameterSet): + """Base adapter for a pyobjcryst 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 + ---------- + scatterer : pyobjcryst.scatterer.Scatterer + The adapted pyobjcryst object. + parent : ParameterSet or None + The ParameterSet this belongs to. + x, y, z : ParameterAdapter + The position of the scatterer in crystal coordinates. + occ : ParameterAdapter + The occupancy of the scatterer on its crystal site. + """ + + def __init__(self, name, scatterer, parent): + """Initialize the scatterer ParameterSet. + + Parameters + ---------- + name : str + The name of the scatterer. + scatterer : pyobjcryst.scatterer.Scatterer + The pyobjcryst scatterer to adapt. + parent : ParameterSet or None + The ParameterSet this belongs to. + """ + ParameterSet.__init__(self, name) + self.scatterer = scatterer + self.parent = parent + + # x, y, z, occ + self.add_parameter(ParameterAdapter("x", self.scatterer, attr="X")) + self.add_parameter(ParameterAdapter("y", self.scatterer, attr="Y")) + self.add_parameter(ParameterAdapter("z", self.scatterer, attr="Z")) + self.add_parameter( + ParameterAdapter("occ", self.scatterer, attr="Occupancy") + ) + return + + def is_dummy(self): + """Return whether this scatterer is a dummy atom. + + Returns + ------- + bool + The flag indicating if this is a dummy atom. Always False for + this class. + """ + return False + + def has_scatterers(self): + """Return whether this scatterer has its own scatterers. + + Returns + ------- + bool + The flag indicating if this scatterer has a ``get_scatterers`` + method. + """ + return hasattr(self, "get_scatterers") + + +# End class ObjCrystScattererParSet + + +class ObjCrystAtomParSet(ObjCrystScattererParSet): + """Adapt a pyobjcryst.atom.Atom to the ParameterSet interface. + + This class derives from ObjCrystScattererParSet. + + Attributes + ---------- + scatterer : pyobjcryst.atom.Atom + The adapted atom. + element : str + The non-refinable name of the element (property). + parent : ObjCrystCrystalParSet + The crystal ParameterSet this belongs to. + occ : ParameterAdapter + The occupancy of the atom on its crystal location. + Biso : ParameterAdapter + The isotropic displacement factor of the atom. + Bij : ParameterAdapter or ParameterProxy + The anisotropic displacement factors B11, B22, B33, B12, B21, B13, + B31, B23 and B32 of the atom. The Bij and Bji parameters are the + same. + """ + + def __init__(self, name, atom, parent): + """Initialize the atom ParameterSet. + + Parameters + ---------- + name : str + The name of the atom. + atom : pyobjcryst.atom.Atom + The atom to adapt. + parent : ObjCrystCrystalParSet + The crystal ParameterSet this belongs to. + """ + ObjCrystScattererParSet.__init__(self, name, atom, parent) + sp = atom.GetScatteringPower() + + # The B-parameters + self.add_parameter(ParameterAdapter("Biso", sp, attr="Biso")) + self.add_parameter(ParameterAdapter("B11", sp, attr="B11")) + self.add_parameter(ParameterAdapter("B22", sp, attr="B22")) + self.add_parameter(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.add_parameter(B12) + self.add_parameter(B21) + self.add_parameter(B13) + self.add_parameter(B31) + self.add_parameter(B23) + self.add_parameter(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.scatterer.GetScatteringPower().GetSymbol() + + element = property(_getelem) + + +# End class ObjCrystAtomParSet + + +class ObjCrystMoleculeParSet(ObjCrystScattererParSet): + """Adapt a pyobjcryst.molecule.Molecule to the ParameterSet + interface. + + This class derives from ObjCrystScattererParSet. Other attributes are + inherited from diffpy.srfit.fitbase.parameterset.ParameterSet. + + Attributes + ---------- + scatterer : pyobjcryst.molecule.Molecule + The adapted molecule. + structure : pyobjcryst.molecule.Molecule + The adapted molecule. + parent : ObjCrystCrystalParSet or None + The crystal ParameterSet this belongs to. This is None when the + ObjCrystMoleculeParSet is used on its own. + atoms : list of ObjCrystMolAtomParSet + The ParameterSets of the atoms in the molecule. + occ : ParameterAdapter + The occupancy of the molecule on its crystal location. + q0, q1, q2, q3 : ParameterAdapter + The orientational quaternion of the molecule. + """ + + def __init__(self, name, molecule, parent=None): + """Initialize the molecule ParameterSet. + + Parameters + ---------- + name : str + The name of the molecule. + molecule : pyobjcryst.molecule.Molecule + The molecule to adapt. + parent : ObjCrystCrystalParSet, optional + The crystal ParameterSet this belongs to (default None). + + Raises + ------ + AttributeError + If a MolAtom in the molecule has no name, or if two MolAtoms + share a name. Give every MolAtom a unique name before + wrapping the molecule. + """ + ObjCrystScattererParSet.__init__(self, name, molecule, parent) + self.structure = molecule + + # Add orientation quaternion + self.add_parameter(ParameterAdapter("q0", self.scatterer, attr="Q0")) + self.add_parameter(ParameterAdapter("q1", self.scatterer, attr="Q1")) + self.add_parameter(ParameterAdapter("q2", self.scatterer, attr="Q2")) + self.add_parameter(ParameterAdapter("q3", self.scatterer, 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, structure): + """Return whether the structure can be adapted by this class. + + Parameters + ---------- + structure : object + The structure object to check. + + Returns + ------- + bool + The flag indicating if `structure` is a pyobjcryst Molecule. + """ + from pyobjcryst.molecule import Molecule + + return isinstance(structure, Molecule) + + # Part of SrRealParSet interface + def use_symmetry(self, use=True): + """Set whether this structure uses symmetry. + + This structure object does not support symmetry, so this does + nothing. + + Parameters + ---------- + use : bool, optional + The flag indicating if symmetry is used (default True). + """ + return + + # Part of SrRealParSet interface + def using_symmetry(self): + """Return whether symmetry is being used. + + Returns + ------- + bool + The flag indicating if symmetry is used. Always False, since + 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.structure + + def get_lattice(self): + """Return a ParameterSet holding a unit cubic lattice. + + A molecule is not periodic, so this returns a new ParameterSet + with a = b = c = 1 and alpha = beta = gamma = 90 degrees. + + Returns + ------- + ParameterSet + The ParameterSet holding the placeholder lattice Parameters. + """ + lattice = ParameterSet("lattice") + lattice.new_parameter("a", 1.0) + lattice.new_parameter("b", 1.0) + lattice.new_parameter("c", 1.0) + lattice.new_parameter("alpha", 90) + lattice.new_parameter("beta", 90) + lattice.new_parameter("gamma", 90) + lattice.angle_units = "deg" + return lattice + + def get_scatterers(self): + """Return the list of ParameterSets that represent the + scatterers. + + Returns + ------- + list of ObjCrystMolAtomParSet + The atom ParameterSets of the molecule. + """ + return self.atoms + + def wrap_restraints(self): + """Wrap the restraints implicit to the molecule. + + This wraps the MolBonds, MolBondAngles and MolDihedralAngles of + the Molecule as ObjCrystMoleculeRestraint objects. Restraints + wrapped this way cannot be modified from within this class. + """ + # Wrap restraints. Restraints wrapped in this way cannot be modified + # from within this class. + for b in self.scatterer.GetBondList(): + restraint = ObjCrystMoleculeRestraint(b) + self._restraints.add(restraint) + + for ba in self.scatterer.GetBondAngleList(): + restraint = ObjCrystMoleculeRestraint(ba) + self._restraints.add(restraint) + + for da in self.scatterer.GetDihedralAngleList(): + restraint = ObjCrystMoleculeRestraint(da) + self._restraints.add(restraint) + + return + + def wrap_stretch_mode_parameters(self): + """Wrap the stretch modes implicit to the Molecule as + Parameters. + + This wraps the StretchModeBondLengths and StretchModeBondAngles + of the Molecule as Parameters. The MolBondAtoms in the Molecule + must have unique names. Torsion angles are not wrapped, as there + is not enough information to determine each MolAtom in the + angle. + + Each Parameter is named after its constituent atoms, as + "bl_aname1_aname2" for bond lengths and + "ba_aname1_aname2_aname3" for bond angles. + """ + for mode in self.scatterer.GetStretchModeBondLengthList(): + name1 = mode.mpAtom0.GetName() + name2 = mode.mpAtom1.GetName() + + name = "bl_" + "_".join((name1, name2)) + + atom1 = getattr(self, name1) + atom2 = getattr(self, name2) + + parameter = ObjCrystBondLengthParameter( + name, atom1, atom2, mode=mode + ) + + atoms = [] + for a in mode.GetAtoms(): + name = a.GetName() + atoms.append(getattr(self, name)) + + parameter.AddAtoms(atoms) + + self.add_parameter(parameter) + + for mode in self.scatterer.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) + + parameter = ObjCrystBondAngleParameter( + name, atom1, atom2, atom3, mode=mode + ) + + atoms = [] + for a in mode.GetAtoms(): + name = a.GetName() + atoms.append(getattr(self, name)) + parameter.AddAtoms(atoms) + + self.add_parameter(parameter) + + return + + def restrain_bond_length( + self, atom1, atom2, length, sigma, delta, scaled=False + ): + """Add a bond length restraint between two atoms. + + This creates an ObjCrystBondLengthRestraint and adds it to the + ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second atom in the bond. + length : float + The length of the bond in Angstroms. + sigma : float + The uncertainty of the bond length in Angstroms. + delta : float + The width of the bond in Angstroms. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondLengthRestraint + The restraint, for use with the ``unrestrain`` method. + """ + restraint = ObjCrystBondLengthRestraint( + atom1, atom2, length, sigma, delta, scaled + ) + self._restraints.add(restraint) + + return restraint + + def restrain_bond_length_parameter( + self, parameter, length, sigma, delta, scaled=False + ): + """Add a bond length restraint on a bond length Parameter. + + This creates an ObjCrystBondLengthRestraint between the atoms of + `parameter` and adds it to the ObjCrystMoleculeParSet. + + Parameters + ---------- + parameter : ObjCrystBondLengthParameter + The bond length Parameter to restrain (see + add_bond_length_parameter). + length : float + The length of the bond in Angstroms. + sigma : float + The uncertainty of the bond length in Angstroms. + delta : float + The width of the bond in Angstroms. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondLengthRestraint + The restraint, for use with the ``unrestrain`` method. + """ + return self.restrain_bond_length( + parameter.atom1, parameter.atom2, length, sigma, delta, scaled + ) + + def restrain_bond_angle( + self, atom1, atom2, atom3, angle, sigma, delta, scaled=False + ): + """Add a bond angle restraint between three atoms. + + This creates an ObjCrystBondAngleRestraint and adds it to the + ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the bond angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the bond angle. + atom3 : ObjCrystMolAtomParSet + The third atom in the bond angle. + angle : float + The bond angle in radians. + sigma : float + The uncertainty of the bond angle in radians. + delta : float + The width of the bond angle in radians. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondAngleRestraint + The restraint, for use with the ``unrestrain`` method. + """ + restraint = ObjCrystBondAngleRestraint( + atom1, atom2, atom3, angle, sigma, delta, scaled + ) + self._restraints.add(restraint) + + return restraint + + def restrain_bond_angle_parameter( + self, parameter, angle, sigma, delta, scaled=False + ): + """Add a bond angle restraint on a bond angle Parameter. + + This creates an ObjCrystBondAngleRestraint between the atoms of + `parameter` and adds it to the ObjCrystMoleculeParSet. + + Parameters + ---------- + parameter : ObjCrystBondAngleParameter + The bond angle Parameter to restrain (see + add_bond_angle_parameter). + angle : float + The bond angle in radians. + sigma : float + The uncertainty of the bond angle in radians. + delta : float + The width of the bond angle in radians. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystBondAngleRestraint + The restraint, for use with the ``unrestrain`` method. + """ + return self.restrain_bond_angle( + parameter.atom1, + parameter.atom2, + parameter.atom3, + angle, + sigma, + delta, + scaled, + ) + + def restrain_dihedral_angle( + self, atom1, atom2, atom3, atom4, angle, sigma, delta, scaled=False + ): + """Add a dihedral angle restraint between four atoms. + + This creates an ObjCrystDihedralAngleRestraint and adds it to the + ObjCrystMoleculeParSet. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the angle. + atom4 : ObjCrystMolAtomParSet + The fourth atom in the angle. + angle : float + The dihedral angle in radians. + sigma : float + The uncertainty of the dihedral angle in radians. + delta : float + The width of the dihedral angle in radians. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystDihedralAngleRestraint + The restraint, for use with the ``unrestrain`` method. + """ + restraint = ObjCrystDihedralAngleRestraint( + atom1, atom2, atom3, atom4, angle, sigma, delta, scaled + ) + self._restraints.add(restraint) + + return restraint + + def restrain_dihedral_angle_parameter( + self, parameter, angle, sigma, delta, scaled=False + ): + """Add a dihedral angle restraint on a dihedral angle Parameter. + + This creates an ObjCrystDihedralAngleRestraint between the atoms of + `parameter` and adds it to the ObjCrystMoleculeParSet. + + Parameters + ---------- + parameter : ObjCrystDihedralAngleParameter + The dihedral angle Parameter to restrain (see + add_dihedral_angle_parameter). + angle : float + The dihedral angle in radians. + sigma : float + The uncertainty of the dihedral angle in radians. + delta : float + The width of the dihedral angle in radians. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + ObjCrystDihedralAngleRestraint + The restraint, for use with the ``unrestrain`` method. + """ + return self.restrain_dihedral_angle( + parameter.atom1, + parameter.atom2, + parameter.atom3, + parameter.atom4, + angle, + sigma, + delta, + scaled, + ) + + def add_bond_length_parameter( + self, name, atom1, atom2, value=None, const=False + ): + """Add a refinable bond length to the Molecule. + + This adds an ObjCrystBondLengthParameter to the + ObjCrystMoleculeParSet that can be adjusted during the fit. + + Parameters + ---------- + name : str + The name of the new Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second (mutated) atom in the bond. + value : float, optional + The initial bond length. If None (default), the current + distance between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + + Returns + ------- + ObjCrystBondLengthParameter + The new bond length Parameter. + """ + parameter = ObjCrystBondLengthParameter( + name, atom1, atom2, value, const + ) + self.add_parameter(parameter) + + return parameter + + def add_bond_angle_parameter( + self, name, atom1, atom2, atom3, value=None, const=False + ): + """Add a refinable bond angle to the Molecule. + + This adds an ObjCrystBondAngleParameter to the + ObjCrystMoleculeParSet that can be adjusted during the fit. + + Parameters + ---------- + name : str + The name of the new Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the bond angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the bond angle. + atom3 : ObjCrystMolAtomParSet + The third (mutated) atom in the bond angle. + value : float, optional + The initial bond angle in radians. If None (default), the + current bond angle between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + + Returns + ------- + ObjCrystBondAngleParameter + The new bond angle Parameter. + """ + parameter = ObjCrystBondAngleParameter( + name, atom1, atom2, atom3, value, const + ) + self.add_parameter(parameter) + + return parameter + + def add_dihedral_angle_parameter( + self, name, atom1, atom2, atom3, atom4, value=None, const=False + ): + """Add a refinable dihedral angle to the Molecule. + + This adds an ObjCrystDihedralAngleParameter to the + ObjCrystMoleculeParSet that can be adjusted during the fit. + + Parameters + ---------- + name : str + The name of the new Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the dihedral angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the dihedral angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the dihedral angle. + atom4 : ObjCrystMolAtomParSet + The fourth (mutated) atom in the dihedral angle. + value : float, optional + The initial dihedral angle in radians. If None (default), the + current dihedral angle between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + + Returns + ------- + ObjCrystDihedralAngleParameter + The new dihedral angle Parameter. + """ + parameter = ObjCrystDihedralAngleParameter( + name, atom1, atom2, atom3, atom4, value, const + ) + self.add_parameter(parameter) + + return parameter + + +# End class ObjCrystMoleculeParSet + + +class ObjCrystMolAtomParSet(ObjCrystScattererParSet): + """Adapt a pyobjcryst.molecule.MolAtom to the ParameterSet + interface. + + This class derives from ObjCrystScattererParSet. MolAtom does not + derive from Scatterer, but the relevant interface is the same within + pyobjcryst. See the ParameterSet class for base attributes. + + Attributes + ---------- + scatterer : pyobjcryst.molecule.MolAtom + The adapted MolAtom. + parent : ObjCrystMoleculeParSet + The molecule ParameterSet this belongs to. + element : str + The non-refinable name of the element, or "dummy" for a dummy + atom (property). + occ : ParameterAdapter + The occupancy of the atom on its crystal location. + Biso : ParameterAdapter + The isotropic displacement factor of the atom. This does not + exist for dummy atoms; see the ``is_dummy`` method. + Bij : ParameterAdapter or ParameterProxy + The anisotropic displacement factors B11, B22, B33, B12, B21, B13, + B31, B23 and B32 of the atom. The Bij and Bji parameters are the + same. These do not exist for dummy atoms. + """ + + def __init__(self, name, scatterer, parent): + """Initialize the MolAtom ParameterSet. + + Parameters + ---------- + name : str + The name of the atom. + scatterer : pyobjcryst.molecule.MolAtom + The MolAtom to adapt. + parent : ObjCrystMoleculeParSet + The molecule ParameterSet this belongs to. + """ + ObjCrystScattererParSet.__init__(self, name, scatterer, parent) + sp = scatterer.GetScatteringPower() + + # Only wrap this if there is a scattering power + if sp is not None: + self.add_parameter(ParameterAdapter("Biso", sp, attr="Biso")) + self.add_parameter(ParameterAdapter("B11", sp, attr="B11")) + self.add_parameter(ParameterAdapter("B22", sp, attr="B22")) + self.add_parameter(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.add_parameter(B12) + self.add_parameter(B21) + self.add_parameter(B13) + self.add_parameter(B31) + self.add_parameter(B23) + self.add_parameter(B32) + + return + + def _getelem(self): + """Getter for the element type.""" + sp = self.scatterer.GetScatteringPower() + if sp: + return sp.GetSymbol() + else: + return "dummy" + + element = property(_getelem) + + def is_dummy(self): + """Return whether this atom is a dummy atom. + + Returns + ------- + bool + The flag indicating if this is a dummy atom. + """ + return self.scatterer.IsDummy() + + +# End class ObjCrystMolAtomParSet + + +class ObjCrystMoleculeRestraint(object): + """Base class for adapting pyobjcryst Molecule restraints to srfit. + + This implements the ``penalty`` method of + diffpy.srfit.fitbase.restraint.Restraint by calling + ``GetLogLikelihood`` of the pyobjcryst restraint. The ``restrain`` + method is not needed or implemented. + + Attributes + ---------- + restraint : object + The pyobjcryst Molecule restraint. + scaled : bool + The flag indicating if the restraint is scaled (multiplied) by + the unrestrained point-average chi^2 (chi^2/numpoints) (default + False). + """ + + def __init__(self, restraint, scaled=False): + """Wrap a pyobjcryst Molecule restraint as a Restraint. + + Parameters + ---------- + restraint : object + The pyobjcryst Molecule restraint. + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + """ + self.restraint = restraint + self.scaled = scaled + return + + def penalty(self, w=1.0): + """Calculate the penalty of the restraint. + + Parameters + ---------- + w : float, optional + The point-average chi^2 which is optionally used to scale the + penalty (default 1.0). + + Returns + ------- + float + The log-likelihood of the pyobjcryst restraint, optionally + scaled by `w`. + """ + penalty = self.restraint.GetLogLikelihood() + if self.scaled: + penalty *= w + return penalty + + +# End class ObjCrystMoleculeRestraint + + +class ObjCrystBondLengthRestraint(ObjCrystMoleculeRestraint): + """Restrain the distance between two atoms. + + Attributes + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second atom in the bond. + length : float + The length of the bond in Angstroms. + sigma : float + The uncertainty of the bond length in Angstroms. + delta : float + The width of the bond in Angstroms. + restraint : pyobjcryst.molecule.MolBond + The pyobjcryst bond length restraint. + scaled : bool + The 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): + """Initialize the bond length restraint. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second atom in the bond. + length : float + The length of the bond in Angstroms. + sigma : float + The uncertainty of the bond length in Angstroms. + delta : float + The width of the bond in Angstroms. + scaled : bool, optional + The 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.scatterer.GetMolecule() + restraint = m.AddBond( + atom1.scatterer, atom2.scatterer, length, sigma, delta + ) + + ObjCrystMoleculeRestraint.__init__(self, restraint, scaled) + return + + # Give access to the parameters of the restraint + length = property( + lambda self: self.restraint.GetLength0(), + lambda self, value: self.restraint.SetLength0(value), + ) + sigma = property( + lambda self: self.restraint.GetLengthSigma(), + lambda self, value: self.restraint.SetLengthSigma(value), + ) + delta = property( + lambda self: self.restraint.GetLengthDelta(), + lambda self, value: self.restraint.SetLengthDelta(value), + ) + + +# End class ObjCrystBondLengthRestraint + + +class ObjCrystBondAngleRestraint(ObjCrystMoleculeRestraint): + """Restrain the angle defined by three atoms. + + Attributes + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the angle. + atom3 : ObjCrystMolAtomParSet + The third atom in the angle. + angle : float + The bond angle in radians. + sigma : float + The uncertainty of the bond angle in radians. + delta : float + The width of the bond angle in radians. + restraint : pyobjcryst.molecule.MolBondAngle + The pyobjcryst bond angle restraint. + scaled : bool + The 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): + """Initialize the bond angle restraint. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the bond angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the bond angle. + atom3 : ObjCrystMolAtomParSet + The third atom in the bond angle. + angle : float + The bond angle in radians. + sigma : float + The uncertainty of the bond angle in radians. + delta : float + The width of the bond angle in radians. + scaled : bool, optional + The 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.scatterer.GetMolecule() + restraint = m.AddBondAngle( + atom1.scatterer, + atom2.scatterer, + atom3.scatterer, + angle, + sigma, + delta, + ) + + ObjCrystMoleculeRestraint.__init__(self, restraint, scaled) + return + + # Give access to the parameters of the restraint + angle = property( + lambda self: self.restraint.GetAngle0(), + lambda self, value: self.restraint.SetAngle0(value), + ) + sigma = property( + lambda self: self.restraint.GetAngleSigma(), + lambda self, value: self.restraint.SetAngleSigma(value), + ) + delta = property( + lambda self: self.restraint.GetAngleDelta(), + lambda self, value: self.restraint.SetAngleDelta(value), + ) + + +# End class ObjCrystBondAngleRestraint + + +class ObjCrystDihedralAngleRestraint(ObjCrystMoleculeRestraint): + """Restrain the dihedral (torsion) angle defined by four atoms. + + Attributes + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the angle. + atom4 : ObjCrystMolAtomParSet + The fourth atom in the angle. + angle : float + The dihedral angle in radians. + sigma : float + The uncertainty of the dihedral angle in radians. + delta : float + The width of the dihedral angle in radians. + restraint : pyobjcryst.molecule.MolDihedralAngle + The pyobjcryst dihedral angle restraint. + scaled : bool + The 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 + ): + """Initialize the dihedral angle restraint. + + Parameters + ---------- + atom1 : ObjCrystMolAtomParSet + The first atom in the angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the angle. + atom4 : ObjCrystMolAtomParSet + The fourth atom in the angle. + angle : float + The dihedral angle in radians. + sigma : float + The uncertainty of the dihedral angle in radians. + delta : float + The width of the dihedral angle in radians. + scaled : bool, optional + The 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.scatterer.GetMolecule() + restraint = m.AddDihedralAngle( + atom1.scatterer, + atom2.scatterer, + atom3.scatterer, + atom4.scatterer, + angle, + sigma, + delta, + ) + + ObjCrystMoleculeRestraint.__init__(self, restraint, scaled) + return + + # Give access to the parameters of the restraint + angle = property( + lambda self: self.restraint.GetAngle0(), + lambda self, value: self.restraint.SetAngle0(value), + ) + sigma = property( + lambda self: self.restraint.GetAngleSigma(), + lambda self, value: self.restraint.SetAngleSigma(value), + ) + delta = property( + lambda self: self.restraint.GetAngleDelta(), + lambda self, value: self.restraint.SetAngleDelta(value), + ) + + +# End class ObjCrystDihedralAngleRestraint + + +class StretchModeParameter(Parameter): + """Partial Parameter class encapsulating pyobjcryst stretch modes. + + This class relies upon attributes that subclasses must set before + calling ``StretchModeParameter.__init__``. Do not instantiate this + class directly. + + Attributes + ---------- + mutated_atoms : set of ObjCrystMolAtomParSet + The set of all mutated atoms. Set by the subclass. + molecule : ObjCrystMoleculeParSet + The molecule the atoms belong to. Set by the subclass. + mode : pyobjcryst.molecule.StretchMode + The stretch mode used to change atomic positions. Set by the + subclass. + keepcenter : bool + The flag indicating whether to keep the center of mass of the + molecule stationary within the crystal when changing the value + of the parameter (default True). + """ + + def __init__(self, name, value=None, const=False): + """Initialize the stretch mode Parameter. + + Parameters + ---------- + name : str + The name of this Parameter. It must be a valid attribute + identifier. + value : float, optional + The initial value of this Parameter (default None). + const : bool, optional + The flag indicating whether the Parameter is a constant + (default False). + + Raises + ------ + ValueError + If `name` is not a valid attribute identifier. + """ + Parameter.__init__(self, name, value, const) + self.keepcenter = True + + def set_value(self, value): + """Set the value of the Parameter by stretching the molecule. + + The stretch mode moves the mutated atoms by the change in value. + + Parameters + ---------- + value : float + The new value of the Parameter. + + Returns + ------- + StretchModeParameter + Return self so that mutators can be chained. + """ + curval = self.get_value() + value = float(value) + + if value == curval: + return self + + # The StretchMode expects the change in mutated value. + delta = value - curval + self.mode.Stretch(delta, self.keepcenter) + + # Let Parameter take care of the general details + Parameter.set_value(self, value) + + return self + + def add_atoms(self, atomlist): + """Associate additional atoms with the Parameter. + + The added atoms are mutated in exactly the same way as the + primary mutated atom. This is useful when a group of atoms should + move rigidly in response to a change in a bond property. + + Parameters + ---------- + atomlist : ObjCrystMolAtomParSet or list of ObjCrystMolAtomParSet + The atom or atoms to associate with the Parameter. + + Returns + ------- + StretchModeParameter + Return self so that mutators can be chained. + """ + if not hasattr(atomlist, "__iter__"): + atomlist = [atomlist] + # Record the added atoms in the Parameter + self.mutated_atoms.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.scatterer 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 observe this Parameter while this + Parameter also observes them. Observable does not allow both, so + the mutated Parameters are notified directly. + + Parameters + ---------- + other : tuple, optional + The objects that have already been notified (default empty). + """ + noneother = () + # Notify the atoms that have moved + for a in self.mutated_atoms: + 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): + """Represent a bond length in a Molecule as 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 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second (mutated) atom in the bond. + mutated_atoms : set of ObjCrystMolAtomParSet + The set of all mutated atoms. + molecule : ObjCrystMoleculeParSet + The molecule the atoms belong to. + mode : pyobjcryst.molecule.StretchModeBondLength + The stretch mode for the bond. + name : str + The name of this Parameter (inherited). + const : bool + The flag indicating whether this is considered a constant + (inherited). + value : float + The property for ``get_value`` and ``set_value`` (inherited). + constraint : callable or None + The callable that calculates the value of this Parameter. If + None, the Parameter is responsible for its own value + (inherited). + bounds : list of float + The lower and upper bounds on the Parameter, which some + optimizers use when the Parameter is varied (inherited). + """ + + def __init__(self, name, atom1, atom2, value=None, const=False, mode=None): + """Initialize the bond length Parameter. + + Parameters + ---------- + name : str + The name of the Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the bond. + atom2 : ObjCrystMolAtomParSet + The second (mutated) atom in the bond. + value : float, optional + The initial bond length. If None (default), the current + distance between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + mode : pyobjcryst.molecule.StretchModeBondLength, optional + The existing stretch mode to use. If None (default), a new + StretchModeBondLength is built. + """ + # Create the mode + self.mode = mode + if mode is None: + self.mode = StretchModeBondLength( + atom1.scatterer, atom2.scatterer, None + ) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom2.scatterer) + self.mutated_atoms = 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.scatterer, atom2.scatterer) + 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 : bool, optional + The flag indicating if the Parameter is constant (default + True). + value : float, optional + The value to set the Parameter to (default None). If this is + not None, the Parameter gets a new value, constant or + otherwise. + + Returns + ------- + StretchModeParameter + Return 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 get_value(self): + """Return the bond length, recalculating it if needed. + + The atoms underlying the bond may have moved, so the bond length + is recalculated whenever the cached value has been cleared. + + Returns + ------- + float + The bond length in Angstroms. + """ + if self._value is None: + value = GetBondLength(self.atom1.scatterer, self.atom2.scatterer) + Parameter.set_value(self, value) + + return self._value + + +# End class ObjCrystBondLengthParameter + + +class ObjCrystBondAngleParameter(StretchModeParameter): + """Represent a bond angle in a Molecule as 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 : ObjCrystMolAtomParSet + The first atom in the bond angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the bond angle. + atom3 : ObjCrystMolAtomParSet + The third (mutated) atom in the bond angle. + mutated_atoms : set of ObjCrystMolAtomParSet + The set of all mutated atoms. + molecule : ObjCrystMoleculeParSet + The molecule the atoms belong to. + mode : pyobjcryst.molecule.StretchModeBondAngle + The stretch mode for the bond angle. + name : str + The name of this Parameter (inherited). + const : bool + The flag indicating whether this is considered a constant + (inherited). + value : float + The property for ``get_value`` and ``set_value`` (inherited). + constraint : callable or None + The callable that calculates the value of this Parameter. If + None, the Parameter is responsible for its own value + (inherited). + bounds : list of float + The lower and upper bounds on the Parameter, which some + optimizers use when the Parameter is varied (inherited). + """ + + def __init__( + self, name, atom1, atom2, atom3, value=None, const=False, mode=None + ): + """Initialize the bond angle Parameter. + + Parameters + ---------- + name : str + The name of the Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the bond angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the bond angle. + atom3 : ObjCrystMolAtomParSet + The third (mutated) atom in the bond angle. + value : float, optional + The initial bond angle in radians. If None (default), the + current bond angle between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + mode : pyobjcryst.molecule.StretchModeBondAngle, optional + The existing stretch mode to use. If None (default), a new + StretchModeBondAngle is built. + """ + # Create the stretch mode + self.mode = mode + if mode is None: + self.mode = StretchModeBondAngle( + atom1.scatterer, atom2.scatterer, atom3.scatterer, None + ) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom3.scatterer) + self.mutated_atoms = 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.scatterer, atom2.scatterer, atom3.scatterer + ) + 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 : bool, optional + The flag indicating if the Parameter is constant (default + True). + value : float, optional + The value to set the Parameter to (default None). If this is + not None, the Parameter gets a new value, constant or + otherwise. + + Returns + ------- + StretchModeParameter + Return 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 get_value(self): + """Return the bond angle, recalculating it if needed. + + The atoms underlying the bond angle may have moved, so the angle + is recalculated whenever the cached value has been cleared. + + Returns + ------- + float + The bond angle in radians. + """ + if self._value is None: + value = GetBondAngle( + self.atom1.scatterer, + self.atom2.scatterer, + self.atom3.scatterer, + ) + Parameter.set_value(self, value) + + return self._value + + +# End class ObjCrystBondAngleParameter + + +class ObjCrystDihedralAngleParameter(StretchModeParameter): + """Represent a dihedral angle in a Molecule as 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 : ObjCrystMolAtomParSet + The first atom in the dihedral angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the dihedral angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the dihedral angle. + atom4 : ObjCrystMolAtomParSet + The fourth (mutated) atom in the dihedral angle. + mutated_atoms : set of ObjCrystMolAtomParSet + The set of all mutated atoms. + molecule : ObjCrystMoleculeParSet + The molecule the atoms belong to. + mode : pyobjcryst.molecule.StretchModeTorsion + The stretch mode for the dihedral angle. + name : str + The name of this Parameter (inherited). + const : bool + The flag indicating whether this is considered a constant + (inherited). + value : float + The property for ``get_value`` and ``set_value`` (inherited). + constraint : callable or None + The callable that calculates the value of this Parameter. If + None, the Parameter is responsible for its own value + (inherited). + bounds : list of float + The lower and upper bounds on the Parameter, which some + optimizers use when the Parameter is varied (inherited). + """ + + def __init__( + self, + name, + atom1, + atom2, + atom3, + atom4, + value=None, + const=False, + mode=None, + ): + """Initialize the dihedral angle Parameter. + + Parameters + ---------- + name : str + The name of the Parameter. + atom1 : ObjCrystMolAtomParSet + The first atom in the dihedral angle. + atom2 : ObjCrystMolAtomParSet + The second (central) atom in the dihedral angle. + atom3 : ObjCrystMolAtomParSet + The third (central) atom in the dihedral angle. + atom4 : ObjCrystMolAtomParSet + The fourth (mutated) atom in the dihedral angle. + value : float, optional + The initial dihedral angle in radians. If None (default), the + current dihedral angle between the atoms is used. + const : bool, optional + The flag indicating whether the Parameter is constant + (default False). + mode : pyobjcryst.molecule.StretchModeTorsion, optional + The existing stretch mode to use. If None (default), a new + StretchModeTorsion is built. + """ + # Create the stretch mode + self.mode = mode + if mode is None: + self.mode = StretchModeTorsion( + atom2.scatterer, atom3.scatterer, None + ) + # We only add the last atom. This is the one that will move + self.mode.AddAtom(atom4.scatterer) + self.mutated_atoms = 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.scatterer, + atom2.scatterer, + atom3.scatterer, + atom4.scatterer, + ) + 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 : bool, optional + The flag indicating if the Parameter is constant (default + True). + value : float, optional + The value to set the Parameter to (default None). If this is + not None, the Parameter gets a new value, constant or + otherwise. + + Returns + ------- + StretchModeParameter + Return 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 get_value(self): + """Return the dihedral angle, recalculating it if needed. + + The atoms underlying the dihedral angle may have been moved by + another Parameter, so the angle is recalculated whenever the + cached value has been cleared. + + Returns + ------- + float + The dihedral angle in radians. + """ + if self._value is None: + value = GetDihedralAngle( + self.atom1.scatterer, + self.atom2.scatterer, + self.atom3.scatterer, + self.atom4.scatterer, + ) + Parameter.set_value(self, value) + + return self._value + + +# End class ObjCrystDihedralAngleParameter + + +class ObjCrystCrystalParSet(SrRealParSet): + """Adapt a pyobjcryst.crystal.Crystal to the ParameterSet interface. + + This class derives from SrRealParSet. See that class for base + attributes. + + Attributes + ---------- + structure : pyobjcryst.crystal.Crystal + The adapted crystal. + scatterers : list of ObjCrystAtomParSet or ObjCrystMoleculeParSet + The scatterer ParameterSets, provided for convenience. + space_group_parameters : SpaceGroupParameters + The free structure Parameters after applying the crystal's space + group constraints, created when first accessed. See the + diffpy.cmistructure.sgconstraints module. + angle_units : str + The units of the lattice angles, always "rad". + a, b, c, alpha, beta, gamma : ParameterAdapter + The unit cell parameters. + """ + + def __init__(self, name, crystal): + """Initialize the crystal ParameterSet. + + Parameters + ---------- + name : str + The name of this ParameterSet. + crystal : pyobjcryst.crystal.Crystal + The crystal to adapt. + + Raises + ------ + ValueError + If a scatterer in the crystal has no name, or if two + scatterers share a name. Give every scatterer a unique name + before wrapping the crystal. + TypeError + If the crystal contains a scatterer that is neither an Atom + nor a Molecule. + """ + SrRealParSet.__init__(self, name) + self.angle_units = "rad" + self.structure = crystal + self._space_group_parameters = None + + self.add_parameter(ParameterAdapter("a", self.structure, attr="a")) + self.add_parameter(ParameterAdapter("b", self.structure, attr="b")) + self.add_parameter(ParameterAdapter("c", self.structure, attr="c")) + self.add_parameter( + ParameterAdapter("alpha", self.structure, attr="alpha") + ) + self.add_parameter( + ParameterAdapter("beta", self.structure, attr="beta") + ) + self.add_parameter( + ParameterAdapter("gamma", self.structure, attr="gamma") + ) + + # Now we must loop over the scatterers and create parameter sets from + # them. + self.scatterers = [] + snames = [] + + for j in range(self.structure.GetNbScatterer()): + s = self.structure.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": + parameter_set = ObjCrystAtomParSet(name, s, self) + elif cname == "Molecule": + parameter_set = ObjCrystMoleculeParSet(name, s, self) + else: + raise TypeError("Unrecognized scatterer '%s'" % cname) + + self.add_parameter_set(parameter_set) + self.scatterers.append(parameter_set) + snames.append(name) + + return + + def _constrain_space_group(self): + """Constrain the space group.""" + if self._space_group_parameters is not None: + return self._space_group_parameters + space_group = self._create_space_group(self.structure.GetSpaceGroup()) + from diffpy.cmistructure.sgconstraints import ( + _constrain_as_space_group, + ) + + adpsymbols = ["B11", "B22", "B33", "B12", "B13", "B23"] + isosymbol = "Biso" + sgoffset = [0, 0, 0] + self._space_group_parameters = _constrain_as_space_group( + self, + space_group, + self.scatterers, + sgoffset, + adpsymbols=adpsymbols, + isosymbol=isosymbol, + ) + return self._space_group_parameters + + space_group_parameters = 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 SymOp, get_space_group + + 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 + # get_space_group, as this function messes with symop_list. + space_group = copy.copy(get_space_group(name)) + + # Replace the symmetry operations to guarantee that we get it right. + symops = sgobjcryst.GetSymmetryOperations() + tranops = sgobjcryst.GetTranslationVectors() + space_group.symop_list = [] + + for trans in tranops: + for shift, rot in symops: + tv = trans + shift + tv -= numpy.floor(tv) + space_group.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) + space_group.symop_list.append(SymOp(-rot, tv)) + + return space_group + + @classmethod + def can_adapt(self, structure): + """Return whether the structure can be adapted by this class. + + Parameters + ---------- + structure : object + The structure object to check. + + Returns + ------- + bool + The flag indicating if `structure` is a pyobjcryst Crystal. + """ + from pyobjcryst.crystal import Crystal + + return isinstance(structure, Crystal) + + def get_lattice(self): + """Return the ParameterSet containing the lattice Parameters. + + Returns + ------- + ObjCrystCrystalParSet + This ParameterSet, which holds the lattice Parameters + directly. + """ + return self + + def get_scatterers(self): + """Return the list of ParameterSets that represent the + scatterers. + + Returns + ------- + list of ObjCrystAtomParSet or ObjCrystMoleculeParSet + The scatterer ParameterSets of the crystal. + """ + return self.scatterers + + +# End class ObjCrystCrystalParSet diff --git a/src/diffpy/cmistructure/sgconstraints.py b/src/diffpy/cmistructure/sgconstraints.py new file mode 100644 index 0000000..38b01b8 --- /dev/null +++ b/src/diffpy/cmistructure/sgconstraints.py @@ -0,0 +1,832 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.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 a P1 structure to a space group. + + This applies space group constraints to a structure ParameterSet with + P1 symmetry. The passed scatterers are explicitly constrained to the + specified space group, and the ADPs and lattice may be constrained as + well. New Parameters used in the constraints are created within the + returned SpaceGroupParameters object. Constraints are created in the + ParameterSet that contains the constrained Parameter. This erases any + constraints or constant flags on the scatterers, lattice or ADPs that + are to be constrained. + + Parameters + ---------- + phase : BaseStructureParSet + The structure ParameterSet to constrain. + spacegroup : int, str or diffpy.structure.spacegroups.SpaceGroup + The space group number, symbol or SpaceGroup instance. + scatterers : list of ParameterSet, optional + The scatterer ParameterSets to constrain. If None (default), all + scatterers returned by ``phase.get_scatterers()`` are constrained. + sgoffset : list of float, optional + The offset of the space group origin (default [0, 0, 0]). + constrainlat : bool, optional + The flag indicating whether to constrain the lattice (default + True). + constrainadps : bool, optional + The flag indicating whether to constrain the ADPs (default True). + adpsymbols : list of str, optional + The ADP names. By default this is + diffpy.structure.symmetryutilities.stdUsymbols (U11, U22, etc.). + The names must be given in the same order as stdUsymbols. + isosymbol : str, optional + The name of the isotropic ADP (default "Uiso"). If None, + isotropic ADPs are constrained via the anisotropic ADPs. + + Returns + ------- + SpaceGroupParameters + The free Parameters of the structure that remain after applying + the space group constraints. + + Notes + ----- + The lattice constraints are applied as follows. + + 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, and 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 SpaceGroup, get_space_group + + space_group = spacegroup + if not isinstance(spacegroup, SpaceGroup): + space_group = get_space_group(spacegroup) + sgp = _constrain_as_space_group( + phase, + space_group, + scatterers, + sgoffset, + constrainlat, + constrainadps, + adpsymbols, + isosymbol, + ) + + return sgp + + +def _constrain_as_space_group( + phase, + space_group, + 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, + space_group, + 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 stores the variable Parameters of a structure, leaving out + those that are constrained or fixed by the space group. It has the same + Parameter attribute access as a ParameterSet, which makes it easy to + access the free variables of a structure when scripting. + + Attributes + ---------- + name : str + The name of this container (default "sgpars"). + """ + + def __init__(self, name="sgpars"): + """Initialize the space group Parameter container. + + Parameters + ---------- + name : str, optional + The name of this container (default "sgpars"). + """ + RecipeContainer.__init__(self, name) + return + + def add_parameter(self, parameter, check=True): + """Store a Parameter. + + Parameters + ---------- + parameter : Parameter + The Parameter to be stored. + check : bool, optional + The flag indicating whether to check for an existing Parameter + of the same name (default True). + + Raises + ------ + ValueError + If the Parameter has no name, or if `check` is True and a + Parameter of the same name has already been stored. + """ + # Store the Parameter + RecipeContainer._add_object(self, parameter, self._parameters, check) + return + + +# End class BaseSpaceGroupParameters + + +class SpaceGroupParameters(BaseSpaceGroupParameters): + """Create and hold the free Parameters of a space group constraint. + + This class stores the variable Parameters of a structure, leaving out + those that are constrained or fixed by the space group, and does the + work of constrain_as_space_group. It has the same Parameter attribute + access as a ParameterSet. + + Attributes + ---------- + name : str + The name of this container, always "sgpars". + phase : BaseStructureParSet + The constrained structure ParameterSet. + space_group : diffpy.structure.spacegroups.SpaceGroup + The space group of the constraints. + sgoffset : list of float + The offset of the space group origin. + scatterers : list of ParameterSet + The constrained scatterer ParameterSets. + constrainlat : bool + The flag indicating whether the lattice is constrained. + constrainadps : bool + The flag indicating whether the ADPs are constrained. + adpsymbols : list of str + The ADP names. + isosymbol : str or None + The name of the isotropic ADP. + xyz_parameters : BaseSpaceGroupParameters + The free xyz Parameters, created on first access. + lattice_parameters : BaseSpaceGroupParameters + The free lattice Parameters, created on first access. + adp_parameters : BaseSpaceGroupParameters + The free ADP Parameters, created on first access. + """ + + def __init__( + self, + phase, + space_group, + scatterers, + sgoffset, + constrainlat, + constrainadps, + adpsymbols, + isosymbol, + ): + """Initialize the space group Parameters. + + The constraints are not applied until the Parameters are first + accessed. + + Parameters + ---------- + phase : BaseStructureParSet + The structure ParameterSet to be constrained. + space_group : diffpy.structure.spacegroups.SpaceGroup + The space group of the constraints. + scatterers : list of ParameterSet + The scatterer ParameterSets to constrain. + sgoffset : list of float + The offset of the space group origin. + constrainlat : bool + The flag indicating whether to constrain the lattice. + constrainadps : bool + The flag indicating whether to constrain the ADPs. + adpsymbols : list of str + The ADP names, in the same order as + diffpy.structure.symmetryutilities.stdUsymbols. + isosymbol : str or None + The name of the isotropic ADP. If None, isotropic ADPs are + constrained via the anisotropic ADPs. + """ + BaseSpaceGroupParameters.__init__(self) + self._lattice_parameters = None + self._xyz_parameters = None + self._adp_parameters = None + + self._parsets = {} + self._manage(self._parsets) + + self.phase = phase + self.space_group = space_group + 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._lattice_parameters is None + or self._xyz_parameters is None + or self._adp_parameters is None + ): + self._make_constraints() + return RecipeContainer.__iter__(self) + + lattice_parameters = property(lambda self: self._get_lat_pars()) + + def _get_lat_pars(self): + """Accessor for _lattice_parameters.""" + if self._lattice_parameters is None: + self._constrain_lattice() + return self._lattice_parameters + + xyz_parameters = property(lambda self: self._get_xyz_pars()) + + def _get_xyz_pars(self): + """Accessor for _xyz_parameters.""" + positions = [] + for scatterer in self.scatterers: + xyz = [scatterer.x, scatterer.y, scatterer.z] + positions.append([p.value for p in xyz]) + if self._xyz_parameters is None: + self._constrain_xyzs(positions) + return self._xyz_parameters + + adp_parameters = property(lambda self: self._get_adp_pars()) + + def _get_adp_pars(self): + """Accessor for _adp_parameters.""" + positions = [] + for scatterer in self.scatterers: + xyz = [scatterer.x, scatterer.y, scatterer.z] + positions.append([p.value for p in xyz]) + if self._adp_parameters is None: + self._constrain_adps(positions) + return self._adp_parameters + + 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 parameter in [scatterer.x, scatterer.y, scatterer.z]: + if scatterer.is_constrained(parameter): + scatterer.remove_constraint(parameter) + parameter.set_constant(False) + + # Clear the lattice + if self.constrainlat: + + lattice = phase.get_lattice() + lattice_parameters = [ + lattice.a, + lattice.b, + lattice.c, + lattice.alpha, + lattice.beta, + lattice.gamma, + ] + for parameter in lattice_parameters: + if lattice.is_constrained(parameter): + lattice.remove_constraint(parameter) + parameter.set_constant(False) + + # Clear ADPs + if self.constrainadps: + for scatterer in scatterers: + if isosymbol: + parameter = scatterer.get(isosymbol) + if parameter is not None: + if scatterer.is_constrained(parameter): + scatterer.remove_constraint(parameter) + parameter.set_constant(False) + + for pname in adpsymbols: + parameter = scatterer.get(pname) + if parameter is not None: + if scatterer.is_constrained(parameter): + scatterer.remove_constraint(parameter) + parameter.set_constant(False) + + return + + def _constrain_lattice(self): + """Constrain the lattice parameters.""" + if not self.constrainlat: + return + + phase = self.phase + space_group = self.space_group + + lattice = phase.get_lattice() + system = space_group.crystal_system + if not system: + system = "Triclinic" + system = system.title() + # This makes the constraints + f = _constraint_map[system] + f(lattice) + + # Now get the unconstrained, non-constant lattice pars and store them. + self._lattice_parameters = BaseSpaceGroupParameters( + "lattice_parameters" + ) + lattice_parameters = [ + lattice.a, + lattice.b, + lattice.c, + lattice.alpha, + lattice.beta, + lattice.gamma, + ] + pars = [ + p for p in lattice_parameters if not p.const and not p.constrained + ] + for parameter in pars: + # FIXME - the original parameter will still appear as + # constrained. + newpar = self.__add_par(parameter.name, parameter) + self._lattice_parameters.add_parameter(newpar) + + return + + def _constrain_xyzs(self, positions): + """Constrain the positions. + + Parameters + ---------- + positions + The coordinates of the scatterers. + """ + from diffpy.structure.symmetryutilities import SymmetryConstraints + + space_group = self.space_group + sgoffset = self.sgoffset + + # We do this without ADPs here so we can skip much complication. See + # the _constrain_adps method for details. + g = SymmetryConstraints(space_group, positions, sgoffset=sgoffset) + + scatterers = self.scatterers + self._xyz_parameters = BaseSpaceGroupParameters("xyz_parameters") + + # Make proxies to the free xyz parameters + xyznames = [name[:1] + "_" + name[1:] for name, value in g.pospars] + for pname in xyznames: + name, index = pname.rsplit("_", 1) + index = int(index) + parameter = scatterers[index].get(name) + newpar = self.__add_par(pname, parameter) + self._xyz_parameters.add_parameter(newpar) + + # Constrain non-free xyz parameters + fpos = g.position_formulas(xyznames) + for index, 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, index, 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 + + space_group = self.space_group + sgoffset = self.sgoffset + scatterers = self.scatterers + isosymbol = self.isosymbol + adpsymbols = self.adpsymbols + adpmap = dict(zip(stdUsymbols, adpsymbols)) + self._adp_parameters = BaseSpaceGroupParameters("adp_parameters") + + # 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 index, parameter in enumerate(pars): + i, j = _idxtoij[index] + Uij[i, j] = Uij[j, i] = parameter.get_value() + + Uijs.append(Uij) + + # Discard any positions for the nonadps + positions = list(positions) + nonadps.reverse() + [positions.pop(index) for index in nonadps] + + # Now we can create symmetry constraints without having to worry about + # the nonadps + g = SymmetryConstraints( + space_group, positions, Uijs, sgoffset=sgoffset + ) + + adpnames = [ + adpmap[name[:3]] + "_" + name[3:] for name, value 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, index = pname.rsplit("_", 1) + index = int(index) + # Check for isotropic ADPs + scatterer = scatterers[index] + if isosymbol and g.Uisotropy[index] and index not in isoidx: + isoidx.append(index) + parameter = scatterer.get(isosymbol) + if parameter is not None: + parname = "%s_%i" % (isosymbol, index) + newpar = self.__add_par(parname, parameter) + self._adp_parameters.add_parameter(newpar) + isonames.append(newpar.name) + else: + parameter = scatterer.get(name) + if parameter is not None: + newpar = self.__add_par(pname, parameter) + self._adp_parameters.add_parameter(newpar) + + # Constrain dependent isotropics + for index, isoname in zip(isoidx[:], isonames): + for j in g.coremap[index]: + if j == index: + continue + isoidx.append(j) + scatterer = scatterers[j] + scatterer.add_constraint( + isosymbol, isoname, params=self._parameters + ) + + fadp = g.u_formulas(adpnames) + + # Constrain dependent anisotropics. We use the fact that an + # anisotropic cannot be dependent on an isotropic. + for index, tmp in enumerate(zip(scatterers, fadp)): + if index 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, index, self._parameters + ) + + def __add_par(self, parname, parameter): + """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, parameter) + self.add_parameter(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.angle_units == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + lattice.alpha.set_constant(True, ang90) + beta = lattice.beta.get_value() + gamma = lattice.gamma.get_value() + + 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.angle_units == "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.angle_units == "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.angle_units == "rad": + afactor = deg2rad + ang90 = 90.0 * afactor + ang120 = 120.0 * afactor + if lattice.gamma.get_value() == 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.angle_units == "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.angle_units == "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. +_constraint_map = { + "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, index, 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. + """ + parameter = scatterer.get(parname) + + if parameter is None: + return + + compname = "%s_%i" % (parname, index) + + # Check to see if this parameter is free + pat = r"%s *([+-] *\d+)?$" % compname + if re.match(pat, formula): + return parameter + + # Check to see if it is a constant + fval = _get_float(formula) + if fval is not None: + parameter.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(parameter, 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/cmistructure/srrealparset.py b/src/diffpy/cmistructure/srrealparset.py new file mode 100644 index 0000000..d9ed71c --- /dev/null +++ b/src/diffpy/cmistructure/srrealparset.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +############################################################################## +"""Structure wrapper class for structures compatible with SrReal.""" + +__all__ = ["SrRealParSet"] + +from diffpy.cmistructure.basestructureparset import BaseStructureParSet +from diffpy.cmistructure.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 + ---------- + structure : object + The adapted structure object. + _usesymmetry : bool + The flag indicating if SrReal calculators that operate on + this object should use symmetry (default True). + """ + + def __init__(self, *args, **kw): + BaseStructureParSet.__init__(self, *args, **kw) + self._usesymmetry = True + self.structure = 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 : float, optional + The uncertainty on the BVS (default 1). + scaled : bool, optional + The flag indicating if the restraint is scaled (multiplied) + by the unrestrained point-average chi^2 (chi^2/numpoints) + (default False). + + Returns + ------- + BVSRestraint + The restraint object, for use with the ``unrestrain`` method. + """ + # Create the Restraint object + restraint = BVSRestraint(self, sig, scaled) + # Add it to the _restraints set + self._restraints.add(restraint) + # Our configuration changed. Notify observers. + self._update_configuration() + # Return the Restraint object + return restraint + + def use_symmetry(self, use=True): + """Set whether this structure uses symmetry. + + This determines how the structure is treated by SrReal + calculators. + + Parameters + ---------- + use : bool, optional + The flag indicating if symmetry is used (default True). + """ + self._usesymmetry = bool(use) + return + + def using_symmetry(self): + """Return whether symmetry is being used. + + Returns + ------- + bool + The flag indicating if symmetry is 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.structure + return nosymmetry(self.structure) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/conftest.py b/tests/conftest.py index e3b6313..8b363f0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,3 +1,4 @@ +import importlib.resources import json from pathlib import Path @@ -17,3 +18,14 @@ def user_filesystem(tmp_path): json.dump(home_config_data, f) yield tmp_path + + +@pytest.fixture(scope="session") +def datafile(): + """Fixture to load a test data file from the testdata package + directory.""" + + def _datafile(filename): + return importlib.resources.files("tests.testdata").joinpath(filename) + + return _datafile diff --git a/tests/test_diffpyparset.py b/tests/test_diffpyparset.py new file mode 100644 index 0000000..4c38ba6 --- /dev/null +++ b/tests/test_diffpyparset.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2010 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Pavol Juhas and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +############################################################################## +"""Tests for diffpy.cmistructure package.""" + +import pickle +import unittest + +import numpy as np + +from diffpy.cmistructure.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) + + actual_name = s.name + expected_name = "CuAg" + assert actual_name == expected_name + + def _testAtoms(): + # Check the atoms thoroughly + actual_atoms = { + "Cu0": { + "element": s.Cu0.element, + "Uiso": s.Cu0.Uiso.get_value(), + "Biso": s.Cu0.Biso.get_value(), + "xyz": [ + s.Cu0.x.get_value(), + s.Cu0.y.get_value(), + s.Cu0.z.get_value(), + ], + }, + "Ag0": { + "element": s.Ag0.element, + "Uiso": s.Ag0.Uiso.get_value(), + "Biso": s.Ag0.Biso.get_value(), + }, + } + expected_atoms = { + "Cu0": { + "element": a1.element, + "Uiso": a1.Uisoequiv, + "Biso": a1.Bisoequiv, + "xyz": [a1.xyz[0], a1.xyz[1], a1.xyz[2]], + }, + "Ag0": { + "element": a2.element, + "Uiso": a2.Uisoequiv, + "Biso": a2.Bisoequiv, + }, + } + assert actual_atoms == expected_atoms + + # The Uij and Uji (Bij and Bji) Parameters both read the + # structure's Uij (Bij). + actual_anisotropic = {} + expected_anisotropic = {} + for i in range(1, 4): + for j in range(i, 4): + for prefix in "UB": + ij = "%s%i%i" % (prefix, i, j) + ji = "%s%i%i" % (prefix, j, i) + actual_anisotropic[ij] = getattr(s.Cu0, ij).get_value() + actual_anisotropic[ji] = getattr(s.Cu0, ji).get_value() + expected_anisotropic[ij] = getattr(a1, ij) + expected_anisotropic[ji] = getattr(a1, ij) + assert actual_anisotropic == expected_anisotropic + return + + def _testLattice(): + # Test the lattice + lattice_names = ["a", "b", "c", "alpha", "beta", "gamma"] + actual_lattice = [ + getattr(s.lattice, name).get_value() for name in lattice_names + ] + expected_lattice = [ + getattr(dsstru.lattice, name) for name in lattice_names + ] + assert actual_lattice == expected_lattice + + # C1: The ParameterSet has just been created from the structure. + # Expected: The Parameters match the atoms and lattice. + _testAtoms() + _testLattice() + + # C2: The diffpy Structure is changed directly. + # Expected: The Parameters follow the changes. + a1.xyz[1] = 0.123 + a1.U11 = 0.321 + a1.B32 = 0.111 + dsstru.lattice.set_latt_parms(a=3.0, gamma=121) + _testAtoms() + _testLattice() + + # C3: The Parameters of the DiffpyStructureParSet are changed. + # Expected: The structure follows the changes, so the distance + # between the atoms changes. + 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() + actual_distance_changed = d != dsstru.lattice.dist(a1.xyz, a2.xyz) + expected_distance_changed = True + assert actual_distance_changed == expected_distance_changed + 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]) + structure = Structure([atom], lattice=lat) + dsps = DiffpyStructureParSet("dsps", structure) + # C1: The structure, lattice and atom ParameterSets are printed. + # Expected: Each repr matches the repr of the adapted object. + actual_reprs = [repr(dsps), repr(dsps.lattice), repr(dsps.atoms[0])] + expected_reprs = [repr(structure), repr(lat), repr(atom)] + assert actual_reprs == expected_reprs + return + + +def test_pickling(): + """Test pickling of DiffpyStructureParSet.""" + from diffpy.structure import Atom, Structure + + structure = Structure([Atom("C", [0, 0.2, 0.5])]) + dsps = DiffpyStructureParSet("dsps", structure) + data = pickle.dumps(dsps) + dsps2 = pickle.loads(data) + # C1: A DiffpyStructureParSet is pickled and unpickled. + # Expected: The copy keeps its single atom and the atom's position. + actual_atom_count = len(dsps2.atoms) + expected_atom_count = 1 + assert actual_atom_count == expected_atom_count + actual_y = dsps2.atoms[0].y.value + expected_y = 0.2 + assert actual_y == expected_y + return + + +# End of class TestParameterAdapter + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_functions.py b/tests/test_functions.py deleted file mode 100644 index 1ec72dd..0000000 --- a/tests/test_functions.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -import pytest - -from diffpy.cmistructure import functions # noqa - - -def test_dot_product_2D_list(): - a = [1, 2] - b = [3, 4] - expected = 11.0 - actual = functions.dot_product(a, b) - assert actual == expected - - -def test_dot_product_3D_list(): - a = [1, 2, 3] - b = [4, 5, 6] - expected = 32.0 - actual = functions.dot_product(a, b) - assert actual == expected - - -@pytest.mark.parametrize( - "a, b, expected", - [ - # Test whether the dot product function works with 2D and 3D vectors - # C1: lists, expect correct float output - ([1, 2], [3, 4], 11.0), - ([1, 2, 3], [4, 5, 6], 32.0), - # C2: tuples, expect correct float output - ((1, 2), (3, 4), 11.0), - ((1, 2, 3), (4, 5, 6), 32.0), - # C3: numpy arrays, expect correct float output - (np.array([1, 2]), np.array([3, 4]), 11.0), - (np.array([1, 2, 3]), np.array([4, 5, 6]), 32.0), - ], -) -def test_dot_product(a, b, expected): - actual = functions.dot_product(a, b) - assert actual == expected diff --git a/tests/test_objcrystparset.py b/tests/test_objcrystparset.py new file mode 100644 index 0000000..f236938 --- /dev/null +++ b/tests/test_objcrystparset.py @@ -0,0 +1,807 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2010 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Pavol Juhas and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +############################################################################## +"""Tests for diffpy.cmistructure 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): + # shared setup + 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.cmistructure.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 + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + def testObjCrystParSet(self): + """Test the structure conversion.""" + occryst = self.occryst + ocmol = self.ocmol + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + actual_name = crystal.name + expected_name = "bucky" + assert actual_name == expected_name + + def _testCrystal(): + # Test the lattice + actual_lattice = [ + crystal.a.value, + crystal.b.get_value(), + crystal.c.get_value(), + crystal.alpha.get_value(), + crystal.beta.get_value(), + crystal.gamma.get_value(), + ] + expected_lattice = [ + occryst.a, + occryst.b, + occryst.c, + occryst.alpha, + occryst.beta, + occryst.gamma, + ] + assert actual_lattice == pytest.approx(expected_lattice) + return + + def _testMolecule(): + # Test position, occupancy and orientation + actual_molecule = [ + m.x.get_value(), + m.y.get_value(), + m.z.get_value(), + m.occ.get_value(), + m.q0.get_value(), + m.q1.get_value(), + m.q2.get_value(), + m.q3.get_value(), + ] + expected_molecule = [ + ocmol.X, + ocmol.Y, + ocmol.Z, + ocmol.Occupancy, + ocmol.Q0, + ocmol.Q1, + ocmol.Q2, + ocmol.Q3, + ] + assert actual_molecule == pytest.approx(expected_molecule) + + # Check the atoms thoroughly + actual_elements = [a.element for a in m.atoms] + expected_elements = [ + ocmol[i].GetScatteringPower().GetSymbol() + for i in range(len(ocmol)) + ] + assert actual_elements == expected_elements + actual_atoms = [ + [ + a.x.get_value(), + a.y.get_value(), + a.z.get_value(), + a.occ.get_value(), + a.Biso.get_value(), + ] + for a in m.atoms + ] + expected_atoms = [ + pytest.approx( + [ + ocmol[i].X, + ocmol[i].Y, + ocmol[i].Z, + ocmol[i].Occupancy, + ocmol[i].GetScatteringPower().Biso, + ] + ) + for i in range(len(ocmol)) + ] + assert actual_atoms == expected_atoms + return + + # C1: The ParameterSet has just been created from the crystal. + # Expected: The Parameters match the pyobjcryst values. + _testCrystal() + _testMolecule() + + # C2: Values are changed through pyobjcryst. + # Expected: The Parameters follow the changes. + 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() + + # C3: Values are changed through the ParameterSet. + # Expected: The pyobjcryst objects follow the changes. + crystal.c60.C44.x.set_value(1.1) + crystal.c60.C44.occ.set_value(1.1) + crystal.c60.C44.Biso.set_value(1.1) + crystal.c60.q3.set_value(1.1) + crystal.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 + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + bonds = ocmol.GetBondList() + expected_penalties = set( + [bonds[0].GetLogLikelihood(), bonds[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + 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 + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + 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 + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.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) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + bonds = ocmol.GetBondList() + actual_bond_count = len(bonds) + expected_bond_count = 2 + assert actual_bond_count == expected_bond_count + actual_penalties = [res0.penalty(), res1.penalty()] + expected_penalties = [b.GetLogLikelihood() for b in bonds] + assert actual_penalties == expected_penalties + + 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 + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.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 + ) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + def testExplicitDihedralAngleRestraints(self): + """Test the structure with explicit dihedral angles.""" + occryst = self.occryst + ocmol = self.ocmol + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.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 + ) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + def testExplicitBondLengthParameter(self): + """Test adding bond length parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.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.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7 = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + + dd = xyz0 - xyz7 + d0 = numpy.dot(dd, dd) ** 0.5 + # C1: A bond length Parameter is added. + # Expected: Its value is the current bond length. + actual_length = p1.get_value() + expected_length = d0 + assert actual_length == pytest.approx(expected_length, abs=1e-6) + + # Record the unit direction of change for later + u = dd / d0 + + # C2: The bond length Parameter is stretched by 5%. + # Expected: The Parameter and the measured bond length both take + # the new value, the first atom stays put, and the second and + # tag-along atoms move along the bond. + scale = 1.05 + p1.set_value(scale * d0) + + actual_length = p1.get_value() + expected_length = scale * d0 + assert actual_length == pytest.approx(expected_length, abs=1e-6) + + xyz0a = numpy.array( + [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7a = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + + dda = xyz0a - xyz7a + d1 = numpy.dot(dda, dda) ** 0.5 + + actual_measured_length = d1 + expected_measured_length = scale * d0 + assert actual_measured_length == pytest.approx( + expected_measured_length, abs=1e-6 + ) + + actual_xyz0 = xyz0a.tolist() + expected_xyz0 = xyz0.tolist() + assert actual_xyz0 == expected_xyz0 + + actual_xyz7 = xyz7a + expected_xyz7 = xyz7 + (1 - scale) * d0 * u + assert actual_xyz7 == pytest.approx(expected_xyz7, abs=1e-5) + + actual_xyz20 = xyz20a + expected_xyz20 = xyz20 + (1 - scale) * d0 * u + assert actual_xyz20 == pytest.approx(expected_xyz20, abs=1e-6) + + return + + def testExplicitBondAngleParameter(self): + """Test adding bond angle parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + a25 = m.atoms[25] + + xyz0 = numpy.array( + [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7 = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25 = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + + 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]) + + # C1: A bond angle Parameter is added. + # Expected: Its value is the current bond angle. + actual_angle = p1.get_value() + expected_angle = angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + # C2: The bond angle Parameter is stretched by 5%. + # Expected: The Parameter and the measured angle both take the new + # value, and only the third and tag-along atoms move. + scale = 1.05 + p1.set_value(scale * angle0) + + actual_angle = p1.get_value() + expected_angle = scale * angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + xyz0a = numpy.array( + [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7a = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25a = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + + 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)) + + actual_measured_angle = angle1 + expected_measured_angle = scale * angle0 + assert actual_measured_angle == pytest.approx( + expected_measured_angle, abs=1e-6 + ) + + actual_moved = { + "C0": not numpy.array_equal(xyz0, xyz0a), + "C7": not numpy.array_equal(xyz7, xyz7a), + "C20": not numpy.array_equal(xyz20, xyz20a), + "C25": not numpy.array_equal(xyz25, xyz25a), + } + expected_moved = {"C0": False, "C7": False, "C20": True, "C25": True} + assert actual_moved == expected_moved + + return + + def testExplicitDihedralAngleParameter(self): + """Test adding dihedral angle parameters to the molecule.""" + occryst = self.occryst + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.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.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7 = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25 = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + xyz33 = numpy.array( + [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] + ) + + 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]) + + # C1: A dihedral angle Parameter is added. + # Expected: Its value is the current dihedral angle. + actual_angle = p1.get_value() + expected_angle = angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + # C2: The dihedral angle Parameter is stretched by 5%. + # Expected: The Parameter and the measured angle both take the new + # value, and only the fourth and tag-along atoms move. + scale = 1.05 + p1.set_value(scale * angle0) + + actual_angle = p1.get_value() + expected_angle = scale * angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + xyz0a = numpy.array( + [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] + ) + xyz7a = numpy.array( + [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] + ) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25a = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + xyz33a = numpy.array( + [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] + ) + + 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)) + actual_measured_angle = angle1 + expected_measured_angle = scale * angle0 + assert actual_measured_angle == pytest.approx( + expected_measured_angle, abs=1e-6 + ) + + actual_moved = { + "C0": not numpy.array_equal(xyz0, xyz0a), + "C7": not numpy.array_equal(xyz7, xyz7a), + "C20": not numpy.array_equal(xyz20, xyz20a), + "C25": not numpy.array_equal(xyz25, xyz25a), + "C33": not numpy.array_equal(xyz33, xyz33a), + } + expected_moved = { + "C0": False, + "C7": False, + "C20": False, + "C25": True, + "C33": True, + } + assert actual_moved == expected_moved + + 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): + # shared setup + global ObjCrystCrystalParSet, spacegroups + from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet + from diffpy.structure import spacegroups + + @staticmethod + def getObjCrystParSetSpaceGroup(space_group): + """Make an ObjCrystCrystalParSet with the proper space group.""" + from pyobjcryst.spacegroup import SpaceGroup + + sgobjcryst = SpaceGroup(space_group.short_name) + sgnew = ObjCrystCrystalParSet._create_space_group(sgobjcryst) + return sgnew + + @staticmethod + def hashDiffPySpaceGroup(space_group): + lines = [str(space_group.number % 1000)] + sorted( + map(str, space_group.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.""" + from cctbx import sgtbx + + for smbls in sgtbx.space_group_symbol_iterator(): + shn = smbls.hermann_mauguin() + short_name = shn.replace(" ", "") + if spacegroups.is_space_group_identifier(short_name): + space_group = spacegroups.get_space_group(shn) + sgnew = self.getObjCrystParSetSpaceGroup(space_group) + actual_equivalent = self.sgsEquivalent(space_group, sgnew) + expected_equivalent = True + assert actual_equivalent == expected_equivalent + 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..15ce753 --- /dev/null +++ b/tests/test_sgconstraints.py @@ -0,0 +1,300 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2010 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Pavol Juhas and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +############################################################################## +"""Tests space group constraints.""" + +import unittest + +import numpy +import pytest + +# ---------------------------------------------------------------------------- + + +def test_ObjCryst_constrain_space_group(): + """Make sure that all Parameters are constrained properly. + + This tests constrainSpaceGroup from + diffpy.cmistructure.sgconstraints, which is performed automatically + when an ObjCrystCrystalParSet is created. + """ + from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet + + pi = numpy.pi + + occryst = makeLaMnO3() + structure = ObjCrystCrystalParSet(occryst.GetName(), occryst) + # Make sure we actually create the constraints + structure._constrain_space_group() + # Make the space group parameters individually + structure.space_group_parameters.lattice_parameters + structure.space_group_parameters.xyz_parameters + structure.space_group_parameters.adp_parameters + + # C1: The orthorhombic lattice of LaMnO3 in P b n m. + # Expected: The angles are fixed at pi / 2, the lengths are free, and + # no constraint equations are needed. + lattice = structure.get_lattice() + lattice_names = ["a", "b", "c", "alpha", "beta", "gamma"] + actual_lattice_const = { + name: getattr(lattice, name).const for name in lattice_names + } + expected_lattice_const = { + "a": False, + "b": False, + "c": False, + "alpha": True, + "beta": True, + "gamma": True, + } + assert actual_lattice_const == expected_lattice_const + actual_angles = [ + lattice.alpha.get_value(), + lattice.beta.get_value(), + lattice.gamma.get_value(), + ] + expected_angles = [pi / 2, pi / 2, pi / 2] + assert actual_angles == expected_angles + actual_lattice_constraint_count = len(lattice._constraints) + expected_lattice_constraint_count = 0 + assert actual_lattice_constraint_count == expected_lattice_constraint_count + + # C2: The scatterers of LaMnO3 on their P b n m sites. + # Expected: Coordinates on special positions are fixed, the rest are + # free, and no constraint equations are needed. + scatterers = structure.get_scatterers() + la, mn, o1, o2 = scatterers + actual_xyz_const = { + "La1": [la.x.const, la.y.const, la.z.const], + "Mn1": [mn.x.const, mn.y.const, mn.z.const], + "O1": [o1.x.const, o1.y.const, o1.z.const], + "O2": [o2.x.const, o2.y.const, o2.z.const], + } + expected_xyz_const = { + "La1": [False, False, True], + "Mn1": [True, True, True], + "O1": [False, False, True], + "O2": [False, False, False], + } + assert actual_xyz_const == expected_xyz_const + actual_constraint_counts = [len(s._constraints) for s in scatterers] + expected_constraint_counts = [0, 0, 0, 0] + assert actual_constraint_counts == expected_constraint_counts + + # C3: Fixed coordinates are constrained or made into variables. + # Expected: A ValueError is raised. + 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): + """Test the constrain_as_space_group function.""" + from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet + from diffpy.cmistructure.sgconstraints import constrain_as_space_group + + structure = makeLaMnO3_P1(datafile) + parameter_set = DiffpyStructureParSet("LaMnO3", structure) + + space_group_parameters = constrain_as_space_group( + parameter_set, + "P b n m", + scatterers=parameter_set.get_scatterers()[::2], + constrainadps=True, + ) + + # C1: The space group Parameters are created. + # Expected: Every Parameter exists and has a value. + actual_unset_parameters = [ + parameter + for parameter in space_group_parameters + if parameter is None or parameter.get_value() is None + ] + expected_unset_parameters = [] + assert actual_unset_parameters == expected_unset_parameters + + # C2: Scatterers that were not passed to constrain_as_space_group. + # Expected: Their positions and ADPs are free and unconstrained. + unconstrained = parameter_set.get_scatterers()[1::2] + actual_free_const = { + scatterer.name: [ + scatterer.x.const, + scatterer.y.const, + scatterer.z.const, + scatterer.U11.const, + scatterer.U22.const, + scatterer.U33.const, + scatterer.U12.const, + scatterer.U13.const, + scatterer.U23.const, + ] + for scatterer in unconstrained + } + expected_free_const = { + scatterer.name: [False] * 9 for scatterer in unconstrained + } + assert actual_free_const == expected_free_const + actual_free_constraint_counts = { + scatterer.name: len(scatterer._constraints) + for scatterer in unconstrained + } + expected_free_constraint_counts = { + scatterer.name: 0 for scatterer in unconstrained + } + assert actual_free_constraint_counts == expected_free_constraint_counts + + proxied = [p.par for p in space_group_parameters] + + def _consttest(parameter): + return parameter.const + + def _constrainedtest(parameter): + return parameter.constrained + + def _proxytest(parameter): + return parameter in proxied + + def _alltests(parameter): + return ( + _consttest(parameter) + or _constrainedtest(parameter) + or _proxytest(parameter) + ) + + # C3: Scatterers that were passed to constrain_as_space_group. + # Expected: At least one position and one ADP Parameter of each is + # fixed, constrained or proxied by a space group Parameter. + constrained = parameter_set.get_scatterers()[::2] + actual_restricted = { + scatterer.name: [ + any( + _alltests(parameter) + for parameter in [scatterer.x, scatterer.y, scatterer.z] + ), + any( + _alltests(parameter) + for parameter in [ + scatterer.U11, + scatterer.U22, + scatterer.U33, + scatterer.U12, + scatterer.U13, + scatterer.U23, + ] + ), + ] + for scatterer in constrained + } + expected_restricted = { + scatterer.name: [True, True] for scatterer in constrained + } + assert actual_restricted == expected_restricted + + return + + +def test_constrain_as_space_group_args(datafile): + """Test the arguments processing of constrain_as_space_group + function.""" + from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet + from diffpy.cmistructure.sgconstraints import constrain_as_space_group + from diffpy.structure.spacegroups import get_space_group + + # C1: The space group is given as a symbol or as a SpaceGroup object. + # Expected: Both create the same space group Parameters. + structure = makeLaMnO3_P1(datafile) + parameter_set = DiffpyStructureParSet("LaMnO3", structure) + symbol_parameters = constrain_as_space_group(parameter_set, "P b n m") + space_group = get_space_group("P b n m") + object_parameter_set = DiffpyStructureParSet( + "LMO", makeLaMnO3_P1(datafile) + ) + object_parameters = constrain_as_space_group( + object_parameter_set, space_group + ) + list(symbol_parameters) + list(object_parameters) + actual_names = symbol_parameters.names + expected_names = object_parameters.names + assert actual_names == expected_names + return + + +def makeLaMnO3_P1(datafile): + from diffpy.structure import Structure + + structure = Structure() + structure.read(datafile("LaMnO3.stru")) + return structure + + +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