diff --git a/examples/ch4.py b/examples/ch4.py new file mode 100644 index 0000000..c4a2c9e --- /dev/null +++ b/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() \ No newline at end of file diff --git a/sphecerix/charactertables/td.json b/sphecerix/charactertables/td.json new file mode 100644 index 0000000..30f96cc --- /dev/null +++ b/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] + } + ] +} \ No newline at end of file diff --git a/sphecerix/symmetry_operations.py b/sphecerix/symmetry_operations.py index edfd78c..4d4e1ae 100644 --- a/sphecerix/symmetry_operations.py +++ b/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 \ No newline at end of file + return M @ R + + def get_wigner_matrix(self, l): + return tesseral_wigner_D_improper(l, self.robj) \ No newline at end of file diff --git a/tests/test_methane_symmetry_operations.py b/tests/test_methane_symmetry_operations.py new file mode 100644 index 0000000..af86c11 --- /dev/null +++ b/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()