From bba2708eb3d6b88f77af8ffdc1826ac935de3411 Mon Sep 17 00:00:00 2001 From: James Ball Date: Fri, 30 Oct 2026 15:47:26 +0100 Subject: [PATCH] unitcell: read cell_sg, find R axes, warn when the cell lacks the symmetry - unitcell(..., sgname=None): an xfab space group name, cell_sg in a .par file. It overrides symmetry for the space group and the hkls, so a name like R-3mr (rhombohedral axes) is used for the rings, not only by transformer. A lattice letter or number that does not match it is logged. unitcell_from_parameters reads cell_sg, to_par_dict writes it. - With a space group number, an R group whose cell is on rhombohedral axes now uses them. Before, 166 on a 5,5,5,70,70,70 cell gave 350 reflections with d* < 1 instead of 438. - With any space group, a warning is logged when the Laue rotations change the cell's metric (e.g. P21/c with unique axis c). - tensor_map.sym_ops_from_phase uses unitcell.cell_choice. Co-Authored-By: Claude Opus 5.5 --- ImageD11/sinograms/tensor_map.py | 12 ++--- ImageD11/unitcell.py | 74 +++++++++++++++++++++++++-- test/test_cell_sg.py | 88 ++++++++++++++++++++++++++++++++ 3 files changed, 161 insertions(+), 13 deletions(-) create mode 100644 test/test_cell_sg.py diff --git a/ImageD11/sinograms/tensor_map.py b/ImageD11/sinograms/tensor_map.py index a305dbd2..343f46a4 100644 --- a/ImageD11/sinograms/tensor_map.py +++ b/ImageD11/sinograms/tensor_map.py @@ -558,9 +558,8 @@ def sym_ops_from_phase(ucell): to U directly. The cell must be in the standard setting of its spacegroup (monoclinic: - unique axis b). R spacegroups may be on hexagonal or rhombohedral axes; - a=b=c, alpha=beta=gamma != 90 is taken as rhombohedral. Pass sym_ops= - explicitly for any other setting. + unique axis b). R spacegroups may be on hexagonal or rhombohedral axes, + as ucell.cell_choice says. Pass sym_ops= explicitly for any other setting. """ try: spacegroup = ucell.spacegroup @@ -571,12 +570,7 @@ def sym_ops_from_phase(ucell): "Bravais letter. Pass sym_ops= explicitly instead." ) - a, b, c, alpha, beta, gamma = ucell.lattice_parameters - rhombohedral_axes = (np.isclose(a, b) and np.isclose(b, c) - and np.isclose(alpha, beta) and np.isclose(beta, gamma) - and not np.isclose(alpha, 90.0)) - cell_choice = 'rhombohedral' if rhombohedral_axes else 'standard' - ops = sym_u.cartesian(sym_u.laue_rotations(spacegroup, cell_choice), ucell.B) + ops = sym_u.cartesian(sym_u.laue_rotations(spacegroup, ucell.cell_choice), ucell.B) if not np.allclose(np.matmul(ops, np.transpose(ops, (0, 2, 1))), np.eye(3), atol=1e-6): raise ValueError( "Cell %s is not in the standard setting of spacegroup %s. " diff --git a/ImageD11/unitcell.py b/ImageD11/unitcell.py index e46d99e4..f908066a 100644 --- a/ImageD11/unitcell.py +++ b/ImageD11/unitcell.py @@ -206,15 +206,35 @@ def get_unitcells(self): # set the name of the unitcell from the phase name self.unitcells[phase_name].name = phase_name +# R space groups can be on hexagonal (standard) or rhombohedral axes +RHOMBOHEDRAL_SPACEGROUPS = (146, 148, 155, 160, 161, 166, 167) +# relative change of the metric tensor accepted as symmetric +METRIC_TOL = 1e-5 + + +def laue_misfit(g, spacegroup, cell_choice="standard"): + """ + Largest change of the real space metric tensor g under the Laue class + rotations of spacegroup, relative to g. Zero when the cell has the symmetry. + """ + from ImageD11 import sym_u + ops = sym_u.laue_rotations(spacegroup, cell_choice).group + return max(np.abs(np.dot(np.dot(o, g), o.T) - g).max() for o in ops) / np.abs(g).max() + + class unitcell: # Unit cell stuff # Generate a list of peaks from a unit cell - def __init__(self, lattice_parameters, symmetry="P", verbose=0, name=None): + def __init__(self, lattice_parameters, symmetry="P", verbose=0, name=None, sgname=None): """ Unit cell class supply a list (tuple etc) of a,b,c,alpha,beta,gamma optionally a lattice symmetry, one of "P","A","B","C","I","F","R" or : a space group number (integer) + optionally sgname : an xfab space group name (cell_sg in a .par file) + e.g. "Fm-3m" or "R-3mr" (rhombohedral axes). It overrides symmetry + for the space group and the hkls. + With a space group, a warning is logged if the cell does not have its symmetry. """ self.lattice_parameters = np.array(lattice_parameters) if self.lattice_parameters.shape[0] != 6: @@ -299,7 +319,41 @@ def __init__(self, lattice_parameters, symmetry="P", verbose=0, name=None): # orix stuff self._orix_phase = None - + + self.sgname = sgname + self.cell_choice = "standard" + if sgname is not None: + self._set_sgname(sgname) + elif isinstance(self.symmetry, int): + if (self.spacegroup in RHOMBOHEDRAL_SPACEGROUPS and + laue_misfit(self.g, self.spacegroup) > METRIC_TOL and + laue_misfit(self.g, self.spacegroup, "rhombohedral") <= METRIC_TOL): + self.cell_choice = "rhombohedral" + if hasattr(self, "spacegroup"): + misfit = laue_misfit(self.g, self.spacegroup, self.cell_choice) + if misfit > METRIC_TOL: + logging.warning("Unit cell %s does not have the symmetry of space group %s (%s axes): " + "the metric changes by %.2g under its rotations. Check the lattice " + "parameters, or the setting (monoclinic: unique axis b)." % ( + str(self.lattice_parameters), str(self.spacegroup), + self.cell_choice, misfit)) + + def _set_sgname(self, sgname): + """ Space group from an xfab name, checked against self.symmetry """ + from xfab import sg + spg = sg.sg(sgname=sgname) + self.cell_choice = spg.cell_choice + if isinstance(self.symmetry, int): + if self.symmetry != spg.no: + logging.warning("Space group %d does not match %s (number %d): using %s" % ( + self.symmetry, sgname, spg.no, sgname)) + else: + centring = "P" if spg.cell_choice == "rhombohedral" else spg.name[0] + if self.symmetry != centring: + logging.warning("Lattice %s does not match %s (lattice %s): using %s" % ( + self.symmetry, sgname, centring, sgname)) + self.spacegroup = spg.no + def __repr__(self): if (not hasattr(self, 'name')) or (self.name is None): return "Unitcell" + " | " + str(self.lattice_parameters) + " | " + str(self.symmetry) @@ -463,10 +517,13 @@ def gethkls_xfab(self, dsmax, spg=None): sgnum = self.symmetry else: sgnum = None + # a name carries its own axes (e.g. 'R-3mr'), a number uses those found for the cell + cell_choice = self.cell_choice if spg is None else "standard" raw_peaks = tools.genhkl_all(self.lattice_parameters, 0, stl_max, sgname=spg, sgno=sgnum, + cell_choice=cell_choice, output_stl=True) peaks = [] for i in range(len(raw_peaks)): @@ -488,6 +545,8 @@ def gethkls(self, dsmax): """ if dsmax == self.limit and self.peaks is not None: return self.peaks + if self.sgname is not None: + return self.gethkls_xfab(dsmax, spg=self.sgname) if isinstance(self.symmetry, int): return self.gethkls_xfab(dsmax, spg=None) h = k = 0 @@ -764,7 +823,9 @@ def to_par_dict(self): pars_dict[parname] = self.lattice_parameters[inc] else: pars_dict[parname] = self.symmetry - + if self.sgname is not None: + pars_dict['cell_sg'] = self.sgname + if hasattr(self, 'name') and self.name is not None: pars_dict['phase_name'] = self.name return pars_dict @@ -896,8 +957,13 @@ def ubi_equiv(ubilist, ublist, tol=1e-8): def unitcell_from_parameters(pars): parnames = "_a _b _c alpha beta gamma".split() + try: + sgname = pars.get("cell_sg") # optional + except KeyError: + sgname = None cell = unitcell([pars.get("cell_%s" % (s)) for s in parnames], - pars.get("cell_lattice_[P,A,B,C,I,F,R]")) + pars.get("cell_lattice_[P,A,B,C,I,F,R]"), + sgname=sgname) return cell diff --git a/test/test_cell_sg.py b/test/test_cell_sg.py new file mode 100644 index 00000000..330ff843 --- /dev/null +++ b/test/test_cell_sg.py @@ -0,0 +1,88 @@ + +import logging +import unittest +from ImageD11 import unitcell, parameters + +RHOMB = [5., 5., 5., 70., 70., 70.] +HEX = [4.9, 4.9, 13., 90., 90., 120.] +CUBIC = [4.05, 4.05, 4.05, 90., 90., 90.] + + +class Catch(logging.Handler): + def __init__(self): + logging.Handler.__init__(self, logging.WARNING) + self.messages = [] + + def emit(self, record): + self.messages.append(record.getMessage()) + + +def warnings_from(*args, **kwds): + """ unitcell(*args, **kwds) and the warnings it logged """ + catch = Catch() + logging.getLogger().addHandler(catch) + try: + cell = unitcell.unitcell(*args, **kwds) + finally: + logging.getLogger().removeHandler(catch) + return cell, catch.messages + + +def hkls(cell, dsmax=1.0): + return set(h for d, h in cell.gethkls(dsmax)) + + +class test_cell_sg(unittest.TestCase): + + def test_rhombohedral_axes_from_number(self): + # 166 on rhombohedral axes : no centring absences + cell, w = warnings_from(RHOMB, 166) + self.assertEqual(cell.cell_choice, "rhombohedral") + self.assertEqual(w, []) + self.assertEqual(hkls(cell), hkls(unitcell.unitcell(RHOMB, "P", sgname="R-3mr"))) + + def test_hexagonal_axes_from_number(self): + cell, w = warnings_from(HEX, 166) + self.assertEqual(cell.cell_choice, "standard") + self.assertEqual(w, []) + self.assertEqual(hkls(cell), hkls(unitcell.unitcell(HEX, "R"))) + + def test_more_specific_name(self): + cell, w = warnings_from(CUBIC, "F", sgname="Fm-3m") + self.assertEqual(w, []) + self.assertEqual(cell.spacegroup, 225) + self.assertEqual(cell.symmetry, "F") + self.assertEqual(hkls(cell), hkls(unitcell.unitcell(CUBIC, 225))) + + def test_incompatible_name(self): + for symmetry in ("I", 229): + cell, w = warnings_from(CUBIC, symmetry, sgname="Fm-3m") + self.assertEqual(len(w), 1, symmetry) + self.assertEqual(cell.spacegroup, 225) # the name wins + self.assertEqual(hkls(cell), hkls(unitcell.unitcell(CUBIC, 225))) + + def test_cell_without_the_symmetry(self): + for cell, sg in [([3., 4., 5., 90., 90., 100.], 14), # unique axis c + ([4.05, 4.06, 4.05, 90., 90., 90.], 225), + ([4.9, 4.9, 13., 90., 90., 90.], 166)]: + self.assertEqual(len(warnings_from(cell, sg)[1]), 1, sg) + for cell, sg in [([3., 4., 5., 90., 100., 90.], 14), (CUBIC, 225), (HEX, 194)]: + self.assertEqual(warnings_from(cell, sg)[1], [], sg) + self.assertEqual(warnings_from([3., 4., 5., 90., 90., 100.], "P")[1], []) + + def test_par_round_trip(self): + pars = unitcell.unitcell(RHOMB, "P", sgname="R-3mr").to_par_dict() + self.assertEqual(pars["cell_sg"], "R-3mr") + for p in (pars, parameters.parameters(**pars)): + cell = unitcell.unitcell_from_parameters(p) + self.assertEqual(cell.sgname, "R-3mr") + self.assertEqual(cell.cell_choice, "rhombohedral") + # old files without cell_sg + pars = unitcell.unitcell(CUBIC, 225).to_par_dict() + self.assertNotIn("cell_sg", pars) + for p in (pars, parameters.parameters(**pars)): + self.assertEqual(unitcell.unitcell_from_parameters(p).sgname, None) + + +if __name__ == "__main__": + unittest.main()