Adding methane and Td group

examples/ch4.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+
+import sys
+import os
+import numpy as np
+
+# add a reference to load the Sphecerix library
+sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
+
+from sphecerix import Molecule, BasisFunction, SymmetryOperations,\
+                      visualize_matrices, CharacterTable, ProjectionOperator
+
+def main():
+    mol = Molecule()
+    mol.add_atom('C', 0.0,  0.0, 0.0, unit='angstrom')
+    mol.add_atom('H', 1,1,1, unit='angstrom')
+    mol.add_atom('H', 1,-1,-1, unit='angstrom')
+    mol.add_atom('H', -1,1,-1, unit='angstrom')
+    mol.add_atom('H', -1,-1,1, unit='angstrom')
+
+    molset = {
+        'C': [BasisFunction(1,0,0),
+              BasisFunction(2,0,0),
+              BasisFunction(2,1,1),
+              BasisFunction(2,1,-1),
+              BasisFunction(2,1,0)],
+        'H': [BasisFunction(1,0,0)]
+    }
+    mol.build_basis(molset)
+
+    symops = SymmetryOperations(mol)
+    symops.add('identity')
+
+    # add C3 rotations
+    for i in range(0,4):
+        axis = mol.atoms[i+1][1]
+        axis /= np.linalg.norm(axis) # normalize axis
+        for j in range(0,2):
+            symops.add('rotation', '3,%i' % (i*2+j+1), axis, (-1)**j * 2.0 * np.pi / 3)
+
+    # C2 rotations
+    for i in range(0,3):
+        axis = np.zeros(3)
+        axis[i] = 1.0
+        symops.add('rotation', '2,%i' % (i+1), axis, np.pi)
+
+    # S4 rotations
+    for i in range(0,3):
+        axis = np.zeros(3)
+        axis[i] = 1.0
+        for j in range(0,2):
+            symops.add('improper', '4,%i' % (i+1), axis, (-1)**j * np.pi/2)
+
+    # sigma_d mirror planes
+    ctr = 0
+    for i in range(0,4):
+        axis1 = mol.atoms[i+1][1]
+        for j in range(i+1,4):
+            axis2 = mol.atoms[j+1][1]
+            normal = np.cross(axis1, axis2)
+            normal /= np.linalg.norm(normal)
+            ctr += 1
+            symops.add('mirror', ',d(%i)' % (ctr), normal)
+
+    symops.run()
+
+    visualize_matrices(symops.operation_matrices, 
+                       [op.name for op in  symops.operations],
+                       [bf.name for bf in symops.mol.basis], 
+                       xlabelrot=90, figsize=(18,12), numcols=6)
+
+    # print result LOT
+    ct = CharacterTable('td')
+    print(ct.lot(np.trace(symops.operation_matrices, axis1=1, axis2=2)))
+
+    # # apply projection operator
+    po = ProjectionOperator(ct, symops)
+    mos = po.build_mos()
+    newmats = [mos @ m @ mos.transpose() for m in symops.operation_matrices]
+
+    visualize_matrices(newmats, 
+                        [op.name for op in  symops.operations],
+                        ['$\phi_{%i}$' % (i+1) for i in range(len(symops.mol.basis))],
+                        xlabelrot=90, figsize=(18,12), numcols=6)
+
+if __name__ == '__main__':
+    main()

sphecerix/charactertables/td.json

@@ -0,0 +1,47 @@
+{
+   "name":"c3v",
+   "classes": [
+       {
+           "symbol": "E",
+           "multiplicity": 1
+       },
+       {
+           "symbol": "C3",
+           "multiplicity": 8
+       },
+       {
+           "symbol": "C2",
+           "multiplicity": 3
+       },
+       {
+           "symbol": "S4",
+           "multiplicity": 6
+       },
+       {
+           "symbol": "sigma_d",
+           "multiplicity": 6
+       }
+   ],
+   "symmetry_groups": [
+       {
+           "symbol": "A1",
+           "characters": [1,1,1,1,1]
+       },
+       {
+           "symbol": "A2",
+           "characters": [1,1,1,-1,-1]
+       },
+       {
+           "symbol": "E",
+           "characters": [2,-1,2,0,0]
+       },
+       {
+           "symbol": "T1",
+           "characters": [3,0,-1,1,-1]
+       },
+       {
+           "symbol": "T2",
+           "characters": [3,0,-1,-1,1]
+       }
+   ]
+}

