Probing symmetry operations (#6)

* Adding methane to test suite

* Adding ammonia and adding matrix generation for improper rotations

* Adding improper rotations to tests

* Expanding documentation

* Expanding documentation

* Initial commit applying matrix operations to molecules

* Establishing result of symmetry operations

* Adding result sets

* Adding dependency

* Fixing deprecated np.float

* Fixing setuptools

* Update meta.yaml

* Update meta.yaml

* Update meta.yaml

* Update meta.yaml

* Update setup.py

* Update meta.yaml

.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

.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

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()

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()

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "sphecerix"
-  version: "0.4.0"
+  version: "0.5.0"
 
 source:
   path: .
@@ -13,27 +13,31 @@ requirements:
   build:
     - numpy>=1.21
     - python>=3.9
+    - scipy
 
   host:
     - pip
     - python>=3.9
-    - setuptools
+    - setuptools<=58.2.0
     - numpy>=1.21
     - scipy
 
   run:
     - 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"
 

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__

sphecerix/_version.py

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

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]

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)
+

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)
+