From 8d9bd30804b8d0c4d806bad86890b424b6c04930 Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Thu, 19 Oct 2023 14:45:21 -0400
Subject: [PATCH] add access for generating SymOp matrix representation

---
 python/src/xtal.cpp        | 40 ++++++++++++++++++++++++++++++++++++++
 python/tests/test_symop.py | 12 ++++++++++++
 2 files changed, 52 insertions(+)

diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp
index ce1aff1..9058871 100644
--- a/python/src/xtal.cpp
+++ b/python/src/xtal.cpp
@@ -2008,6 +2008,46 @@ PYBIND11_MODULE(_xtal, m) {
           )pbdoc")
       .def("matrix", &xtal::get_matrix,
            "Returns the transformation matrix value.")
+      .def(
+          "matrix_rep",
+          [](xtal::SymOp const &op, std::string key) -> Eigen::MatrixXd {
+            Eigen::MatrixXd M;
+            try {
+              AnisoValTraits traits(key);
+              M = traits.symop_to_matrix(get_matrix(op), get_translation(op),
+                                         get_time_reversal(op));
+            } catch (std::exception &e) {
+              std::stringstream msg;
+              msg << "Error getting matrix rep: CASM does not know how to "
+                     "transform the property '"
+                  << key << "'.";
+              throw std::runtime_error(msg.str());
+            }
+            return M;
+          },
+          R"pbdoc(
+          Returns the matrix representation of a symmetry operation for transforming \
+          properties
+
+          Parameters
+          ----------
+          key: str
+              The name of the CASM-supported property to be transformed.
+
+          Returns
+          -------
+          matrix_rep : numpy.ndarray[numpy.float64[m, m]]
+              The matrix representation for transforming properties. The matrix is
+              square, with dimension equal to the standard dimension of the specified
+              property. For example, `m=3` for `key="disp"`, and `m=6` for
+              `key="Hstrain"`. Local properties, such as `"disp"`, stored as columns of
+              array `local_values`, can then be transformed using
+              ``matrix_rep @ local_values``. Global properties, such as `"Hstrain"`,
+              stored as array `global_values` with a single column, can similarly be
+              transformed using  ``matrix_rep @ global_values``.
+
+          )pbdoc",
+          py::arg("key"))
       .def("translation", &xtal::get_translation,
            "Returns the translation value.")
       .def("time_reversal", &xtal::get_time_reversal,
diff --git a/python/tests/test_symop.py b/python/tests/test_symop.py
index 625a88b..1cf459b 100644
--- a/python/tests/test_symop.py
+++ b/python/tests/test_symop.py
@@ -127,6 +127,10 @@ def test_SymOp_mul_properties():
     assert "disp" in transformed_properties
     assert transformed_properties["disp"].shape == (3, 4)
 
+    matrix_rep = op.matrix_rep("disp")
+    transformed_disp = matrix_rep @ disp
+    assert np.allclose(transformed_disp, transformed_properties["disp"])
+
     # check 1d array - accepted, but returned as (6,1) array
     Hstrain = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6])
     assert Hstrain.shape == (6,)
@@ -138,6 +142,10 @@ def test_SymOp_mul_properties():
     assert "Hstrain" in transformed_properties
     assert transformed_properties["Hstrain"].shape == (6, 1)
 
+    matrix_rep = op.matrix_rep("Hstrain")
+    transformed_Hstrain = matrix_rep @ Hstrain.reshape(-1, 1)
+    assert np.allclose(transformed_Hstrain, transformed_properties["Hstrain"])
+
     # check 2d column array - accepted, stays as (6,1) array
     Hstrain = np.array([[0.1, 0.2, 0.3, 0.4, 0.5, 0.6]]).transpose()
     assert Hstrain.shape == (6, 1)
@@ -149,6 +157,10 @@ def test_SymOp_mul_properties():
     assert "Hstrain" in transformed_properties
     assert transformed_properties["Hstrain"].shape == (6, 1)
 
+    matrix_rep = op.matrix_rep("Hstrain")
+    transformed_Hstrain = matrix_rep @ Hstrain
+    assert np.allclose(transformed_Hstrain, transformed_properties["Hstrain"])
+
 
 def test_SymOp_mul_lattice(tetragonal_lattice):
     R = np.array(