sphecerix/symmetry_operations.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.spatial.transform import Rotation as R
-from . import tesseral_wigner_D, tesseral_wigner_D_mirror
+from . import tesseral_wigner_D, tesseral_wigner_D_mirror, tesseral_wigner_D_improper
 
 class SymmetryOperations:
     """
@@ -143,12 +143,15 @@ class ImproperRotation(Operation):
     Rotation operation "S"
     """
     def __init__(self, label, axis, angle):
-        super().__init__('sigma' + label)
+        super().__init__('S' + label)
         self.axis = axis / np.linalg.norm(axis)
         self.angle = angle
         self.robj = R.from_rotvec(self.axis * self.angle)
 
     def get_matrix(self):
-        M = np.identity(3) - 2.0 * self.normal @ self.normal.transpose()
+        M = np.identity(3) - 2 * np.outer(self.axis, self.axis)
         R = self.robj.as_matrix()
-        return M @ R
+        return M @ R
+
+    def get_wigner_matrix(self, l):
+        return tesseral_wigner_D_improper(l, self.robj)

tests/test_methane_symmetry_operations.py

@@ -0,0 +1,78 @@
+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 Molecule, BasisFunction, SymmetryOperations,\
+                      CharacterTable
+
+class TestEthyleneSymmetryOperations(unittest.TestCase):
+    """
+    Test all the symmetry operations under ethylene (D2h)
+    """
+
+    def test_ethylene(self):
+        mol = Molecule()
+        mol.add_atom('C', 0.0,  0.0, 0.0, unit='angstrom')
+        mol.add_atom('H', 1,1,1, unit='angstrom')
+        mol.add_atom('H', 1,-1,-1, unit='angstrom')
+        mol.add_atom('H', -1,1,-1, unit='angstrom')
+        mol.add_atom('H', -1,-1,1, unit='angstrom')
+
+        molset = {
+            'C': [BasisFunction(1,0,0),
+                  BasisFunction(2,0,0),
+                  BasisFunction(2,1,1),
+                  BasisFunction(2,1,-1),
+                  BasisFunction(2,1,0)],
+            'H': [BasisFunction(1,0,0)]
+        }
+        mol.build_basis(molset)
+
+        symops = SymmetryOperations(mol)
+        symops.add('identity')
+
+        # add C3 rotations
+        for i in range(0,4):
+            axis = mol.atoms[i+1][1]
+            axis /= np.linalg.norm(axis) # normalize axis
+            for j in range(0,2):
+                symops.add('rotation', '3,%i' % (i*2+j+1), axis, (-1)**j * 2.0 * np.pi / 3)
+
+        # C2 rotations
+        for i in range(0,3):
+            axis = np.zeros(3)
+            axis[i] = 1.0
+            symops.add('rotation', '2,%i' % (i+1), axis, np.pi)
+
+        # S4 rotations
+        for i in range(0,3):
+            axis = np.zeros(3)
+            axis[i] = 1.0
+            for j in range(0,2):
+                symops.add('improper', '4,%i' % (i+1), axis, (-1)**j * np.pi/2)
+
+        # sigma_d mirror planes
+        ctr = 0
+        for i in range(0,4):
+            axis1 = mol.atoms[i+1][1]
+            for j in range(i+1,4):
+                axis2 = mol.atoms[j+1][1]
+                normal = np.cross(axis1, axis2)
+                normal /= np.linalg.norm(normal)
+                ctr += 1
+                symops.add('mirror', ',d(%i)' % (ctr), normal)
+
+        symops.run()
+
+        # print result LOT
+        ct = CharacterTable('td')
+        irreps = ct.lot(np.trace(symops.operation_matrices, axis1=1, axis2=2))
+        np.testing.assert_almost_equal(irreps, [3,0,0,0,2])
+
+if __name__ == '__main__':
+    unittest.main()