✨ Add fractional coordinates and displacement parameter commands

README.md

@@ -120,6 +120,9 @@ The following instructions are currently supported by picometer:
   - fit `line` to the current atom / centroid selection;
   - fit `plane` to the currect atom / centroid selection;
 - **Evaluation instructions**
+  - write out fractional `coordinates` of currently selected centroids or atoms.
+  - write out `displacement` parameters of currently selected centroids or atoms
+    (note: currently does not correctly handle symmetry transformations).
   - measure `distance` between 2 selected objects; if the selection includes
     groups of atoms, measure closes distance to the group of atoms.
   - measure `angle` between 2–3 selected objects: planes, lines, or (ordered) atoms.

picometer/atom.py

@@ -2,17 +2,22 @@
 import logging
 from typing import Dict, NamedTuple, List, Sequence
 
-import hikari.symmetry
 from hikari.dataframes import BaseFrame, CifFrame
-import numpy as np
 from numpy.linalg import norm
+import numpy as np
 import pandas as pd
 
 from picometer.shapes import (are_synparallel, degrees_between, Line,
                               Plane, Shape, Vector3)
 from picometer.utility import ustr2float
 
 
+try:
+    from hikari.symmetry import Operation
+except ImportError:  # hikari version < 0.3.0
+    from hikari.symmetry import SymmOp as Operation
+
+
 logger = logging.getLogger(__name__)
 
 
@@ -85,6 +90,19 @@
                 'fract_y': [ustr2float(v) for v in cb['_atom_site_fract_y']],
                 'fract_z': [ustr2float(v) for v in cb['_atom_site_fract_z']],
             }
