Symmetry (`qcip_tools.symmetry`)

Handle symmetry.

Introduction

What’s a group anyway ?

A group is a set \(G\), together with a binary operation \(\star\) (group law), form a group \((G, \star)\) if it satisfies 4 requirements:

Closure \(\forall a, b \in \mathcal{G}: c=a\star b \land c\in\mathcal{G}\) ;
Associativity: \(\forall a, b, c \in \mathcal{G}: (a\star b)\star c = a\star (b\star c)\)
(unique) Identity: \(\exists! e\in\mathcal{G},\forall a\in\mathcal{G}: a\star e = a = e\star a\);
(unique) Inverse: \(\forall a\in\mathcal{G}:\exists! a'\in\mathcal{G}: a\star a' = e\).

In the code, the binary operation is represented by BinaryOperation, which works on a given codomain and uses a surjective function:

from qcip_tools import symmetry

f = symmetry.BinaryOperation(
    symmetry.Set(range(4)),  # domain is {0, 1, 2, 3}
    lambda e: (e[0] + e[1]) % 4 # function is modulo 4
)

g = symmetry.Group(f)  # cyclic group of order 4

e = g.identity()  # there is an identity element (here 0)
inv1 = g.inverse(1)  # the inverse of "1" is "3"

# and the inverse of 1 times 1 is identity
assert symmetry.GroupElement(1, g) * inv1 == e

This implementation also compute the different classes (which includes, at least, an element and its inverse) and the Cayley table (multiplication table).

Point group

By using symmetry operation as group element (and matrix multiplication as the binary operation), one can create a point group, which apply on finite objects (in oposition to space group, which aply to infinite objects).

Symmetry operation are of two types, basically:

Proper rotation of \(k\times \frac{2\pi}{n}\) around an axis, noted \(C_n^k\).
Improper rotation of \(k\times \frac{2\pi}{n}\) around an axis, noted \(S_n^k\) (a rotation of \(\frac{2\pi}{n}\) followed by a reflexion, repeated \(k\) times).

From that, one can generate:

Identity, \(E=C_n^n\) ;
Reflexion, \(\sigma=S_1\) ;
Inversion, \(i = S_2^1\).

Those operation are generated trough function of the Operation class (internally, it uses quaternions). One can generate simply any point group by using the methods of PointGroup class.

C_3v = symmetry.PointGroup.C_nv(3)  # generate C3v

self.assertEqual(C_3v.conjugacy_classes[0], {C_3v.e})  # first class is always identity
self.assertEqual(  # second class contains the proper rotation
    {e.element for e in C_3v.conjugacy_classes[1]}, {symmetry.Operation.C(3), symmetry.Operation.C(3, 2)})

# third class contains the reflexion planes (and, among other, the one that sets x)
self.assertIn(symmetry.Operation.sigma(axis=numpy.array([0, 1., 0])), C_3v.conjugacy_classes[2])

Character table

An useful tools when working with group theory is the character table, which contains irreducible representations of character for each class of the group.

C_2v = symmetry.PointGroup.C_nv(2)
t = C_2v.character_table()

A1 = t.irreducible_representations[0]  # A1 (trivial) representation
A2 = t.irreducible_representations[1]  # A2
r = A1 + A2  # reducible representation contains A1 + A2
assert A2 in A1 * A2 # direct product gives A2

One can also do print(t) to get the formated output:

C_2v (h=4)
------------------------------------------------------
       E           C_2         σ_v         σ_v
------------------------------------------------------
A1     1.00        1.00        1.00        1.00
A2     1.00        1.00       -1.00       -1.00
B1     1.00       -1.00        1.00       -1.00
B2     1.00       -1.00       -1.00        1.00
------------------------------------------------------

Discover symmetry

Use SymmetryFinder. In particular, the documentation of find_symmetry() describes the way the point group is found.

Warning

This was not tested for all point groups. In particular, it detects icosahedral point groups correctly, but outputs a wrong orientation.

In the case of molecule, see qcip_tools.molecule.MolecularSymmetryFinder.

