From 6f356b8c75a01a2952498027d864dc81a4c993a0 Mon Sep 17 00:00:00 2001 From: ifilot Date: Fri, 30 Jun 2023 21:49:39 +0200 Subject: [PATCH] Initial commit dodecahedrane --- examples/dodecahedrane.py | 119 ++++++++++++++++++++++++++++++ sphecerix/charactertables/ih.json | 87 ++++++++++++++++++++++ sphecerix/molecule.py | 15 ++++ 3 files changed, 221 insertions(+) create mode 100644 examples/dodecahedrane.py create mode 100644 sphecerix/charactertables/ih.json diff --git a/examples/dodecahedrane.py b/examples/dodecahedrane.py new file mode 100644 index 0000000..e2a5e64 --- /dev/null +++ b/examples/dodecahedrane.py @@ -0,0 +1,119 @@ +# -*- coding: utf-8 -*- + +import sys +import os +import numpy as np +from itertools import combinations + +# 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.from_file('molecules/dodecahedrane.xyz') + + 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') + + # Find C5 axes; for this specific implementation of the molecule, it is + # given that indices [0-4] form a pentagon. The center of all pentagons lie + # at the same distance with respec to the origin + c5_axes = [] + d = np.linalg.norm(np.average(np.take([at[1] for at in mol.atoms], range(5), axis=0), axis=0)) + for k in combinations(range(20), 5): + dd = np.linalg.norm(np.average(np.take([at[1] for at in mol.atoms], k, axis=0), axis=0)) + if np.abs(dd - d) < 1e-5: + c5_axes.append(np.average(np.take([at[1] for at in mol.atoms], k, axis=0), axis=0)) + c5_axes[-1] /= np.linalg.norm(c5_axes[-1]) + + for i,ax in enumerate(c5_axes): + symops.add('rotation', '5,%i' % (i+1), ax, 2.0 * np.pi / 5) + + # create c5^2 + for i,ax in enumerate(c5_axes): + symops.add('rotation', '5^2,%i' % (i+1), ax, 4.0 * np.pi / 5) + + # create c3 axes + for i,at in enumerate(mol.atoms[0:20]): + axis = at[1] / np.linalg.norm(at[1]) + symops.add('rotation', '3,%i' % (i+1), axis, 2.0 * np.pi / 3) + + # Find C2 axes; all C2 axes lie at the vertices of two adjacent atoms. We can + # adopt the same strategy as for the C5 axes, though it will yield exactly + # double duplicates as similar rotational axes lie at opposite sites of the + # icosphere. Note that we can also immediately use this procedure to find the + # mirror planes + d = np.linalg.norm(np.average(np.take([at[1] for at in mol.atoms], range(2), axis=0), axis=0)) + c2_axes = [] + mirror_normals = [] + for k in combinations(range(20), 2): + dd = np.linalg.norm(np.average(np.take([at[1] for at in mol.atoms], k, axis=0), axis=0)) + if np.abs(dd - d) < 1e-5: + axis = np.average(np.take([at[1] for at in mol.atoms], k, axis=0), axis=0) + if axis[0] < 0: # prune duplicates + continue + c2_axes.append(axis) + c2_axes[-1] /= np.linalg.norm(c2_axes[-1]) + + mirror_normals.append(np.cross(mol.atoms[k[0]][1], mol.atoms[k[1]][1])) + mirror_normals[-1] /= np.linalg.norm(mirror_normals[-1]) + + # C2 axes + for i,ax in enumerate(c2_axes): + symops.add('rotation', '2,%i' % (i+1), ax, np.pi) + + symops.add('inversion') + + # S10 improper rotations + for i,ax in enumerate(c5_axes): + symops.add('improper', '10,%i' % (i+1), ax, 2.0 * np.pi / 10) + + # S10^3 improper rotations + for i,ax in enumerate(c5_axes): + symops.add('improper', '10^3,%i' % (i+1), ax, 6.0 * np.pi / 10) + + # S6 operations + for i,at in enumerate(mol.atoms[0:20]): + axis = at[1] / np.linalg.norm(at[1]) + symops.add('improper', '6,%i' % (i+1), axis, 2.0 * np.pi / 6) + + # sigma mirror planes + for i,n in enumerate(mirror_normals): + symops.add('mirror', 'd,%i' % (i+1), n) + + symops.run() + + # print result LOT + ct = CharacterTable('ih') + print(ct.lot(np.trace(symops.operation_matrices, axis1=1, axis2=2))) + + # visualize_matrices(symops.operation_matrices[0:10], + # [op.name for op in symops.operations[0:10]], + # [bf.name for bf in symops.mol.basis], + # xlabelrot=90, figsize=(20,10), numcols=5) + + # # 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))], + # figsize=(18,10), numcols=4) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/sphecerix/charactertables/ih.json b/sphecerix/charactertables/ih.json new file mode 100644 index 0000000..85c8aff --- /dev/null +++ b/sphecerix/charactertables/ih.json @@ -0,0 +1,87 @@ +{ + "name":"Ih", + "classes": [ + { + "symbol": "E", + "multiplicity": 1 + }, + { + "symbol": "C5", + "multiplicity": 12 + }, + { + "symbol": "C5^2", + "multiplicity": 12 + }, + { + "symbol": "C3", + "multiplicity": 20 + }, + { + "symbol": "C2", + "multiplicity": 15 + }, + { + "symbol": "i", + "multiplicity": 1 + }, + { + "symbol": "S10", + "multiplicity": 12 + }, + { + "symbol": "S10^3", + "multiplicity": 12 + }, + { + "symbol": "S6", + "multiplicity": 20 + }, + { + "symbol": "sigma", + "multiplicity": 15 + } + ], + "symmetry_groups": [ + { + "symbol": "Ag", + "characters": [1,1,1,1,1,1,1,1,1,1] + }, + { + "symbol": "T1g", + "characters": [3,1.61803398875,-0.61803398875,0,-1,3,-0.61803398875,1.61803398875,0,-1] + }, + { + "symbol": "T2g", + "characters": [3,-0.61803398875,1.61803398875,0,-1,3,1.61803398875,-0.61803398875,0,-1] + }, + { + "symbol": "Gg", + "characters": [4,-1,-1,1,0,4,-1,-1,1,0] + }, + { + "symbol": "Hg", + "characters": [5,0,0,-1,1,5,0,0,-1,1] + }, + { + "symbol": "Au", + "characters": [1,1,1,1,1,-1,-1,-1,-1,-1] + }, + { + "symbol": "T1u", + "characters": [3,1.61803398875,-0.61803398875,0,-1,-3,0.61803398875,-1.61803398875,0,1] + }, + { + "symbol": "T2u", + "characters": [3,-0.61803398875,1.61803398875,0,-1,-3,-1.61803398875,0.61803398875,0,1] + }, + { + "symbol": "Gu", + "characters": [4,-1,-1,1,0,-4,1,1,-1,0] + }, + { + "symbol": "Hu", + "characters": [5,0,0,-1,1,-5,0,0,1,-1] + } + ] +} \ No newline at end of file diff --git a/sphecerix/molecule.py b/sphecerix/molecule.py index babd05b..db1544e 100644 --- a/sphecerix/molecule.py +++ b/sphecerix/molecule.py @@ -13,6 +13,21 @@ def __init__(self, _name='unknown'): self.name = _name self.basis = None + def from_file(self, path, molname=None): + """ + Build molecule from file and return it + """ + self.name = molname + + with open(path, 'r') as f: + lines = f.readlines() + + nratoms = int(lines[0].strip()) + + for line in lines[2:2+nratoms]: + pieces = line.split() + self.add_atom(pieces[0], float(pieces[1]), float(pieces[2]), float(pieces[3]), unit='angstrom') + def __str__(self): res = "Molecule: %s\n" % self.name for atom in self.atoms: