Adding improper rotations (#5)

* Adding methane to test suite

* Adding ammonia and adding matrix generation for improper rotations

* Adding improper rotations to tests

* Expanding documentation

* Expanding documentation

docs/background.rst

@@ -95,11 +95,11 @@ the tesseral spherical harmonics are associated to the values for :math:`m`.
      - :math:`-3`
      - :math:`-2`
      - :math:`-1`
-     - :math:`-0`
-     - :math:`-1`
-     - :math:`-2`
-     - :math:`-3`
-     - :math:`-4`
+     - :math:`0`
+     - :math:`1`
+     - :math:`2`
+     - :math:`3`
+     - :math:`4`
    * - :math:`-`
      - :math:`-`
      - :math:`-`

docs/wignerd.rst

@@ -9,5 +9,5 @@ wignerd module
 .. _wigner_d:
 
 .. automodule:: sphecerix.wignerd
-    :members: wigner_D, tesseral_wigner_D, tesseral_wigner_D_mirror
+    :members: wigner_D, tesseral_wigner_D, tesseral_wigner_D_mirror, tesseral_wigner_D_improper
     :show-inheritance:

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "sphecerix"
-  version: "0.3.0"
+  version: "0.4.0"
 
 source:
   path: .

sphecerix/__init__.py

@@ -1,4 +1,5 @@
-from .wignerd import tesseral_wigner_D, wigner_D, tesseral_wigner_D_mirror
+from .wignerd import tesseral_wigner_D, wigner_D, tesseral_wigner_D_mirror,\
+                     tesseral_wigner_D_improper
 from .tesseral import tesseral_transformation, permutation_sh_car
 from .atomic_wave_functions import wfcart, wf, wffield, wffield_l
 

sphecerix/_version.py

@@ -1 +1 @@
-__version__ = "0.3.0"
+__version__ = "0.4.0"

sphecerix/wignerd.py

@@ -8,7 +8,7 @@
 
 def tesseral_wigner_D(l, Robj):
     """
-    Produce the Wigner D-matrix for tesseral spherical harmonics
+    Produce the Wigner D-matrix for tesseral spherical harmonics for a rotation
 
     Parameters
     ----------
@@ -123,6 +123,68 @@ def tesseral_wigner_D_mirror(l, normal):
 
     return inv * np.real(T @ D @ T.conjugate().transpose())
 
+def tesseral_wigner_D_improper(l, Robj):
+    """
+    Produce the Wigner D-matrix for tesseral spherical harmonics under an
+    improper rotation
+
+    Parameters
+    ----------
+    l : int
+        Order of the spherical harmonics
+    Robj : scipy.spatial.transform.Rotation
+        Rotation in :math:`\mathbb{R}^{3}`
+
+    Returns
+    -------
+    D : numpy.ndarray
+        Real-valued Wigner-D matrix with dimensions :math:`(2l+1) \\times (2l+1)`
+
+    Raises
+    ------
+    TypeError
+        If the Robj object is not of type scipy.spatial.transform.R.
+        
+    Examples
+    --------
+    >>> from sphecerix import tesseral_wigner_D_improper
+    ... from scipy.spatial.transform import Rotation as R
+    ... import numpy as np
+    ... 
+    ... # construct (improper) rotation vector
+    ... axis = np.array([1,0,0])
+    ... Robj = R.from_rotvec(axis * np.pi / 2)
+    ... 
+    ... # construct wigner D matrix
+    ... D = tesseral_wigner_D_improper(1, Robj)
+    ... 
+    ... print(D)
+    [[ 7.49879891e-33 -1.00000000e+00  1.83697020e-16]
+     [ 1.00000000e+00 -2.24963967e-32 -1.83697020e-16]
+     [-1.83697020e-16 -1.83697020e-16 -1.00000000e+00]]
+
+    Construct the Wigner-D matrix for the tesseral p-orbitals for an improper
+    rotation by 90 degrees around the cartesian x-axis.
+    """
+    # verify that Robj is a rotation object
+    if not isinstance(Robj, R):
+        raise TypeError('Second argument Robj should be of type scipy.spatial.transform.R')
+
+    T = tesseral_transformation(l)
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', r'Gimbal lock detected. Setting third angle to zero since it is not possible to uniquely determine all angles.')
+        alpha, beta, gamma = Robj.as_euler('zyz', degrees=False)
+
+    D = wigner_D(l, Robj)
+
+    # extract rotation vector from rotation object and use it to construct
+    # a mirror matrix
+    normal = Robj.as_rotvec()
+    normal /= np.linalg.norm(normal)
+    M = tesseral_wigner_D_mirror(l, normal)
+
+    return M @ np.real(T @ D @ T.conjugate().transpose())
+
 def wigner_D(l, Robj):
     """
     Produce Wigner D-matrix for canonical spherical harmonics

tests/test_ammonia.py

@@ -0,0 +1,99 @@
+import unittest
+import numpy as np
+import sys
+import os
+
+# add a reference to load the Sphecerix library
+sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
+
+# import functions
+from sphecerix import tesseral_wigner_D, tesseral_wigner_D_mirror
+from scipy.spatial.transform import Rotation as R
+
+class TestAmmonia(unittest.TestCase):
+    """
+    Test the construction of tesseral transformation matrices for the ammonia
+    molecule (C3v symmetry)
+    """
+
+    def test_rotation_c3(self):
+        """
+        The p-orbitals on N in ammonia belong to the A1 (z) and E (x,y) groups, 
+        thus the character under C3 should be 0.
+        """
+        axes = [[0,0,1], [0,0,-1]]
+        angle = np.radians(120)  # C3 rotation
+
+        for axis in axes:
+            axis /= np.linalg.norm(axis)
+            Robj = R.from_rotvec(np.array(axis) * angle)
+            T = tesseral_wigner_D(1, Robj)
+
+            # test that determinant of the transformation matrix is 1.0
+            # because the operation contains *no* reflection
+            np.testing.assert_almost_equal(np.linalg.det(T), 1.0)
+
+            # test the the matrix multiplied with its complex conjugate
+            # equals the identity matrix
+            np.testing.assert_almost_equal(T @ T.transpose(), np.identity(3))
+
+            # test that trace equals the character under C3
+            np.testing.assert_almost_equal(np.trace(T), 0)
+
+            # assert that z-component gives a 1
+            np.testing.assert_almost_equal(T[1,1], 1)
+
+            # assert that x,y-components give -1
+            np.testing.assert_almost_equal(T[0,0] + T[2,2], -1)
+
+    def test_mirror_sigma_v(self):
+        """
+        The p-orbitals on N in ammonia belong to the A1 (z) and E (x,y) groups, 
+        thus the character under sigma_v should be +1.
+        """
+        # construct the normal vectors for the sigma_v mirror planes. Given the
+        # orientation of the molecule, one normal vector of a mirror plane lies
+        # alongside the +x axis. The other normal vectors can be found using
+        # 120 degree rotations of that vector
+        normals = []
+        for i in range(0,3):
+            normals.append([np.cos(2.0 / 3.0 * i * np.pi), 
+                            np.sin(2.0 / 3.0 * i * np.pi), 
+                            0.0])
+
+        for normal in normals:
+            normal /= np.linalg.norm(normal)
+            M = tesseral_wigner_D_mirror(1, normal)
+
+            # test that determinant of the transformation matrix is -1.0
+            # because the operation contains a reflection
+            np.testing.assert_almost_equal(np.linalg.det(M), -1.0)
+
+            # test the the matrix multiplied with its complex conjugate
+            # equals the identity matrix
+            np.testing.assert_almost_equal(M @ M.transpose(), np.identity(3))
+
+            # test that trace equals the character under sigma_d
+            np.testing.assert_almost_equal(np.trace(M), 1)
+
+            # assert that z yields a 1 under sigma_d
+            np.testing.assert_almost_equal(M[1,1], 1)
+
+            # assert that x,y-components give 0 under sigma_d
+            np.testing.assert_almost_equal(M[0,0] + M[2,2], 0)
+
+    def build_coordinates_ammonia(self):
+        """
+        Construct fictitious ammonia coordinates
+        """
+        coords = [
+            [ 0.00000000,    -0.00000000,    -0.06931370],
+            [ 0.00000000,     0.94311105,     0.32106944],
+            [-0.81675813,    -0.47155553,     0.32106944],
+            [ 0.81675813,    -0.47155553,     0.32106944]
+        ]
+
+        return coords
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_improper_rotation.py

@@ -0,0 +1,39 @@
+import unittest
+import numpy as np
+import sys
+import os
+from scipy.spatial.transform import Rotation as R
+
+# add a reference to load the Sphecerix library
+sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
+
+# import functions
+from sphecerix import tesseral_wigner_D_improper
+
+class TestImproperRotation(unittest.TestCase):
+    """
+    Test improper rotation matrices
+    """
+
+    def test_improper_s4(self):
+        """
+        Test the result under an improper rotation
+        """
+        # construct rotation vector
+        axis = np.array([1,0,0])
+        Robj = R.from_rotvec(axis * np.pi / 2)
+
+        # construct wigner D matrix
+        D = tesseral_wigner_D_improper(1, Robj)
+
+        # assert that determinant is -1
+        np.testing.assert_almost_equal(np.linalg.det(D), -1)
+
+        # rotate a point at +x,+y,+z under S4 to -x,-y,+z
+        # note that the ordering in the vector (using increasing value of m)
+        # is [y,z,x]
+        np.testing.assert_almost_equal(D @ np.array([1,1,1]), np.array([-1,1,-1]))
+        print(D)
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_methane.py

@@ -0,0 +1,134 @@
+import unittest
+import numpy as np
+import sys
+import os
+
+# add a reference to load the Sphecerix library
+sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
+
+# import functions
+from sphecerix import tesseral_wigner_D, tesseral_wigner_D_mirror, \
+                      tesseral_wigner_D_improper
+from scipy.spatial.transform import Rotation as R
+
+class TestMethane(unittest.TestCase):
+    """
+    Test the construction of tesseral transformation matrices for the
+    p-orbitals on the C atom of methane (lying at the origin) 
+    
+    Methane has Td symmetry and the p-orbitals belong to the T2 symmetry group,
+    thus we expect the traces for the transformation matrices under
+    (C3, C2, S4, sigma_d) to correspond to (0,-1,-1,+1)
+    """
+
+    def test_rotation_c3(self):
+        """
+        Assert that the characters under C3 correspond to 0
+        """
+        coords = self.build_coordinates_methane()
+        angle = np.radians(120)  # C3 rotation
+
+        for axis in coords[1:]:
+            axis /= np.linalg.norm(axis)
+            Robj = R.from_rotvec(np.array(axis) * angle)
+            T = tesseral_wigner_D(1, Robj)
+
+            # test that determinant of the transformation matrix is 1.0
+            # because the operation contains *no* reflection
+            np.testing.assert_almost_equal(np.linalg.det(T), 1.0)
+
+            # test the the matrix multiplied with its complex conjugate
+            # equals the identity matrix
+            np.testing.assert_almost_equal(T @ T.transpose(), np.identity(3))
+
+            # test that trace equals the character under C3
+            np.testing.assert_almost_equal(np.trace(T), 0)
+
+    def test_rotation_c2(self):
+        """
+        Assert that the characters under C2 correspond to -1
+        """
+        axes = [
+            [1,0,0],
+            [0,1,0],
+            [0,0,1]
+        ]
+        angle = np.radians(180) # C2 rotation
+
+        for axis in axes:
+            axis /= np.linalg.norm(axis)
+            Robj = R.from_rotvec(np.array(axis) * angle)
+            T = tesseral_wigner_D(1, Robj)
+
+            # test that determinant of the transformation matrix is 1.0
+            # because the operation contains *no* reflection
+            np.testing.assert_almost_equal(np.linalg.det(T), 1.0)
+
+            # test the the matrix multiplied with its complex conjugate
+            # equals the identity matrix
+            np.testing.assert_almost_equal(T @ T.transpose(), np.identity(3))
+
+            # test that trace equals the character under C2
+            np.testing.assert_almost_equal(np.trace(T), -1)
+
+    def test_improper_rotation_s4(self):
+        """
+        Assert that the characters under S4 correspond to -1
+        """
+        coords = self.build_coordinates_methane()
+        angle = np.radians(90)  # S4 rotation
+
+        for axis in coords[1:]:
+            axis /= np.linalg.norm(axis)
+            Robj = R.from_rotvec(np.array(axis) * angle)
+            T = tesseral_wigner_D_improper(1, Robj)
+
+            # test that determinant of the transformation matrix is -1.0
+            # because the operation contains a reflection
+            np.testing.assert_almost_equal(np.linalg.det(T), -1.0)
+
+            # test the the matrix multiplied with its complex conjugate
+            # equals the identity matrix
+            np.testing.assert_almost_equal(T @ T.transpose(), np.identity(3))
+
+            # test that trace equals the character under S4
+            np.testing.assert_almost_equal(np.trace(T), -1)
+
+    def test_mirror_sigma_d(self):
+        """
+        Assert that the characters under S4 correspond to 1
+        """
+        coords = self.build_coordinates_methane()
+
+        for i in range(1,5):
+            for j in range(i+1,5):
+                axis = (coords[i] + coords[j]) / 2
+                M = tesseral_wigner_D_mirror(1, axis)
+
+                # test that determinant of the transformation matrix is -1.0
+                # because the operation contains a reflection
+                np.testing.assert_almost_equal(np.linalg.det(M), -1.0)
+
+                # test the the matrix multiplied with its complex conjugate
+                # equals the identity matrix
+                np.testing.assert_almost_equal(M @ M.transpose(), np.identity(3))
+
+                # test that trace equals the character under sigma_d
+                np.testing.assert_almost_equal(np.trace(M), 1)
+
+    def build_coordinates_methane(self):
+        """
+        Construct fictitious methane coordinates
+        """
+        coords = [
+            [0,0,0],
+            [1.0, 1.0, 1.0],
+            [1.0, -1.0, -1.0],
+            [-1.0, 1.0, -1.0],
+            [-1.0, -1.0, 1.0],
+        ]
+
+        return np.array(coords)
+
+if __name__ == '__main__':
+    unittest.main()