API documentation

Symmetry handling.

See https://en.wikipedia.org/wiki/Group_(mathematics)

Base implementation inspired by https://github.com/naftaliharris/Abstract-Algebra/

and https://pdfs.semanticscholar.org/4edd/1ac673ea4cab251bb0b64bf0f5b65f18cc9d.pdf

class qcip_tools.symmetry.BinaryOperation(codomain, func)

Define a binary operation \(f\) such that,

\[f:S\times S \rightarrow S.\]

Thus, we only define a binary operation trough its codomain.

See https://en.wikipedia.org/wiki/Binary_operation

check_associativity()

Check if the binary operation is associative for all the element of the domain

Return type:: bool

check_closure()

Check that the relation gives elements that are in the codomain

Return type:: bool

check_surjectivity()

Check surjectivity (image equals codomain)

Return type:: bool

classmethod generate(generators, func)

Create a complete set out of the generators.

Algorithm from https://physics.stackexchange.com/a/351400j, simplified version of http://pymatgen.org/_modules/pymatgen/symmetry/analyzer.html

:param generators; the generators of the set :type generators: list :param func: function :type func: function :rtype: BinaryOperation

image()

Get the image set (which is part of the codomain)

Return type:: Set

class qcip_tools.symmetry.CharacterTable(group, table)

List of irreducible representations

LIMIT = 1e-05

property size

exception qcip_tools.symmetry.CharacterTableError

class qcip_tools.symmetry.Group(binary_operation)

Abstract group structure.

A set \(G\), together with an operation \(\star\) (group law), form a group \((G, \star)\) if it satisfies 4 requirements:

Closure (not checked here) ;
Associativity (not checked here) ;
(unique) Identity ;
(unique) Inverse.

character_table(table=None)

Compute the character table

Parameters:: table (numpy.ndarray) – character table
Return type:: CharacterTable

class_matrices()

Get all the class matrices

Return type:: list

class_matrix(class_number)

Compute class matrix: given the \(k\) classes \(C_i\) (\(0\leq i < k\)),

For each \(x\in C_r\), compute \(y=x^{-1}z\), with zin C_t and \(y\in C_s\). Add 1 to the \((s,t)\) component of \(M^r\), a \(k\times k\) matrix. Here, class_number gives the value of \(r\).

This results in a matrix where each component is the number of solution \((y,z)\) for which \(xy=z\).

Parameters:: class_number (int) – class id
Return type:: numpy.ndarray

conjugacy_class(g)

Get the elements of the conjugacy class of g

Parameters:: g (GroupElement) – a group element

gen_character_table(tweak_matrices_func=<function Group.<lambda>>)

Generate a character table for the group.

Parameters:: tweak_matrices_func (function) – function to alter the list of class matrices (if needed): takes a list of matrices as input and returns the new list. Useful in some cases (especially if the classes are sorted).
Returns:: the character table, with the column ordered as the different class, and the line in a random order
Type:: numpy.ndarray

identity()

Get identity

Returns:: an element of G

inverse(element)

Return an inverse of a given element of G

Parameters:: element – the element
Returns:: another element of G

class qcip_tools.symmetry.GroupElement(element, group): Abstract group element: syntaxic sugar

class qcip_tools.symmetry.IrreducibleRepresentation(characters, table, label=None)

An irreducible representation (irrep.)

property degree

dot(other)

Compute the scalar product between two representations, \(\chi^\lambda\) (self) and \(\chi^\mu\) (other), as

