✨ New feature: set to get displacement_get_cartesian_eigenvalues

picometer/atom.py

@@ -185,6 +185,21 @@ def transform(self, symm_op_code: str) -> 'AtomSet':
                 data.loc[mask, 'U23'] = uij_rot[:, 1, 2]
         return self.__class__(self.base, data)
 
+    @property
+    def u_cartesian_eigenvalues(self):
+        u_columns = ['U11', 'U12', 'U13', 'U22', 'U23', 'U33']
+        eigenvalues = np.full((len(self), 3), np.nan)
+        if not set(u_columns).issubset(self.table.keys()):
+            return eigenvalues
+        u_fract = self.fract_uij  # Nx3n3 stack of abc-normalized U_cif tensors
+        a_mat = self.base.A_d.T  # fractional to orthogonal cartesian metric
+        n_mat = np.diag([self.base.a_r, self.base.b_r, self.base.c_r])
+        u_star = (n_mat @ u_fract) @ n_mat  # eq. 4b @ S0021889802008580
+        u_cart = (a_mat @ u_star) @ a_mat.T  # eq. 3a @ S0021889802008580
+        mask = ~self.table[u_columns].isna().any(axis=1)
+        eigenvalues[mask, :] = np.linalg.eigh(u_cart[mask, :, :])[0]
+        return eigenvalues
+
     @property
     def centroid(self) -> np.ndarray:
         """A 3-vector with average atom position."""

picometer/instructions.py

@@ -225,7 +225,8 @@ def _load_model_state(self, cif_path, block_name):
                 if 'Uiso' not in atoms.table.columns:
                     atoms.table['Uiso'] = pd.NA
                 u_equiv = atoms.table[['U11', 'U22', 'U33']].mean(axis=1)
-                atoms.table['Uiso'].fillna(u_equiv, inplace=True)
+                mask = atoms.table['Uiso'].isna()
+                atoms.table.loc[mask, 'Uiso'] = u_equiv[mask]
 
         if self.processor.settings['complete_umatrix_from_uiso']:
             u_columns = ['U11', 'U12', 'U13', 'U22', 'U23', 'U33']
@@ -359,6 +360,11 @@ def handle_one(self, instruction: Instruction, ms_key: str, ms: ModelState) -> N
                 value = getattr(displacements, suffix, None)
                 if value is not None:
                     self.processor.evaluation_table.loc[ms_key, label_] = value
+        if self.processor.settings['displacement_get_cartesian_eigenvalues']:
+            for label, es in zip(focus.table.index, focus.u_cartesian_eigenvalues):
+                for e, suffix in zip(es, ['Uce1', 'Uce2', 'Uce3']):
+                    label_ = label + '_' + suffix
+                    self.processor.evaluation_table.loc[ms_key, label_] = e
         logger.info(f'Noted displacement for current selection in model state {ms_key}')
 
 

picometer/settings.py

@@ -19,6 +19,7 @@ class DefaultSettings:
     """Store default values of all settings. Use `AnyValue` if no default."""
     auto_write_unit_cell: bool = True
     clear_selection_after_use: bool = True
+    displacement_get_cartesian_eigenvalues: bool = False
     complete_uiso_from_umatrix: bool = False
     complete_umatrix_from_uiso: bool = False
 

picometer/settings.yaml

@@ -1,5 +1,6 @@
 settings:
   auto_write_unit_cell: True
   clear_selection_after_use: True
+  displacement_get_cartesian_eigenvalues: False
   complete_uiso_from_umatrix: False
   complete_umatrix_from_uiso: False

tests/test_instructions.py

@@ -316,7 +316,20 @@ def test_displacement_complete_umatrix_from_Uiso(self):
         p = process(Routine.from_string(r))
         results = p.evaluation_table['C(11)_U13'].to_numpy()
         self.assertAlmostEqual(results[0], 0.010286, places=6)
-        self.assertAlmostEqual(results[0], 0.010286, places=6)
+        self.assertAlmostEqual(results[1], 0.010286, places=6)
+        np.testing.assert_equal(results[2], np.nan)
+        np.testing.assert_equal(results[3], np.nan)
+        np.testing.assert_equal(results[4], np.nan)
+        np.testing.assert_equal(results[5], np.nan)
+
+    def test_displacement_displacement_get_cartesian_eigenvalues(self):
+        r = ('settings:  \n  complete_umatrix_from_uiso: True\n' +
+             '  displacement_get_cartesian_eigenvalues: True\n' + self.routine_text)
+        r += '  - select: C(11)\n  - displacement\n'
+        p = process(Routine.from_string(r))
+        results = p.evaluation_table['C(11)_Uce1'].to_numpy()
+        self.assertAlmostEqual(results[0], 0.020000, places=6)
+        self.assertAlmostEqual(results[1], 0.019980, places=6)
         np.testing.assert_equal(results[2], np.nan)
         np.testing.assert_equal(results[3], np.nan)
         np.testing.assert_equal(results[4], np.nan)