Skip to content

Commit

Permalink
Adding projection operator for NH3 and ethylene
Browse files Browse the repository at this point in the history
  • Loading branch information
ifilot committed Jun 30, 2023
1 parent f51ddbc commit 03b23e5
Show file tree
Hide file tree
Showing 8 changed files with 230 additions and 20 deletions.
21 changes: 19 additions & 2 deletions examples/ethylene.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
# 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
from sphecerix import Molecule, BasisFunction, SymmetryOperations,\
visualize_matrices, CharacterTable, ProjectionOperator

def main():
mol = Molecule()
Expand Down Expand Up @@ -40,8 +41,24 @@ def main():

symops.run()

visualize_matrices(symops, xlabelrot=90, figsize=(18,10), numcols=4)
visualize_matrices(symops.operation_matrices,
[op.name for op in symops.operations],
[bf.name for bf in symops.mol.basis],
xlabelrot=90, figsize=(18,10), numcols=4)

# print result LOT
ct = CharacterTable('d2h')
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))],
figsize=(18,10), numcols=4)

if __name__ == '__main__':
main()
19 changes: 16 additions & 3 deletions examples/nh3.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
# 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
from sphecerix import Molecule, BasisFunction, SymmetryOperations,\
visualize_matrices, CharacterTable, ProjectionOperator

def main():
mol = Molecule()
Expand Down Expand Up @@ -39,11 +39,24 @@ def main():

symops.run()

visualize_matrices(symops, xlabelrot=90, figsize=(9,6))
visualize_matrices(symops.operation_matrices,
[op.name for op in symops.operations],
[bf.name for bf in symops.mol.basis],
xlabelrot=90, figsize=(9,6))

# print result LOT
ct = CharacterTable('c3v')
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))],
figsize=(9,6))

if __name__ == '__main__':
main()
1 change: 1 addition & 0 deletions sphecerix/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,5 +7,6 @@
from .symmetry_operations import *
from .matrixplot import plot_matrix, visualize_matrices
from .character_table import CharacterTable
from .projection_operator import ProjectionOperator

from ._version import __version__
12 changes: 8 additions & 4 deletions sphecerix/character_table.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,14 +18,18 @@ def __init__(self, name):

# build expanded table
self.expandedtable = np.zeros((self.nrgroups, self.order), dtype=np.int64)
ctr = 0
for i,g in enumerate(self.chartablelib['symmetry_groups']):
idx = 0
for i,g in enumerate(self.chartablelib['symmetry_groups']): # loop over irreps
idx = 0 # loop over operations
for j,c in enumerate(self.chartablelib['classes']):
for k in range(c['multiplicity']):
self.expandedtable[i,idx] = g['characters'][j]
idx += 1

def lot(self, traces):
return np.round(self.expandedtable @ traces / self.order)


def get_label_irrep(self, irrep_idx):
return self.chartablelib['symmetry_groups'][irrep_idx]['symbol']

def get_character(self, irrep_idx, op_idx):
return self.expandedtable[irrep_idx, op_idx]
71 changes: 71 additions & 0 deletions sphecerix/charactertables/d2h.json
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
{
"name":"d2h",
"classes": [
{
"symbol": "E",
"multiplicity": 1
},
{
"symbol": "C2(z)",
"multiplicity": 1
},
{
"symbol": "C2(y)",
"multiplicity": 1
},
{
"symbol": "C2(x)",
"multiplicity": 1
},
{
"symbol": "i",
"multiplicity": 1
},
{
"symbol": "sigma(xy)",
"multiplicity": 1
},
{
"symbol": "sigma(xz)",
"multiplicity": 1
},
{
"symbol": "sigma(yz)",
"multiplicity": 1
}
],
"symmetry_groups": [
{
"symbol": "Ag",
"characters": [1,1,1,1,1,1,1,1]
},
{
"symbol": "B1g",
"characters": [1,1,-1,-1,1,1,-1,-1]
},
{
"symbol": "B2g",
"characters": [1,-1,1,-1,1,-1,1,-1]
},
{
"symbol": "B3g",
"characters": [1,-1,-1,1,1,-1,-1,1]
},
{
"symbol": "Au",
"characters": [1,1,1,1,-1,-1,-1,-1]
},
{
"symbol": "B1u",
"characters": [1,1,-1,-1,-1,-1,1,1]
},
{
"symbol": "B2u",
"characters": [1,-1,1,-1,-1,1,-1,1]
},
{
"symbol": "B3u",
"characters": [1,-1,-1,1,-1,1,1,-1]
}
]
}
15 changes: 5 additions & 10 deletions sphecerix/matrixplot.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,23 +4,19 @@
import numpy as np
import matplotlib.patches as patches