\[\frac{1}{\#\mathcal G}\sum_{r\in \mathcal G} [\chi^\lambda(r)]^\star \chi^\mu(r).\]

Parameters:: other (IrreducibleRepresentation) – other irreducible representation
Returns:: scalar product
Return type:: complex

exception qcip_tools.symmetry.NotInGroup

class qcip_tools.symmetry.Operation(q, improper=False, description=None)

Define a symmetry element (fully based on quaternion).

According Mamone et al.,

\[S^k_n = \hat{R}\left(\pi + \frac{2\,k\,\pi}{n}\right),\]

where \(\hat{R}\) is the normal rotation operator. To use this operator on a point, the opposite of the sandwich product must be taken, so that the conjugate of the quaternion that represent the point is actually used.

classmethod C(n, k=1, axis=array([0., 0., 1.]))

Generate the representation of a proper rotation

Parameters:

n (int) – order of the axis
k (int) – number of step
axis (numpy.ndarray) – to which axis this rotation axis should be parallel ?

Return type:

Operation

classmethod E()

Generate the identity

Return type:: Operation

classmethod S(n, k=1, axis=array([0., 0., 1.]))

Generate the representation of a improper rotation

Parameters:

n (int) – order of the axis
k (int) – number of step
axis (numpy.ndarray) – to which axis this rotation axis should be parallel ?

Return type:

Operation

apply(pin)

Apply the operation on a point

Parameters:: pin (numpy.ndarray) – the point
Return type:: numpy.ndarray

classmethod from_axangle(axis, angle=6.283185307179586, improper=False, in_degree=False, description=None)

Create from axe and angle

Parameters:

axis (numpy.ndarray) – the rotation axe
angle (float) – the rotation angle
improper (bool) – is it an improper rotation ?
in_degree (bool) – is the angle in degree
description (OperationDescription) – description of the operation

Return type:

Operation

get_description(force=False)

Try to recognize the transformation:

Extract axis and angle from quaternion (using quat2axangle()) ;
Treat angle (\(\theta\)):
- if the axis have a negative z, take the opposite and set \(\theta = 2\pi-\theta\),
- if this is an improper rotation, set \(\theta=\theta-\pi\),
- Cast between 0 and \(4\pi\),
- divide by \(2\pi\) (so that the result is clamped between 0 and 2).
Try to extract a simple fraction (thus \(k\) and \(n\)) out of this number ;
Further reduce k to be smaller than \(n\) if proper rotation or odd \(n\) ;
Recognize the specific \(n\) cases (1 or 2).

Parameters:: force (bool) – use self.description if False, or try a new recognition otherwise
Return type:: OperationDescription

classmethod i()

Generate an inversion center

Return type:: Operation

matrix_representation(): Get a 3x3 matrix representation of the symmetry operation

classmethod sigma(axis=array([0., 0., 1.]))

Generate a symmetry plane

Parameters:: axis (numpy.ndarray) – normal axis to the plane
Return type:: Operation

class qcip_tools.symmetry.OperationDescription(symbol, axis=array([0., 0., 0.]), n=1, k=1)

Store the description of a symmetry operation

static find_textual_axis(axis)

Try to find the main axis:

Find sole axis (x, y, z) ;
Find the axis at 45° with respect to two axis ;
Find the axis at 45° with respect to three axis (the 8 apexes of a cube).

Some example of the different axis:

    (x'yz') +------*-------+ (xyz')
           /              /|
    (x'y) *      o (y)   * |
         /              /  * (xz')
 (x'yz) +------*-------+   |
        |              | o |
        |              |   + (xy',z')    y
  (x'z) *      o (z)   *  /              |
        |              | * (y'z)         +--- x
        |              |/               /
(x'y'z) +------*-------+ (xy'z)        z

with

o : single axis,
* : double axis,
+ : triple axis.

Parameters:: axis (numpy.ndarray) – the axis
Return type:: str|None

class qcip_tools.symmetry.OperationType(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)

Symbol of the symmetry element, in Schoenflies notation.

identity = 'E'

improper_rotation = 'S'

inversion = 'i'

proper_rotation = 'C'

reflexion_plane = 'σ'

class qcip_tools.symmetry.PointGroup(func, description=None)

Concrete implementation of a point group

classmethod C_n(n=1)

Create a cyclic group of order n

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod C_nh(n=1)

Create a reflexion point group of order n:

\[C_{nv} = C_n \times \sigma_h.\]

Note that \(C_{1h}\) is usually noted \(C_s\).

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod C_nv(n=2)

Create a pyramidal group of order n:

\[C_{nv} = C_n \times \sigma_v(y).\]

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod D_n(n=2)

Create a dihedral point group of order n:

\[D_{n} = C_n(z) \times C_2(x).\]

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod D_nd(n=2)

Create an anti prismatic point group of order n:

\[D_{nd} = S_{2n}(z) \times \sigma_d(y).\]

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod D_nh(n=2)

Create a prismatic point group of order n:

\[\begin{split}D_{nh} &= D_n \times \sigma_h \\ &= C_n(z) \times C_2(x) \times\sigma_h.\end{split}\]

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod I()

Create the chiral octahedral group

\[I = C_5(1, 0, \phi) \times C_3(1, 1, 1) \times C_2(x), \text{ with } \phi = \frac{1+\sqrt 5}{2}.\]

Return type:: PointGroup

classmethod I_h()

Create the chiral octahedral group

\[I_h = C_5(1, 0, \phi) \times C_3(1, 1, 1) \times C_2(x) \times \sigma_h, \text{ with } \phi = \frac{1+\sqrt 5}{2}.\]

Return type:: PointGroup

classmethod O()

Create the chiral octahedral group

\[O = C_4(z) \times C_3(1, 1, 1).\]

Return type:: PointGroup

classmethod O_h()

Create the achiral octahedral group

\[O_h = C_4 \times \sigma_h \times C_3(1, 1, 1).\]

Return type:: PointGroup

classmethod S_n(n=1)

Create a reflexion (odd n, \(C_{(2k+1)h}\)) or improper rotation (even n, \(S_{2k}\)) point group of order n.

Note that \(S_{2}\) is usually noted \(C_i\).

Parameters:: n (int) – order of the group
Return type:: PointGroup

classmethod T()

Create the chiral tetrahedral group

\[T = C_2(z) \times C_3(1, 1, 1).\]

Return type:: PointGroup

classmethod T_d()

Create the achiral tetrahedral group

\[T_d = S_4 \times C_3(1, 1, 1).\]

Return type:: PointGroup

classmethod T_h()

Create the pyritohedral group

\[T_h = C_2 \times \sigma_h \times C_3(1, 1, 1).\]

Return type:: PointGroup

character_table(table=None)

Generate a character table for the group.

Warning

This uses the simultaneous_diagonalization() function, but for unknown reasons, it “fails” for \(D_{4h}\), \(O_{h}\), \(T_{h}\) and \(I_{h}\) (it generates eigenvectors that are not representations, but are probably orthogonal with each other).

Thus, with the tweak_matrices_func parameter,

For \(D_{4h}\) it uses list(reversed(matrices)) ;
For \(T_{h}\) it uses list(reversed(matrices)) ;
For \(O_{h}\) it uses list(reversed(matrices)) with del mat[0] ;
For \(I_{h}\) it uses list(reversed(matrices)) with del mat[0] and del mat[1].

That’s the best I can do for the moment ;)

Return type:: CharacterTable

classmethod generate(generators, description=None)

static product(e)

class qcip_tools.symmetry.PointGroupDescription(symbol, n=0)

gen_point_group(lowers_infinite=-1)

Generate the point group out of the description.

Not able to generate infinite groups \(D_{\infty h}\) or \(C_{\infty v}\), except if lowers_infinite is set.

Parameters:: lowers_infinite (int) – axis order used instead of \(\infty\).
Return type:: PointGroup

exception qcip_tools.symmetry.PointGroupError

class qcip_tools.symmetry.PointGroupType(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)

antiprismatic = 'D_nd'

cyclic = 'C_n'

dihedral = 'D_n'

icosahedral_achiral = 'I_h'

icosahedral_chiral = 'I'

improper_rotation = 'S_n'

octahedral_achiral = 'O_h'

octahedral_chiral = 'O'

prismatic = 'D_nh'

pyramidal = 'C_nv'

pyritohedral = 'T_h'

reflexion = 'C_nh'

tetrahedral_achiral = 'T_d'

tetrahedral_chiral = 'T'

class qcip_tools.symmetry.ReducibleRepresentation(table, irreducible_representations=None, characters=None)

Sum of irreducible representation

Reduction of \(\chi^{(red)}\) into a sum of irreducible representation of \(n^\lambda\) times a irreducible representation \(\chi^\lambda\) is performed using

\[n^\lambda = \frac{1}{\#\mathcal G} \sum_{r\in\mathcal G} [\chi^{(red)}(r)]^\star \chi^\lambda(r)\]

property characters

property size

exception qcip_tools.symmetry.RepresentationError

class qcip_tools.symmetry.Set: Define a (finite) set of (unique) elements

class qcip_tools.symmetry.SymmetryFinder(points, tol=1e-05)

Find the symmetry

Inspired by https://github.com/sunqm/pyscf/blob/master/pyscf/symm/geom.py and http://pymatgen.org/_modules/pymatgen/symmetry/analyzer.html

Parameters:

points (numpy.ndarray) – an Nx4 array of N points with (Z, x, y, z) for each point. Z is a label for points that are equivalent.
tol (float) – tolerance threshold

static cartesian_tensor(tensor, n)

Compute multipole expansion of cartesian tensor of order n (what you find in the right column of a character table) and get its principal components.

https://en.wikipedia.org/wiki/Multipole_expansion

Parameters:

tensor (numpy.ndarray) – tensor
n (int) – order

Returns:

eigenvalue and eigenvectors of the tensor

static degeneracies(e, decimals=-5)

find_best_x_axis(principal_axis, mirrors_or_C2, are_mirrors=False)

Find best x axis among the list of mirrors or C2: the one that goes by the largest number of atoms

Parameters:

principal_axis (numpy.ndarray) – main axis
mirrors_or_C2 (list) – axes to check
are_mirrors (bool) – is axis the one of a mirror ?

find_c_highest(probable_cn)

Among all rotation axis, find highest.

Parameters:: probable_cn (list(tuple)) – list of rotation axis (order, axis)

find_probable_parallel_mirrors(parallel_to)

Find possible mirrors parallel to a given axis (so that the normal is perpendicular).

For each combination of points, find if the corresponding normal is perpendicular to the axis (dot product is null), then remove duplicates.

Returns a list of normals.

Parameters:: parallel_to (numpy.ndarray|list) – axis to which the plane should be parallel to

find_probable_rotations(parallel_to=None)

Find rotations

Inspired by https://github.com/sunqm/pyscf/blob/master/pyscf/symm/geom.py !

The idea is to find, for each set of equidistant points, the ones that are at equal distance from the first point, which means that a rotation axis may pass by those points.

Parameters:: parallel_to (numpy.ndarray|list) – get only the rotation axis parallel to a given axis (should be normalized)

find_symmetry()

(Try to) find symmetry.

Algorithm is inspired by http://pymatgen.org/_modules/pymatgen/symmetry/analyzer.html (and actually any book on group theory):

Find inertia tensor, its eigenvalues and eigenvectors (axes)
Based on the eigenvalues:

Linear molecules have one zero egeinvalue, and belongs to \(C_{\infty v}\) or \(D_{\infty h}\) ;

Spherical top molecules have three equal (nonzero) eigenvalues, which corresponds to \(T\), \(I\) or \(O\) groups ;

Asymmetric top molecules have all different (non zero) eigenvalues, which corresponds to groups with at most rotation axis of order 2 (\(D_{2d}\) or subgroups), and here they are treated like symmetric top (since it is essentially the same thing) ;

Symmetric top molecules have two equals eigenvalues, thus corresponding to axial point group with (maybe) a unique principal rotation axis. To discriminate a little bit more,

Either there is perpendicular \(C_2\) axes, and the group is a dihedral one (depending on the presence or not of mirror planes, \(D_{nh}\), \(D_{nd}\) or \(D_n\)) ;

There is a main axis: depending on the presence of mirror plane, one can discriminate between the different cyclic groups (\(C_{nh}\), \(C_{nv}\), \(S_{2n}\) or \(C_n\)) ;

There is no main axis, and the group can be either \(C_s\), \(C_i\) or \(C_1\).

Find best orientation: for (a)symmetric, tetrahedral and octahedral molecules,

\(\vec{e}_z\) is the main axis for (a)symmetric top molecules, as well as for octahedral ones. For tetrahedral ones, the main axis is any \(C_2\).

The \(\vec{e}_x\) axis is either parallel to a reflexion plane, or a symmetry axis perpendicular to the main one if there is no reflexion plane in the case. If there is none of those, the inertia tensor is used. Note: for octahedral molecules, the \(C_4\) are the axes, against any logic.

\(\vec{e}_y = \vec{e}_z \times \vec{e}_x\).

In the other case, the orientation is taken from the inertia tensor.

Returns:: the point group, the center and the rotational matrix
Return type:: tuple(PointGroupDescription, numpy.ndarray, numpy.ndarray)

get_mirrors(parallel_to, other_axes)

Find the different mirror types (if they actually exists)

Parameters:

parallel_to (numpy.ndarray) – principal axis
other_axes (list) – the other rotation axes

Return type:

list

get_perpendicular_C2(perpendicular_to, probable_cn)

Get perpendicular \(C_2\) (if they actually exist)

Parameters:

probable_cn (list(tuple)) – list of rotation axis (order, axis)
perpendicular_to (numpy.ndarray) – main axis

Return type:

list

get_symmetry(force=False)

Find symmetry and store it, so that it is not computed next time (except if force)

Parameters:: force (bool) – force new search
Return type:: tuple(PointGroupDescription, numpy.ndarray, numpy.ndarray)

static group_points(points, decimals)

Group points per label and distances wrt center (within decimals)

Parameters:

points (numpy.ndarray) – the points, with the first axis being the label
decimals (int) – precision

Return type:

list

has_improper_rotation(axis, n)

has_inversion()

has_mirror(axis)

has_rotation(axis, n)

static inertia_tensor_moments(tensor)

Get inertia egeinvalue and egeinvectors, sorted (smallest first, so the main axis is the last eigenvector, if any).

Parameters:: tensor (numpy.ndarray) – carthesian tensor
Return type:: tuple

pass_trough(principal_axis, axis, is_mirror=False)

Determine by how much atoms the plane or the plane formed by the rotation axis pass trough

Parameters:

principal_axis (numpy.ndarray) – main axis
axis (numpy.ndarray) – axis to check
is_mirror (bool) – is axis the one of a mirror ?

Return type:

int

symmetric_for(op)

Check if symmetric for a given operation

Parameters:: op (Operation) – the operation

static vec_in_vecs(vec, vecs, tol)

Check if vec is in the set of vecs

Parameters:

vec – vector to find
vec – numpy.ndarray
vecs (numpy.ndarray|list) – set of vectors
tol (float) – tolerence threshold

Return type:

bool

exception qcip_tools.symmetry.SymmetryFinderError

qcip_tools.symmetry.simultaneous_diagonalization(matrices, tol=1e-14, in_vecs=None)

Simultaneous diagonalization of (communting Hermitian) matrices.

Seems to works by computing a basis for the first matrix, then by diagonalizing the different blocks (organized by the degeneracy of the eigenvalues).

Original source: http://qutip.org/docs/latest/modules/qutip/simdiag.html (I did a little bit of rewriting). Also check https://ocw.mit.edu/courses/physics/8-05-quantum-physics-ii-fall-2013/lecture-notes/MIT8_05F13_Chap_05.pdf (which is the closest explanation I found of this implementation, which is not common).

Parameters:

matrices (list) – list of matrices
tol (float) – tolerance on the eigenvalue (to check if they are degenerate or not)
in_vecs (numpy.ndarray) – previously computed eigenvectors

Symmetry (qcip_tools.symmetry)

Introduction

What’s a group anyway ?

Point group

Character table

Discover symmetry

API documentation

Symmetry (`qcip_tools.symmetry`)