diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2989d56723784d816e4ca2720dcba8a92a663a31..29bdf782bad45db5164c873ec610201eb3c71201 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -5,12 +5,27 @@ variables: import kwant.tests.test_qsymm; print(os.path.abspath(kwant.tests.test_qsymm.__file__)); +.only-default: &only-default + only: + - branches + - merge_requests + - tags + image: quantumtinkerer/research ## Documentation for the format of this file can be found here: ## https://docs.gitlab.com/ce/ci/yaml/README.html#configuration-of-your-builds-with-gitlab-ci-yml +test if changelog entry present: + script: + - git diff --name-only origin/$CI_MERGE_REQUEST_TARGET_BRANCH_NAME | grep "CHANGELOG" + stage: test + allow_failure: true + only: + - merge_requests + test minimal requirements: + <<: *only-default script: - conda env create -f environment-minimal.yml - source activate qsymm-minimal @@ -18,6 +33,7 @@ test minimal requirements: stage: test test latest requirements: + <<: *only-default script: - conda env create -f environment-latest.yml - source activate qsymm-latest @@ -25,6 +41,7 @@ test latest requirements: stage: test test kwant against latest qsymm: + <<: *only-default script: - conda env create -f environment-latest.yml - source activate qsymm-latest diff --git a/CHANGELOG.md b/CHANGELOG.md index c7b9a0b514c3155f4903f59b02dc002b6d70ebcb..38fd0da623f09c430ce08bccdaa0bd9d159eda9c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,10 +8,23 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -+ integration test: run the tests for Kwant's qsymm module against the current version - of qsymm. +- Check that a Changelog entry has been added for each merge request. +## [1.2.7] - 2019-12-04 + +### Added +- Added integration test: run the tests for Kwant's qsymm module against the + current version of qsymm. +- Added tests for fail cases in symmetry finder and Model with locals. + +### Fixed +- Bug when Model is initialized with numpy arrays in `locals`. +- Bug in symmetry finding when transformed Hamiltonian has very small coefficients, + this could lead to not finding symmetries in continuum models. +- Bug in Model.copy with sparse arrays did not return a copy of the stored arrays. +- Move docstrings to class level from __init__ functions and extend docstrings. + ## [1.2.6] - 2019-11-12 ### Fixed diff --git a/qsymm/groups.py b/qsymm/groups.py index eed8e6a9138517b139f147194d43e1e70fd37379..6c2a49ce3c0eb961845cd76d5c9786028a819e29 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -20,6 +20,13 @@ from .model import Model @lru_cache(maxsize=10000) def _mul(R1, R2): # Cached multiplication of spatial parts. + if R1 is None and R2 is None: + return None + elif R1 is None: + return R2 + elif R2 is None: + return R1 + if is_sympy_matrix(R1) and is_sympy_matrix(R2): # If spatial parts are sympy matrices, use cached multiplication. R = R1 * R2 @@ -41,7 +48,9 @@ def _mul(R1, R2): @lru_cache(maxsize=1000) def _inv(R): - if isinstance(R, ta.ndarray_int): + if R is None: + Rinv = None + elif isinstance(R, ta.ndarray_int): Rinv = ta.array(_inv_int(R), int) elif isinstance(R, ta.ndarray_float): Rinv = ta.array(la.inv(R)) @@ -54,6 +63,10 @@ def _inv(R): @lru_cache(maxsize=10000) def _eq(R1, R2): + if _is_identity(R1) and _is_identity(R2): + return True + elif (R1 is None) ^ (R2 is None): + return False if type(R1) != type(R2): R1 = ta.array(np.array(R1).astype(float), float) R2 = ta.array(np.array(R2).astype(float), float) @@ -78,6 +91,14 @@ def _make_int(R): return R +@lru_cache(maxsize=1000) +def _is_identity(R): + if R is None or R.shape[0] == 0: + return True + prop, coeff = prop_to_id(np.array(R).astype(float)) + return prop and np.isclose(coeff, 1) + + def _is_hermitian(a): return allclose(a, a.conjugate().transpose()) @@ -105,7 +126,12 @@ class PointGroupElement: Whether the operator flips the sign of the Hamiltonian (antisymmetry) U : ndarray (optional) The unitary action on the Hilbert space. - May be None, to be able to treat symmetry candidates + May be None, to be able to treat symmetry candidates, or as a shorthand + to act as the identity operator. + transpose: boolean (default False) + Whether the operator also includes transposition. Useful in non-Hermitian + systems to declare Hermiticity-like constraints. Similarly to conjugation + applied inside the unitary part: g(H) = U H(- R^−1 k)^T U^−1. _strict_eq : boolean (default False) Whether to test the equality of the unitary parts when comparing with other PointGroupElements. By default the unitary parts are ignored. @@ -129,10 +155,12 @@ class PointGroupElement: this works with floating point rotations. """ - __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq') + __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq', 'transpose') - def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False): - if isinstance(R, sympy.ImmutableMatrix): + def __init__(self, R=None, conjugate=False, antisymmetry=False, U=None, transpose=False, _strict_eq=False): + if R is None: + pass + elif isinstance(R, sympy.ImmutableMatrix): # If it is integer, recast to integer tinyarray R = _make_int(R) elif isinstance(R, ta.ndarray_int): @@ -149,15 +177,18 @@ 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 + self.transpose = transpose # Calculating sympy inverse is slow, remember it self._strict_eq = _strict_eq def __repr__(self): - return ('\nPointGroupElement(\nR = {},\nconjugate = {},\nantisymmetry = {},\nU = {})' + return ('\nPointGroupElement(\nR = {},\nconjugate = {},\nantisymmetry = {},\nU = {},\ntranspose = {})' .format(repr(self.R).replace('\n', '\n '), self.conjugate, self.antisymmetry, - repr(self.U).replace('\n', '\n ') if self.U is not None else 'None')) + repr(self.U).replace('\n', '\n ') if self.U is not None else 'None', + self.transpose) + ) def __str__(self): return pretty_print_pge(self, full=True) @@ -170,56 +201,78 @@ 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)) + basic_eq = R_eq and ((self.conjugate, self.antisymmetry, self.transpose) + == (other.conjugate, other.antisymmetry, other.transpose)) # 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): 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 + elif self.U is None: + U_eq = prop_to_id(other.U)[0] + elif other.U is None: + U_eq = prop_to_id(self.U)[0] else: - prop, coeff = prop_to_id((self.inv() * other).U) - U_eq = (prop and np.isclose(abs(coeff), 1)) + U_eq = prop_to_id((self.inv() * other).U)[0] else: U_eq = True return basic_eq and U_eq def __lt__(self, other): # Sort group elements: - # First by conjugate and a, then R = identity, then the rest + # First by conjugate and antisymmetry, then R = identity, then the rest # lexicographically + if not ((self.conjugate, self.antisymmetry, self.transpose) + == (other.conjugate, other.antisymmetry, other.transpose)): + return ((self.conjugate, self.antisymmetry, self.transpose) + < (other.conjugate, other.antisymmetry, other.transpose)) + elif (self.R is None) ^ (other.R is None): + return self.R is None + elif (self.R is None) and (other.R is None): + return False + Rs = ta.array(np.array(self.R).astype(float), float) Ro = ta.array(np.array(other.R).astype(float), float) identity = ta.array(np.eye(Rs.shape[0], dtype=int)) - - if not (self.conjugate, self.antisymmetry) == (other.conjugate, other.antisymmetry): - return (self.conjugate, self.antisymmetry) < (other.conjugate, other.antisymmetry) - elif (Rs == identity) ^ (Ro == identity): + if (Rs == identity) ^ (Ro == identity): return Rs == identity else: return Rs < Ro def __hash__(self): - # U is not hashed, if R is floating point it is also not hashed - R, c, a = self.R, self.conjugate, self.antisymmetry - if isinstance(R, ta.ndarray_float): - return hash((c, a)) + # U is not hashed, if R is close to identity, it is hashed as None, + # if R is floating point it is not hashed + R, c, a, t = self.R, self.conjugate, self.antisymmetry, self.transpose + if _is_identity(R): + return hash((None, c, a, t)) + + if is_sympy_matrix(R): + return hash((R, c, a, t)) + elif isinstance(R, ta.ndarray_float): + return hash((c, a, t)) else: - return hash((R, c, a)) + return hash((R, c, a, t)) def __mul__(self, g2): g1 = self - R1, c1, a1, U1 = g1.R, g1.conjugate, g1.antisymmetry, g1.U - R2, c2, a2, U2 = g2.R, g2.conjugate, g2.antisymmetry, g2.U + R1, c1, a1, t1, U1 = g1.R, g1.conjugate, g1.antisymmetry, g1.transpose, g1.U + R2, c2, a2, t2, U2 = g2.R, g2.conjugate, g2.antisymmetry, g2.transpose, g2.U - if (U1 is None) or (U2 is None): + # Take care of the transformation from conjugation on the second operator. + # Transpose behaves the same way as conjugation, hermitian adjoint behaves normally + if U2 is not None and c1^t1: + U2 = U2.conj() + + if (U1 is None) and (U2 is None): U = None - elif c1: - U = U1.dot(U2.conj()) + elif U1 is None: + U = U2 + elif U2 is None: + U = U1 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, transpose=t1^t2, _strict_eq=(self._strict_eq or g2._strict_eq)) def __pow__(self, n): result = self.identity() @@ -230,10 +283,11 @@ class PointGroupElement: def inv(self): """Invert PointGroupElement""" - R, c, a, U = self.R, self.conjugate, self.antisymmetry, self.U + R, c, a, U, t = self.R, self.conjugate, self.antisymmetry, self.U, self.transpose if U is None: Uinv = None - elif c: + elif c^t: + # transpose behaves the same way as conjugation, hermitian adjoint behaves normally Uinv = la.inv(U).conj() else: Uinv = la.inv(U) @@ -246,7 +300,8 @@ class PointGroupElement: # Stricter equality, testing the unitary parts to be approx. equal if (self.U is None) or (other.U is None): return False - return ((self.R, self.conjugate, self.antisymmetry) == (other.R, other.conjugate, other.antisymmetry) and + return ((self.R, self.conjugate, self.antisymmetry, self.transpose) + == (other.R, other.conjugate, other.antisymmetry, other.transpose) and allclose(self.U, other.U)) def apply(self, model): @@ -254,19 +309,22 @@ class PointGroupElement: if unitary: (+/-) U H(inv(R) k) U^dagger if antiunitary: (+/-) U H(- inv(R) k).conj() U^dagger + if transpose: (+/-) U H(- inv(R) k).T() U^dagger (+/-) stands for (symmetry / antisymmetry) - If self.U is None, U is taken as the identity. + If self.R or self.U is None, it is taken as the identity. """ R, antiunitary, antisymmetry, U = self.R, self.conjugate, self.antisymmetry, self.U + transpose = self.transpose R = _inv(R) - R = R * (-1 if antiunitary else 1) - result = model.rotate_momenta(R) + result = model.rotate(R, antiunitary, transpose) if antiunitary: result = result.conj() if antisymmetry: result = -result + if transpose: + result = result.T() if U is not None: result = U @ result @ U.T.conj() @@ -274,8 +332,10 @@ class PointGroupElement: def identity(self): """Return identity element with the same structure as self.""" - dim = self.R.shape[0] - R = ta.identity(dim, int) + if self.R is not None: + R = ta.identity(self.R.shape[0], int) + else: + R = None if self.U is not None: U = np.eye(self.U.shape[0]) else: @@ -284,7 +344,7 @@ class PointGroupElement: ## Factory functions for point group elements -def identity(dim, shape=None): +def identity(dim=None, shape=None): """Return identity operator with appropriate shape. Parameters @@ -299,7 +359,10 @@ def identity(dim, shape=None): ------- id : PointGroupElement """ - R = ta.identity(dim, int) + if dim is not None: + R = ta.identity(dim, int) + else: + R = None if shape is not None: U = np.eye(shape) else: @@ -307,7 +370,7 @@ def identity(dim, shape=None): return PointGroupElement(R, False, False, U) -def time_reversal(realspace_dim, U=None, spin=None): +def time_reversal(realspace_dim=None, U=None, spin=None): """Return a time-reversal symmetry operator parameters @@ -330,15 +393,18 @@ def time_reversal(realspace_dim, U=None, spin=None): ------- T : PointGroupElement """ + if realspace_dim is not None: + R = ta.identity(realspace_dim, int) + else: + R = None if U is not None and spin is not None: raise ValueError('Only one of `U` and `spin` may be provided.') if spin is not None: U = spin_rotation(np.pi * np.array([0, 1, 0]), spin) - R = ta.identity(realspace_dim, int) return PointGroupElement(R, conjugate=True, antisymmetry=False, U=U) -def particle_hole(realspace_dim, U=None): +def particle_hole(realspace_dim=None, U=None): """Return a particle-hole symmetry operator parameters @@ -353,11 +419,14 @@ def particle_hole(realspace_dim, U=None): ------- P : PointGroupElement """ - R = ta.identity(realspace_dim, int) + if realspace_dim is not None: + R = ta.identity(realspace_dim, int) + else: + R = None return PointGroupElement(R, conjugate=True, antisymmetry=True, U=U) -def chiral(realspace_dim, U=None): +def chiral(realspace_dim=None, U=None): """Return a chiral symmetry operator parameters @@ -372,7 +441,10 @@ def chiral(realspace_dim, U=None): ------- P : PointGroupElement """ - R = ta.identity(realspace_dim, int) + if realspace_dim is not None: + R = ta.identity(realspace_dim, int) + else: + R = None return PointGroupElement(R, conjugate=False, antisymmetry=True, U=U) @@ -492,6 +564,32 @@ def mirror(axis, U=None, spin=None): return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) + +def hermitian_adjoint(dim, shape=None): + """Return hermitian adjoint operator with appropriate shape. + Hermitian adjoint takes the adjoint of the matrix and reverses + all `hopping_vector` coordinates. + + Parameters + ---------- + dim : int + Dimension of real space. + shape : int (optional) + Size of the unitary part of the operator. + If not provided, U is set to None. + + Returns + ------- + herm : PointGroupElement + """ + R = ta.identity(dim, int) + if shape is not None: + U = np.eye(shape) + else: + U = None + return PointGroupElement(R, True, False, U, transpose=True) + + ## Continuous symmetry generators (conserved quantities) class ContinuousGroupGenerator: @@ -545,15 +643,12 @@ class ContinuousGroupGenerator: momenta = model.momenta R_nonzero = not (R is None or allclose(R, 0)) U_nonzero = not (U is None or allclose(U, 0)) + if R_nonzero: - dim = R.shape[0] - assert len(momenta) == dim - result = model.zeros_like() - if R_nonzero: - def trf(key): - return sum([sympy.diff(key, momenta[i]) * R[i, j] * momenta[j] - for i, j in product(range(dim), repeat=2)]) - result += 1j * model.transform_symbolic(trf) + result = model.rotate(R, infinitesimal=True) + else: + result = model.zeros_like() + if U_nonzero: result += model @ (1j*U) + (-1j*U) @ model return result @@ -900,13 +995,19 @@ def pretty_print_pge(g, full=False, latex=False): "" if den == 1 else ("/" + str(den))) return angle - R = np.array(g.R).astype(float) - if R.shape[0] == 1: + if g.R is None: + rot_name = '1' + dim = None + else: + R = np.array(g.R).astype(float) + dim = R.shape[0] + + if dim == 1: if R[0, 0] == 1: rot_name = '1' else: rot_name = 'I' - elif R.shape[0] == 2: + elif dim == 2: if np.isclose(la.det(R), 1): # pure rotation theta = np.arctan2(R[1, 0], R[0, 0]) @@ -926,7 +1027,7 @@ def pretty_print_pge(g, full=False, latex=False): rot_name = r'M\left({}\right)'.format(_round_axis(n)) else: rot_name = 'M({})'.format(_round_axis(n)) - elif R.shape[0] == 3: + elif dim == 3: if np.isclose(la.det(R), 1): # pure rotation n, theta = rotation_to_angle(R) @@ -961,12 +1062,21 @@ def pretty_print_pge(g, full=False, latex=False): .format(name_angle(theta, latex), _round_axis(n))) if full: + if g.conjugate: + if g.transpose: + conj_tr = "^{\dagger}" if latex else "^+" + else: + conj_tr = "^*" + elif g.transpose: + conj_tr = "^T" + else: + conj_tr = "" if latex: name = r'U H(\mathbf{{k}}){} U^{{-1}} = {}H({}R\mathbf{{k}}) \\' - name = name.format("^*" if g.conjugate else "", "-" if g.antisymmetry else "", + name = name.format(conj_tr, "-" if g.antisymmetry else "", "-" if g.conjugate else "") else: - name = '\nU⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", + name = '\nU⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format(conj_tr, "-" if g.antisymmetry else "", "-" if g.conjugate else "") name += 'R = {}'.format(rot_name) + (r'\\' if latex else '\n') @@ -1026,7 +1136,7 @@ def pretty_print_cgg(g, latex=False): if latex: rot_name = r'R_{{\phi}} = R\left(\phi, {}\right)\\'.format(_round_axis(n)) else: - rot_name = '\nR_ϕ = R(ϕ, {})'.format(_round_axis(n)) + rot_name = '\nR_ϕ = R(ϕ, {})\n'.format(_round_axis(n)) elif R is not None and R.shape[0] == 2: rot_name = r'R_{{\phi}} = R\left(\phi\right)\\' if latex else 'R_ϕ = R(ϕ)\n' else: @@ -1036,7 +1146,9 @@ def pretty_print_cgg(g, latex=False): if latex: L_name = r'L = {} \\'.format(_array_to_latex(np.round(L, 3))) else: - L_name = '\nL = {}'.format(str(np.round(L, 3)).replace('\n', '\n ')) + L_name = 'L = {}\n'.format(str(np.round(L, 3)).replace('\n', '\n ')) + else: + L_name = '' if latex: if R is not None: @@ -1047,11 +1159,11 @@ def pretty_print_cgg(g, latex=False): name = r'\left[ H(\mathbf{{k}}), L \right] = 0 \\' else: if R is not None: - name = '\n' + r'{}H(k){} = H(R_ϕ⋅k)' + name = r'{}H(k){} = H(R_ϕ⋅k)' + '\n' name = name.format(r'exp(-i ϕ L)⋅' if L is not None else '', r'⋅exp(i ϕ L)' if L is not None else '') else: - name = '\n[H(k), L] = 0' + name = '[H(k), L] = 0\n' name += rot_name + L_name return '$' + name + '$' if latex else name diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py index 031547698c7f353ede272f80e6db34ba3597e881..013df340b01e52b1b2133a173d9242f0056ad3b6 100644 --- a/qsymm/hamiltonian_generator.py +++ b/qsymm/hamiltonian_generator.py @@ -7,35 +7,46 @@ from copy import deepcopy import tinyarray as ta from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref, allclose -from .model import Model, BlochModel, BlochCoeff, _commutative_momenta, e, I +from .model import Model, BlochModel, BlochCoeff, _commutative_momenta, e, I, _symbol_normalizer from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group from . import kwant_continuum -def continuum_hamiltonian(symmetries, dim, total_power, all_powers=None, - momenta=_commutative_momenta, sparse_linalg=False, - prettify=False, num_digits=10): - """Generate a family of continuum Hamiltonians that satisfy symmetry constraints. +def continuum_hamiltonian(symmetries, dim=None, total_power=2, all_powers=None, + momenta=_commutative_momenta, + coordinates=None, + hopping_vector=None, + sparse_linalg=False, + prettify=False, num_digits=10, hermitian=True): + """Generate a family of continuum Hamiltonians `H(k, x, d)` as matrix-valued + polynomials of vector quantities `momenta` (`k`), `coordinates` (`x`) and + `hopping_vector` (`d`) satisfying specified symmetries. Parameters ---------- - symmetries: iterable of PointGroupElement objects. - An iterable of PointGroupElement objects, each describing a symmetry - that the family should possess. - dim: integer - The number of spatial dimensions along which the Hamiltonian family is - translationally invariant. Only the first `dim` entries in `all_powers` and - `momenta` are used. - total_power: integer or list of integers - Allowed total powers of the momentum variables in the continuum Hamiltonian. + symmetries: iterable of PointGroupElement and ContinuousGroupGenerator objects. + List of all symmetries that the family should possess. Sufficient to provide + a generator set, the result possesses all symmetries generated as combinations + of the symmetries. The size of the matrix is deduced from the size of the + unitary transformations in `symmetries`. + dim: integer (optional) + The number of spatial dimensions. If symmetries are provided, the + dimensionality is deduced from the size of rotation matrices. Only the first + `dim` entries in `all_powers`, `momenta`, `coordinates` and `hopping_vector` + are used. + total_power: integer or list of integers, default 2 + Allowed total powers of the vector variables in the continuum Hamiltonian. If an integer is given, all powers below it are included as well. - all_powers: list of integer or list of list of integers + all_powers: list of vector or list of list of integers Allowed powers of the momentum variables in the continuum Hamiltonian listed for each spatial direction. If an integer is given, all powers below it are included as well. If a list of integers is given, only these powers are used. - momenta : iterable of strings or Sympy symbols - Names of momentum variables, default ``['k_x', 'k_y', 'k_z']`` or - corresponding sympy symbols. + momenta : iterable of strings or Sympy symbols or None + Names of momentum variables, default ``('k_x', 'k_y', 'k_z')``. + coordinates : iterable of strings or Sympy symbols or None (default) + Names of coordinate variables. + hopping_vector : iterable of strings or Sympy symbols or None (default) + Names of hopping vector variables. sparse_linalg : bool Whether to use sparse linear algebra. Using sparse solver can result in performance increase for large, highly constrained families, @@ -46,6 +57,9 @@ def continuum_hamiltonian(symmetries, dim, total_power, all_powers=None, num_digits: integer, default 10 Number of significant digits to which the basis is rounded using prettify. This argument is ignored if prettify = False. + hermitian: bool, default True + Whether to generate only Hermitian models. If `hopping_vector` is provided, + the appropriate hermiticity constraint is used: `H(k, x, d)^† = H(k, x, -d)`. Returns ------- @@ -59,18 +73,52 @@ def continuum_hamiltonian(symmetries, dim, total_power, all_powers=None, max_power = total_power total_power = range(max_power + 1) + for symmetry in symmetries: + if symmetry.R is not None and dim is None: + dim = symmetry.R.shape[0] + elif symmetry.R is not None and symmetry.R.shape[0] != dim: + raise ValueError("Symmetry operator {} has incompatible shaped rotation " + "matrix R = {} with dim = {}.".format(symmetry, symmetry.R, dim)) + if dim is None: + raise ValueError("Could not deduce real space dimensionality, please provide `dim`.") + + var = {k: (None if v is None else v[:dim]) + for k, v in zip(["momenta", "coordinates", "hopping_vector"], + [momenta, coordinates, hopping_vector]) + } + + if len([v for v in var.values() if v is not None]) > 1: + raise NotImplementedError("Only one of momenta, coordinates, hopping_vector " + "may be provided.") + + # Dimension of Hamiltonian matrix + N = list(symmetries)[0].U.shape[0] + + # Hermiticity needs special treatment with hopping vector + if hopping_vector is not None: + hermitian_bas = False + if hermitian: + herm = hermitian_adjoint(dim, shape=N) + symmetries = list(symmetries) + [herm] + else: + hermitian_bas = hermitian + # Generate a Hamiltonian family - N = list(symmetries)[0].U.shape[0] # Dimension of Hamiltonian matrix - momenta = momenta[:dim] + keys = [v for v in var.values() if v is not None][0] # Symmetries do not mix the total degree of momenta. We can thus work separately at each # fixed degree. family = [] for degree in total_power: # Make all momentum variables given the constraints on dimensions and degrees in the family - momentum_variables = continuum_variables(dim, degree, all_powers=all_powers, momenta=momenta) - degree_family = [Model({momentum_variable: matrix}, momenta=momenta) + momentum_variables = continuum_variables(dim, degree, all_powers=all_powers, momenta=keys) + # Make matrix basis depending on whether hermitian or not + if hermitian_bas: + mat_bas = matrix_basis(N) + else: + mat_bas = it.chain(matrix_basis(N), matrix_basis(N, antihermitian=True)) + degree_family = [Model({momentum_variable: matrix}, **var) for momentum_variable, matrix - in it.product(momentum_variables, matrix_basis(N))] + in it.product(momentum_variables, mat_bas)] degree_family = constrain_family(symmetries, degree_family, sparse_linalg=sparse_linalg) if prettify: family_size = len(degree_family) @@ -287,7 +335,9 @@ def hamiltonian_from_family(family, coeffs=None, nsimplify=True, tosympy=True): if coeffs is None: coeffs = list(sympy.symbols('c0:%d'%len(family))) else: - assert len(coeffs) == len(family), 'Length of family and coeffs do not match.' + # sympify strings + coeffs = [(_symbol_normalizer(s) if isinstance(s, str) else s) for s in coeffs] + assert len(coeffs) == len(family), 'Length of family and coeffs do not match.' # The order of multiplication is important here, so that __mul__ of 'term' # gets used. 'c' is a sympy symbol, which multiplies 'term' incorrectly. ham = sum(term * c for c, term in zip(coeffs, family)) @@ -392,16 +442,6 @@ def constrain_family(symmetries, family, sparse_linalg=False): # Fix ordering family = list(family) symmetries = list(symmetries) - # Check compatibility of family members and symmetries - shape = family[0].shape - momenta = family[0].momenta - for term in family: - assert term.shape == shape - assert term.momenta == momenta - for symmetry in symmetries: - assert symmetry.U.shape == shape - if symmetry.R is not None: - assert symmetry.R.shape[0] == len(momenta) # Need all the linearly independent variables before and after # rotation to make the matrix of linear constraints. @@ -573,7 +613,7 @@ def symmetrize_monomial(monomial, symmetries): def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, momenta=_commutative_momenta, symmetrize=True, prettify=True, num_digits=10, - bloch_model=False): + bloch_model=False, hermitian=True): """Generate a family of symmetric Bloch-Hamiltonians. Parameters @@ -610,6 +650,8 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, a list of Model objects. If True, returns a list of BlochModel objects. BlochModel objects are more suitable than Model objects if the hopping vectors include floating point numbers. + hermitian: bool, default True + Whether to generate only Hermitian models. Returns ------- @@ -690,7 +732,8 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, term = BlochModel({bloch_coeff: matrix}, momenta=momenta[:dim]) else: term = Model({factor: matrix}, momenta=momenta[:dim]) - term = term + term.T().conj() + if hermitian: + term = term + term.T().conj() hopfamily.append(term) # If there are conserved quantities, constrain the hopping, it is assumed that # conserved quantities do not mix different sites diff --git a/qsymm/kwant_continuum.py b/qsymm/kwant_continuum.py index 9a720522ad97cf806c9f30e2c0cfa45efe645706..06b8336ccf39f6789466f57345eff52067c0106b 100644 --- a/qsymm/kwant_continuum.py +++ b/qsymm/kwant_continuum.py @@ -206,7 +206,10 @@ def sympify(expr, locals=None): "identifiers and may not be keywords".format(repr(k))) # sympify values of locals before updating it with extra_ns - locals = {k: sympify(v) for k, v in locals.items()} + # Cast numpy array values in locals to sympy matrices to make sure they have + # correct format + locals = {k: (sympy.Matrix(v) if isinstance(v, np.ndarray) else sympify(v)) + for k, v in locals.items()} for k, v in extra_ns.items(): locals.setdefault(k, v) try: diff --git a/qsymm/linalg.py b/qsymm/linalg.py index 7f925b68fe398ea094cd4c58433fcd187a59b678..16e43f169b80ad90b55883125b11d9777c8d9ad2 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -4,7 +4,7 @@ import scipy.linalg as la import scipy.sparse.linalg as sla import scipy from numbers import Number -from scipy.optimize import minimize +from scipy.optimize import minimize, linear_sum_assignment from scipy.spatial.distance import cdist from scipy.sparse.csgraph import connected_components import itertools as it @@ -61,6 +61,29 @@ def allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): return np.allclose(diff, 0, rtol=0, atol=atol, equal_nan=equal_nan) +def spectra_allclose(a, b, rtol=1e-05, atol=1e-08): + """Decide whether the spectra a and b (iterables of complex numbers + with the same length) are identical up to permutations and numerical + precision. Useful for spectra of non-Hermitian matrices, as these are + typically returned in arbitrary order. Relative tolerance is relative + to the maximum absolute value entry. Returns True if the spectra are + close.""" + tol = atol + rtol * max(np.max(np.abs(a)), np.max(np.abs(b))) + # Pre-test the sums + if not np.isclose(np.sum(a), np.sum(b), atol=tol * len(a)): + return False + # Make the distance matrix + dist = np.abs(a[:, np.newaxis] - b[np.newaxis, :]) + # set smaller than tolerance distances to 0 + dist[dist < tol] = 0 + row_ind, col_ind = linear_sum_assignment(dist) + # check if the minimal assignment is 0 + if np.all(dist[row_ind, col_ind] == 0): + return True + else: + return False + + def mtm(a, B, c): # matrix-tensor-matrix multiplication for dimensions 2, 3, 2. # Equivalent to 'np.einsum('ij,ajk,kl->ail', a, B, c)' diff --git a/qsymm/model.py b/qsymm/model.py index 857ee54cbcf29530da80d6a5db940ed3dba287d8..36e6e14d4d73ced75da27d12cf19fb40216194cc 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -3,7 +3,7 @@ import scipy import tinyarray as ta import scipy.linalg as la from itertools import product -from copy import copy +import copy as copy_module from numbers import Number from warnings import warn from functools import lru_cache @@ -22,6 +22,17 @@ e = kwant_continuum.sympify('e') I = kwant_continuum.sympify('I') +# Scipy sparse matrices defined a 'copy' method (which does +# a deep-copy of their data) but no '__copy__' method, to work +# correctly with the 'copy' module. We therefore make our own +# wrapper here that does the right thing +def copy(a): + if callable(getattr(a, 'copy', None)): + return a.copy() + else: + return copy_module.copy(a) + + def substitute_exponents(expr): """Substitute trignometric functions with exp. @@ -52,12 +63,12 @@ def substitute_exponents(expr): class BlochCoeff(tuple): + """ + Container for Bloch coefficient in ``BlochModel``, in the form of + ``(hop, coeff)``, equivalent to ``coeff * exp(I * hop.dot(k))``. + """ def __new__(cls, hop, coeff): - """ - Container for Bloch coefficient in ``BlochModel``, in the form of - ``(hop, coeff)``, equivalent to ``coeff * exp(I * hop.dot(k))``. - """ if not (isinstance(hop, np.ndarray) and isinstance(coeff, sympy.Expr)): raise ValueError('`hop` must be a 1D numpy array and `coeff` a sympy expression.') if isinstance(coeff, sympy.add.Add): @@ -112,6 +123,89 @@ class BlochCoeff(tuple): class Model(UserDict): + """ + Symbolic matrix-valued function that depends on momenta and other parameters. + + Implements the algebra of matrix valued functions. + Implements many sympy and numpy methods and overrides arithmetic operators. + Internally it represents ``sum(symbol * value)``, where ``symbol`` is a symbolic + expression, and ``value`` can be scalar, array (both dense and sparse) + or LinearOperator. This is accessible as a dict ``{symbol: value}``. + + Parameters + ---------- + hamiltonian : str, SymPy expression, dict or None (default) + Symbolic representation of a Hamiltonian. If a string, it is + first converted to a SymPy expression using `kwant_continuum.sympify`. + If a dict is provided, it should have the form + ``{symbol: array}`` with all arrays the same size (dense or sparse). + ``symbol`` by default is passed through sympy.sympify, and should + consist purely of a product of symbolic coefficients, no constant + factors other than 1, except if ``normalize=True``. ``None`` initializes + a zero ``Model``. + locals : dict or ``None`` (default) + Additional namespace entries for `~kwant_continuum.sympify`. May be + used to simplify input of matrices or modify input before proceeding + further. For example: + ``locals={'k': 'k_x + I * k_y'}`` or + ``locals={'sigma_plus': [[0, 2], [0, 0]]}``. + keep : iterable of expressions (optional) + Set of symbolic coefficients that are kept, anything that does not + appear here is discarded. Useful for perturbative calculations where + only terms to a given order are needed. By default all keys are kept. + momenta : iterable of strings or Sympy symbols or None + Names of momentum variables, as strings or sympy symbols, default + ``('k_x', 'k_y', 'k_z')``. Momenta are treated the same as other keys + for the purpose of `keep`. Momenta are transformed under spatial rotations + and reverse sing under complex conjugation and transposition in + PointGroupElement. + coordinates : iterable of strings or Sympy symbols or None (default) + Names of coordinate variables, as strings or sympy symbols. If both provided, + must have the same length as ``momenta``. Coordinates are treated the + same as other keys for the purpose of `keep`. Coordinates are + transformed under spatial rotations and do not reverse sign under complex + conjugation and transposition in PointGroupElement. + hopping_vector : iterable of strings or Sympy symbols or None (default) + Names of coordinate variables describing the hopping vector as strings + or sympy symbols, if the model represents a single hopping matrix. + If both provided, must have the same length as ``momenta``. Treated the + same as other keys for the purpose of `keep`. hopping_vector is + transformed under spatial rotations and reverse sign under transposition + in PointGroupElement. + symbol_normalizer : callable (optional) + Function applied to symbols when initializing the internal dict. By default the + keys are passed through ``sympy.sympify`` and ``sympy.expand_power_exp``. + Keys when accessing a term and keys in ``keep`` are also passed through + ``symbol_normalizer``. + normalize : bool, default False + Whether to clean input dict by splitting summands in symbols, + moving numerical factors in the symbols to values, removing entries + with values allclose to zero. Ignored if hamiltonian is not a dict. + shape : tuple or None (default) + Shape of the Model, must match the shape of all the values. If not + provided, it is automatically found based on the shape of the input. + Must be provided if ``hamiltonian`` is ``None`` or ``{}``. Empty tuple + corresponds to scalar values. + format : class or None (default) + Type of the values in the model. Supported types are + ``np.complex128``, ``scipy.sparse.linalg.LinearOperator``, ``np.ndarray``, + and subclasses of ``scipy.sparse.spmatrix`` . If ``hamiltonian`` is + provided as a dict, all values must be of this type, except for + scalar values, which are recast to ``np.complex128``. If ``format`` is + not provided, it is inferred from the type of the values. Must be + provided if ``hamiltonian`` is `None` or ``{}``. If ``hamiltonian`` is + not a dictionary, ``format`` is ignored an set to ``np.ndarray``. + + Notes + ----- + Sympy symbols are immutable and references to the same symbols is + stored in different Models. Be warned that setting any assumptions + for symbols (such as ``real``) will result in an identically named, + but different symbol, and these are not handled properly. Model assumes + that all sympy symbols are real and commutative without any assumptions + explicitly set. Non-commutative symbols are not implemented, momenta and + coordinates also commute. + """ # Make it work with numpy arrays __array_ufunc__ = None @@ -121,78 +215,37 @@ class Model(UserDict): hamiltonian=None, locals=None, momenta=('k_x', 'k_y', 'k_z'), + coordinates=None, + hopping_vector=None, keep=None, symbol_normalizer=None, normalize=False, shape=None, format=None ): - """ - Symbolic matrix-valued function that depends on momenta and other parameters. - - Implements the algebra of matrix valued functions. - Implements many sympy and numpy methods and overrides arithmetic operators. - Internally it represents ``sum(symbol * value)``, where ``symbol`` is a symbolic - expression, and ``value`` can be scalar, array (both dense and sparse) - or LinearOperator. This is accessible as a dict ``{symbol: value}``. - - Parameters - ---------- - hamiltonian : str, SymPy expression, dict or None (default) - Symbolic representation of a Hamiltonian. If a string, it is - first converted to a SymPy expression using `kwant_continuum.sympify`. - If a dict is provided, it should have the form - ``{symbol: array}`` with all arrays the same size (dense or sparse). - ``symbol`` by default is passed through sympy.sympify, and should - consist purely of a product of symbolic coefficients, no constant - factors other than 1, except if ``normalize=True``. ``None`` initializes - a zero ``Model``. - locals : dict or ``None`` (default) - Additional namespace entries for `~kwant_continuum.sympify`. May be - used to simplify input of matrices or modify input before proceeding - further. For example: - ``locals={'k': 'k_x + I * k_y'}`` or - ``locals={'sigma_plus': [[0, 2], [0, 0]]}``. - keep : iterable of expressions (optional) - Set of symbolic coefficients that are kept, anything that does not - appear here is discarded. Useful for perturbative calculations where - only terms to a given order are needed. By default all keys are kept. - momenta : iterable of strings or Sympy symbols - Names of momentum variables, default ``['k_x', 'k_y', 'k_z']`` or - corresponding sympy symbols. Momenta are treated the same as other - keys for the purpose of `keep`. - symbol_normalizer : callable (optional) - Function applied to symbols when initializing the internal dict. By default the - keys are passed through ``sympy.sympify`` and ``sympy.expand_power_exp``. - Keys when accessing a term and keys in ``keep`` are also passed through - ``symbol_normalizer``. - normalize : bool, default False - Whether to clean input dict by splitting summands in symbols, - moving numerical factors in the symbols to values, removing entries - with values allclose to zero. Ignored if hamiltonian is not a dict. - shape : tuple or None (default) - Shape of the Model, must match the shape of all the values. If not - provided, it is automatically found based on the shape of the input. - Must be provided if ``hamiltonian`` is ``None`` or ``{}``. Empty tuple - corresponds to scalar values. - format : class or None (default) - Type of the values in the model. Supported types are - ``np.complex128``, ``scipy.sparse.linalg.LinearOperator``, ``np.ndarray``, - and subclasses of ``scipy.sparse.spmatrix`` . If ``hamiltonian`` is - provided as a dict, all values must be of this type, except for - scalar values, which are recast to ``np.complex128``. If ``format`` is - not provided, it is inferred from the type of the values. Must be - provided if ``hamiltonian`` is `None` or ``{}``. If ``hamiltonian`` is - not a dictionary, ``format`` is ignored an set to ``np.ndarray``. - """ if hamiltonian is None: hamiltonian = {} + if symbol_normalizer is None: symbol_normalizer = _symbol_normalizer - self.momenta = _find_momenta(tuple(momenta)) + + if momenta is not None: + self.momenta = _find_momenta(tuple(momenta)) + else: + self.momenta = None if keep is not None: self.keep = {symbol_normalizer(k) for k in keep} else: self.keep = set() + if coordinates is not None: + self.coordinates = tuple(symbol_normalizer(x) for x in coordinates) + else: + self.coordinates = None + + if hopping_vector is not None: + self.hopping_vector = tuple(symbol_normalizer(x) for x in hopping_vector) + else: + self.hopping_vector = None + if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): # Initialize as dict sympifying the keys self.data = {symbol_normalizer(k): v for k, v in hamiltonian.items() @@ -262,6 +315,10 @@ class Model(UserDict): # remove zero entries, apply symbol_normalizer self.data = {symbol_normalizer(k): v for k, v in self.items() if not allclose(v, 0)} + @property + def _special_symbols(self): + return {k: getattr(self, k) for k in ('momenta', 'coordinates', 'hopping_vector') if hasattr(self, k)} + # Make sure values have the correct format def __setitem__(self, key, item): if (isinstance(item, self.format) and self.shape == item.shape): @@ -316,8 +373,9 @@ class Model(UserDict): if not (self.format is other.format and self.shape == other.shape): raise ValueError('Addition is only possible for Models with the same shape and data type.') # other is not empty, so the result is not empty - if self.momenta != other.momenta: - raise ValueError("Can only add Models with the same momenta") + if self._special_symbols != other._special_symbols: + raise ValueError("Can only add Models with the same momenta, " + "coordinates and hopping_vector") result = self.zeros_like() for key in self.keys() & other.keys(): result[key] = self[key] + other[key] @@ -366,7 +424,8 @@ class Model(UserDict): keep = self.keep result = sum((type(self)({key * other: copy(val)}, keep=keep, - momenta=self.momenta) + **self._special_symbols, + ) for key, val in self.items() if (key * other in keep or not keep)), self.zeros_like()) @@ -376,12 +435,14 @@ class Model(UserDict): raise ValueError('Elementwise multiplication only allowed for scalar ' 'and ndarra data types. With sparse matrices use `@` ' 'for matrix multiplication.') - if self.momenta != other.momenta: - raise ValueError("Can only multiply Models with the same momenta") + if self._special_symbols != other._special_symbols: + raise ValueError("Can only multiply Models with the same momenta, " + "coordinates and hopping_vector") keep = self.keep | other.keep result = sum(type(self)({k1 * k2: v1 * v2}, keep=keep, - momenta=self.momenta) + **self._special_symbols, + ) for (k1, v1), (k2, v2) in product(self.items(), other.items()) if (k1 * k2 in keep or not keep)) # Find out the shape of the result even if it is empty @@ -406,7 +467,8 @@ class Model(UserDict): # is correct as long as the symbols in 'key' and 'other' commute. result = sum((type(self)({key * other: copy(val)}, keep=keep, - momenta=self.momenta) + **self._special_symbols, + ) for key, val in self.items() if (key * other in keep or not keep)), self.zeros_like()) @@ -420,12 +482,14 @@ class Model(UserDict): def __matmul__(self, other): # Multiplication by arrays and Model if isinstance(other, Model): - if self.momenta != other.momenta: - raise ValueError("Can only multiply Models with the same momenta") + if self._special_symbols != other._special_symbols: + raise ValueError("Can only multiply Models with the same momenta, " + "coordinates and hopping_vector") keep = self.keep | other.keep result = sum(type(self)({k1 * k2: v1 @ v2}, keep=keep, - momenta = self.momenta) + **self._special_symbols, + ) for (k1, v1), (k2, v2) in product(self.items(), other.items()) if (k1 * k2 in keep or not keep)) # Find out the shape of the result even if it is empty @@ -468,9 +532,9 @@ class Model(UserDict): def zeros_like(self): """Return an empty model object that inherits the other properties""" result = type(self)(shape=self.shape, - format=self.format) - result.keep=self.keep.copy() - result.momenta=self.momenta + format=self.format, + **self._special_symbols) + result.keep = self.keep.copy() return result def transform_symbolic(self, func): @@ -479,23 +543,70 @@ class Model(UserDict): # Add possible duplicate keys that only differ in constant factors result = sum((type(self)({func(key): copy(val)}, normalize=True, - momenta=self.momenta) - for key, val in self.items()), + keep=self.keep, + **self._special_symbols, + ) + for key, val in self.items() + if func(key) in self.keep or not self.keep), self.zeros_like()) return result - def rotate_momenta(self, R): - """Rotate momenta with rotation matrix R.""" - momenta = self.momenta - assert len(momenta) == R.shape[0], (momenta, R) + def rotate(self, R, conjugate=False, transpose=False, infinitesimal=False): + r""" + Rotate momenta and real space coordinates with rotation matrix ``R``: + + .. math:: H(k, x, d) → H(±R k, R x, ±R d) + + with the transformation preformed for all vector-like quantities the + model depends on. - k_prime = R @ sympy.Matrix(momenta) - rotated_subs = {k: k_prime for k, k_prime in zip(momenta, k_prime)} + Momenta ``k`` are inverted by conjugate xor transpose, + hopping_vector ``d`` is inverted by transpose. + The matrix part of the model is not conjugated or transposed! - def trf(key): - return key.subs(rotated_subs, simultaneous=True) + If ``infinitesimal=True``, perform infinitesimal spatial rotation + generated by ``exp(i λ R)``: - return self.transform_symbolic(trf) + .. math:: H(k) → sum_{ij} 1j * dH(k)/dk_i R_{ij} k_j + + plus the same transformation for ``x`` and ``d``. + """ + if infinitesimal and (conjugate or transpose): + raise ValueError("If `infinitesimal=True` `conjugate` and `transpose` " + "must be False.") + for key, vec in self._special_symbols.items(): + if not (vec is None or R is None or len(vec) == R.shape[0]): + raise ValueError("Shape of rotation matrix {} incompatible with " + "number of {} variables {}.".format(R, key, vec)) + + if infinitesimal: + def trf(key): + return sum([sympy.diff(key, var[i]) * R[i, j] * var[j] + for i, j in product(range(R.shape[0]), repeat=2) + for var in self._special_symbols.values() if var is not None + ]) + return 1j * self.transform_symbolic(trf) + + else: + signs = {'momenta' : (-1 if (conjugate^transpose) else 1), + 'coordinates' : 1, + 'hopping_vector' : (-1 if transpose else 1)} + + subs = dict() + for key, vec in self._special_symbols.items(): + if vec is None: + continue + if R is None: + vec_prime = sympy.Matrix(vec) + else: + vec_prime = R @ sympy.Matrix(vec) + # change sign as appropriate + vec_prime = signs[key] * vec_prime + subs.update(dict(zip(vec, vec_prime))) + + def trf(key): + return key.subs(subs, simultaneous=True) + return self.transform_symbolic(trf) def subs(self, *args, **kwargs): """Substitute symbolic expressions. See documentation of @@ -519,6 +630,7 @@ class Model(UserDict): args = ([(key, value) for key, value in args[0].items()], ) momenta = self.momenta + #### TODO implement subs for coordinates and hoppin_coordinates for (old, new) in args[0]: # Substitution of a momentum variable with a symbol # is a renaming of the momentum. @@ -704,62 +816,83 @@ class Model(UserDict): return all(allclose(self[key], other[key], rtol, atol, equal_nan) for key in self.keys() | other.keys()) + def eliminate_zeros(self, rtol=1e-05, atol=1e-08): + """Return a model with small terms removed. Tolerances are + in Frobenius matrix norm, relative tolerance compares to the + value with largest norm.""" + if not issubclass(self.format, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Operation only supported for Models with ' + '`np.ndarray` or `scipy.sparse.spmatrix` data type.') + result = self.zeros_like() + # For empty model do nothing + if self.data == {}: + return result + # Write it explicitely so it works with sparse arrays + norm = lambda mat: np.sqrt(np.sum(np.abs(mat)**2)) + max_norm = np.max([norm(val) for val in self.values()]) + tol = max(atol, max_norm * rtol) + result.data = {key: copy(val) for key, val in self.items() if not norm(val) < tol} + return result + class BlochModel(Model): + """ + A ``Model`` where coefficients are periodic functions of momenta. + + Internally it is a ``sum(BlochCoeff * value)``, where ``BlochCoeff`` is + a symbolic representation of coefficients and a periodic function of ``k``. + ``value`` can be scalar, array (both dense and sparse) or LinearOperator. + This is accessible as a dict ``{BlochCoeff: value}``. + + Parameters + ---------- + hamiltonian : Model, str, SymPy expression, dict or None (default) + Symbolic representation of a Hamiltonian. If a string, it is + converted to a SymPy expression using ``kwant_continuum.sympify``. + If a dict is provided, it should have the form + ``{symbol: array}`` with all arrays the same size (dense or sparse). + If symbol is not a BlochCoeff, it is passed through sympy.sympify, + and should consist purely of a product of symbolic coefficients, + no constant factors other than 1. `symbol` is then converted to BlochCoeff. + `None` initializes a zero ``BlochModel``. + locals : dict or ``None`` (default) + Additional namespace entries for `~kwant_continuum.sympify`. May be + used to simplify input of matrices or modify input before proceeding + further. For example: + ``locals={'k': 'k_x + I * k_y'}`` or + ``locals={'sigma_plus': [[0, 2], [0, 0]]}``. + momenta : iterable of strings or Sympy symbols + Names of momentum variables, default ``('k_x', 'k_y', 'k_z')`` or + corresponding sympy symbols. Momenta are treated the same as other + keys for the purpose of `keep`. Ignored when initialized with Model. + keep : iterable of BlochCoeff (optional) + Set of symbolic coefficients that are kept, anything that does not + appear here is discarded. Useful for perturbative calculations where + only terms to a given order are needed. By default all keys are kept. + Ignored when initialized with Model. + shape : tuple or None (default) + Shape of the Model, must match the shape of all the values. If not + provided, it is automatically found based on the shape of the input. + Must be provided is ``hamiltonian`` is `None` or ``{}``. Empty tuple + corresponds to scalar values. Ignored when initialized with Model. + format : class or None (default) + Type of the values in the model. Supported types are `np.complex128`, + ``np.ndarray``, ``scipy.sparse.spmatrix`` and ``scipy.sparse.linalg.LinearOperator``. + If ``hamiltonian`` is provided as a dict, all values must be of this type, + except for scalar values, which are recast to ``np.complex128``. + If ``format`` is not provided, it is inferred from the type of the values. + If ``hamiltonian`` is not a dictionary, ``format`` is ignored and set to + ``np.ndarray`` or ``hamiltonian.format`` if it is a ``Model``. + """ def __init__(self, hamiltonian=None, locals=None, momenta=('k_x', 'k_y', 'k_z'), keep=None, shape=None, format=None): - """ - A ``Model`` where coefficients are periodic functions of momenta. - - Internally it is a ``sum(BlochCoeff * value)``, where ``BlochCoeff`` is - a symbolic representation of coefficients and a periodic function of ``k``. - ``value`` can be scalar, array (both dense and sparse) or LinearOperator. - This is accessible as a dict ``{BlochCoeff: value}``. - - Parameters - ---------- - hamiltonian : Model, str, SymPy expression, dict or None (default) - Symbolic representation of a Hamiltonian. If a string, it is - converted to a SymPy expression using ``kwant_continuum.sympify``. - If a dict is provided, it should have the form - ``{symbol: array}`` with all arrays the same size (dense or sparse). - If symbol is not a BlochCoeff, it is passed through sympy.sympify, - and should consist purely of a product of symbolic coefficients, - no constant factors other than 1. `symbol` is then converted to BlochCoeff. - `None` initializes a zero ``BlochModel``. - locals : dict or ``None`` (default) - Additional namespace entries for `~kwant_continuum.sympify`. May be - used to simplify input of matrices or modify input before proceeding - further. For example: - ``locals={'k': 'k_x + I * k_y'}`` or - ``locals={'sigma_plus': [[0, 2], [0, 0]]}``. - momenta : iterable of strings or Sympy symbols - Names of momentum variables, default ``['k_x', 'k_y', 'k_z']`` or - corresponding sympy symbols. Momenta are treated the same as other - keys for the purpose of `keep`. Ignored when initialized with Model. - keep : iterable of BlochCoeff (optional) - Set of symbolic coefficients that are kept, anything that does not - appear here is discarded. Useful for perturbative calculations where - only terms to a given order are needed. By default all keys are kept. - Ignored when initialized with Model. - shape : tuple or None (default) - Shape of the Model, must match the shape of all the values. If not - provided, it is automatically found based on the shape of the input. - Must be provided is ``hamiltonian`` is `None` or ``{}``. Empty tuple - corresponds to scalar values. Ignored when initialized with Model. - format : class or None (default) - Type of the values in the model. Supported types are `np.complex128`, - ``np.ndarray``, ``scipy.sparse.spmatrix`` and ``scipy.sparse.linalg.LinearOperator``. - If ``hamiltonian`` is provided as a dict, all values must be of this type, - except for scalar values, which are recast to ``np.complex128``. - If ``format`` is not provided, it is inferred from the type of the values. - If ``hamiltonian`` is not a dictionary, ``format`` is ignored and set to - ``np.ndarray`` or ``hamiltonian.format`` if it is a ``Model``. - """ momenta = tuple(momenta) if hamiltonian is None: hamiltonian = {} if isinstance(hamiltonian, Model): + if hamiltonian.coordinates is not None or hamiltonian.hopping_vector is not None: + raise NotImplementedError("`coordinates` and `hopping_vector` attributes " + "are not supported in `BlochModel`.") # Use Model's init, only need to recast keys to BlochCoeff super().__init__(hamiltonian=hamiltonian.data, locals=locals, @@ -805,6 +938,9 @@ class BlochModel(Model): keep=keep, shape=shape, format=format)) + @property + def _special_symbols(self): + return {k: getattr(self, k) for k in ('momenta',) if hasattr(self, k)} # Allow getting values using text keys def __getitem__(self, key): @@ -818,14 +954,21 @@ class BlochModel(Model): def transform_symbolic(self, func): raise NotImplementedError('`transform_symbolic` is not implemented for `BlochModel`') - def rotate_momenta(self, R): - """Rotate momenta with rotation matrix R.""" + def rotate(self, R, conjugate=False, transpose=False): + """Rotate momenta with rotation matrix R. + Momenta ``k`` are inverted by conjugate xor transpose. + The matrix part of the model is not conjugated or transposed!""" momenta = self.momenta - assert len(momenta) == R.shape[0], (momenta, R) - # do rotation on hopping vectors with transpose matrix - R_T = np.array(R).astype(float).T result = self.zeros_like() - result.data = {BlochCoeff(R_T @ hop, coeff): copy(val) for (hop, coeff), val in self.items()} + # if conjugated or transposed, need to reverse k + sign = (-1 if (conjugate^transpose) else 1) + if R is not None: + assert len(momenta) == R.shape[0], (momenta, R) + # do rotation on hopping vectors with transpose matrix + R_T = np.array(R).astype(float).T + result.data = {BlochCoeff((sign * R_T) @ hop, coeff): copy(val) for (hop, coeff), val in self.items()} + else: + result.data = {BlochCoeff(sign * hop, coeff): copy(val) for (hop, coeff), val in self.items()} return result def conj(self): @@ -961,8 +1104,6 @@ def _to_bloch_coeff(key, momenta): def _find_momenta(momenta): if any(isinstance(i, int) for i in momenta): raise TypeError('Momenta should be strings or sympy symbols.') - elif all(m in _commutative_momenta for m in momenta): - return tuple(momenta) else: _momenta = [kwant_continuum.sympify(k) for k in momenta] return tuple(kwant_continuum.make_commutative(k, k) diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index b1a67b0c1e9fcf1e5fa3f0a9cb1476b9d069e809..44cd5bc94c3b02652f24f51deedfc39bb4195581 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -6,7 +6,7 @@ from copy import deepcopy from .linalg import matrix_basis, nullspace, split_list, simult_diag, commutator, \ prop_to_id, sparse_basis, mtm, family_to_vectors, solve_mat_eqn, \ - allclose + allclose, spectra_allclose from .model import Model, BlochModel, BlochCoeff from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, \ set_multiply, L_matrices, spin_rotation, time_reversal, \ @@ -19,7 +19,7 @@ from .kwant_linalg_lll import lll, cvp, voronoi def symmetries(model, candidates=None, continuous_rotations=False, generators=False, prettify=False, num_digits=8, verbose=False, - sparse_linalg=False): + sparse_linalg=False, _check_result=False): """ Find symmetries of the Hamiltonian family described by model. @@ -57,6 +57,8 @@ def symmetries(model, candidates=None, continuous_rotations=False, sparse_linalg : bool Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems. + _check_result : bool + Whether to perform consistency checks during the calculation. Returns ------- @@ -75,7 +77,9 @@ def symmetries(model, candidates=None, continuous_rotations=False, '`np.ndarray` or `scipy.sparse.spmatrix` data type.') # Find onsite conserved quantites and projectors onto blocks. - Ps = _reduce_hamiltonian(np.array(list(model.values())), sparse=sparse_linalg) + Ps = _reduce_hamiltonian(np.array(list(model.values())), + sparse=sparse_linalg, + check=_check_result) cont_sym = conserved_quantities(Ps, prettify=prettify, num_digits=num_digits) # Find continuous rotations @@ -92,7 +96,8 @@ def symmetries(model, candidates=None, continuous_rotations=False, if candidates: disc_sym, _ = discrete_symmetries(model, set(candidates), Ps=Ps, generators=generators, verbose=verbose, - sparse_linalg=sparse_linalg) + sparse_linalg=sparse_linalg, + check=_check_result) disc_sym = sorted(list(disc_sym)) else: disc_sym = [] @@ -267,7 +272,7 @@ def symmetry_adapted_sun(gens, check=False): ### Continuous onsite symmetry finding -def _reduce_hamiltonian(H, sparse=False): +def _reduce_hamiltonian(H, sparse=False, check=False): """ Find the unitary symmetries and the symmetry adapted basis of a family of matrices H. In the symmetry adapted basis the matrices are @@ -282,6 +287,8 @@ def _reduce_hamiltonian(H, sparse=False): sparse : bool Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems. + check : bool + Whether to check the final result. Returns ------- @@ -317,7 +324,7 @@ def _reduce_hamiltonian(H, sparse=False): if len(gensr) == 0: P2s = [np.eye(Hr[0].shape[-1])] else: - P2s = symmetry_adapted_sun(gensr) + P2s = symmetry_adapted_sun(gensr, check=check) Preduced += (np.array([np.dot(P, P2i) for P2i in P2s]),) return Preduced @@ -369,7 +376,7 @@ def conserved_quantities(Ps, prettify=False, num_digits=3): ### Point group symmetry finding def discrete_symmetries(model, candidates, Ps=None, generators=False, - verbose=False, sparse_linalg=False): + verbose=False, sparse_linalg=False, check=False): """Find point group symmetries of Hamiltonians family. Optimized version to reduce number of tests, uses sympy exact rotation matrices @@ -390,6 +397,9 @@ def discrete_symmetries(model, candidates, Ps=None, generators=False, Whether to use sparse linear algebra in the calculation. Can give large performance gain in large systems. verbose : bool + Whether to print statistics during run. + check : bool + Whether to check the final result. Returns ------- @@ -406,7 +416,7 @@ def discrete_symmetries(model, candidates, Ps=None, generators=False, m = len(symmetry_candidates) # Reduce Hamiltonian by onsite unitaries if not Ps: - Ps = _reduce_hamiltonian(list(model.values()), sparse=sparse_linalg) + Ps = _reduce_hamiltonian(list(model.values()), sparse=sparse_linalg, check=check) # After every step, symlist is guaranteed to form a group, start with the trivial group e = next(iter(candidates)).identity() e.U = np.eye(Ps[0].shape[1]) @@ -423,10 +433,10 @@ def discrete_symmetries(model, candidates, Ps=None, generators=False, g = min(symmetry_candidates) symmetry_candidates -= {g} # Find unitary part - gr = _find_unitary(model, Ps, g, sparse=sparse_linalg) + gr = _find_unitary(model, Ps, g, sparse=sparse_linalg, check=check) if gr.U is not None: # Check that it's indeed a symmetry - assert gr.apply(model) == model, (n, gr) + assert gr.apply(model).allclose(model, atol=1e-6), (n, gr) genset.add(gr) # Needless to test anything in the group generated by the # symmetries found already, they are symmetries for sure. @@ -456,7 +466,7 @@ def discrete_symmetries(model, candidates, Ps=None, generators=False, return symset, Ps -def _find_unitary(model, Ps, g, sparse=False, checks=False): +def _find_unitary(model, Ps, g, sparse=False, check=False): """Test if the candidate k-space symmetry is (anti)unitary (anti)symmetry, if not, return 'None', if yes, return the unitary part of the transformation 'U'. Checked condition if unitary: U H(inv(R) k) = (+/-) H(k) U @@ -472,7 +482,7 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False): as returned by '_reduce_hamiltonian' g : PointGroupElement Standard representation of symmetry operator. g.U must be None. - checks : bool + check : bool Whether to perform checks. Returns @@ -484,72 +494,57 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False): if not issubclass(model.format, (np.ndarray, scipy.sparse.spmatrix)): raise ValueError('Symmetry finding is only supported for Models with ' '`np.ndarray` or `scipy.sparse.spmatrix` data type.') + # Remove potential small terms + model = model.eliminate_zeros() if g.U is not None: raise ValueError('g.U must be None.') - Rmodel = g.apply(model) - if set(model) != set(Rmodel): + # After eliminating small entries, if g is a symmetry only the same keys should appear + if model.keys() != g.apply(model).eliminate_zeros().keys(): return g - HR, HL = [], [] - # Only test eigenvalues if all matrices are Hermitian - ev_test = True - for key, matL in model.items(): - HR.append(matL) - matR = Rmodel[key] - HL.append(matR) - ev_test = ev_test and allclose(matL, matL.T.conj()) and allclose(matR, matR.T.conj()) - # Need to carry conjugation on left side through P - if g.conjugate: - PsL = [P.conj() for P in Ps] - else: - PsL = Ps - HRs = [[P[0].T.conj() @ hR @ P[0] for hR in HR] for P in Ps] - HLs = [[P[0].T.conj() @ hL @ P[0] for hL in HL] for P in PsL] + # Transform to symmetry adapted basis, then apply the symmetry. This makes sure that the + # conjugation gets applied properly to the basis transformation matrix as well. + reduced_models = _reduced_model(model, Ps, check=check) + reduced_Rmodels = [g.apply(model).eliminate_zeros() for model in reduced_models] squares_to_1 = g * g == g.identity() - block_dict = _find_unitary_blocks(HLs, HRs, Ps, conjugate=g.conjugate, ev_test=ev_test, + # Find candidate blocks of the unitary part + block_dict = _find_unitary_blocks(reduced_Rmodels, reduced_models, + Ps, conjugate=g.conjugate, + transpose=g.transpose, squares_to_1=squares_to_1, sparse=sparse) - S = _construct_unitary(block_dict, Ps, conjugate=g.conjugate, squares_to_1=squares_to_1) + # Try to construct the unitary from the blocks + S = _construct_unitary(block_dict, Ps, conjugate=g.conjugate, + squares_to_1=squares_to_1, transpose=g.transpose) - if checks: - HR, HL = np.array(HR), np.array(HL) - for i, j in it.product(range(len(Ps)), range(len(Ps))): - for a, b in it.product(range(len(Ps[i])), range(len(Ps[j]))): - if i !=j or a!=b: - assert np.allclose(mtm(PsL[i][a].T.conj(), HL, PsL[j][b]), 0) - assert np.allclose(mtm(Ps[i][a].T.conj(), HR, Ps[j][b]), 0) - else: - assert allclose(mtm(PsL[i][a].T.conj(), HL, PsL[j][b]), HLs[i]) - assert allclose(mtm(Ps[i][a].T.conj(), HR, Ps[j][b]), HRs[i]) - if (not g.conjugate) and (not g.antisymmetry) and (S is not None): - assert allclose(S @ HL, HR @ S) + g_new = PointGroupElement(g.R, g.conjugate, g.antisymmetry, S, transpose=g.transpose) + + if check and S is not None: + assert model.allclose(g_new.apply(model)), g - return PointGroupElement(g.R, g.conjugate, g.antisymmetry, S) + return g_new -def _find_unitary_blocks(HLs, HRs, projectors, squares_to_1=True, conjugate=False, - ev_test=True, sparse=False): +def _find_unitary_blocks(modelsL, modelsR, projectors, squares_to_1=True, conjugate=False, + sparse=False, transpose=False): """Find candidate symmetry linear spaces in all blocks. - HLs and HRs are lists of reduced Hamiltonians (families) that go to left and right side + modelsL and modelsR are lists of reduced Hamiltonian Models that go to left and right side of the equations. Returns a dictionary {(i, j): Uij} of all symmetry candidate blocks that have a - nonzero solution of Uij @ HLs[j] = HRs[i] @ Uij for Uij. + nonzero solution of Uij @ modelsL[j] = modelsR[i] @ Uij for Uij. If squares_to_1=True, it is assumed that the operators square is proportional to 1 in every block. The search is limited to j <= i, the diagonal blocks have a phase choice and the off-diagonal blocks with j > i are constructed to ensure squaring to +-1. Otherwise the blocks Uij and Uji are calculated independently. - - If ev_test=True the eigenvalues of the matrices are tested first """ # Only need to find symmetries in half of each block of the Hamiltonian. # We take blocks in the lower triangular half and on the diagonal. block_dict = {} ind = range(len(projectors)) - # Pretest eigenvalues - if ev_test: - evsL = [[la.eigvalsh(h) for h in HLs[i]] for i in ind] - evsR = [[la.eigvalsh(h) for h in HRs[i]] for i in ind] + # Precalculate eigenvalues. + evsL = [{key: la.eigvals(h) for key, h in modelL.items()} for modelL in modelsL] + evsR = [{key: la.eigvals(h) for key, h in modelR.items()} for modelR in modelsR] for (i, j) in it.product(ind, ind): # Only do j <= i if squares to 1 if squares_to_1 and j>i: @@ -557,12 +552,19 @@ def _find_unitary_blocks(HLs, HRs, projectors, squares_to_1=True, conjugate=Fals # Only allowed between blocks of identical shape if projectors[i].shape != projectors[j].shape: continue + # Test that the two reduced models have the same keys + if not modelsL[j].keys() == modelsR[i].keys(): + continue + keys = list(modelsL[j].keys()) # Pretest eigenvalues - if ev_test: - if not allclose(evsL[j], evsR[i]): - continue + if not all(spectra_allclose(evsL[j][key], evsR[i][key]) for key in keys): + continue + # Convert models to lists of 3 index arrays + # need to make sure keys are ordered the same + HR = np.array([modelsR[i][key] for key in keys]) + HL = np.array([modelsL[j][key] for key in keys]) # Find block ij of the symmetry operator - block_dsymm = solve_mat_eqn(HLs[j], HRs[i], hermitian=False, traceless=False, sparse=sparse, k_max=2) + block_dsymm = solve_mat_eqn(HL, HR, hermitian=False, traceless=False, sparse=sparse, k_max=2) # Normalize block_dsymm such that it is close to unitary. The matrix # returned by solve_mat_eqn is normalized such that Tr(X.T.conj() @ X) is close to 1. block_dsymm = np.sqrt(block_dsymm.shape[-1]) * block_dsymm @@ -572,8 +574,7 @@ def _find_unitary_blocks(HLs, HRs, projectors, squares_to_1=True, conjugate=Fals if len(block_dsymm) > 1 or np.isclose(la.det(block_dsymm[0]), 0): raise ValueError('Hamiltonian blocks have residual symmetry.') block_dsymm = block_dsymm[0] - assert allclose(block_dsymm @ HLs[j], - HRs[i] @ block_dsymm) + assert allclose(block_dsymm @ HL, HR @ block_dsymm) # The block must be proportional to a unitary prop_to_I, coeff = prop_to_id(block_dsymm.dot(block_dsymm.T.conj())) assert prop_to_I and np.isclose(np.imag(coeff), 0) and np.real(coeff)>0 @@ -582,16 +583,17 @@ def _find_unitary_blocks(HLs, HRs, projectors, squares_to_1=True, conjugate=Fals block_dict[(i, j)] = block_dsymm # If squares to 1, fill out the lower triangle if squares_to_1: - block_dict[(i, j)], block_dict[(j, i)] = _nice_square(block_dsymm, (i == j), conjugate) + block_dict[(i, j)], block_dict[(j, i)] = _nice_square(block_dsymm, (i == j), + conjugate, transpose=transpose) return block_dict -def _nice_square(block_dsymm, diagonal, conjugate): +def _nice_square(block_dsymm, diagonal, conjugate, transpose=False): # Make sure blocks square to +-1 # Diagonal blocks need proper phase choice if unitary, # nothing to be done if antiunitary, must square to +-1 if diagonal: - if conjugate: + if conjugate^transpose: prop_to_I, coeff = prop_to_id(block_dsymm.dot(block_dsymm.conj())) assert prop_to_I and (np.isclose(coeff, 1) or np.isclose(coeff, -1)) else: @@ -601,13 +603,14 @@ def _nice_square(block_dsymm, diagonal, conjugate): return block_dsymm, block_dsymm # Off-diagonal blocks are chosen such that it squares to +1 else: - if conjugate: + if conjugate^transpose: return block_dsymm, block_dsymm.T else: return block_dsymm, block_dsymm.T.conj() -def _construct_unitary(block_dict, projectors, conjugate=False, squares_to_1=True): +def _construct_unitary(block_dict, projectors, + conjugate=False, squares_to_1=True, transpose=False): """Search for combinations of blocks of the symmetry operator that when combined give a symmetry in canonical form, i.e. with only one nonzero block per row and column, and attempt to construct the operator. """ @@ -634,8 +637,8 @@ def _construct_unitary(block_dict, projectors, conjugate=False, squares_to_1=Tru block = block_dict[(i, j)] # Rebuild full operator using projectors pi, pj = projectors[i], projectors[j] - # Use conjugate projector if antiunitary - if conjugate: + # Use conjugate projector if antiunitary or transposes + if conjugate^transpose: # 'aij, jk, alk -> il' S += np.tensordot(pi @ block, pj, axes=((0, 2), (0, 2))) else: @@ -698,7 +701,9 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li blockdim = rham.shape[0] L = None # Generate all reduced hamiltonians transformed by spatial part - trf_hams = [ContinuousGroupGenerator(1j * R, L).apply(rham) for R in Rs()] + # ignoring small terms + trf_hams = [ContinuousGroupGenerator(1j * R, L).apply(rham).eliminate_zeros() + for R in Rs()] # Generate all keys that appear as transformed reduced hamiltonians keys = rham.keys() for th in trf_hams: @@ -710,10 +715,16 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li Rblocks.append(Rblock) # Iterate over all reduced hamiltonians transformed by only hilbert space part + # ignoring small terms R = None Ls = matrix_basis(blockdim, traceless=True) - trf_hams = [ContinuousGroupGenerator(R, L).apply(rham) for L in Ls] - Lblock = family_to_vectors(trf_hams, all_keys=keys).T + trf_hams = [ContinuousGroupGenerator(R, L).apply(rham).eliminate_zeros() + for L in Ls] + # Need to treat the 1x1 case separately where Ls has zero length + if len(trf_hams) == 0: + Lblock = np.empty((rham['1'].reshape(-1).shape[0] * 2 * len(keys), 0)) + else: + Lblock = family_to_vectors(trf_hams, all_keys=keys).T Lblocks.append(Lblock) # Build constraint matrix: first the spatial blocks in a column, then the @@ -748,9 +759,9 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li return PrettyList(symmetries) -def _reduced_model(model, Ps=None): +def _reduced_model(model, Ps=None, check=False): """ - Construct reduced Hamiltonians in monomial form. + Construct reduced Hamiltonians in Model form. Parameters ---------- @@ -759,6 +770,8 @@ def _reduced_model(model, Ps=None): Ps : list fo 3 index ndarrays projectors 'Ps' returned '_reduce_hamiltonian()' Optional, can be provided to speed up the calculation. + check : bool + Whether to perform checks. Returns ------- @@ -773,8 +786,15 @@ def _reduced_model(model, Ps=None): Ps = _reduce_hamiltonian(np.array([H for H in model.values()])) reduced_hamiltonians = [] for P in Ps: - Hr = P[0].T.conj() @ model @ P[0] + Hr = (P[0].T.conj() @ model @ P[0]).eliminate_zeros() reduced_hamiltonians.append(Hr) + if check: + for i, j in it.product(range(len(Ps)), range(len(Ps))): + for a, b in it.product(range(len(Ps[i])), range(len(Ps[j]))): + if i !=j or a!=b: + assert (Ps[i][a].T.conj() @ model @ Ps[j][b]).allclose(0) + else: + assert (Ps[i][a].T.conj() @ model @ Ps[j][b]).allclose(reduced_hamiltonians[i]) return reduced_hamiltonians ### Bravais lattice symmetry finding diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index 8c86d7693c5c95b11720ed05d87722a05607afe6..7770c4d5881e23f5e27856bfd81517bb8e8d67f6 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -293,6 +293,15 @@ def test_Model(): assert m2[k_x] == m1[k_x].conj() assert m2[e**(-I*k_y)] == m1[e**(I*k_y)].conj() +def test_Model_locals(): + # Test that locals are treated properly + ham1 = Model('alpha * sigma_z') + ham2 = Model('alpha * sz', locals=dict(sz=np.diag([1, -1]))) + assert ham1 == ham2 + ham3 = Model('alpha * sz', locals=dict(sz='[[1, 0], [0, -1]]')) + assert ham2 == ham3 + ham4 = Model('Hz', locals=dict(Hz='[[alpha, 0], [0, -alpha]]')) + assert ham3 == ham4 def test_BlochModel(): diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index 3a57d3a0f8ef6e7f6653ce16f8b9b08d73cb857c..37d81998a2518416521759f6f41cb3ef38814a51 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -5,6 +5,7 @@ import sympy from .. import kwant_rmt from ..symmetry_finder import * from ..symmetry_finder import _reduced_model, _reduce_hamiltonian, bravais_point_group +from ..groups import hexagonal, cubic from ..linalg import * from .. import kwant_continuum @@ -29,7 +30,7 @@ def sumrep(*args): def test_cont_finder(): # Test symmetry adapted basis - gens = sumrep((*2*[0.5*sigma[[3, 1, 2]]])) + gens = sumrep(*2*[0.5*sigma[[3, 1, 2]]]) U2 = kwant_rmt.circular(len(gens[0])) gens2 = np.einsum('ij,ajk,kl->ail',(U2.T).conjugate(),gens,U2) U = symmetry_adapted_sun(gens, check=True) @@ -104,7 +105,7 @@ def test_disc_finder(verbose = False): momenta=[]) # Find symmetry operators sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) if verbose: print('sym', sg) T, P, C = kwant_rmt.t(sym), kwant_rmt.p(sym), kwant_rmt.c(sym) @@ -173,7 +174,7 @@ def test_disc_finder(verbose = False): momenta=[]) # Find symmetry operators sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) if verbose: print('sym', sg) T, P, C = kwant_rmt.t(sym), kwant_rmt.p(sym), kwant_rmt.c(sym) @@ -242,7 +243,7 @@ def test_disc_finder(verbose = False): momenta=[]) # Find symmetry operators sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) if verbose: print('sym', sg) T, P, C = kwant_rmt.t(sym), kwant_rmt.p(sym), kwant_rmt.c(sym) @@ -314,7 +315,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) assert [P.shape for P in Ps] == [(1, 2*dim, dim), (1, 2*dim, dim)] assert len(sg) == 2 tr = list(sg)[0] @@ -352,7 +353,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) assert [P.shape for P in Ps] == [(1, 3*dim, dim), (1, 3*dim, dim), (1, 3*dim, dim)] assert len(sg) == 2 tr = list(sg)[0] @@ -388,7 +389,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) assert [P.shape for P in Ps] == [(2, 2*dim, dim)] assert len(sg) == 2 tr = list(sg)[0] @@ -422,7 +423,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) sg = {time_reversal(0), particle_hole(0), chiral(0)} - sg, Ps = discrete_symmetries(Hs, sg) + sg, Ps = discrete_symmetries(Hs, sg, check=True) assert sorted([P.shape for P in Ps]) == sorted([(1, 4*dim, dim), (1, 4*dim, dim), (2, 4*dim, dim)]) assert len(sg) == 2 tr = list(sg)[0] @@ -474,16 +475,20 @@ def test_continuum(): "alpha * sigma_x * k_x + alpha * sigma_y * k_y + alpha * sigma_z * k_z") # Convert to standard model form H1 = Model(ham1) - sg, Ps = discrete_symmetries(H1, cubic_group) + sg, Ps = discrete_symmetries(H1, cubic_group, check=True) assert [P.shape for P in Ps] == [(1, 2, 2)] assert len(sg) == 48 assert sg == generate_group({C4, C3, TR}) + # Test continuous rotations + sg, cg = symmetries(H1, cubic_group, continuous_rotations=True) + assert set(sg) == generate_group({C4, C3, TR}), len(sg) + assert len(cg) == 3, len(cg) # Add a degeneracy ham2 = "kron(eye(2), " + ham1 + ")" # Convert to standard model form H2 = Model(ham2) - sg, Ps = discrete_symmetries(H2, cubic_group) + sg, Ps = discrete_symmetries(H2, cubic_group, check=True) assert [P.shape for P in Ps] == [(2, 4, 2)] assert len(sg) == 48 assert sg == generate_group({C4, C3, TR}) @@ -492,7 +497,7 @@ def test_continuum(): ham2 = "kron(sigma_z, " + ham1 + ")" # Convert to standard model form H3 = Model(ham2) - sg, Ps = discrete_symmetries(H3, cubic_group) + sg, Ps = discrete_symmetries(H3, cubic_group, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 96 assert sg == generate_group({C4, C3, TR, PH}) @@ -503,22 +508,41 @@ def test_continuum(): # Test sparse H3sp = H3.tocsr() - sg, Ps = discrete_symmetries(H3, cubic_group, sparse_linalg=True) + sg, Ps = discrete_symmetries(H3, cubic_group, sparse_linalg=True, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 96 assert sg == generate_group({C4, C3, TR, PH}) assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 - sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=True) + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=True, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 96 assert sg == generate_group({C4, C3, TR, PH}) assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 - sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=False) + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=False, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 96 assert sg == generate_group({C4, C3, TR, PH}) assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + # Test continuous rotations on constant 1x1 Hamiltonian + H4 = Model('1') + sg, cg = symmetries(H4, cubic_group, continuous_rotations=True) + assert set(sg) == generate_group({I, C4, C3, TR}), len(sg) + # 3 real space rotation generators + assert len(cg) == 3, len(cg) + # Test continuous rotations on constant 2x2 Hamiltonian + H4 = Model('eye(2)') + sg, cg = symmetries(H4, cubic_group, continuous_rotations=True) + assert set(sg) == generate_group({I, C4, C3, TR}), len(sg) + # 3 real space rotation generators + 3 onsite spin rotations + assert len(cg) == 6, len(cg) + # Test continuous rotations on constant 2x2 Hamiltonian with PH + H4 = Model('sigma_z') + sg, cg = symmetries(H4, cubic_group, continuous_rotations=True) + assert set(sg) == cubic_group, len(sg) + # 3 real space rotation generators + 1 onsite spin rotation + assert len(cg) == 4, len(cg) + def test_bloch(): # Simple tests for Bloch models @@ -540,7 +564,7 @@ def test_bloch(): # First example 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))' H6 = Model(ham6, momenta=['k_x', 'k_y']) - sg, Ps = discrete_symmetries(H6, hex_group_2D) + sg, Ps = discrete_symmetries(H6, hex_group_2D, check=True) assert [P.shape for P in Ps] == [(1, 1, 1)] assert len(sg) == 24 assert sg == generate_group({Mx, C6, TR}) @@ -550,7 +574,7 @@ def test_bloch(): 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))' H62 = Model(ham62, momenta=['k_x', 'k_y']) - sg, Ps = discrete_symmetries(H62, hex_group_2D) + sg, Ps = discrete_symmetries(H62, hex_group_2D, check=True) assert [P.shape for P in Ps] == [(1, 2, 2)] assert len(sg) == 24 assert sg == generate_group({Mx, C6, TR}) @@ -558,7 +582,7 @@ def test_bloch(): # Add degeneracy ham63 = 'kron(eye(2), ' + ham62 + ')' H63 = Model(ham63, momenta=['k_x', 'k_y']) - sg, Ps = discrete_symmetries(H63, hex_group_2D) + sg, Ps = discrete_symmetries(H63, hex_group_2D, check=True) assert [P.shape for P in Ps] == [(2, 4, 2)] assert len(sg) == 24 assert sg == generate_group({Mx, C6, TR}) @@ -566,23 +590,23 @@ def test_bloch(): # Add PH states ham64 = 'kron(sigma_z, ' + ham62 + ')' H64 = Model(ham64, momenta=['k_x', 'k_y']) - sg, Ps = discrete_symmetries(H64, hex_group_2D) + sg, Ps = discrete_symmetries(H64, hex_group_2D, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 48 assert sg == generate_group({Mx, C6, TR, PH}) # Test sparse H64 = Model(ham64, momenta=['k_x', 'k_y']) - sg, Ps = discrete_symmetries(H64, hex_group_2D, sparse_linalg=True) + sg, Ps = discrete_symmetries(H64, hex_group_2D, sparse_linalg=True, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 48 assert sg == generate_group({Mx, C6, TR, PH}) Hcsr = H64.tocsr() - sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=True) + sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=True, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 48 assert sg == generate_group({Mx, C6, TR, PH}) - sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=False) + sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=False, check=True) assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] assert len(sg) == 48 assert sg == generate_group({Mx, C6, TR, PH}) @@ -634,3 +658,269 @@ def test_bravais_symmetry(): assert len(group) == n, (name, periods, group, n) group = bravais_point_group(periods @ R, tr=False, ph=False) assert len(group) == n, (name, periods, group, n) + + +def test_3d_dirac(): + # Define 3D Dirac Hamiltonian + ham = """ + (C + D1*k_z**2 + D2*(k_x**2+k_y**2))*eye(4) + + kron(sigma_x,sigma_x)*A2*k_y + + kron(sigma_x,sigma_y)*(-A2*k_x) + + kron(sigma_y,eye(2))*A1*k_z + + kron(sigma_z, eye(2))*(M+B1*k_z**2+B2*(k_x**2+k_y**2)) + + m5*kron(sigma_x,sigma_z) + """ + + H = Model(ham, momenta=['k_x', 'k_y','k_z']) + candidates = hexagonal(dim=3) + sg, cg = symmetries(H, candidates, continuous_rotations=True, prettify=True) + assert len(sg) == 24 + assert len(cg) == 1 + + # Find symmetries with m5 set to 0 + sg, cg = symmetries(H.subs('m5', 0), candidates, continuous_rotations=True, prettify=True) + assert len(sg) == 48 + assert len(cg) == 1 + + +def test_nonhermitian_reduce(): + # Test solve_mat_eqn + n = np.random.randint(2, 5) + h1 = np.random.random((2, n, n)) + 1j * np.random.random((2, n, n)) + h2 = np.random.random((2, n, n)) + 1j * np.random.random((2, n, n)) + + H = np.array([la.block_diag(h1[i], h1[i]) for i in range(len(h1))]) + H = H + (np.einsum('ijk->ikj', H)).conjugate() + assert solve_mat_eqn(H, H, hermitian=True, traceless=True).shape == (3, 2*n, 2*n) + assert solve_mat_eqn(H, H, hermitian=False, traceless=True).shape == (3, 2*n, 2*n) + assert solve_mat_eqn(H, H, hermitian=True, traceless=False).shape == (4, 2*n, 2*n) + + H = np.array([la.block_diag(h1[i], h1[i], h1[i]) for i in range(len(h1))]) + H = H + (np.einsum('ijk->ikj', H)).conjugate() + assert solve_mat_eqn(H, H, hermitian=True, traceless=True).shape == (8, 3*n, 3*n) + assert solve_mat_eqn(H, H, hermitian=False, traceless=True).shape == (8, 3*n, 3*n) + assert solve_mat_eqn(H, H, hermitian=True, traceless=False).shape == (9, 3*n, 3*n) + + H = np.array([la.block_diag(h1[i], h2[i]) for i in range(len(h1))]) + H = H + (np.einsum('ijk->ikj', H)).conjugate() + assert solve_mat_eqn(H, H, hermitian=True, traceless=True).shape == (1, 2*n, 2*n) + assert solve_mat_eqn(H, H, hermitian=False, traceless=True).shape == (1, 2*n, 2*n) + assert solve_mat_eqn(H, H, hermitian=True, traceless=False).shape == (2, 2*n, 2*n) + + # Test _reduce_hamiltonian + Hs = [h1, np.array([la.block_diag(h1[i], h1[i]) for i in range(len(h1))]), + np.array([la.block_diag(h1[i], h2[i]) for i in range(len(h1))])] + for H in Hs: + H = H + (np.einsum('ijk->ikj', H)).conjugate() + Ps = _reduce_hamiltonian(H) + # Check it is proportional to the identity in every block + for P in Ps: + Hr = mtm(P[0].T.conjugate(), H, P[0]) + for i, j in it.product(range(len(P)), repeat=2): + if i ==j: + assert np.allclose(mtm(P[i].T.conjugate(), H, P[j]), Hr) + else: + assert np.allclose(mtm(P[i].T.conjugate(), H, P[j]), 0) + # Check that it vanishes between every block + for P1, P2 in it.combinations(Ps, 2): + P1, P2 = np.hstack(P1), np.hstack(P2) + assert np.allclose(mtm(P1.T.conjugate(), H, P2), 0) + +def test_nonhermitian_continuum(): + # Simple tests for nonhermitian continuum models + + # Cubic point group with PH, TR and hermiticity + eye = np.array(np.eye(3), int) + E = PointGroupElement(eye, False, False, None) + I = PointGroupElement(-eye, False, False, None) + C4 = PointGroupElement(np.array([[1, 0, 0], + [0, 0, 1], + [0, -1, 0]], int), False, False, None) + C3 = PointGroupElement(np.array([[0, 0, 1], + [1, 0, 0], + [0, 1, 0]], int), False, False, None) + TR = time_reversal(realspace_dim=3) + PH = particle_hole(realspace_dim=3) + herm = PointGroupElement(eye, True, False, None, transpose=True) + cubic_gens = {I, C4, C3, TR, PH, herm} + cubic_group = generate_group(cubic_gens) + assert len(cubic_group) == 384 + + ## Check that hermiticity is foundd in hermitian models + # First model + 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 model form + H1 = Model(ham1) + sg, Ps = discrete_symmetries(H1, cubic_group, check=True) + assert [P.shape for P in Ps] == [(1, 2, 2)] + assert len(sg) == 96 + assert sg == generate_group({C4, C3, TR, herm}) + + # Add a degeneracy + ham2 = "kron(eye(2), " + ham1 + ")" + # Convert to standard model form + H2 = Model(ham2) + sg, Ps = discrete_symmetries(H2, cubic_group, check=True) + assert [P.shape for P in Ps] == [(2, 4, 2)] + assert len(sg) == 96 + assert sg == generate_group({C4, C3, TR, herm}) + + # Add hole degrees of freedom + ham2 = "kron(sigma_z, " + ham1 + ")" + # Convert to standard model form + H3 = Model(ham2) + sg, Ps = discrete_symmetries(H3, cubic_group, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 192 + assert sg == generate_group({C4, C3, TR, PH, herm}) + + # Continuous rotation symmetry + for H in [H1, H2, H3]: + assert len(continuous_symmetries(H)) == 3 + + # Test sparse + H3sp = H3.tocsr() + sg, Ps = discrete_symmetries(H3, cubic_group, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 192 + assert sg == generate_group({C4, C3, TR, PH, herm}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 192 + assert sg == generate_group({C4, C3, TR, PH, herm}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=False, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 192 + assert sg == generate_group({C4, C3, TR, PH, herm}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + + # Check that unitary symmetries work in nonhermitian models + # Rotation by a generic complex phase should remove non-unitary + # symmetries + # Don't use herm, as these have some weird transposition symmetries + cubic_gens = {I, C4, C3, TR, PH} + cubic_group = generate_group(cubic_gens) + assert len(cubic_group) == 192 + + # First model + 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 model form + # twist it with a complex number + H1 = (2 + 1j) * Model(ham1) + sg, Ps = discrete_symmetries(H1, cubic_group, check=True) + assert [P.shape for P in Ps] == [(1, 2, 2)] + assert len(sg) == 24 + assert sg == generate_group({C4, C3}) + + # Add a degeneracy + ham2 = "kron(eye(2), " + ham1 + ")" + # Convert to standard model form + H2 = (2 + 1j) * Model(ham2) + sg, Ps = discrete_symmetries(H2, cubic_group, check=True) + assert [P.shape for P in Ps] == [(2, 4, 2)] + assert len(sg) == 24 + assert sg == generate_group({C4, C3}) + + # Add hole degrees of freedom, adds a chiral symmetry with sigma_x + ham2 = "kron(sigma_z, " + ham1 + ")" + # Convert to standard model form + H3 = (2 + 1j) * Model(ham2) + sg, Ps = discrete_symmetries(H3, cubic_group, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 48, len(sg) + assert sg == generate_group({C4, C3, PH * TR}) + + # Continuous rotation symmetry + for H in [H1, H2, H3]: + assert len(continuous_symmetries(H)) == 3 + + # Test sparse + H3sp = H3.tocsr() + sg, Ps = discrete_symmetries(H3, cubic_group, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 48 + assert sg == generate_group({C4, C3, PH * TR}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 48 + assert sg == generate_group({C4, C3, PH * TR}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + sg, Ps = discrete_symmetries(H3sp, cubic_group, sparse_linalg=False, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 48 + assert sg == generate_group({C4, C3, PH * TR}) + assert len(continuous_symmetries(H, sparse_linalg=True)) == 3 + + +def test_nonhermitian_bloch(): + # Simple tests for Bloch models with only unitary symmetries + + # Hexagonal point group + eyesym = sympy.ImmutableMatrix(sympy.eye(2)) + Mx = PointGroupElement(sympy.ImmutableMatrix([[-1, 0], + [0, 1]]), + False, False, None) + C6 = PointGroupElement(sympy.ImmutableMatrix([[sympy.Rational(1, 2), sympy.sqrt(3)/2], + [-sympy.sqrt(3)/2, sympy.Rational(1, 2)]]), + False, False, None) + TR = time_reversal(realspace_dim=2) + PH = particle_hole(realspace_dim=2) + gens_hex_2D ={Mx, C6, TR, PH} + hex_group_2D = generate_group(gens_hex_2D) + assert len(hex_group_2D) == 48 + + # First example + 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))' + # with complex twist + H6 = (3 - 2j) * Model(ham6, momenta=['k_x', 'k_y']) + sg, Ps = discrete_symmetries(H6, hex_group_2D, check=True) + assert [P.shape for P in Ps] == [(1, 1, 1)] + assert len(sg) == 12 + assert sg == generate_group({Mx, C6}) + + # extend model to add SOC + 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))' + H62 = (3 - 2j) * Model(ham62, momenta=['k_x', 'k_y']) + sg, Ps = discrete_symmetries(H62, hex_group_2D, check=True) + assert [P.shape for P in Ps] == [(1, 2, 2)] + assert len(sg) == 12 + assert sg == generate_group({Mx, C6}) + + # Add degeneracy + ham63 = 'kron(eye(2), ' + ham62 + ')' + H63 = (3 - 2j) * Model(ham63, momenta=['k_x', 'k_y']) + sg, Ps = discrete_symmetries(H63, hex_group_2D, check=True) + assert [P.shape for P in Ps] == [(2, 4, 2)] + assert len(sg) == 12 + assert sg == generate_group({Mx, C6}) + + # Add PH states + ham64 = 'kron(sigma_z, ' + ham62 + ')' + H64 = (3 - 2j) * Model(ham64, momenta=['k_x', 'k_y']) + sg, Ps = discrete_symmetries(H64, hex_group_2D, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 24 + assert sg == generate_group({Mx, C6, PH * TR}), sg + + # Test sparse + H64 = (3 - 2j) * Model(ham64, momenta=['k_x', 'k_y']) + sg, Ps = discrete_symmetries(H64, hex_group_2D, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 24 + assert sg == generate_group({Mx, C6, PH * TR}) + Hcsr = H64.tocsr() + sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=True, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 24 + assert sg == generate_group({Mx, C6, PH * TR}) + sg, Ps = discrete_symmetries(Hcsr, hex_group_2D, sparse_linalg=False, check=True) + assert [P.shape for P in Ps] == [(1, 4, 2), (1, 4, 2)] + assert len(sg) == 24 + assert sg == generate_group({Mx, C6, PH * TR}) \ No newline at end of file diff --git a/qsymm/tests/test_util.py b/qsymm/tests/test_util.py index 97dc776f1cb2d16995588b23de5fd0e577eada21..a1dcddee1cb1b9b8f9d280abc8618af0b462f1f9 100644 --- a/qsymm/tests/test_util.py +++ b/qsymm/tests/test_util.py @@ -34,6 +34,8 @@ def test_spatial_types(): np.eye(3)) S3 = PointGroupElement(np.eye(2), False, False, 1j * np.eye(3)) + S4 = PointGroupElement(None, False, False, + np.eye(3)) C6s = PointGroupElement(sympy.ImmutableMatrix( [[sympy.Rational(1, 2), sympy.sqrt(3)/2], [-sympy.sqrt(3)/2, sympy.Rational(1, 2)]] @@ -43,8 +45,11 @@ def test_spatial_types(): [-np.sqrt(3)/2, 1/2]] )) + assert S1 == S4 assert S2**2 == S1 + assert S2**2 == S4 assert not S1 == S2 + assert not S4 == S2 assert S1 == S3 assert C6s == C6f # Mixing sympy with other types raises an error