Skip to content

Commit

Permalink
Polishing algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@ifilot
ifilot committed Jul 1, 2023
1 parent 592c571 commit a4a68c9
Show file tree
Hide file tree
Showing 7 changed files with 170 additions and 29 deletions.
97 changes: 97 additions & 0 deletions examples/ccl4.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
# -*- coding: utf-8 -*-

import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# 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, \
plot_matrix

def main():
mol = Molecule()
mol.add_atom('C', 0.0, 0.0, 0.0, unit='angstrom')
mol.add_atom('Cl', 1,1,1, unit='angstrom')
mol.add_atom('Cl', 1,-1,-1, unit='angstrom')
mol.add_atom('Cl', -1,1,-1, unit='angstrom')
mol.add_atom('Cl', -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)],
'Cl': [BasisFunction(1,0,0),
BasisFunction(2,1,1),
BasisFunction(2,1,-1),
BasisFunction(2,1,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=(36,24), 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(verbose=True)

fig, ax = plt.subplots(1,1,dpi=144,figsize=(20,20))
plot_matrix(ax, mos,[bf.name for bf in symops.mol.basis],None,0)

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=(36,24), numcols=6,
highlight_groups=po.get_blocks())

if __name__ == '__main__':
main()
2 changes: 1 addition & 1 deletion examples/ch4.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ def main():
[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,
highlight_groups=po.get_block_sizes())
highlight_groups=po.get_blocks())

if __name__ == '__main__':
main()
15 changes: 7 additions & 8 deletions examples/dodecahedrane.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,11 @@ def main():
mol.from_file('molecules/dodecahedrane.xyz')

molset = {
'C': [BasisFunction(1,0,0),
'C': [#BasisFunction(1,0,0),
#BasisFunction(2,0,0),
#BasisFunction(2,1,1),
#BasisFunction(2,1,-1),
#BasisFunction(2,1,0)],
],
BasisFunction(2,1,1),
BasisFunction(2,1,-1),
BasisFunction(2,1,0)],
'H': [BasisFunction(1,0,0)]
}
mol.build_basis(molset)
Expand Down Expand Up @@ -110,21 +109,21 @@ def main():

# apply projection operator
po = ProjectionOperator(ct, symops)
mos = po.build_mos()
mos = po.build_mos(verbose=True)
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)

fig, ax = plt.subplots(1,2,dpi=144,figsize=(25,15))
fig, ax = plt.subplots(1,2,dpi=144,figsize=(45,45))

plot_matrix(ax[0], symops.operation_matrices[5],[bf.name for bf in symops.mol.basis],None,0)

plot_matrix(ax[1], newmats[5],
['$\phi_{%i}$' % (i+1) for i in range(len(symops.mol.basis))],
None,0,highlight_groups=po.get_block_sizes())
None,0,highlight_groups=po.get_blocks())

plt.tight_layout()

Expand Down
3 changes: 2 additions & 1 deletion examples/nh3.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,8 @@ def main():
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))
figsize=(9,6),
highlight_groups=po.get_blocks())

if __name__ == '__main__':
main()
14 changes: 12 additions & 2 deletions sphecerix/character_table.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,15 +16,25 @@ def __init__(self, name):
self.nrclasses = len(self.chartablelib['classes'])
self.nrgroups = len(self.chartablelib['symmetry_groups'])

# build expanded table
self.expandedtable = np.zeros((self.nrgroups, self.order), dtype=np.int64)
# build condensed and expanded table
self.table = np.zeros((self.nrgroups, self.nrgroups), dtype=np.float64)
self.expandedtable = np.zeros((self.nrgroups, self.order), dtype=np.float64)
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']):
self.table[i,j] = g['characters'][j]
for k in range(c['multiplicity']):
self.expandedtable[i,idx] = g['characters'][j]
idx += 1

# check that the character table is correct
try:
np.testing.assert_almost_equal(self.expandedtable @ self.expandedtable.transpose() / self.order,
np.identity(self.nrgroups))
except AssertionError as e:
print('Invalid character table encountered in: %s' % filename)
raise e

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

