unitcell: read cell_sg, find R axes, warn when the cell lacks the symmetry

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. "

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
 

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()