+            try:
+                atoms_dict['Uiso'] = [ustr2float(v) for v in cb['_atom_site_U_iso_or_equiv']]
+            except KeyError:
+                pass
+            try:
+                atoms_dict['U11'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_11']]
+                atoms_dict['U22'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_22']]
+                atoms_dict['U33'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_33']]
+                atoms_dict['U23'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_23']]
+                atoms_dict['U13'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_13']]
+                atoms_dict['U12'] = [ustr2float(v) for v in cb['_atom_site_aniso_U_12']]
+            except KeyError:
+                pass
             atoms = pd.DataFrame.from_records(atoms_dict).set_index('label')
         except KeyError:
             atoms = pd.DataFrame()
@@ -131,7 +149,7 @@
         return self.__class__(self.base, deepcopy(self.table[mask]))
 
     def transform(self, symm_op_code: str) -> 'AtomSet':
-        symm_op = hikari.symmetry.SymmOp.from_code(symm_op_code)
+        symm_op = Operation.from_code(symm_op_code)
         fract_xyz = symm_op.transform(self.fract_xyz.T)
         data = deepcopy(self.table)
         data['fract_x'] = fract_xyz[:, 0]

picometer/instructions.py

@@ -307,6 +307,34 @@ def handle_one(self, instruction: Instruction, ms_key: str, ms: ModelState) -> N
         logger.info(f'Defined plane {label}: {plane} for model state {ms_key}')
 
 
+class CoordinatesInstructionHandler(SerialInstructionHandler):
+    name = 'coordinates'
+    kwargs = None
+
+    def handle_one(self, instruction: Instruction, ms_key: str, ms: ModelState) -> None:
+        focus = ms.nodes.locate(self.processor.selection)
+        for label, coords in focus.table.iterrows():
+            self.processor.evaluation_table.loc[ms_key, label + '_x'] = coords.fract_x
+            self.processor.evaluation_table.loc[ms_key, label + '_y'] = coords.fract_y
+            self.processor.evaluation_table.loc[ms_key, label + '_z'] = coords.fract_z
+        logger.info(f'Noted coordinates for current selection in model state {ms_key}')
+
+
+class DisplacementInstructionHandler(SerialInstructionHandler):
+    name = 'displacement'
+    kwargs = None
+
+    def handle_one(self, instruction: Instruction, ms_key: str, ms: ModelState) -> None:
+        focus = ms.nodes.locate(self.processor.selection)
+        for label, displacements in focus.table.iterrows():
+            for suffix in 'Uiso U11 U22 U33 U23 U13 U12'.split():
+                label_ = label + '_' + suffix
+                value = getattr(displacements, suffix, None)
+                if value is not None:
+                    self.processor.evaluation_table.loc[ms_key, label_] = value
+        logger.info(f'Noted displacement for current selection in model state {ms_key}')
+
+
 class DistanceInstructionHandler(SerialInstructionHandler):
     name = 'distance'
     kwargs = dict(label=str)

tests/ferrocene1.cif

@@ -52,14 +52,15 @@ _atom_site_label
 _atom_site_fract_x
 _atom_site_fract_y
 _atom_site_fract_z
-Fe .0 .0 .0
-C(11) .0227(4) .2623(4) -.0218(10)
-C(12) .0395(3) .1790(5) -.2170(6)
-C(13) .1593(4) .0635(4) -.0900(8)
-C(14) .2165(4) .0753(5) .1836(7)
-C(15) .1322(5) .1982(5) .2258(7)
-H(11) -.0624 .3567 -.0567
-H(12) -.0304 .1975 -.4297
-H(13) .1985 -.0233 -.1871
-H(14) .3079 -.0006 .3359
-H(15) .1467 .2342 .4165
+_atom_site_U_iso_or_equiv
+Fe .0 .0 .0 0.01
+C(11) .0227(4) .2623(4) -.0218(10) .02
+C(12) .0395(3) .1790(5) -.2170(6) .02
+C(13) .1593(4) .0635(4) -.0900(8) .02
+C(14) .2165(4) .0753(5) .1836(7) .02
+C(15) .1322(5) .1982(5) .2258(7) .02
+H(11) -.0624 .3567 -.0567 .03
+H(12) -.0304 .1975 -.4297 .03
+H(13) .1985 -.0233 -.1871 .03
+H(14) .3079 -.0006 .3359 .03
+H(15) .1467 .2342 .4165 .03

tests/ferrocene2.cif

@@ -63,3 +63,22 @@ H(12) -.0236 .1961 -.4072
 H(13) .1983 -.0170 -.1710
 H(14) .3019 .0053 .3297
 H(15) .1440 .2321 .4048
+loop_
+_atom_site_aniso_label
+_atom_site_aniso_U_11
+_atom_site_aniso_U_22
+_atom_site_aniso_U_33
+_atom_site_aniso_U_23
+_atom_site_aniso_U_13
+_atom_site_aniso_U_12
+Fe .01 .01 .01 .0 .0 .0
+C(11) .02 .02 .02 .0 .0 .0
+C(12) .02 .02 .02 .0 .0 .0
+C(13) .02 .02 .02 .0 .0 .0
+C(14) .02 .02 .02 .0 .0 .0
+C(15) .02 .02 .02 .0 .0 .0
+H(11) .03 .03 .03 .0 .0 .0
+H(12) .03 .03 .03 .0 .0 .0
+H(13) .03 .03 .03 .0 .0 .0
+H(14) .03 .03 .03 .0 .0 .0
+H(15) .03 .03 .03 .0 .0 .0

tests/test_instructions.py

@@ -7,6 +7,7 @@
 import unittest
 
 import numpy as np
+import pandas
 import pandas as pd
 from pandas.testing import assert_frame_equal
 
@@ -175,7 +176,7 @@ def test_transformed_group(self):
         for _, ms in p.model_states.items():
             carbons_a = ms.atoms.locate([Locator('cp_A')]).table
             carbons_b = ms.atoms.locate([Locator('cp_B')]).table
-            for key in carbons_a.keys():
+            for key in [k for k in carbons_a.keys() if str(k).startswith('fract')]:
                 self.assertEqual(carbons_a[key].iloc[0],
                                  -carbons_b[key].iloc[0])
 
@@ -263,6 +264,41 @@ class TestMeasuringInstructions(unittest.TestCase):
     def setUp(self) -> None:
         self.routine_text = self.routine_prefix
 
+    def test_coordinates(self):
+        self.routine_text += '  - select: cp_A\n'
+        self.routine_text += '  - coordinates\n'
+        p = process(Routine.from_string(self.routine_text))
+        results = p.evaluation_table['C(11)_y'].to_numpy()
+        correct = np.array([0.2623, 0.2612, 0.2662, 0.2622, 0.2624, 0.2615])
+        self.assertTrue(np.allclose(results, correct))
+        results = p.evaluation_table['C(21)_y'].to_numpy()
+        correct = np.array([0.2576, 0.2583, 0.2654, 0.258])
+        np.testing.assert_equal(results[0], np.nan)
+        self.assertTrue(np.allclose(results[2:], correct))
+
+    @unittest.expectedFailure
+    def test_coordinates2(self):
+        """Known failure: coordinates of atoms with the same name are overwritten"""
+        self.routine_text += '  - select: cp_A\n'
+        self.routine_text += '  - coordinates\n'
+        t1 = process(Routine.from_string(self.routine_text)).evaluation_table
+        self.routine_text += '  - select: cp_B\n'
+        self.routine_text += '  - coordinates\n'
+        t2 = process(Routine.from_string(self.routine_text)).evaluation_table
+        self.assertEqual(t1.shape, t2.shape)
+        self.assertTrue(t1.equals(t2[t1.keys()]))
+
+    def test_displacement(self):
+        self.routine_text += '  - select: cp_A\n'
+        self.routine_text += '  - displacement\n'
+        p = process(Routine.from_string(self.routine_text))
+        results = p.evaluation_table['C(11)_Uiso'].to_numpy()
+        self.assertEqual(results[0], 0.02)
+        np.testing.assert_equal(results[1], np.nan)
+        results = p.evaluation_table['C(11)_U11'].to_numpy()
+        self.assertEqual(results[1], 0.02)
+        np.testing.assert_equal(results[0], np.nan)
+
     def test_distance_plane_plane(self):
         self.routine_text += '  - select: cp_A_plane\n'
         self.routine_text += '  - select: cp_B_plane\n'