def visualize_matrices(symops, numcols = 3,
def visualize_matrices(matrices, opnames, labels, 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,
fig, ax = plt.subplots(len(matrices) // numcols,
numcols, dpi=144, figsize=figsize)

for i,(op,mat) in enumerate(zip(operations,matrices)):
for i,(opname,mat) in enumerate(zip(opnames,matrices)):
axh = ax[i//numcols, i%numcols]
plot_matrix(axh, mat, bfs, title=op.name, xlabelrot = xlabelrot)
plot_matrix(axh, mat, labels, title=opname, xlabelrot = xlabelrot)

if highlight_groups:
plot_highlight_groups(axh, highlight_groups, mat)
Expand Down Expand Up @@ -64,7 +60,7 @@ def plot_highlight_groups(axh, groups, mat):
cum += g


def plot_matrix(ax, mat, bfs, title = None, xlabelrot = 0):
def plot_matrix(ax, mat, labels, title = None, xlabelrot = 0):
"""
Produce plot of matrix
"""
Expand All @@ -81,7 +77,6 @@ def plot_matrix(ax, mat, bfs, title = None, xlabelrot = 0):
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]))
Expand Down
101 changes: 101 additions & 0 deletions sphecerix/projection_operator.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
# -*- coding: utf-8 -*-

import numpy as np
from collections import Counter

class ProjectionOperator:

def __init__(self, ct, so):
#
# ToDo: assert that the objects have the right type
#

self.ct = ct # character table
self.so = so # symmetry operations
self.groups = None
self.irreps = None

def collect(self):
"""
Construct groups of basis functions that can transform among each
other and figure out which irreps these span
"""
# create one logical matrix that shows the interaction
groupmatrix = np.zeros_like(self.so.operation_matrices[0], dtype=np.bool8)
for m in self.so.operation_matrices:
bm = np.abs(m) > 1e-9
groupmatrix = np.logical_or(groupmatrix, bm)

#print(groupmatrix)

# collect the groups
self.groups = []
for i,row in enumerate(groupmatrix):
if np.sum(row) > 1:
idx = np.where(row)[0]
self.groups.append(np.array(idx, dtype=np.int64))
else:
self.groups.append(np.array([i], dtype=np.int64))

# create unique list
self.groups = np.array(self.groups, dtype=object)
res = Counter(map(tuple, self.groups))
self.groups = [r for r in res]

# determine irreps per lips
self.irreps = []
for g in self.groups:
chars = [np.sum(np.take(np.diagonal(m),g)) for m in self.so.operation_matrices]
self.irreps.append(self.ct.lot(chars))

def build_mos(self):
# check if groups have been collected
if self.groups is None:
self.collect()

# build molecular orbitals
self.mos = np.zeros((len(self.so.mol.basis), len(self.so.mol.basis)),
dtype=np.float64)
mo_idx = 0
for g,irreplist in zip(self.groups,self.irreps):
for j,irrep in enumerate(irreplist): # loop over irreps

# early exit if there are no irreps of this type
if irrep == 0:
continue

irrep_e = self.ct.get_character(j,0) # get dimensionality of irrep

irrepres = np.zeros((int(irrep * irrep_e), len(self.so.mol.basis)))
for k in range(int(irrep * irrep_e)): # loop over times the irrep is present
#label = self.ct.get_label_irrep(j)
res = self.apply_projection_operator(g[k], j)
irrepres[k,:] = res

# perform orthogonalization if the irrep space is larger than 1
if len(irrepres) > 1:
S = irrepres @ irrepres.transpose()
e,v = np.linalg.eigh(S)
X = v @ np.diag(1.0 / np.sqrt(e)) @ v.transpose()
irrepres = X @ irrepres

for mo in irrepres:
self.mos[mo_idx] = mo
mo_idx += 1
else:
self.mos[mo_idx] = irrepres[0,:]
mo_idx += 1

return self.mos


def apply_projection_operator(self, bf_idx, irrep_idx):
"""
Apply the projection operator to a basis function
"""
res = np.zeros(len(self.so.mol.basis))
for i,m in enumerate(self.so.operation_matrices):
res += self.ct.get_character(irrep_idx, i) * m[bf_idx,:]

# return result and normalize it
return np.array(res, dtype=np.float64) / np.linalg.norm(res)
10 changes: 9 additions & 1 deletion tests/test_ammonia_symmetry_operations.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
sys.path.append(os.path.join(os.path.dirname(__file__), '..'))

# import functions
from sphecerix import Molecule, BasisFunction, SymmetryOperations
from sphecerix import Molecule, BasisFunction, SymmetryOperations, \
visualize_matrices, CharacterTable

class TestAmmoniaSymmetryOperations(unittest.TestCase):
"""
Expand Down Expand Up @@ -46,6 +47,13 @@ def test_nh3(self):
result = np.load(os.path.join(os.path.dirname(__file__), 'results', 'nh3.npy'))

np.testing.assert_almost_equal(symops.operation_matrices, result)

# assess LOT results
ct = CharacterTable('c3v')
irreps = ct.lot(np.trace(symops.operation_matrices, axis1=1, axis2=2))
np.testing.assert_almost_equal(irreps, [4,0,2])



if __name__ == '__main__':
unittest.main()

0 comments on commit 03b23e5

Please sign in to comment.