Commit f020bca9 authored by Dániel Varjas's avatar Dániel Varjas
Browse files

start work on point group representations

parent ca4a5fc6
Pipeline #18150 passed with stage
in 8 minutes and 18 seconds
......@@ -91,6 +91,17 @@ def is_sympy_matrix(R):
return isinstance(R, (sympy.ImmutableMatrix, sympy.matrices.MatrixBase))
@ft.lru_cache(maxsize=10000)
def _U_phase_eq(U1, U2):
prop, coeff = prop_to_id(np.dot(la.inv(U1), U2))
return (prop and np.isclose(abs(coeff), 1))
@ft.lru_cache(maxsize=10000)
def _U_strict_eq(U1, U2):
return allclose(U1, U2)
class PointGroupElement:
"""An element of a point group.
......@@ -106,11 +117,10 @@ class PointGroupElement:
U : ndarray (optional)
The unitary action on the Hilbert space.
May be None, to be able to treat symmetry candidates
_strict_eq : boolean (default False)
Whether to test the equality of the unitary parts when comparing with
_U_eq : callable or None (default)
Function to test the equality of the unitary parts when comparing with
other PointGroupElements. By default the unitary parts are ignored.
If True, PointGroupElements are considered equal, if the unitary parts
are proportional, an overall phase difference is still allowed.
Should take two tinyarrays and return a boolean.
Notes
-----
......@@ -129,9 +139,9 @@ class PointGroupElement:
this works with floating point rotations.
"""
__slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq')
__slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_U_eq')
def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False):
def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _U_eq=None):
if isinstance(R, sympy.ImmutableMatrix):
# If it is integer, recast to integer tinyarray
R = _make_int(R)
......@@ -149,8 +159,8 @@ class PointGroupElement:
else:
raise ValueError('Real space rotation must be provided as a sympy matrix or an array.')
self.R, self.conjugate, self.antisymmetry, self.U = R, conjugate, antisymmetry, U
# Calculating sympy inverse is slow, remember it
self._strict_eq = _strict_eq
self._U_eq = _U_eq
def __repr__(self):
return ('\nPointGroupElement(\nR = {},\nconjugate = {},\nantisymmetry = {},\nU = {})'
......@@ -171,15 +181,14 @@ class PointGroupElement:
def __eq__(self, other):
R_eq = _eq(self.R, other.R)
basic_eq = R_eq and ((self.conjugate, self.antisymmetry) == (other.conjugate, other.antisymmetry))
# Equality is only checked for U if _strict_eq is True and basic_eq is True.
if basic_eq and (self._strict_eq is True or other._strict_eq is True):
# Equality is only checked for U if basic_eq is True.
if basic_eq and self._U_eq is not None:
if (self.U is None) and (other.U is None):
U_eq = True
elif (self.U is None) ^ (other.U is None):
U_eq = False
else:
prop, coeff = prop_to_id((self.inv() * other).U)
U_eq = (prop and np.isclose(abs(coeff), 1))
U_eq = self._U_eq(ta.array(self.U), ta.array(other.U))
else:
U_eq = True
return basic_eq and U_eq
......@@ -219,7 +228,7 @@ class PointGroupElement:
else:
U = U1.dot(U2)
R = _mul(R1, R2)
return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq))
return PointGroupElement(R, c1^c2, a1^a2, U, _U_eq=self._U_eq)
def __pow__(self, n):
result = self.identity()
......@@ -239,7 +248,7 @@ class PointGroupElement:
Uinv = la.inv(U)
# Check if inverse is stored, if not, calculate it
Rinv = _inv(R)
result = PointGroupElement(Rinv, c, a, Uinv, _strict_eq=self._strict_eq)
result = PointGroupElement(Rinv, c, a, Uinv, _U_eq=self._U_eq)
return result
def _strictereq(self, other):
......
......@@ -8,7 +8,7 @@ import tinyarray as ta
from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref, allclose
from .model import Model, BlochModel, _commutative_momenta, e, I
from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group
from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, _U_phase_eq
from . import kwant_continuum
......@@ -683,7 +683,7 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True,
if symmetrize:
# Make sure that group is generated while keeping track of unitary part.
for g in pg:
g._strict_eq = True
g._U_eq = _U_phase_eq
pg = generate_group(set(pg))
# Symmetrize every term and remove linearly dependent or zero ones
family2 = []
......
%% Cell type:code id: tags:
``` python
import numpy as np
import sympy
import qsymm
sympy.init_printing(print_builtin=True)
np.set_printoptions(precision=2, suppress=True, linewidth=150)
```
%% Cell type:markdown id: tags:
# Find symmetries
%% Cell type:markdown id: tags:
Define Rashba-Hamiltonian
%% Cell type:code id: tags:
``` python
ham1 = ("hbar^2 / (2 * m) * (k_x**2 + k_y**2 + k_z**2) * eye(2) +" +
"alpha * sigma_x * k_x + alpha * sigma_y * k_y + alpha * sigma_z * k_z")
# Convert to standard monomials form
H1 = qsymm.Model(ham1)
H1.tosympy()
```
%% Cell type:markdown id: tags:
Define `cubic_group` as the set of candidates for point group symmetries.
%% Cell type:code id: tags:
``` python
cubic_group = qsymm.groups.cubic()
```
%% Cell type:markdown id: tags:
Find point group symmetries, and discrete onsite and continuous symmetries.
%% Cell type:code id: tags:
``` python
sg, cg = qsymm.symmetries(H1, cubic_group)
print(len(sg))
print(cg)
```
%% Cell type:markdown id: tags:
Print the names of the group elements. It includes time reversal, but not inversion.
%% Cell type:code id: tags:
``` python
[display(g) for g in sg];
```
%% Cell type:markdown id: tags:
Detailed printout in LateX format.
%% Cell type:code id: tags:
``` python
from IPython.display import display, Math
[display(Math(qsymm.groups.pretty_print_pge(g, latex=True, full=True))) for g in sg];
```
%% Cell type:markdown id: tags:
Check that this is the same as the full cubic group without inversion plus TR.
%% Cell type:code id: tags:
``` python
C4 = qsymm.PointGroupElement(np.array([[1, 0, 0],
[0, 0, 1],
[0, -1, 0]]), False, False, None)
C3 = qsymm.PointGroupElement(np.array([[0, 0, 1],
[1, 0, 0],
[0, 1, 0]]), False, False, None)
TR = qsymm.PointGroupElement(np.eye(3), True, False, None)
set(sg) == qsymm.groups.generate_group({C4, C3, TR})
```
%% Cell type:markdown id: tags:
Add degeneracy, adds an `SU(2)` continuous symmetry group.
%% Cell type:code id: tags:
``` python
ham2 = "kron(eye(2), " + ham1 + ")"
print(ham2)
# Convert to standard monomials form
H2 = qsymm.Model(ham2)
display(H2.tosympy())
sg, cg = qsymm.symmetries(H2, cubic_group)
print(len(sg))
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
Add hole degrees of freedom instead, only adds `U(1)` charge conservation, as there is no pairing.
%% Cell type:code id: tags:
``` python
ham3 = "kron(sigma_z, " + ham1 + ")"
print(ham3)
# Convert to standard monomials form
H3 = qsymm.Model(ham3)
display(H3.tosympy())
sg, cg = qsymm.symmetries(H3, cubic_group)
print(len(sg))
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
Define Rashba-Hamiltonian with J = 3/2
%% Cell type:code id: tags:
``` python
J_x, J_y, J_z = qsymm.groups.spin_matrices(3/2)
ham32 = ("hbar^2 / (2 * m) * (k_x**2 + k_y**2 + k_z**2) * eye(4) +" +
"alpha * J_x * k_x + alpha * J_y * k_y + alpha * J_z * k_z")
# Convert to standard monomials form
H1 = qsymm.Model(ham32, locals=dict(J_x=J_x, J_y=J_y, J_z=J_z))
H1.tosympy(nsimplify=True)
```
%% Cell type:markdown id: tags:
Define `cubic_group` as the set of candidates for point group symmetries.
%% Cell type:code id: tags:
``` python
sg, cg = qsymm.symmetries(H1, cubic_group)
print(len(sg))
```
%% Cell type:markdown id: tags:
Check taht this is the same as the full cubic group without inversion plus TR.
%% Cell type:code id: tags:
``` python
set(sg) == qsymm.groups.generate_group({C4, C3, TR})
```
%% Cell type:markdown id: tags:
# Bloch Hamiltonian
%% Cell type:markdown id: tags:
Define hexagonal point group in 2D
%% Cell type:code id: tags:
``` python
hex_group_2D = qsymm.groups.hexagonal()
print(len(hex_group_2D))
```
%% Cell type:markdown id: tags:
Single band nearest neighbor hopping in triangular lattice
%% Cell type:code id: tags:
``` python
ham6 = 'm * (cos(k_x) + cos(1/2*k_x + sqrt(3)/2*k_y) + cos(-1/2*k_x + sqrt(3)/2*k_y))'
display(qsymm.sympify(ham6))
H6 = qsymm.Model(ham6, momenta=[0, 1])
sg, cg = qsymm.symmetries(H6, hex_group_2D)
print(len(sg))
print(set(sg) == qsymm.groups.hexagonal(ph=False))
print(cg)
```
%% Cell type:markdown id: tags:
Add spin and Rashba term
%% Cell type:code id: tags:
``` python
ham62 = 'eye(2) * (' + ham6 + ') +'
ham62 += 'alpha * (sin(k_x) * sigma_x + sin(1/2*k_x + sqrt(3)/2*k_y) * (1/2*sigma_x + sqrt(3)/2*sigma_y) +'
ham62 += 'sin(-1/2*k_x + sqrt(3)/2*k_y) * (-1/2*sigma_x + sqrt(3)/2*sigma_y))'
print(ham62)
display(qsymm.sympify(ham62))
H62 = qsymm.Model(ham62, momenta=[0, 1])
sg, cg = qsymm.symmetries(H62, hex_group_2D)
print(len(sg))
print(set(sg) == qsymm.groups.hexagonal(ph=False))
print(cg)
```
%% Cell type:markdown id: tags:
Add degeneracy
%% Cell type:code id: tags:
``` python
ham63 = 'kron(eye(2), ' + ham62 + ')'
display(qsymm.sympify(ham63))
H63 = qsymm.Model(ham63, momenta=[0, 1])
sg, cg = qsymm.symmetries(H63, hex_group_2D)
print(len(sg))
print(set(sg) == qsymm.groups.hexagonal(ph=False))
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
Add hole degrees of freedom
%% Cell type:code id: tags:
``` python
ham64 = 'kron(sigma_z, ' + ham62 + ')'
print(ham64)
display(qsymm.sympify(ham64))
H64 = qsymm.Model(ham64, momenta=[0, 1])
sg, cg = qsymm.symmetries(H64, hex_group_2D)
print(len(sg))
print(set(sg) == hex_group_2D)
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
# Continuous Rotation Symmetry
%% Cell type:markdown id: tags:
Rashba Hamiltonian actually has full rotation invariance.
%% Cell type:code id: tags:
``` python
display(qsymm.sympify(ham1))
pg, cg = qsymm.symmetries(H1, continuous_rotations=True, prettify=True)
```
%% Cell type:code id: tags:
``` python
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
# More complicated examples
## k.p Hamiltonian
%% Cell type:markdown id: tags:
Use an 8x8 k.p Hamiltonian of zinc blende semiconductors.
%% Cell type:code id: tags:
``` python
serialized_kp = 'Matrix([[hbar**2*k_x*gamma_0*k_x/(2*m_0)+hbar**2*k_y*gamma_0*k_y/(2*m_0)+hbar**2*k_z*gamma_0*k_z/(2*m_0)+E_0+E_v,0,-sqrt(2)*P*k_x/2-sqrt(2)*I*P*k_y/2,sqrt(6)*P*k_z/3,sqrt(6)*P*k_x/6-sqrt(6)*I*P*k_y/6,0,-sqrt(3)*P*k_z/3,-sqrt(3)*P*k_x/3+sqrt(3)*I*P*k_y/3],[0,hbar**2*k_x*gamma_0*k_x/(2*m_0)+hbar**2*k_y*gamma_0*k_y/(2*m_0)+hbar**2*k_z*gamma_0*k_z/(2*m_0)+E_0+E_v,0,-sqrt(6)*P*k_x/6-sqrt(6)*I*P*k_y/6,sqrt(6)*P*k_z/3,sqrt(2)*P*k_x/2-sqrt(2)*I*P*k_y/2,-sqrt(3)*P*k_x/3-sqrt(3)*I*P*k_y/3,sqrt(3)*P*k_z/3],[-sqrt(2)*k_x*P/2+sqrt(2)*I*k_y*P/2,0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)-hbar**2*k_x*gamma_2*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(2*m_0)-3*I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)-hbar**2*k_y*gamma_2*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(2*m_0)+3*I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)+hbar**2*k_z*gamma_2*k_z/m_0+E_v,sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)+sqrt(3)*hbar**2*k_x*k_z/(6*m_0)+sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)+sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)-sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),0,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(6)*hbar**2*k_x*k_z/(12*m_0)-sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(6)*hbar**2*k_z*k_x/(12*m_0)-sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)+sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0)],[sqrt(6)*k_z*P/3,-sqrt(6)*k_x*P/6+sqrt(6)*I*k_y*P/6,sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)-sqrt(3)*hbar**2*k_x*k_z/(6*m_0)-sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)+sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)+sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)+sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)+hbar**2*k_x*gamma_2*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(6*m_0)-I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)+hbar**2*k_y*gamma_2*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(6*m_0)+I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)-hbar**2*k_z*gamma_2*k_z/m_0+E_v,hbar**2*k_x*k_z/(3*m_0)+hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0-hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)-hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)-sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)+sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(2)*hbar**2*k_x*k_z/(12*m_0)-sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)-3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(2)*hbar**2*k_z*k_x/(12*m_0)-sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0)],[sqrt(6)*k_x*P/6+sqrt(6)*I*k_y*P/6,sqrt(6)*k_z*P/3,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-hbar**2*k_x*k_z/(3*m_0)-hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0+hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)+hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)+hbar**2*k_x*gamma_2*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(6*m_0)+I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)+hbar**2*k_y*gamma_2*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(6*m_0)-I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)-hbar**2*k_z*gamma_2*k_z/m_0+E_v,-sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)+sqrt(3)*hbar**2*k_x*k_z/(6*m_0)+sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)-sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)-sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(2)*hbar**2*k_x*k_z/(12*m_0)-sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)+3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(2)*hbar**2*k_z*k_x/(12*m_0)+sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)+sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)-sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0],[0,sqrt(2)*k_x*P/2+sqrt(2)*I*k_y*P/2,0,sqrt(3)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(3)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(3)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(3)*hbar**2*k_x*gamma_3*k_z/(2*m_0)-sqrt(3)*hbar**2*k_x*k_z/(6*m_0)-sqrt(3)*hbar**2*k_x*kappa*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*gamma_3*k_z/(2*m_0)-sqrt(3)*I*hbar**2*k_y*k_z/(6*m_0)-sqrt(3)*I*hbar**2*k_y*kappa*k_z/(2*m_0)-sqrt(3)*hbar**2*k_z*gamma_3*k_x/(2*m_0)-sqrt(3)*I*hbar**2*k_z*gamma_3*k_y/(2*m_0)+sqrt(3)*hbar**2*k_z*k_x/(6*m_0)+sqrt(3)*I*hbar**2*k_z*k_y/(6*m_0)+sqrt(3)*hbar**2*k_z*kappa*k_x/(2*m_0)+sqrt(3)*I*hbar**2*k_z*kappa*k_y/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)-hbar**2*k_x*gamma_2*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(2*m_0)+3*I*hbar**2*k_x*kappa*k_y/(2*m_0)-hbar**2*k_y*gamma_1*k_y/(2*m_0)-hbar**2*k_y*gamma_2*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(2*m_0)-3*I*hbar**2*k_y*kappa*k_x/(2*m_0)-hbar**2*k_z*gamma_1*k_z/(2*m_0)+hbar**2*k_z*gamma_2*k_z/m_0+E_v,sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)+sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)-sqrt(6)*hbar**2*k_x*k_z/(12*m_0)-sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)+sqrt(6)*hbar**2*k_z*k_x/(12*m_0)+sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)+sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0)],[-sqrt(3)*k_z*P/3,-sqrt(3)*k_x*P/3+sqrt(3)*I*k_y*P/3,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(6)*hbar**2*k_x*k_z/(12*m_0)+sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(6)*hbar**2*k_z*k_x/(12*m_0)-sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)-sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)+sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(2)*hbar**2*k_x*k_z/(12*m_0)+sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)-3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)-3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(2)*hbar**2*k_z*k_x/(12*m_0)+sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)-sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),-hbar**2*k_x*gamma_1*k_x/(2*m_0)-I*hbar**2*k_x*k_y/(3*m_0)-I*hbar**2*k_x*kappa*k_y/m_0-hbar**2*k_y*gamma_1*k_y/(2*m_0)+I*hbar**2*k_y*k_x/(3*m_0)+I*hbar**2*k_y*kappa*k_x/m_0-hbar**2*k_z*gamma_1*k_z/(2*m_0)-Delta+E_v,hbar**2*k_x*k_z/(3*m_0)+hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0-hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)-hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0],[-sqrt(3)*k_x*P/3-sqrt(3)*I*k_y*P/3,sqrt(3)*k_z*P/3,-sqrt(6)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(6)*I*hbar**2*k_x*gamma_3*k_y/(2*m_0)+sqrt(6)*hbar**2*k_y*gamma_2*k_y/(2*m_0)-sqrt(6)*I*hbar**2*k_y*gamma_3*k_x/(2*m_0),3*sqrt(2)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(2)*hbar**2*k_x*k_z/(12*m_0)+sqrt(2)*hbar**2*k_x*kappa*k_z/(4*m_0)+3*sqrt(2)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)+sqrt(2)*I*hbar**2*k_y*k_z/(12*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_z/(4*m_0)+3*sqrt(2)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+3*sqrt(2)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(2)*hbar**2*k_z*k_x/(12*m_0)-sqrt(2)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(2)*hbar**2*k_z*kappa*k_x/(4*m_0)-sqrt(2)*I*hbar**2*k_z*kappa*k_y/(4*m_0),sqrt(2)*hbar**2*k_x*gamma_2*k_x/(2*m_0)-sqrt(2)*I*hbar**2*k_x*k_y/(6*m_0)-sqrt(2)*I*hbar**2*k_x*kappa*k_y/(2*m_0)+sqrt(2)*hbar**2*k_y*gamma_2*k_y/(2*m_0)+sqrt(2)*I*hbar**2*k_y*k_x/(6*m_0)+sqrt(2)*I*hbar**2*k_y*kappa*k_x/(2*m_0)-sqrt(2)*hbar**2*k_z*gamma_2*k_z/m_0,-sqrt(6)*hbar**2*k_x*gamma_3*k_z/(4*m_0)+sqrt(6)*hbar**2*k_x*k_z/(12*m_0)+sqrt(6)*hbar**2*k_x*kappa*k_z/(4*m_0)+sqrt(6)*I*hbar**2*k_y*gamma_3*k_z/(4*m_0)-sqrt(6)*I*hbar**2*k_y*k_z/(12*m_0)-sqrt(6)*I*hbar**2*k_y*kappa*k_z/(4*m_0)-sqrt(6)*hbar**2*k_z*gamma_3*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*gamma_3*k_y/(4*m_0)-sqrt(6)*hbar**2*k_z*k_x/(12*m_0)+sqrt(6)*I*hbar**2*k_z*k_y/(12*m_0)-sqrt(6)*hbar**2*k_z*kappa*k_x/(4*m_0)+sqrt(6)*I*hbar**2*k_z*kappa*k_y/(4*m_0),-hbar**2*k_x*k_z/(3*m_0)-hbar**2*k_x*kappa*k_z/m_0-I*hbar**2*k_y*k_z/(3*m_0)-I*hbar**2*k_y*kappa*k_z/m_0+hbar**2*k_z*k_x/(3*m_0)+I*hbar**2*k_z*k_y/(3*m_0)+hbar**2*k_z*kappa*k_x/m_0+I*hbar**2*k_z*kappa*k_y/m_0,-hbar**2*k_x*gamma_1*k_x/(2*m_0)+I*hbar**2*k_x*k_y/(3*m_0)+I*hbar**2*k_x*kappa*k_y/m_0-hbar**2*k_y*gamma_1*k_y/(2*m_0)-I*hbar**2*k_y*k_x/(3*m_0)-I*hbar**2*k_y*kappa*k_x/m_0-hbar**2*k_z*gamma_1*k_z/(2*m_0)-Delta+E_v]])'
display(qsymm.sympify(serialized_kp))
Hkp = qsymm.Model(serialized_kp)
```
%% Cell type:markdown id: tags:
Symmetry is full cubic symmetry and time reversal.
%% Cell type:code id: tags:
``` python
sg, cg = qsymm.symmetries(Hkp, cubic_group)
print(len(sg))
[display(g) for g in cg];
```
%% Cell type:markdown id: tags:
Examples of how restricting parameters changes the symmetry group
%% Cell type:code id: tags:
``` python
kp_examples = [
{'hbar': 1},
{'hbar': 1, 'P': 0, 'Delta': 0},
{'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0},
{'hbar': 1, 'gamma_2': 'gamma_1', 'gamma_3': 'gamma_1', 'k_z': 0},
{'hbar': 1, 'P': 0, 'Delta': 0, 'gamma_3': 0},
{'hbar': 1, 'Delta': 0, 'k_z': 0},
]
for subs in kp_examples:
print(subs)
Hkp = qsymm.Model(serialized_kp, locals=subs)
%time sg, cg = qsymm.symmetries(Hkp, cubic_group)
print(len(sg), len(cg))
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
import itertools as it
import scipy.linalg as la
from qsymm.linalg import prop_to_id, allclose
from qsymm.groups import PointGroupElement, set_multiply, generate_group, _U_strict_eq
```
%% Cell type:code id: tags:
``` python
def identity_power(g, Ps=None, double_group=False):
# Append on-site U(1) symmetries to PG element to ensure that
# n'th power is exactly the identity.
# If double_group=True, if the n'th power is 2\pi rotation,
# -identity is chosen
def find_order(r):
dim = r.shape[0]
rpower = r
order = 1
while not np.allclose(rpower, np.eye(dim)):
rpower = r.dot(rpower)
order += 1
return order
R = np.array(g.R).astype(float)
n = find_order(R)
sign = 1
if double_group and n>1:
if np.isclose(la.det(R), 1):
# If not rotoinversion, its n'th power is 2pi rotation
sign = -1
elif find_order(-R) == n:
# here we assume R is 3d rotoinversion
# R^n is 2pi rotation
sign = -1
# if -R has a lower order (half of n) then
# R^n is a 4pi rotation or identity
Un = np.linalg.matrix_power(g.U, n)