#!/usr/bin/env python
# Copyright 2020-2022 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Xing Zhang <zhangxing.nju@gmail.com>
#
from abc import ABC, abstractmethod
import numpy as np
from pyscf import lib
from pyscf.pbc.symm import geom
from pyscf.pbc.symm.tables import SchoenfliesNotation
def _round_zero(a, tol=1e-9):
a[np.where(abs(a) < tol)] = 0
return a
[docs]
class GroupElement(ABC):
'''
The abstract class for group elements.
'''
def __call__(self, other):
return self.__matmul__(other)
@abstractmethod
def __matmul__(self, other):
pass
def __mul__(self, other):
assert isinstance(other, self.__class__)
return self @ other
@abstractmethod
def __hash__(self):
pass
[docs]
@abstractmethod
def inv(self):
'''
Inverse of the group element.
'''
pass
[docs]
class PGElement(GroupElement):
'''
The class for crystallographic point group elements.
The group elements are rotation matrices represented
in lattice translation vector basis.
Attributes:
matrix : (d,d) array of ints
Rotation matrix in lattice translation vector basis.
dimension : int
Dimension of the space: `d`.
'''
def __init__(self, matrix):
self.matrix = matrix
self.dimension = matrix.shape[0]
def __matmul__(self, other):
if not isinstance(other, PGElement):
raise TypeError(f"{other} is not a point group element.")
return PGElement(np.dot(self.matrix, other.matrix))
def __repr__(self):
return self.matrix.__repr__()
def __hash__(self):
def _id(op):
s = op.flatten() + 1
return lib.inv_base_repr_int(s, op.shape[0])
r = _id(self.rot)
# move identity to the first place
d = self.dimension
r -= _id(np.eye(d, dtype=int))
if r < 0:
r += _id(np.ones((d,d), dtype=int)) + 1
return int(r)
[docs]
@staticmethod
def decrypt_hash(h, dimension=3):
if dimension == 3:
id_eye = int('211121112', 3)
id_max = int('222222222', 3)
elif dimension == 2:
id_eye = int('2112', 3)
id_max = int('2222', 3)
else:
raise NotImplementedError
r = h + id_eye
if r > id_max:
r -= id_max + 1
#s = np.base_repr(r, 3)
#s = '0'*(dimension**2-len(s)) + s
#rot = np.asarray([int(i) for i in s]) - 1
rot = np.asarray(lib.base_repr_int(r, 3, dimension**2)) - 1
rot = rot.reshape(dimension,dimension)
# sanity check
#element = PGElement(rot)
#assert hash(element) == h
return rot
def __lt__(self, other):
if not isinstance(other, PGElement):
raise TypeError(f"{other} is not a point group element.")
return self.__hash__() < other.__hash__()
@property
def matrix(self):
return self._matrix
@matrix.setter
def matrix(self, a):
self._matrix = a
rot = matrix
[docs]
def inv(self):
return PGElement(np.asarray(np.linalg.inv(self.matrix), dtype=np.int32))
[docs]
class FiniteGroup(ABC):
'''
The class for finite groups.
Attributes:
elements : list
Group elements.
order : int
Group order.
'''
def __init__(self, elements, from_hash=False):
if from_hash:
elements = self.__class__.elements_from_hash(elements)
self.elements = np.asarray(elements)
self._order = None
self._hash_table = None
self._inverse_table = None
self._multiplication_table = None
self._conjugacy_table = None
self._conjugacy_mask = None
self._chartab_full = None
self._chartab = None
self._group_name = None
self._group_index = None
self._check_sanity()
[docs]
@staticmethod
@abstractmethod
def elements_from_hash(hashes, **kwargs):
pass
def __len__(self):
return self.order
def __getitem__(self, i):
return self.elements[i]
def __and__(self, other):
if type(self) is not type(other):
raise TypeError(f'{self} and {other} must be the same type')
if self is other:
return self
hi = list(self.hash_table.keys())
hj = list(other.hash_table.keys())
hij = np.intersect1d(hi, hj)
return self.__class__(hij, from_hash=True)
def __or__(self, other):
if type(self) is not type(other):
raise TypeError(f'{self} and {other} must be the same type')
if self is other:
return self
hi = list(self.hash_table.keys())
hj = list(other.hash_table.keys())
hij = np.union1d(hi, hj)
return self.__class__(hij, from_hash=True)
[docs]
def issubset(self, other):
if type(self) is not type(other):
raise TypeError(f'{self} and {other} must be the same type')
if self is other:
return True
hi = list(self.hash_table.keys())
hj = list(other.hash_table.keys())
return set(hi).issubset(set(hj))
@property
def order(self):
if self._order is None:
self._order = len(self.elements)
return self._order
@order.setter
def order(self, n):
self._order = n
@property
def hash_table(self):
'''
Hash table for group elements: {hash : index}.
'''
if self._hash_table is None:
self._hash_table = {hash(g) : i for i, g in enumerate(self.elements)}
return self._hash_table
@hash_table.setter
def hash_table(self, table):
self._hash_table = table
@property
def inverse_table(self):
'''
Table for inverse of the group elements.
Return : (n,) array of ints
The indices of elements.
'''
if self._inverse_table is None:
_table = [self.hash_table[hash(g.inv())] for g in self.elements]
self._inverse_table = np.asarray(_table)
return self._inverse_table
@inverse_table.setter
def inverse_table(self, table):
self._inverse_table = table
@property
def multiplication_table(self):
'''
Multiplication table of the group.
Return : (n, n) array of ints
The indices of elements.
'''
if self._multiplication_table is None:
prod = self.elements[:,None] * self.elements[None,:]
_table = [self.hash_table[hash(gh)] for gh in prod.flatten()]
self._multiplication_table = np.asarray(_table).reshape(prod.shape)
return self._multiplication_table
@multiplication_table.setter
def multiplication_table(self, table):
self._multiplication_table = table
@property
def conjugacy_table(self):
'''
conjugacy_table[`index_g`, `index_x`] returns the index of element `h`,
where :math:`h = x * g * x^{-1}`.
'''
if self._conjugacy_table is None:
prod_table = self.multiplication_table
g_xinv = prod_table[:,self.inverse_table]
self._conjugacy_table = prod_table[np.arange(self.order)[None,:], g_xinv]
return self._conjugacy_table
@conjugacy_table.setter
def conjugacy_table(self, table):
self._conjugacy_table = table
@property
def conjugacy_mask(self):
'''
Boolean mask array indicating whether two elements
are conjugate with each other.
'''
if self._conjugacy_mask is None:
n = self.order
is_conjugate = np.zeros((n,n), dtype=bool)
is_conjugate[np.arange(n)[:,None], self.conjugacy_table] = True
self._conjugacy_mask = is_conjugate
return self._conjugacy_mask
[docs]
def conjugacy_classes(self):
'''
Compute conjugacy classes.
Returns:
classes : (n_irrep,n) boolean array
The indices of `True` correspond to the
indices of elements in this class.
representatives : (n_irrep,) array of ints
Representive elements' indices in each class.
inverse : (n,) array of ints
The indices to reconstruct `conjugacy_mask` from `classes`.
'''
_, idx = np.unique(self.conjugacy_mask, axis=0, return_index=True)
representatives = np.sort(idx.ravel())
classes = self.conjugacy_mask[representatives]
inverse = -np.ones((self.order), dtype=int)
diff = (self.conjugacy_mask[None,:,:]==classes[:,None,:]).all(axis=-1)
for i, a in enumerate(diff):
inverse[np.where(a)[0]] = i
assert (inverse >= 0).all()
assert (classes[inverse] == self.conjugacy_mask).all()
return classes, representatives, inverse
[docs]
def character_table(self, return_full_table=False, recompute=False):
'''
Character table of the group.
Args:
return_full_table : bool
If True, return the characters for all elements.
recompute : bool
Whether to recompute the character table. Default is False,
meaning to use the cached table if possible.
Returns:
chartab : array
Character table for classes.
chartab_full : array, optional
Character table for all elements.
'''
if not recompute:
if not return_full_table and self._chartab is not None:
return self._chartab
if return_full_table and self._chartab_full is not None:
return self._chartab_full
classes, _, inverse = self.conjugacy_classes()
class_sizes = classes.sum(axis=1)
ginv_h = self.multiplication_table[self.inverse_table]
M = classes @ np.random.rand(self.order)[ginv_h] @ classes.T
M /= class_sizes
_, Rchi = np.linalg.eig(M)
chi = Rchi.T / class_sizes
norm = np.sum(np.abs(chi) ** 2 * class_sizes[None,:], axis=1) ** 0.5
chi = chi / norm[:,None] * self.order ** 0.5
chi /= (chi[:, 0] / np.abs(chi[:, 0]))[:,None]
chi = np.round(chi, 9)
chi_copy = chi.copy()
chi_copy[:,1:] *= -1
idx = np.lexsort(np.rot90(chi_copy))
chi = chi[idx]
chi = _round_zero(chi)
self._chartab = chi
self._chartab_full = chi[:, inverse]
if return_full_table:
return self._chartab_full
else:
return self._chartab
[docs]
def project_chi(self, chi, other):
'''
Project characters to another group.
'''
if self is other:
return chi
i_ind, j_ind = self.get_elements_map(other)
chi_j = np.zeros((other.order,))
chi_j[j_ind] = chi[i_ind]
return chi_j
[docs]
def get_elements_map(self, other):
if not (other.issubset(self) or self.issubset(other)):
raise KeyError(f'{self} or {other} must be a subset of the other.')
hi = [hash(g) for g in self.elements]
hj = [hash(g) for g in other.elements]
_, i_ind, j_ind = np.intersect1d(hi, hj, return_indices=True)
return i_ind, j_ind
[docs]
def get_irrep_chi(self, ir):
chartab = self.character_table(True)
return chartab[ir]
def _check_sanity(self):
try:
# check duplication
assert len(self.hash_table) == self.order
# check inversion
assert (abs(np.sort(self.inverse_table)
-np.arange(self.order)).max() == 0)
# check multiplication
assert (abs(np.sort(self.multiplication_table, axis=-1)
-np.arange(self.order)[None,:]).max() == 0)
except (AssertionError, KeyError, ValueError):
raise ValueError('The elements do not form a group.')
[docs]
class PointGroup(FiniteGroup):
'''
The class for crystallographic point groups.
'''
[docs]
def group_name(self, notation='international'):
if self._group_name is not None and notation=='international':
return self._group_name
name = geom.get_crystal_class(None, self.elements)[0]
self._group_name = name
if notation.lower().startswith('scho'): # Schoenflies
name = SchoenfliesNotation[name]
return name
@property
def group_index(self):
if self._group_index is None:
name = self.group_name()
self._group_index = list(SchoenfliesNotation.keys()).index(name)
return self._group_index
[docs]
@staticmethod
def elements_from_hash(hashes, dimension=3):
elements = [PGElement(PGElement.decrypt_hash(h, dimension)) for h in hashes]
return elements
[docs]
class Representation:
'''
Helper class for representation reductions.
Only characters are stored at the moment.
'''
def __init__(self, group, rep=None, chi=None):
self.group = group
self.rep = rep
self.chi = chi
@property
def rep(self):
if self._rep is None:
self._rep = self.chi_to_rep(self.chi)
return self._rep
@rep.setter
def rep(self, value):
self._rep = value
@property
def chi(self):
if self._chi is None:
self._chi = self.rep_to_chi(self.rep)
return self._chi
@chi.setter
def chi(self, value):
self._chi = value
[docs]
def rep_to_chi(self, rep):
chartab = self.group.character_table(True)
chi = np.einsum('ni,n->i', chartab, rep)
return chi
[docs]
def chi_to_rep(self, chi):
group = self.group
chartab = group.character_table(True)
assert len(chi) == chartab.shape[1]
nA = np.einsum('ni,i->n', chartab.conj(), chi) / group.order
assert (abs(nA - nA.round()) < 1e-9).all()
nA = np.rint(nA).astype(int)
return nA
def __matmul__(self, other):
g1, chi1 = self.group, self.chi
g2, chi2 = other.group, other.chi
g12 = g1 & g2
chi1_proj = g1.project_chi(chi1, g12)
chi2_proj = g2.project_chi(chi2, g12)
chi12_proj = chi1_proj * chi2_proj
return self.__class__(g12, chi=chi12_proj)