Expand Down
33 changes: 25 additions & 8 deletions sphecerix/matrixplot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,36 +28,53 @@ def visualize_matrices(matrices, opnames, labels, numcols = 3,
plt.savefig(filename)
plt.close()

def plot_highlight_groups(axh, groups, mat):
def plot_highlight_groups(axh, blocks, mat):
# add semitransparent hash
sumblocks = np.sum(g[0] for g in blocks)
cum = 0
for g in groups:
rect = patches.Rectangle((cum - 0.5, cum - 0.5), g, g,
for g in blocks:
rect = patches.Rectangle((cum - 0.5, cum - 0.5), g[0], g[0],
linewidth=1,
zorder=5,
fill = None,
hatch='///',
alpha=0.5)
axh.add_patch(rect)
cum += g
cum += g[0]

# add red outline
cum = 0
for g in groups:
rect = patches.Rectangle((cum - 0.5, cum - 0.5), g, g,
for g in blocks:
rect = patches.Rectangle((cum - 0.5, cum - 0.5), g[0], g[0],
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])),
axh.text(cum+g[0]/2-0.5, cum+g[0]/2-0.5, '%i' % round(np.trace(mat[cum:cum+g[0],cum:cum+g[0]])),
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
k = (cum+1) > sumblocks/2

if k:
label = r'%i $\times$ %s $\rightarrow$' % (g[1],g[2])
else:
label = r'$\leftarrow$ %i $\times$ %s' % (g[1],g[2])

axh.text(cum+(-0.75 if k else g[0] - 0.25), cum+g[0]/2-0.5,
label,
fontsize=6,
color='red',
horizontalalignment='right' if k else 'left',
verticalalignment='center',
bbox=dict(boxstyle="round", ec=(1., 0.5, 0.5), fc=(1., 0.8, 0.8), ),
zorder=6)

cum += g[0]


def plot_matrix(ax, mat, labels, title = None, xlabelrot = 0,
Expand Down
35 changes: 26 additions & 9 deletions sphecerix/projection_operator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

import numpy as np
from collections import Counter
import random

class ProjectionOperator:

Expand Down Expand Up @@ -49,14 +50,15 @@ def collect(self):
self.irreps = []
self.irreplabels = []
self.block_sizes = []
self.blocks = []
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))

for j,irrepdim in enumerate(self.irreps[-1]):
for k in range(int(irrepdim)):
# get character under E and store this
self.block_sizes.append(self.ct.chartablelib['symmetry_groups'][j]['characters'][0])
if irrepdim > 0:
self.block_sizes.append(int(irrepdim * self.ct.chartablelib['symmetry_groups'][j]['characters'][0]))
self.blocks.append((self.block_sizes[-1], irrepdim, self.ct.get_label_irrep(j)))

def build_mos(self, verbose=False):
# check if groups have been collected
Expand All @@ -69,6 +71,8 @@ def build_mos(self, verbose=False):
mo_idx = 0
for g,irreplist in zip(self.groups,self.irreps): # loop over the groups

bf_idx = 0

if verbose:
print('Group: ', g)
print('Irreps:')
Expand All @@ -83,19 +87,29 @@ def build_mos(self, verbose=False):
label = self.ct.get_label_irrep(j)
print(' - %s: %i' % (label, irrep))


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
zero_modes = True

while zero_modes: # keep on looping until no zero eigenvalues are found
bf_indices = random.sample(g, int(irrep * irrep_e))
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(bf_indices[k], j)
irrepres[k,:] = res
bf_idx += 1

S = irrepres @ irrepres.transpose()
e,v = np.linalg.eigh(S)

#print('Trying: ', e)
zero_modes = np.any(np.abs(e) < 1e-2)

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

# for symmetric orthogonalization:
# X = v @ np.diag(1.0 / np.sqrt(e)) @ v.transpose()
Expand Down Expand Up @@ -135,6 +149,9 @@ def get_block_sizes(self):
"""
return self.block_sizes

def get_blocks(self):
return self.blocks

def check_transformation(self):
return
newmats = [self.mos @ m @ self.mos.transpose() for m in self.so.operation_matrices]
Expand Down

0 comments on commit a4a68c9

Please sign in to comment.