diff --git a/docs/source/tutorial/basics.rst b/docs/source/tutorial/basics.rst
index e63ffc5f75a51eccc383905103e98b0c44a695c7..54bc94206b4f8174b71f5e7d518306010d43146b 100644
--- a/docs/source/tutorial/basics.rst
+++ b/docs/source/tutorial/basics.rst
@@ -114,6 +114,12 @@ Below are a few examples of the sorts of things you can define with Qsymm:
 The documentation page of the `qsymm.groups` module contains an exhaustive list
 of what can be generated.
 
+Alternatively, to define a point group element we initialize it directly by providing its real and orbital space action:
+
+.. jupyter-execute::
+
+    qsymm.PointGroupElement(R=np.diag([1, -1]), U='kron(sigma_y, sigma_z)')
+
 As with other Qsymm objects we can get a readable representation of these
 group elements:
 
diff --git a/qsymm/groups.py b/qsymm/groups.py
index b34de1315528a7963ee7defa87449b8d338150c9..5160d904a9ed31d3f20294d24b68cfc496c0d3f8 100644
--- a/qsymm/groups.py
+++ b/qsymm/groups.py
@@ -14,6 +14,7 @@ from sympy.matrices.matrices import MatrixBase
 
 from .linalg import prop_to_id, _inv_int, allclose
 from .model import Model
+from .kwant_continuum import sympify
 
 
 # Cache the operations that are potentially slow and happen a lot in
@@ -104,7 +105,7 @@ class PointGroupElement:
         Whether the operation includes conplex conjugation (antiunitary operator)
     antisymmetry : boolean (default False)
         Whether the operator flips the sign of the Hamiltonian (antisymmetry)
-    U : ndarray (optional)
+    U : array, str, SymPy expression, or None (default)
         The unitary action on the Hilbert space.
         May be None, to be able to treat symmetry candidates
     _strict_eq : boolean (default False)
@@ -112,6 +113,11 @@ class PointGroupElement:
         other PointGroupElements. By default the unitary parts are ignored.
         If True, PointGroupElements are considered equal, if the unitary parts
         are proportional, an overall phase difference is still allowed.
+    locals : dict or ``None`` (default)
+        Additional namespace entries for `~kwant_continuum.sympify`.  May be
+        used to simplify input of matrices or modify input before proceeding
+        further. For example:
+        ``locals={'sigma_plus': [[0, 2], [0, 0]]}``.
 
     Notes
     -----
@@ -132,7 +138,7 @@ class PointGroupElement:
 
     __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq')
 
-    def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False):
+    def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False, *, locals=None):
         if isinstance(R, sympy.ImmutableMatrix):
             # If it is integer, recast to integer tinyarray
             R = _make_int(R)
@@ -149,6 +155,16 @@ class PointGroupElement:
             R = _make_int(R)
         else:
             raise ValueError('Real space rotation must be provided as a sympy matrix or an array.')
+        # Normalize U
+        if U is None:
+            pass
+        else:
+            try:
+                U = np.atleast_2d(np.array(U, dtype=complex))
+            except (ValueError, TypeError):
+                U = sympify(U, locals=locals)
+                U = np.atleast_2d(np.array(U, dtype=complex))
+
         self.R, self.conjugate, self.antisymmetry, self.U = R, conjugate, antisymmetry, U
         # Calculating sympy inverse is slow, remember it
         self._strict_eq = _strict_eq
@@ -235,9 +251,9 @@ class PointGroupElement:
         if U is None:
             Uinv = None
         elif c:
-            Uinv = la.inv(U).conj()
+            Uinv = U.T
         else:
-            Uinv = la.inv(U)
+            Uinv = U.T.conj()
         # Check if inverse is stored, if not, calculate it
         Rinv = _inv(R)
         result = PointGroupElement(Rinv, c, a, Uinv, _strict_eq=self._strict_eq)
diff --git a/qsymm/model.py b/qsymm/model.py
index 21dbd615484939e1550a4edd83b525a0425cb759..48ad102a8ebdb5ed32b919a71ee454549a975d78 100644
--- a/qsymm/model.py
+++ b/qsymm/model.py
@@ -196,7 +196,7 @@ class Model(UserDict):
         else:
             self.keep = set()
 
-        if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping):
+        if isinstance(hamiltonian, abc.Mapping):
             # Initialize as dict sympifying the keys
             self.data = {symbol_normalizer(k): v for k, v in hamiltonian.items()
                               if symbol_normalizer(k) in self.keep
diff --git a/qsymm/tests/test_groups.py b/qsymm/tests/test_groups.py
new file mode 100644
index 0000000000000000000000000000000000000000..55032f986d9ac8d8a6f29d169820833f6ca53553
--- /dev/null
+++ b/qsymm/tests/test_groups.py
@@ -0,0 +1,23 @@
+import pytest
+import sympy
+import numpy as np
+import tinyarray as ta
+
+from ..groups import PointGroupElement, Model, time_reversal, chiral, rotation
+
+
+@pytest.mark.parametrize(
+    "U,U_expected",
+    [
+        (1, np.eye(1)),
+        ("1", np.eye(1)),
+        ("kron(sigma_0, sigma_0)", np.eye(4)),
+        (ta.array([[0, 1], [1, 0]]), 1 - np.eye(2))
+    ]
+)
+def test_U_normalization(U, U_expected):
+    R = np.eye(2)
+    
+    element = PointGroupElement(R, conjugate=False, antisymmetry=False, U=U)
+    np.testing.assert_equal(element.R, R)
+    np.testing.assert_equal(element.U, U_expected)
\ No newline at end of file