diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 0cf246d..e4c576a 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -18,14 +18,14 @@ jobs: yum -y install mlocate yum -y install gcc ls -alh /opt/python/cp* - /opt/python/cp39-cp39/bin/python -m pip install numpy tqdm cython pytest scipy pylebedev + /opt/python/cp39-cp39/bin/python -m pip install numpy tqdm cython pytest scipy pylebedev matplotlib - name: Checkout repo uses: actions/checkout@v3 - name: Build run: /opt/python/cp39*/bin/python setup.py bdist_wheel - name: Test run: | - /opt/python/cp39-cp39/bin/python -m pip install numpy scipy nose pylebedev + /opt/python/cp39-cp39/bin/python -m pip install numpy scipy nose pylebedev matplotlib /opt/python/cp39-cp39/bin/pip install sphecerix --no-index -f ./dist /opt/python/cp39-cp39/bin/pytest --verbose tests - name: Archive wheels diff --git a/.github/workflows/deploy.yml b/.github/workflows/deploy.yml index a587ce0..0c36f47 100644 --- a/.github/workflows/deploy.yml +++ b/.github/workflows/deploy.yml @@ -17,14 +17,14 @@ jobs: yum -y install mlocate yum -y install gcc ls -alh /opt/python/cp* - /opt/python/cp39-cp39/bin/python -m pip install numpy tqdm cython pytest scipy pylebedev + /opt/python/cp39-cp39/bin/python -m pip install numpy tqdm cython pytest scipy pylebedev matplotlib - name: Checkout repo uses: actions/checkout@v3 - name: Build run: /opt/python/cp39*/bin/python setup.py bdist_wheel - name: Test run: | - /opt/python/cp39-cp39/bin/python -m pip install numpy scipy nose pylebedev + /opt/python/cp39-cp39/bin/python -m pip install numpy scipy nose pylebedev matplotlib /opt/python/cp39-cp39/bin/pip install sphecerix --no-index -f ./dist /opt/python/cp39-cp39/bin/pytest --verbose tests - name: Archive wheels diff --git a/examples/ethylene.py b/examples/ethylene.py new file mode 100644 index 0000000..d6a3b34 --- /dev/null +++ b/examples/ethylene.py @@ -0,0 +1,47 @@ +# -*- 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 + +def main(): + mol = Molecule() + mol.add_atom('C', -0.6530176758, 0.0000000000 ,0.0000000000, unit='angstrom') + mol.add_atom('C', 0.6530176758, 0.0000000000 ,0.0000000000, unit='angstrom') + mol.add_atom('H', -1.2288875372, -0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', -1.2288875372, 0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', 1.2288875372, 0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', 1.2288875372, -0.9156191261 ,0.0000000000, 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') + symops.add('rotation', '2(z)', np.array([0,0,1]), np.pi) + symops.add('rotation', '2(y)', np.array([0,1,0]), np.pi) + symops.add('rotation', '2(x)', np.array([1,0,0]), np.pi) + symops.add('inversion') + symops.add('mirror', 'v(xy)', np.array([0,0,1])) + symops.add('mirror', 'v(xz)', np.array([0,1,0])) + symops.add('mirror', 'v(yz)', np.array([1,0,0])) + + symops.run() + + visualize_matrices(symops, xlabelrot=90, figsize=(18,10), numcols=4) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/examples/nh3.py b/examples/nh3.py new file mode 100644 index 0000000..0ee701f --- /dev/null +++ b/examples/nh3.py @@ -0,0 +1,44 @@ +# -*- 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 + +def main(): + mol = Molecule() + mol.add_atom('N', 0.00000000, 0.00000000, -0.06931370, unit='angstrom') + mol.add_atom('H', 0.00000000, 0.94311105, 0.32106944, unit='angstrom') + mol.add_atom('H', -0.81675813, -0.47155553, 0.32106944, unit='angstrom') + mol.add_atom('H', 0.81675813, -0.47155553, 0.32106944, unit='angstrom') + + molset = { + 'N': [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') + symops.add('rotation', '3+', np.array([0,0,1]), 2.0 * np.pi / 3) + symops.add('rotation', '3-', -np.array([0,0,1]), 2.0 * np.pi / 3) + + for i in range(0,3): + symops.add('mirror', 'v1', np.array([np.cos(i * 2.0 * np.pi / 3), + np.sin(i * 2.0 * np.pi / 3), + 0.0])) + + symops.run() + + visualize_matrices(symops, xlabelrot=90, figsize=(9,6)) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/meta.yaml b/meta.yaml index 2f2b103..1bd5f52 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,6 +1,6 @@ package: name: "sphecerix" - version: "0.4.0" + version: "0.5.0" source: path: . @@ -13,11 +13,12 @@ requirements: build: - numpy>=1.21 - python>=3.9 + - scipy host: - pip - python>=3.9 - - setuptools + - setuptools<=58.2.0 - numpy>=1.21 - scipy @@ -25,15 +26,18 @@ requirements: - python>=3.9 - numpy>=1.21 - scipy + - matplotlib test: requires: - numpy - scipy - - setuptools + - matplotlib + - setuptools<=58.2.0 - nose source_files: - tests/*.py + - tests/results/*.npy commands: - nosetests tests --exclude="test_wave_functions.py" diff --git a/sphecerix/__init__.py b/sphecerix/__init__.py index ee9c375..767e458 100644 --- a/sphecerix/__init__.py +++ b/sphecerix/__init__.py @@ -2,5 +2,9 @@ tesseral_wigner_D_improper from .tesseral import tesseral_transformation, permutation_sh_car from .atomic_wave_functions import wfcart, wf, wffield, wffield_l +from .molecule import Molecule +from .basis_functions import BasisFunction +from .symmetry_operations import * +from .matrixplot import plot_matrix, visualize_matrices from ._version import __version__ \ No newline at end of file diff --git a/sphecerix/_version.py b/sphecerix/_version.py index 49c37a8..5a28280 100644 --- a/sphecerix/_version.py +++ b/sphecerix/_version.py @@ -1 +1 @@ -__version__ = "0.4.0" \ No newline at end of file +__version__ = "0.5.0" \ No newline at end of file diff --git a/sphecerix/basis_functions.py b/sphecerix/basis_functions.py new file mode 100644 index 0000000..c9db4a9 --- /dev/null +++ b/sphecerix/basis_functions.py @@ -0,0 +1,27 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +class BasisFunction: + + def __init__(self, n, l, m): + self.r = np.array([0,0,0]) + self.n = n + self.l = l + self.m = m + self.atomid = None + self.name = None + + self.name = self.__get_name() + + def __get_name(self): + return str(self.n) + self.__get_type() + + def __get_type(self): + results = [ + ['s'], + ['py', 'pz', 'px'], + ['dxy', 'dyz', 'dz2', 'dxz', 'dx2-y2'] + ] + + return results[self.l][self.m + self.l] \ No newline at end of file diff --git a/sphecerix/matrixplot.py b/sphecerix/matrixplot.py new file mode 100644 index 0000000..6f4bae0 --- /dev/null +++ b/sphecerix/matrixplot.py @@ -0,0 +1,94 @@ +# -*- coding: utf-8 -*- + +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.patches as patches + +def visualize_matrices(symops, numcols = 3, + highlight_groups = None, filename = None, + figsize=(7,5), xlabelrot = 0): + """ + Visualize matrix representations of the symmetry operations + """ + # grab data + matrices = symops.operation_matrices + operations = symops.operations + bfs = symops.mol.basis + + fig, ax = plt.subplots(len(operations) // numcols, + numcols, dpi=144, figsize=figsize) + + for i,(op,mat) in enumerate(zip(operations,matrices)): + axh = ax[i//numcols, i%numcols] + plot_matrix(axh, mat, bfs, title=op.name, xlabelrot = xlabelrot) + + if highlight_groups: + plot_highlight_groups(axh, highlight_groups, mat) + + plt.tight_layout() + + if filename: + print('Storing: %s' % filename) + plt.savefig(filename) + plt.close() + +def plot_highlight_groups(axh, groups, mat): + # add semitransparent hash + cum = 0 + for g in groups: + rect = patches.Rectangle((cum - 0.5, cum - 0.5), g, g, + linewidth=1, + zorder=5, + fill = None, + hatch='///', + alpha=0.5) + axh.add_patch(rect) + cum += g + + # add red outline + cum = 0 + for g in groups: + rect = patches.Rectangle((cum - 0.5, cum - 0.5), g, g, + linewidth=1.5, edgecolor='red', + linestyle='solid', + facecolor='none', + zorder=5, + alpha=1.0) + axh.add_patch(rect) + + axh.text(cum+g/2-0.5, cum+g/2-0.5, '%i' % round(np.trace(mat[cum:cum+g,cum:cum+g])), + color='red', horizontalalignment='center', verticalalignment='center', + bbox=dict(boxstyle="round", ec=(1., 0.5, 0.5), fc=(1., 0.8, 0.8), ), + zorder=6) + + cum += g + + +def plot_matrix(ax, mat, bfs, title = None, xlabelrot = 0): + """ + Produce plot of matrix + """ + ax.imshow(mat, vmin=-1, vmax=1, cmap='PiYG') + for i in range(mat.shape[0]): + for j in range(mat.shape[1]): + ax.text(i, j, '%.2f' % mat[j,i], ha='center', va='center', + fontsize=7) + ax.set_xticks([]) + ax.set_yticks([]) + ax.hlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + ax.vlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + + # add basis functions as axes labels + labels = [bf.name for bf in bfs] + ax.set_xticks(np.arange(0, mat.shape[0])) + ax.set_xticklabels(labels, rotation=xlabelrot) + ax.set_yticks(np.arange(0, mat.shape[0])) + ax.set_yticklabels(labels, rotation=0) + ax.tick_params(axis='both', which='major', labelsize=7) + + # add title if supplied + if title: + ax.set_title(title) + \ No newline at end of file diff --git a/sphecerix/molecule.py b/sphecerix/molecule.py new file mode 100644 index 0000000..babd05b --- /dev/null +++ b/sphecerix/molecule.py @@ -0,0 +1,50 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from copy import deepcopy + +class Molecule: + """ + Molecule class + """ + def __init__(self, _name='unknown'): + self.atoms = [] + self.charges = [] + self.name = _name + self.basis = None + + def __str__(self): + res = "Molecule: %s\n" % self.name + for atom in self.atoms: + res += " %s (%f,%f,%f)\n" % (atom[0], atom[1][0], atom[1][1], atom[1][2]) + + return res + + def add_atom(self, atom, x, y, z, unit='bohr'): + ang2bohr = 1.8897259886 + + x = float(x) + y = float(y) + z = float(z) + + if unit == "bohr": + self.atoms.append([atom, np.array([x, y, z])]) + elif unit == "angstrom": + self.atoms.append([atom, np.array([x*ang2bohr, y*ang2bohr, z*ang2bohr])]) + else: + raise RuntimeError("Invalid unit encountered: %s. Accepted units are 'bohr' and 'angstrom'." % unit) + + self.charges.append(0) + + def build_basis(self, molset): + self.basis = [] + + for i,atom in enumerate(self.atoms): + if atom[0] in molset.keys(): + for bf in molset[atom[0]]: + abf = deepcopy(bf) + abf.name = atom[0] + bf.name + abf.r = np.array([atom[1][i] for i in range(3)]) + abf.atomid = i + self.basis.append(abf) + \ No newline at end of file diff --git a/sphecerix/symmetry_operations.py b/sphecerix/symmetry_operations.py new file mode 100644 index 0000000..edfd78c --- /dev/null +++ b/sphecerix/symmetry_operations.py @@ -0,0 +1,154 @@ +import numpy as np +from scipy.spatial.transform import Rotation as R +from . import tesseral_wigner_D, tesseral_wigner_D_mirror + +class SymmetryOperations: + """ + Class containing all symmetry operations as applied to molecule object + """ + def __init__(self, mol): + self.mol = mol + self.operations = [] + + # build atoms matrix + self.positions = np.zeros((len(self.mol.atoms), 3)) + for i,atom in enumerate(self.mol.atoms): + self.positions[i,:] = atom[1] + + def add(self, name, label = None, vec = None, angle = None): + # ensure vector is of float type + if vec is not None: + vec = np.array(vec, dtype=np.float64) + + if name == 'identity': + self.operations.append(Identity()) + elif name == 'rotation': + self.operations.append(Rotation(label, vec, angle)) + elif name == 'mirror': + self.operations.append(Mirror(label, vec)) + elif name == 'improper': + self.operations.append(ImproperRotation(label, vec, angle)) + elif name == 'inversion': + self.operations.append(Inversion()) + else: + raise Exception('Unknown operation: %s' % name) + + def run(self): + N = len(self.mol.atoms) # number of atoms + nbf = len(self.mol.basis) # number of basis functions + self.atomic_transformations = np.zeros((len(self.operations),N), dtype=np.uint64) + self.basis_function_transformations = np.zeros((len(self.mol.basis), N, N)) + + # assert atomic operations + for k,operation in enumerate(self.operations): + tpos = self.positions @ operation.get_matrix() + + for i,p in enumerate(self.positions): + for j,t in enumerate(tpos): + r = t-p + if np.sum(r**2) < 1e-5: + self.atomic_transformations[k,i] = j + + # cache on which atoms basis functions are located + bfindices = -np.ones((N, 5, 3, 7), dtype=np.int64) # encode all possibilities up to 5f-orbitals + for i,bf in enumerate(self.mol.basis): + bfindices[bf.atomid, bf.n, bf.l, bf.m+bf.l] = i + + # assert basis function operations + self.operation_matrices = np.zeros((len(self.operations), nbf, nbf)) + for k,operation in enumerate(self.operations): + for i,bf in enumerate(self.mol.basis): + # establish on which atom id the basis function lands + atomid = self.atomic_transformations[k,bf.atomid] + + # establish tesseral transformation + mvec = np.zeros(2*bf.l+1) + mvec[bf.m + bf.l] = 1 + mres = operation.get_wigner_matrix(bf.l).dot(mvec) + + for m,v in enumerate(mres): + idx = bfindices[atomid, bf.n, bf.l, m] + if idx == -1: + raise Exception('Cannot perform this operation') + self.operation_matrices[k,i,idx] = v + +class Operation: + """ + Base operation class + """ + def __init__(self, name): + self.name = name + + def set_atomic_id(self, idx): + self.atomid = idx + +class Identity(Operation): + """ + Identity operation "E" + """ + def __init__(self): + super().__init__('E') + + def get_matrix(self): + return np.identity(3) + + def get_wigner_matrix(self, l): + return np.identity(2*l+1) + +class Inversion(Operation): + """ + Identity operation "i" + """ + def __init__(self): + super().__init__('i') + + def get_matrix(self): + return -np.identity(3) + + def get_wigner_matrix(self, l): + return np.identity(2*l+1) * (-1)**l + +class Rotation(Operation): + """ + Rotation operation "C" + """ + def __init__(self, label, axis, angle): + super().__init__('C' + label) + self.axis = axis / np.linalg.norm(axis) + self.angle = angle + self.robj = R.from_rotvec(self.axis * self.angle) + + def get_matrix(self): + return self.robj.as_matrix() + + def get_wigner_matrix(self, l): + return tesseral_wigner_D(l, self.robj) + +class Mirror(Operation): + """ + Rotation operation "σ" + """ + def __init__(self, label, normal): + super().__init__('σ' + label) + self.normal = normal + + def get_matrix(self): + return np.identity(3) - 2 * np.outer(self.normal, self.normal) + + def get_wigner_matrix(self, l): + return tesseral_wigner_D_mirror(l, self.normal) + +class ImproperRotation(Operation): + """ + Rotation operation "S" + """ + def __init__(self, label, axis, angle): + super().__init__('sigma' + 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() + R = self.robj.as_matrix() + return M @ R \ No newline at end of file diff --git a/tests/results/ethylene.npy b/tests/results/ethylene.npy new file mode 100644 index 0000000..a0f4ebc Binary files /dev/null and b/tests/results/ethylene.npy differ diff --git a/tests/results/nh3.npy b/tests/results/nh3.npy new file mode 100644 index 0000000..7c4dca0 Binary files /dev/null and b/tests/results/nh3.npy differ diff --git a/tests/test_ammonia_symmetry_operations.py b/tests/test_ammonia_symmetry_operations.py new file mode 100644 index 0000000..63ab681 --- /dev/null +++ b/tests/test_ammonia_symmetry_operations.py @@ -0,0 +1,51 @@ +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 + +class TestAmmoniaSymmetryOperations(unittest.TestCase): + """ + Test all the symmetry operations under ethylene (D2h) + """ + + def test_nh3(self): + mol = Molecule() + mol.add_atom('N', 0.00000000, 0.00000000, -0.06931370, unit='angstrom') + mol.add_atom('H', 0.00000000, 0.94311105, 0.32106944, unit='angstrom') + mol.add_atom('H', -0.81675813, -0.47155553, 0.32106944, unit='angstrom') + mol.add_atom('H', 0.81675813, -0.47155553, 0.32106944, unit='angstrom') + + molset = { + 'N': [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') + symops.add('rotation', '3+', np.array([0,0,1]), 2.0 * np.pi / 3) + symops.add('rotation', '3-', -np.array([0,0,1]), 2.0 * np.pi / 3) + + for i in range(0,3): + symops.add('mirror', 'v1', np.array([np.cos(i * 2.0 * np.pi / 3), + np.sin(i * 2.0 * np.pi / 3), + 0.0])) + + symops.run() + + result = np.load(os.path.join(os.path.dirname(__file__), 'results', 'nh3.npy')) + + np.testing.assert_almost_equal(symops.operation_matrices, result) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_ethylene_symmetry_operations.py b/tests/test_ethylene_symmetry_operations.py new file mode 100644 index 0000000..933607f --- /dev/null +++ b/tests/test_ethylene_symmetry_operations.py @@ -0,0 +1,52 @@ +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 + +class TestEthyleneSymmetryOperations(unittest.TestCase): + """ + Test all the symmetry operations under ethylene (D2h) + """ + + def test_ethylene(self): + mol = Molecule() + mol.add_atom('C', -0.6530176758, 0.0000000000 ,0.0000000000, unit='angstrom') + mol.add_atom('C', 0.6530176758, 0.0000000000 ,0.0000000000, unit='angstrom') + mol.add_atom('H', -1.2288875372, -0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', -1.2288875372, 0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', 1.2288875372, 0.9156191261 ,0.0000000000, unit='angstrom') + mol.add_atom('H', 1.2288875372, -0.9156191261 ,0.0000000000, 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') + symops.add('rotation', '2(z)', np.array([0,0,1]), np.pi) + symops.add('rotation', '2(y)', np.array([0,1,0]), np.pi) + symops.add('rotation', '2(x)', np.array([1,0,0]), np.pi) + symops.add('inversion') + symops.add('mirror', 'v(xy)', np.array([0,0,1])) + symops.add('mirror', 'v(xz)', np.array([0,1,0])) + symops.add('mirror', 'v(yz)', np.array([1,0,0])) + + symops.run() + + result = np.load(os.path.join(os.path.dirname(__file__), 'results', 'ethylene.npy')) + np.testing.assert_almost_equal(symops.operation_matrices, result) + +if __name__ == '__main__': + unittest.main()