diff --git a/picometer/atom.py b/picometer/atom.py index da7ece1..88ccaa1 100644 --- a/picometer/atom.py +++ b/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.""" diff --git a/picometer/instructions.py b/picometer/instructions.py index a337974..850c2a4 100644 --- a/picometer/instructions.py +++ b/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}') diff --git a/picometer/settings.py b/picometer/settings.py index bfb2d83..36398af 100644 --- a/picometer/settings.py +++ b/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 diff --git a/picometer/settings.yaml b/picometer/settings.yaml index a2aa79b..b9dcf90 100644 --- a/picometer/settings.yaml +++ b/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 diff --git a/tests/test_instructions.py b/tests/test_instructions.py index 73687b1..eaea075 100644 --- a/tests/test_instructions.py +++ b/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)