From 79006383b80202d5123d71af71438c27e1ea2603 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 16 May 2019 22:33:34 +0200 Subject: [PATCH 01/41] implement matmul in Model --- qsymm/groups.py | 4 +- qsymm/model.py | 125 +++++++++++++++++++++++++-------------- qsymm/symmetry_finder.py | 2 +- 3 files changed, 83 insertions(+), 48 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 5152b20..e9ab6c6 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -251,7 +251,7 @@ class PointGroupElement: if antisymmetry: result = -result if U is not None: - result = U * result * U.T.conj() + result = U @ result @ U.T.conj() return result @@ -344,7 +344,7 @@ class ContinuousGroupGenerator: for i, j in it.product(range(dim), repeat=2)]) result += 1j * monomials.transform_symbolic(trf) if U_nonzero: - result += monomials * (1j*U) + (-1j*U) * monomials + result += monomials @ (1j*U) + (-1j*U) @ monomials return result diff --git a/qsymm/model.py b/qsymm/model.py index af7a6ba..5a9b607 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -108,7 +108,7 @@ class Model(UserDict): # Make it work with numpy arrays __array_ufunc__ = None - def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2]): + def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], interesting_keys=None): """General class to store Hamiltonian families. Can be used to efficiently store any matrix valued function. Implements many sympy and numpy methods. Arithmetic operators are overloaded, @@ -132,6 +132,10 @@ class Model(UserDict): momenta : list of int or list of Sympy objects Indices of momenta the monomials depend on from 'k_x', 'k_y' and 'k_z' or a list of names for the momentum variables. + interesting_keys : iterable of sympy expressions or None (default) + 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. """ self.momenta = _find_momenta(momenta) if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): @@ -154,6 +158,10 @@ class Model(UserDict): self.data = monomials self.shape = _find_shape(self.data) + if interesting_keys is not None: + self.interesting_keys = set(interesting_keys) + else: + self.interesting_keys = set() # Restructure # Clean internal data by: @@ -204,23 +212,28 @@ class Model(UserDict): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. # Define addition of 0 and {} + result = self.zeros_like() if not other: - result = self.copy() + result.data = self.data.copy() # If self is empty return other elif not self and isinstance(other, type(self)): - result = other.copy() + result = other.zeros_like() + result.data = other.data.copy() elif isinstance(other, type(self)): if self.momenta != other.momenta: raise ValueError("Can only add Models with the same momenta") - result = self.copy() - for key, val in other.items(): - if allclose(result[key], -val): - try: - del result[key] - except KeyError: - pass - else: - result[key] += val + for key in self.keys() & other.keys(): + total = self[key] + other[key] + # If only one is sparse matrix, the result is np.matrix, recast it to np.ndarray + if isinstance(total, np.matrix): + total = total.A + # Only add dense matrix if it is nonzero + if not isinstance(total, np.ndarray) or not allclose(total, 0): + result[key] = total + for key in self.keys() - other.keys(): + result[key] = copy(self[key]) + for key in other.keys() - self.keys(): + result[key] = copy(other[key]) else: raise NotImplementedError('Addition of {} with {} not supported'.format(type(self), type(other))) return result @@ -243,33 +256,27 @@ class Model(UserDict): def __mul__(self, other): # Multiplication by numbers, sympy symbols, arrays and Model + result = self.zeros_like() if isinstance(other, Number): - if np.isclose(other, 0): - result = self.zeros_like() - else: - result = self.zeros_like() - result.data = {key: other * val for key, val in self.items()} + result.data = {key: val * other for key, val in self.items()} elif isinstance(other, Basic): - result = type(self)({key * other: val for key, val in self.items()}) - result.momenta = self.momenta - elif isinstance(other, np.ndarray): - result = self.copy() - for key, val in self.items(): - prod = np.dot(val, other) - if np.allclose(prod, 0): - del result[key] - else: - result[key] = prod - result.shape = _mul_shape(self.shape, other.shape) - elif isinstance(other, type(self)): + result.data = {key * other: val for key, val in self.items() + if (key * other in interesting_keys or not interesting_keys)} + elif isinstance(other, Model): if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") - result = sum(type(self)({k1 * k2: np.dot(v1, v2)}) - for (k1, v1), (k2, v2) in it.product(self.items(), other.items())) + interesting_keys = self.interesting_keys | other.interesting_keys + result = sum(Model({k1 * k2: v1 * v2}, + interesting_keys=interesting_keys) + for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) + if (k1 * k2 in interesting_keys or not interesting_keys)) result.momenta = list(set(self.momenta) | set(other.momenta)) - result.shape = _mul_shape(self.shape, other.shape) + # Find out the shape of the result even if it is empty + result.shape = _find_shape(result.data) if result.data else (self[1] * other[1]).shape else: - raise NotImplementedError('Multiplication of {} with {} not supported'.format(type(self), type(other))) + # Otherwise try to multiply every value with other + result.data = {key: val * other for key, val in self.items()} + result.shape = _find_shape(result.data) if result.data else (self[1] * other).shape return result def __rmul__(self, other): @@ -277,26 +284,46 @@ class Model(UserDict): if isinstance(other, Number): result = self.__mul__(other) elif isinstance(other, Basic): - result = type(self)({other * key: val for key, val in self.items()}) - result.momenta = self.momenta - elif isinstance(other, np.ndarray): result = self.zeros_like() - for key, val in self.items(): - prod = np.dot(other, val) - if np.allclose(prod, 0): - del result[key] - else: - result[key] = prod - result.shape = _mul_shape(other.shape, self.shape) + result.data = {other * key: val for key, val in self.items()} + else: + # Otherwise try to multiply every value with other + result = self.zeros_like() + result.data = {key: other * val for key, val in self.items()} + result.shape = _find_shape(result.data) if result.data else (other * self[1]).shape + return result + + def __matmul__(self, other): + # Multiplication by arrays and PerturbativeModel + if isinstance(other, Model): + if self.momenta != other.momenta: + raise ValueError("Can only multiply Models with the same momenta") + interesting_keys = self.interesting_keys | other.interesting_keys + result = sum(Model({k1 * k2: v1 @ v2}, + interesting_keys=interesting_keys) + for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) + if (k1 * k2 in interesting_keys or not interesting_keys)) + result.momenta = list(set(self.momenta) | set(other.momenta)) + result.shape = _find_shape(result.data) if result.data else (self[1] @ other[1]).shape else: - raise NotImplementedError('Multiplication with type {} not implemented'.format(type(other))) + # Otherwise try to multiply every value with other + result = self.zeros_like() + result.data = {key: val @ other for key, val in self.items()} + result.shape = _find_shape(result.data) if result.data else (self[1] @ other).shape + return result + + def __rmatmul__(self, other): + # Left multiplication by arrays + result = self.zeros_like() + result.data = {key: other @ val for key, val in self.items()} + result.shape = _find_shape(result.data) if result.data else (other @ self[1]).shape return result def __truediv__(self, other): result = self.zeros_like() if isinstance(other, Number): - result.data = {key : val / other for key, val in self.items()} + result.data = {key : val * (1/other) for key, val in self.items()} else: raise TypeError( "unsupported operand type for /: {} and {}".format(type(self), type(other))) @@ -312,6 +339,7 @@ class Model(UserDict): def zeros_like(self): """Return an empty model object that inherits the other properties""" result = type(self)() + result.interesting_keys = self.interesting_keys.copy() result.momenta = self.momenta.copy() result.shape = self.shape return result @@ -586,6 +614,13 @@ class BlochModel(Model): return Model({key.tosympy(self.momenta, nsimplify=nsimplify): val for key, val in self.items()}, momenta=self.momenta) + def zeros_like(self): + """Return an empty model object that inherits the other properties""" + result = type(self)() + result.momenta = self.momenta.copy() + result.shape = self.shape + return result + def _to_bloch_coeff(key, momenta): """Transform sympy expression to BlochCoeff if possible.""" diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index 76b6520..a6e5e27 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -756,7 +756,7 @@ 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] reduced_hamiltonians.append(Hr) return reduced_hamiltonians -- GitLab From 962c80fe223f31fff1332185387f414e674093d9 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 17 May 2019 00:52:20 +0200 Subject: [PATCH 02/41] migrate more functionality to Model --- qsymm/model.py | 136 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 88 insertions(+), 48 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index 5a9b607..8478c63 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -1,4 +1,5 @@ import numpy as np +import scipy import tinyarray as ta import scipy.linalg as la import itertools as it @@ -109,10 +110,17 @@ class Model(UserDict): __array_ufunc__ = None def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], interesting_keys=None): - """General class to store Hamiltonian families. - Can be used to efficiently store any matrix valued function. - Implements many sympy and numpy methods. Arithmetic operators are overloaded, - such that * corresponds to matrix multiplication. + """ + General class to efficiently store any matrix valued function. + The internal structure is a dict with {symbol: value}, where + symbol is a sympy expression, the object representing sum(symbol * value). + The values can be scalars, arrays (both dense and sparse) or LinearOperators. + Implements many sympy and numpy methods and arithmetic operators. + Multiplication is distributed over the sum, * is passed down to + both symbols and values, @ is passed to symbols as * and to values + as @. Assumes that symbols form a commutative group. + Enhances the functionality of Model by allowing interesting_keys to be + specified, symbols not listed there are discarded. Parameters ---------- @@ -120,7 +128,7 @@ class Model(UserDict): Symbolic representation of a Hamiltonian. It is converted to a SymPy expression using kwant_continuum.sympify. If a dict is provided, it should have the form - {expression: np.ndarray} with all arrays the same size. + {expression: array} with all arrays the same size (dense or sparse). expression will be passed through sympy.sympify, and should consist purely of symbolic coefficients, no constant factors other than 1. locals : dict or None (default) @@ -129,18 +137,28 @@ class Model(UserDict): further. For example: locals={'k': 'k_x + I * k_y'} or locals={'sigma_plus': [[0, 2], [0, 0]]}. - momenta : list of int or list of Sympy objects - Indices of momenta the monomials depend on from 'k_x', 'k_y' and 'k_z' - or a list of names for the momentum variables. - interesting_keys : iterable of sympy expressions or None (default) + interesting_keys : iterable of expressions or None (default) 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 : list of int or list of Sympy objects + Indices of momentum variables from ['k_x', 'k_y', 'k_z'] + or a list of names for the momentum variables as sympy symbols. + Momenta are treated the same as other keys for the purpose of + interesting_keys, need to list interesting powers of momenta. + """ self.momenta = _find_momenta(momenta) + + if interesting_keys is not None: + self.interesting_keys = {sympy.sympify(k) for k in interesting_keys} + else: + self.interesting_keys = set() + if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): # Initialize as dict sympifying the keys - super().__init__({sympy.sympify(k): v for k, v in hamiltonian.items()}) + super().__init__({sympy.sympify(k): v for k, v in hamiltonian.items() + if sympy.sympify(k) in self.interesting_keys or not self.interesting_keys}) else: hamiltonian = kwant_continuum.sympify(hamiltonian, locals=locals) if not isinstance(hamiltonian, sympy.matrices.MatrixBase): @@ -157,12 +175,6 @@ class Model(UserDict): if not np.allclose(v, 0)} self.data = monomials - self.shape = _find_shape(self.data) - if interesting_keys is not None: - self.interesting_keys = set(interesting_keys) - else: - self.interesting_keys = set() - # Restructure # Clean internal data by: # * splitting summands in keys @@ -193,18 +205,35 @@ class Model(UserDict): # overwrite internal data self.data = new_data + self.shape = _find_shape(self.data) + + # Keep track of whether this is a dense array + self._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + # Defaultdict functionality def __missing__(self, key): if self.shape is not None: - return np.zeros(self.shape, dtype=complex) + if self.shape == (): + #scalar + return 0 + elif self._isarray: + # Return dense zero array if dense + return np.zeros(self.shape, dtype=complex) + else: + # Otherwise return a csr_matrix + return scipy.sparse.csr_matrix(self.shape, dtype=complex) else: return None def __eq__(self, other): if not set(self) == set(other): return False - for k, v in self.data.items(): - if not allclose(v, other[k]): + if other == {} == self.data: + return True + a = self.toarray() + b = other.toarray() + for key in a.keys() | b.keys(): + if not allclose(a[key], b[key]): return False return True @@ -266,13 +295,14 @@ class Model(UserDict): if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") interesting_keys = self.interesting_keys | other.interesting_keys - result = sum(Model({k1 * k2: v1 * v2}, - interesting_keys=interesting_keys) - for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) - if (k1 * k2 in interesting_keys or not interesting_keys)) - result.momenta = list(set(self.momenta) | set(other.momenta)) + result = sum(type(self)({k1 * k2: v1 * v2}, interesting_keys=interesting_keys) + for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) + if (k1 * k2 in interesting_keys or not interesting_keys)) # Find out the shape of the result even if it is empty - result.shape = _find_shape(result.data) if result.data else (self[1] * other[1]).shape + if result is 0: + result = self.zeros_like() + result.shape = (self[1] * other[1]).shape + result.momenta = list(set(self.momenta) | set(other.momenta)) else: # Otherwise try to multiply every value with other result.data = {key: val * other for key, val in self.items()} @@ -299,12 +329,14 @@ class Model(UserDict): if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") interesting_keys = self.interesting_keys | other.interesting_keys - result = sum(Model({k1 * k2: v1 @ v2}, - interesting_keys=interesting_keys) - for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) - if (k1 * k2 in interesting_keys or not interesting_keys)) + result = sum(type(self)({k1 * k2: v1 @ v2}, interesting_keys=interesting_keys) + for (k1, v1), (k2, v2) in it.product(self.items(), other.items()) + if (k1 * k2 in interesting_keys or not interesting_keys)) + # Find out the shape of the result even if it is empty + if result is 0: + result = self.zeros_like() + result.shape = (self[1] * other[1]).shape result.momenta = list(set(self.momenta) | set(other.momenta)) - result.shape = _find_shape(result.data) if result.data else (self[1] @ other[1]).shape else: # Otherwise try to multiply every value with other result = self.zeros_like() @@ -342,6 +374,7 @@ class Model(UserDict): result.interesting_keys = self.interesting_keys.copy() result.momenta = self.momenta.copy() result.shape = self.shape + result._isarray = self._isarray return result def transform_symbolic(self, func): @@ -465,18 +498,36 @@ class Model(UserDict): # Return sympy representation of the term # If nsimplify=True, attempt to rewrite numerical coefficients as exact formulas if not nsimplify: - result = sympy.sympify(sum(key * val for key, val in self.items())) + result = sympy.sympify(sum(key * val for key, val in self.toarray().items())) else: # Vectorize nsimplify vnsimplify = np.vectorize(sympy.nsimplify, otypes=[object]) result = sympy.MatAdd(*[key * sympy.Matrix(vnsimplify(val)) - for key, val in self.items()]).doit() + for key, val in self.toarray().items()]).doit() if any([isinstance(result, matrix_type) for matrix_type in (sympy.MatrixBase, sympy.ImmutableDenseMatrix, sympy.ImmutableDenseNDimArray)]): result = sympy.Matrix(result).reshape(*result.shape) return result + def evalf(self, subs=None): + return sum(float(key.evalf(subs=subs)) * val for key, val in self.items()) + + def tocsr(self): + result = self.zeros_like() + result.data = {key: scipy.sparse.csr_matrix(val, dtype=complex) + for key, val in self.items()} + result._isarray = False + return result + + def toarray(self): + result = self.zeros_like() + result.data = {key : (val if isinstance(val, np.ndarray) or isinstance(val, Number) + else val.toarray()) + for key, val in self.items()} + result._isarray = True + return result + def copy(self): result = self.zeros_like() # This is faster than deepcopy of the dict @@ -579,7 +630,12 @@ class BlochModel(Model): else: # Use Model to parse input self.__init__(Model(hamiltonian, locals, momenta)) + self.shape = _find_shape(self.data) + self.interesting_keys = set() + + # Keep track of whether this is a dense array + self._isarray = any(isinstance(val, np.ndarray) for val in self.values()) def transform_symbolic(self, func): raise NotImplementedError('transform_symbolic is not implemented for BlochModel') @@ -614,13 +670,6 @@ class BlochModel(Model): return Model({key.tosympy(self.momenta, nsimplify=nsimplify): val for key, val in self.items()}, momenta=self.momenta) - def zeros_like(self): - """Return an empty model object that inherits the other properties""" - result = type(self)() - result.momenta = self.momenta.copy() - result.shape = self.shape - return result - def _to_bloch_coeff(key, momenta): """Transform sympy expression to BlochCoeff if possible.""" @@ -744,12 +793,3 @@ def _find_shape(data): if not all([v.shape == shape for v in data.values()]): raise ValueError('All terms must have the same shape') return shape - -def _mul_shape(shape1, shape2): - # Find the shape of the product - if shape1 == (): - return shape2 - elif shape2 == (): - return shape1 - else: - return (shape1[0], shape2[-1]) -- GitLab From 9931b3155ade0ac74ca9c0f5173f4dc6c26849ec Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 17 May 2019 21:16:34 +0200 Subject: [PATCH 03/41] rework model equality testing and model add --- qsymm/linalg.py | 6 +++--- qsymm/model.py | 30 ++++++++++++++---------------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/qsymm/linalg.py b/qsymm/linalg.py index a516427..f19ed27 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -33,14 +33,14 @@ def allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): """Alternative to numpy.allclose to test that two ndarrays are elementwise close. Unlike the numpy version, the relative tolerance is not elementwise, but it is relative to - the largest entry in the array.""" + the largest entry in the arrays.""" a = np.asarray(a) b = np.asarray(b) # Check if empty arrays, only compare shape if a.size == 0: return a.shape == b.shape - atol = atol + rtol * np.max(np.abs(a)) - return np.allclose(a, b, rtol=0, atol=atol, equal_nan=equal_nan) + atol = atol + rtol * max(np.max(np.abs(a)), np.max(np.abs(b))) + return np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan) def mtm(a, B, c): diff --git a/qsymm/model.py b/qsymm/model.py index 8478c63..845d9a6 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -160,6 +160,7 @@ class Model(UserDict): super().__init__({sympy.sympify(k): v for k, v in hamiltonian.items() if sympy.sympify(k) in self.interesting_keys or not self.interesting_keys}) else: + # Try to parse the input with kwant_continuum.sympify hamiltonian = kwant_continuum.sympify(hamiltonian, locals=locals) if not isinstance(hamiltonian, sympy.matrices.MatrixBase): hamiltonian = sympy.Matrix([[hamiltonian]]) @@ -226,39 +227,36 @@ class Model(UserDict): return None def __eq__(self, other): - if not set(self) == set(other): - return False - if other == {} == self.data: - return True a = self.toarray() - b = other.toarray() - for key in a.keys() | b.keys(): - if not allclose(a[key], b[key]): - return False + if other == {} or other == 0: + if self.data == {}: + return True + else: + return all(allclose(a[key], 0) for key in a.keys()) + else: + b = other.toarray() + return all(allclose(a[key], b[key]) for key in a.keys() | b.keys()) return True def __add__(self, other): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. # Define addition of 0 and {} - result = self.zeros_like() if not other: - result.data = self.data.copy() + result = self.copy() # If self is empty return other elif not self and isinstance(other, type(self)): - result = other.zeros_like() - result.data = other.data.copy() + result = other.copy() elif isinstance(other, type(self)): if self.momenta != other.momenta: raise ValueError("Can only add Models with the same momenta") + result = self.zeros_like() for key in self.keys() & other.keys(): total = self[key] + other[key] # If only one is sparse matrix, the result is np.matrix, recast it to np.ndarray if isinstance(total, np.matrix): total = total.A - # Only add dense matrix if it is nonzero - if not isinstance(total, np.ndarray) or not allclose(total, 0): - result[key] = total + result[key] = total for key in self.keys() - other.keys(): result[key] = copy(self[key]) for key in other.keys() - self.keys(): @@ -315,7 +313,7 @@ class Model(UserDict): result = self.__mul__(other) elif isinstance(other, Basic): result = self.zeros_like() - result.data = {other * key: val for key, val in self.items()} + result.data = {other * key: copy(val) for key, val in self.items()} else: # Otherwise try to multiply every value with other result = self.zeros_like() -- GitLab From a9b55ce10dca5de992127c81eef08b0c787a9f22 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 17 May 2019 22:18:33 +0200 Subject: [PATCH 04/41] add more tests for Model --- qsymm/model.py | 29 +++++++++++++++++------------ qsymm/tests/test_model.py | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 12 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index 845d9a6..90d731e 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -227,16 +227,8 @@ class Model(UserDict): return None def __eq__(self, other): - a = self.toarray() - if other == {} or other == 0: - if self.data == {}: - return True - else: - return all(allclose(a[key], 0) for key in a.keys()) - else: - b = other.toarray() - return all(allclose(a[key], b[key]) for key in a.keys() | b.keys()) - return True + # Call allclose with default tolerances + return self.allclose(other) def __add__(self, other): # Addition of Models. It is assumed that both Models are @@ -300,7 +292,7 @@ class Model(UserDict): if result is 0: result = self.zeros_like() result.shape = (self[1] * other[1]).shape - result.momenta = list(set(self.momenta) | set(other.momenta)) + result.momenta = self.momenta else: # Otherwise try to multiply every value with other result.data = {key: val * other for key, val in self.items()} @@ -334,7 +326,7 @@ class Model(UserDict): if result is 0: result = self.zeros_like() result.shape = (self[1] * other[1]).shape - result.momenta = list(set(self.momenta) | set(other.momenta)) + result.momenta = self.momenta else: # Otherwise try to multiply every value with other result = self.zeros_like() @@ -578,6 +570,19 @@ class Model(UserDict): result.shape = _find_shape(result.data) return result + def allclose(self, other, rtol=1e-05, atol=1e-08, equal_nan=False): + # Test whether two arrays are approximately equal + a = self.toarray() + if other == {} or other == 0: + if self.data == {}: + return True + else: + return all(allclose(val, 0, rtol, atol, equal_nan) for val in a.values()) + else: + b = other.toarray() + return all(allclose(a[key], b[key], rtol, atol, equal_nan) + for key in a.keys() | b.keys()) + class BlochModel(Model): def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2]): diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index a025bbb..545b6db 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -13,6 +13,39 @@ k_x, k_y = _commutative_momenta[:2] momenta = _commutative_momenta[:2] c0, c1 = sympy.symbols('c0 c1', real=True) +def test_dense_algebra(): + # scalar models + a = Model({1: 1, 'x': 3}) + assert a * 0 == 0 * a + assert 0 * a == a.zeros_like() + assert 0 * a == {} + assert a - a == {} + assert a + a == 2 * a + assert a / 2 == a * 0.5 + a += a + assert a == Model({1: 2, 'x': 6}) + a /= 2 + assert a == Model({1: 1, 'x': 3}) + b = Model({1: 2, 'x**2': 3}) + assert a + b == Model({1: 3, 'x': 3, 'x**2': 3}) + assert a * b == Model({1: 2, 'x': 6, 'x**2': 3, 'x**3': 9}) + assert (a - a) * b == {} + assert (a - a) * (b - b) == {} + # vector models + vec = np.ones((3,)) + c = a * vec + mat = np.ones((3, 3)) + d = b * mat + assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) + assert c.T() @ d == 3 * a * b + assert c.T() @ d @ c == 9 * a * a * b + # numpy broadcasting rules apply + assert d + c == (a + b) * mat + assert d * c == a * b * mat + assert d.trace() == 3 * b + assert d.reshape(-1) == b * np.ones((9,)) + + def test_to_bloch_coeff(): key = sympy.sqrt(5)*e**(2*c0)*e**(I*(k_x/2 + np.sqrt(3)*k_y + c1)) -- GitLab From a029cfb5ced6c4e6511b15904c7cfa9650559d8f Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sat, 18 May 2019 13:13:49 +0200 Subject: [PATCH 05/41] finish reworking Model --- qsymm/groups.py | 11 +-- qsymm/hamiltonian_generator.py | 4 +- qsymm/linalg.py | 2 +- qsymm/model.py | 151 ++++++++++++++++++--------------- qsymm/tests/test_model.py | 7 +- 5 files changed, 97 insertions(+), 78 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index e9ab6c6..d849f8f 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -326,25 +326,26 @@ class ContinuousGroupGenerator: def __repr__(self): return 'ContinuousGroupGenerator(\n{},\n{},\n)'.format(self.R, self.U) - def apply(self, monomials): + def apply(self, model): """Return copy of monomials with applied infinitesimal generator. 1j * (H(k) U - U H(k)) + 1j * dH(k)/dk_i R_{ij} k_j """ R, U = self.R, self.U - momenta = monomials.momenta + 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 = monomials.zeros_like() + 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 it.product(range(dim), repeat=2)]) - result += 1j * monomials.transform_symbolic(trf) + result += 1j * model.transform_symbolic(trf) + assert result._isarray == model._isarray if U_nonzero: - result += monomials @ (1j*U) + (-1j*U) @ monomials + result += model @ (1j*U) + (-1j*U) @ model return result diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py index f416dc4..567bde3 100644 --- a/qsymm/hamiltonian_generator.py +++ b/qsymm/hamiltonian_generator.py @@ -658,13 +658,13 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, # Hopping direction in real space # Dot product with momentum vector phase = sum([coordinate * momentum for coordinate, momentum in - zip(vec, _commutative_momenta[:len(vec)])]) + zip(vec, _commutative_momenta[:dim])]) factor = e**(I*phase) hopfamily = [] for mat in block_basis: matrix = np.zeros((N, N), dtype=complex) matrix[ranges[a], ranges[b]] = mat - term = Model({factor: matrix}, momenta=range(len(vec))) + term = Model({factor: matrix}, momenta=range(dim)) term = term + term.T().conj() hopfamily.append(term) # If there are conserved quantities, constrain the hopping, it is assumed that diff --git a/qsymm/linalg.py b/qsymm/linalg.py index f19ed27..1ccb52f 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -40,7 +40,7 @@ def allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): if a.size == 0: return a.shape == b.shape atol = atol + rtol * max(np.max(np.abs(a)), np.max(np.abs(b))) - return np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan) + return np.allclose(a, b, rtol=0, atol=atol, equal_nan=equal_nan) def mtm(a, B, c): diff --git a/qsymm/model.py b/qsymm/model.py index 90d731e..3af25da 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -109,35 +109,35 @@ class Model(UserDict): # Make it work with numpy arrays __array_ufunc__ = None - def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], interesting_keys=None): + def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], interesting_keys=None, + symbol_normalizer=None, restructure_dict=False): """ General class to efficiently store any matrix valued function. - The internal structure is a dict with {symbol: value}, where - symbol is a sympy expression, the object representing sum(symbol * value). - The values can be scalars, arrays (both dense and sparse) or LinearOperators. + The Model represents sum(symbol * value), where symbol is a symbolic + expression, and value can be scalar, array (both dense and sparse) + or LinearOperator. The internal structure is a dict {symbol: value}. Implements many sympy and numpy methods and arithmetic operators. Multiplication is distributed over the sum, * is passed down to both symbols and values, @ is passed to symbols as * and to values - as @. Assumes that symbols form a commutative group. - Enhances the functionality of Model by allowing interesting_keys to be - specified, symbols not listed there are discarded. + as @. By default symbols are sympified and assumed commutative. Parameters ---------- hamiltonian : str, SymPy expression or dict - Symbolic representation of a Hamiltonian. It is + 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 - {expression: array} with all arrays the same size (dense or sparse). - expression will be passed through sympy.sympify, and should consist - purely of symbolic coefficients, no constant factors other than 1. + {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. 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]]}. - interesting_keys : iterable of expressions or None (default) + interesting_keys : 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. @@ -146,7 +146,13 @@ class Model(UserDict): or a list of names for the momentum variables as sympy symbols. Momenta are treated the same as other keys for the purpose of interesting_keys, need to list interesting powers of momenta. - + symbol_normalizer : callable (optional) + Function to apply symbols when initializing with dict. By default the + keys are passed through sympy.sympify and sympy.expand_power_exp. + restructure_dict : 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 np.allclose to zero """ self.momenta = _find_momenta(momenta) @@ -156,9 +162,12 @@ class Model(UserDict): self.interesting_keys = set() if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): + if symbol_normalizer is None: + symbol_normalizer = lambda x: sympy.expand_power_exp(sympy.sympify(x)) # Initialize as dict sympifying the keys - super().__init__({sympy.sympify(k): v for k, v in hamiltonian.items() - if sympy.sympify(k) in self.interesting_keys or not self.interesting_keys}) + super().__init__({symbol_normalizer(k): v for k, v in hamiltonian.items() + if symbol_normalizer(k) in self.interesting_keys + or not self.interesting_keys}) else: # Try to parse the input with kwant_continuum.sympify hamiltonian = kwant_continuum.sympify(hamiltonian, locals=locals) @@ -175,48 +184,47 @@ class Model(UserDict): monomials = {k: v for k, v in monomials.items() if not np.allclose(v, 0)} self.data = monomials - - # Restructure - # Clean internal data by: - # * splitting summands in keys - # * moving numerical factors to values - # * removing entries which values care np.allclose to zero - new_data = defaultdict(lambda: list()) - ### TODO: this loop seems quite inefficient. Maybe add an option - # to skip it? - for key, val in self.data.items(): - for summand in key.expand().powsimp(combine='exp').as_ordered_terms(): - factors = summand.as_ordered_factors() - symbols, numbers = [], [] - for f in factors: - # This catches sqrt(2) and much faster than f.is_constant() - if f.is_number: - numbers.append(f) - else: - symbols.append(f) - new_key = sympy.Mul(*symbols) - new_val = complex(sympy.Mul(*numbers)) * val - new_data[new_key].append(new_val) - # translate list to single values - new_data = {k: np.sum(np.array(v), axis=0) - for k, v in new_data.items()} - # remove zero entries - new_data = {k: v for k, v in new_data.items() - if not allclose(v, 0)} - # overwrite internal data - self.data = new_data - - self.shape = _find_shape(self.data) + restructure_dict = True # Keep track of whether this is a dense array self._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + if restructure_dict: + # Clean internal data by: + # * splitting summands in keys + # * moving numerical factors to values + # * removing entries which values care np.allclose to zero + new_data = defaultdict(lambda: list()) + for key, val in self.data.items(): + for summand in key.expand().powsimp(combine='exp').as_ordered_terms(): + factors = summand.as_ordered_factors() + symbols, numbers = [], [] + for f in factors: + # This catches sqrt(2) and much faster than f.is_constant() + if f.is_number: + numbers.append(f) + else: + symbols.append(f) + new_key = sympy.Mul(*symbols) + new_val = complex(sympy.Mul(*numbers)) * val + new_data[new_key].append(new_val) + # translate list to single values + new_data = {k: np.sum(np.array(v), axis=0) + for k, v in new_data.items()} + # remove zero entries + new_data = {k: v for k, v in new_data.items() + if not allclose(v, 0)} + # overwrite internal data + self.data = new_data + + self.shape = _find_shape(self.data) + # Defaultdict functionality def __missing__(self, key): if self.shape is not None: if self.shape == (): #scalar - return 0 + return np.complex128(0) elif self._isarray: # Return dense zero array if dense return np.zeros(self.shape, dtype=complex) @@ -234,11 +242,8 @@ class Model(UserDict): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. # Define addition of 0 and {} - if not other: + if other.data == {} or other == 0 or other == {}: result = self.copy() - # If self is empty return other - elif not self and isinstance(other, type(self)): - result = other.copy() elif isinstance(other, type(self)): if self.momenta != other.momenta: raise ValueError("Can only add Models with the same momenta") @@ -253,6 +258,7 @@ class Model(UserDict): result[key] = copy(self[key]) for key in other.keys() - self.keys(): result[key] = copy(other[key]) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) else: raise NotImplementedError('Addition of {} with {} not supported'.format(type(self), type(other))) return result @@ -292,11 +298,14 @@ class Model(UserDict): if result is 0: result = self.zeros_like() result.shape = (self[1] * other[1]).shape - result.momenta = self.momenta + else: + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result.momenta = self.momenta.copy() else: # Otherwise try to multiply every value with other result.data = {key: val * other for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (self[1] * other).shape + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) return result def __rmul__(self, other): @@ -311,6 +320,7 @@ class Model(UserDict): result = self.zeros_like() result.data = {key: other * val for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (other * self[1]).shape + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) return result def __matmul__(self, other): @@ -326,12 +336,15 @@ class Model(UserDict): if result is 0: result = self.zeros_like() result.shape = (self[1] * other[1]).shape - result.momenta = self.momenta + else: + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result.momenta = self.momenta.copy() else: # Otherwise try to multiply every value with other result = self.zeros_like() result.data = {key: val @ other for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (self[1] @ other).shape + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) return result def __rmatmul__(self, other): @@ -339,6 +352,7 @@ class Model(UserDict): result = self.zeros_like() result.data = {key: other @ val for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (other @ self[1]).shape + result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) return result def __truediv__(self, other): @@ -370,12 +384,12 @@ class Model(UserDict): def transform_symbolic(self, func): """Transform keys by applying func to all of them. Useful for symbolic substitutions, differentiation, etc.""" - if self == {}: - result = self.zeros_like() - else: - # Add possible duplicate keys that only differ in constant factors - result = sum(type(self)({func(key): val}, momenta=self.momenta) - for key, val in self.items()) + # Add possible duplicate keys that only differ in constant factors + result = sum((type(self)({func(key): copy(val)}, + restructure_dict=True, + momenta=self.momenta.copy()) + for key, val in self.items()), + self.zeros_like()) return result def rotate_momenta(self, R): @@ -412,7 +426,7 @@ class Model(UserDict): elif isinstance(args[0], dict): # Input is a dictionary args = ([(key, value) for key, value in args[0].items()], ) - momenta = self.momenta + momenta = self.momenta.copy() for (old, new) in args[0]: # Substitution of a momentum variable with a symbol # is a renaming of the momentum. @@ -431,7 +445,7 @@ class Model(UserDict): substituted.momenta = momenta # If there are exponentials, evaluate any numerical exponents, # so they can be moved to the matrix valued part of the Model - result = type(substituted)({}, momenta=momenta) + result = substituted.zeros_like() for key, value in substituted.items(): # Expand sums in the exponent to products of exponentials, # find all exponentials. @@ -444,9 +458,9 @@ class Model(UserDict): # Otherwise, leave the exponential unchanged. expos = [expo.subs(e, np.e).evalf() if expo.subs(e, np.e).evalf().is_number else expo for expo in find_expos] - result += type(substituted)({rest * np.prod(expos): value}, momenta=momenta) + result += type(substituted)({rest * np.prod(expos): value}, momenta=momenta, restructure_dict=True) else: - result += type(substituted)({key: value}, momenta=momenta) + result += type(substituted)({key: value}, momenta=momenta, restructure_dict=True) return result def conj(self): @@ -629,7 +643,7 @@ class BlochModel(Model): self.shape = _find_shape(self.data) self.momenta = _find_momenta(momenta) elif symbolic: - self.__init__(Model(hamiltonian, locals, momenta)) + self.__init__(Model(hamiltonian, locals, momenta=momenta)) else: # Use Model to parse input self.__init__(Model(hamiltonian, locals, momenta)) @@ -789,10 +803,13 @@ def _find_shape(data): val = next(iter(data.values())) if isinstance(val, Number): shape = () - if not all([isinstance(v, Number) for v in data.values()]): + if not all(isinstance(v, Number) for v in data.values()): raise ValueError('All terms must have the same shape') + # Recast numbers to numpy complex128 + for key, val in data.items(): + data[key] = np.complex128(val) else: shape = val.shape - if not all([v.shape == shape for v in data.values()]): + if not all(v.shape == shape for v in data.values()): raise ValueError('All terms must have the same shape') return shape diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index 545b6db..b8a51c3 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -100,7 +100,7 @@ def test_Model(): assert allclose(m[keys[1]], 0.5) mat = np.random.rand(2,2) - m = Model({np.sqrt(5)*e**(I*k_x) : mat}) + m = Model({np.sqrt(5)*e**(I*k_x) : mat}, restructure_dict=True) assert allclose(m[list(m.keys())[0]], mat*np.sqrt(5)) m1 = Model({k_x : 2.3}) @@ -138,10 +138,11 @@ def test_Model(): def test_BlochModel(): - m = Model({e**(I*k_y): 3*np.eye(2), k_x : np.eye(2)}) + m = Model({e**(I*k_y): 3*np.eye(2), k_x : np.eye(2)}, restructure_dict=True) with raises(AssertionError, message="All momentum dependence should be in the hopping."): bm = BlochModel(m) - m = Model({c1*e**(I*k_y): 3*np.eye(2), np.sqrt(3)*e**(I*k_x) : np.eye(2)}, momenta=[0, 1]) + m = Model({c1*e**(I*k_y): 3*np.eye(2), np.sqrt(3)*e**(I*k_x) : np.eye(2)}, + momenta=[0, 1], restructure_dict=True) bm = BlochModel(m) keys = [BlochCoeff(np.array([0, 1]), c1), BlochCoeff(np.array([1, 0]), sympy.numbers.One())] assert all([key in keys for key in bm.keys()]) -- GitLab From b8cea826424d179840c4f5fde00e0625f7420a68 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sat, 18 May 2019 13:17:24 +0200 Subject: [PATCH 06/41] one more optimization --- qsymm/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsymm/model.py b/qsymm/model.py index 3af25da..76b54e7 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -242,7 +242,7 @@ class Model(UserDict): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. # Define addition of 0 and {} - if other.data == {} or other == 0 or other == {}: + if other.data == {} or other is 0 or other is {}: result = self.copy() elif isinstance(other, type(self)): if self.momenta != other.momenta: -- GitLab From f5497e22032c9fbe7c084a2ff4aaf538a1bc5103 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sat, 18 May 2019 17:17:55 +0200 Subject: [PATCH 07/41] clean up BlochModel --- qsymm/model.py | 78 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 52 insertions(+), 26 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index 76b54e7..f60474f 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -242,9 +242,11 @@ class Model(UserDict): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. # Define addition of 0 and {} - if other.data == {} or other is 0 or other is {}: + if (not isinstance(other, type(self)) and (other == 0 or other == {}) + or (isinstance(other, type(self)) and other.data=={})): result = self.copy() elif isinstance(other, type(self)): + # 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") result = self.zeros_like() @@ -599,19 +601,29 @@ class Model(UserDict): class BlochModel(Model): - def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2]): - """Class to store Bloch Hamiltonian families. - Can be used to efficiently store any matrix valued function. - Implements many sympy and numpy methods. Arithmetic operators are overloaded, - such that * corresponds to matrix multiplication. + def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], + interesting_keys=None): + """ + Class to efficiently store matrix valued Bloch Hamiltonians. + The BlochModel represents sum(BlochCoeff * value), where BlochCoeff + is a symbolic representation of coefficient and periodic functions. + value can be scalar, array (both dense and sparse) + or LinearOperator. The internal structure is a dict {BlochCoeff: value}. + Implements many sympy and numpy methods and arithmetic operators. + Multiplication is distributed over the sum, * is passed down to + both symbols and values, @ is passed to symbols as * and to values + as @. By default symbols are sympified and assumed commutative. Parameters ---------- - hamiltonian : str, SymPy expression or dict - Symbolic representation of a Hamiltonian. It is + hamiltonian : Model, str, SymPy expression or dict + 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 - {BlochCoeff: np.ndarray} with all arrays the same square size. + {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. 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 @@ -620,30 +632,47 @@ class BlochModel(Model): locals={'sigma_plus': [[0, 2], [0, 0]]}. momenta : list of int or list of Sympy objects Indices of momenta the monomials depend on from 'k_x', 'k_y' and 'k_z' - or a list of names for the momentum variables. + or a list of names for the momentum variables. Ignored when + initialized with Model. + interesting_keys : 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. """ if isinstance(hamiltonian, Model): - # Recast keys into BlochCoeffs - # This works if some keys are different but close, such that BlochCoeff - # is the same - self.__init__({}, momenta=hamiltonian.momenta) - data = defaultdict(lambda: np.zeros(hamiltonian.shape, dtype=complex)) + # First initialize an empty BlochModel, this is the same as init for Model + self.__init__(hamiltonian={}, + locals=locals, + momenta=hamiltonian.momenta, + interesting_keys=hamiltonian.interesting_keys) + # Initialize same as input model, so __missing__ works + self._isarray = hamiltonian._isarray + self.shape = hamiltonian.shape + # Recast keys into BlochCoeffs, if some keys are different but close, + # such that BlochCoeff is the same, collect them. for key, val in hamiltonian.items(): - data[_to_bloch_coeff(key, hamiltonian.momenta)] += val - self.data = data + self[_to_bloch_coeff(key, hamiltonian.momenta)] += val elif isinstance(hamiltonian, abc.Mapping): keys = hamiltonian.keys() - symbolic = all(isinstance(k, Basic) for k in keys) + symbolic = all(not isinstance(k, BlochCoeff) for k in keys) hopping = all(isinstance(k, BlochCoeff) for k in keys) if not (symbolic or hopping): raise ValueError('All keys must have the same type (sympy expression or BlochCoeff).') if hopping or hamiltonian == {}: - # initialize as dict - super(Model, self).__init__(hamiltonian) - self.shape = _find_shape(self.data) - self.momenta = _find_momenta(momenta) + # initialize as Model without any of the preprocessing + super().__init__(hamiltonian, + locals=locals, + momenta=momenta, + interesting_keys=interesting_keys, + symbol_normalizer=lambda x: x, + restructure_dict=False) elif symbolic: - self.__init__(Model(hamiltonian, locals, momenta=momenta)) + # First cast it to model, then try to interpret it as BlochModel + self.__init__(Model(hamiltonian, + locals=locals, + momenta=momenta, + interesting_keys=interesting_keys)) else: # Use Model to parse input self.__init__(Model(hamiltonian, locals, momenta)) @@ -651,9 +680,6 @@ class BlochModel(Model): self.shape = _find_shape(self.data) self.interesting_keys = set() - # Keep track of whether this is a dense array - self._isarray = any(isinstance(val, np.ndarray) for val in self.values()) - def transform_symbolic(self, func): raise NotImplementedError('transform_symbolic is not implemented for BlochModel') -- GitLab From 4a21ed396a19c23156daedf54d7630dd58e2dd51 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sat, 18 May 2019 17:39:18 +0200 Subject: [PATCH 08/41] more cleaning of BlochModel --- qsymm/model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index f60474f..f68d3bd 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -675,10 +675,10 @@ class BlochModel(Model): interesting_keys=interesting_keys)) else: # Use Model to parse input - self.__init__(Model(hamiltonian, locals, momenta)) - - self.shape = _find_shape(self.data) - self.interesting_keys = set() + self.__init__(Model(hamiltonian, + locals=locals, + momenta=momenta, + interesting_keys=interesting_keys)) def transform_symbolic(self, func): raise NotImplementedError('transform_symbolic is not implemented for BlochModel') -- GitLab From 0f552c4f8e3afd49bcfd3710f3135c98c71a67df Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sun, 19 May 2019 22:51:28 +0200 Subject: [PATCH 09/41] add more tests and fix things about Model --- qsymm/model.py | 12 +++++----- qsymm/tests/test_model.py | 47 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+), 6 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index f68d3bd..1e222cf 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -301,13 +301,13 @@ class Model(UserDict): result = self.zeros_like() result.shape = (self[1] * other[1]).shape else: - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) result.momenta = self.momenta.copy() else: # Otherwise try to multiply every value with other result.data = {key: val * other for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (self[1] * other).shape - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) return result def __rmul__(self, other): @@ -322,7 +322,7 @@ class Model(UserDict): result = self.zeros_like() result.data = {key: other * val for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (other * self[1]).shape - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) return result def __matmul__(self, other): @@ -339,14 +339,14 @@ class Model(UserDict): result = self.zeros_like() result.shape = (self[1] * other[1]).shape else: - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) result.momenta = self.momenta.copy() else: # Otherwise try to multiply every value with other result = self.zeros_like() result.data = {key: val @ other for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (self[1] @ other).shape - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) return result def __rmatmul__(self, other): @@ -354,7 +354,7 @@ class Model(UserDict): result = self.zeros_like() result.data = {key: other @ val for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (other @ self[1]).shape - result._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) return result def __truediv__(self, other): diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index b8a51c3..bcad131 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -1,6 +1,8 @@ import pytest import warnings import sympy +from scipy.sparse import csr_matrix +from scipy.sparse.linalg import LinearOperator import numpy as np from pytest import raises @@ -16,6 +18,8 @@ c0, c1 = sympy.symbols('c0 c1', real=True) def test_dense_algebra(): # scalar models a = Model({1: 1, 'x': 3}) + assert a.shape == () + assert a._isarray == False assert a * 0 == 0 * a assert 0 * a == a.zeros_like() assert 0 * a == {} @@ -34,11 +38,18 @@ def test_dense_algebra(): # vector models vec = np.ones((3,)) c = a * vec + assert c.shape == (3,) + assert c._isarray == True + # matrix models mat = np.ones((3, 3)) d = b * mat + assert d.shape == (3, 3) + assert d._isarray == True assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) assert c.T() @ d == 3 * a * b assert c.T() @ d @ c == 9 * a * a * b + assert d @ d == 3 * b * d + assert (d @ d)._isarray == True # numpy broadcasting rules apply assert d + c == (a + b) * mat assert d * c == a * b * mat @@ -46,6 +57,42 @@ def test_dense_algebra(): assert d.reshape(-1) == b * np.ones((9,)) +def test_sparse_algebra(): + # Test sparse matrices + a = Model({1: 1, 'x': 3}) + b = Model({1: 2, 'x**2': 3}) + # sparse vector model + vec = csr_matrix(np.ones((3, 1))) + c = a * vec + assert c.shape == (3, 1) + assert c._isarray == False + mat = csr_matrix(np.ones((3, 3))) + d = b * mat + assert d.shape == (3, 3) + assert d._isarray == False + assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) + assert (d @ d) @ c == 9 * b * b * c + assert (d @ d)._isarray == False + assert c.T() @ d == 3 * a * b + assert c.T() @ d @ c == 9 * a * a * b + assert d.trace() == 3 * b + assert d.reshape((1, 9)) @ np.ones((9,)) == 9 * b + e = d @ np.ones((3,)) + assert e == 3 * b * np.ones((3,)) + assert e._isarray == True + + # Test LinearOperator + d = b * LinearOperator((3, 3), matvec=lambda v: mat @ v) + c = a * np.ones((3, 1)) + assert d.shape == (3, 3) + assert d._isarray == False + assert (d @ d) @ c == 9 * b * b * c + assert (d @ d)._isarray == False + assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) + assert (d @ c)._isarray == True + assert c.T() @ (d @ c) == 9 * a * a * b + + def test_to_bloch_coeff(): key = sympy.sqrt(5)*e**(2*c0)*e**(I*(k_x/2 + np.sqrt(3)*k_y + c1)) -- GitLab From 99b736acff0098091ef82dc15a451718a58ce880 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 20 May 2019 11:56:12 +0200 Subject: [PATCH 10/41] avoid mutable defaults --- qsymm/model.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index 1e222cf..d65f58a 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -109,7 +109,7 @@ class Model(UserDict): # Make it work with numpy arrays __array_ufunc__ = None - def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], interesting_keys=None, + def __init__(self, hamiltonian=None, locals=None, momenta=(0, 1, 2), interesting_keys=None, symbol_normalizer=None, restructure_dict=False): """ General class to efficiently store any matrix valued function. @@ -123,14 +123,14 @@ class Model(UserDict): Parameters ---------- - hamiltonian : str, SymPy expression or dict + hamiltonian : 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). symbol by default is passed through sympy.sympify, and should consist purely of a product of symbolic coefficients, no constant - factors other than 1. + factors other than 1. 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 @@ -141,7 +141,7 @@ class Model(UserDict): 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 : list of int or list of Sympy objects + momenta : iterable of int or list of Sympy objects Indices of momentum variables from ['k_x', 'k_y', 'k_z'] or a list of names for the momentum variables as sympy symbols. Momenta are treated the same as other keys for the purpose of @@ -161,6 +161,9 @@ class Model(UserDict): else: self.interesting_keys = set() + if hamiltonian is None: + hamiltonian = {} + if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): if symbol_normalizer is None: symbol_normalizer = lambda x: sympy.expand_power_exp(sympy.sympify(x)) @@ -601,7 +604,7 @@ class Model(UserDict): class BlochModel(Model): - def __init__(self, hamiltonian={}, locals=None, momenta=[0, 1, 2], + def __init__(self, hamiltonian=None, locals=None, momenta=(0, 1, 2), interesting_keys=None): """ Class to efficiently store matrix valued Bloch Hamiltonians. @@ -616,7 +619,7 @@ class BlochModel(Model): Parameters ---------- - hamiltonian : Model, str, SymPy expression or dict + 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 @@ -624,13 +627,14 @@ class BlochModel(Model): 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 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]]}. - momenta : list of int or list of Sympy objects + momenta : iterable of int or list of Sympy objects Indices of momenta the monomials depend on from 'k_x', 'k_y' and 'k_z' or a list of names for the momentum variables. Ignored when initialized with Model. @@ -640,6 +644,8 @@ class BlochModel(Model): only terms to a given order are needed. By default all keys are kept. Ignored when initialized with Model. """ + if hamiltonian is None: + hamiltonian = {} if isinstance(hamiltonian, Model): # First initialize an empty BlochModel, this is the same as init for Model self.__init__(hamiltonian={}, -- GitLab From 6ea3923ce0b1a207b03e838f9e12acc8f0bcd8bf Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 20 May 2019 11:57:48 +0200 Subject: [PATCH 11/41] fix docstring --- qsymm/groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index d849f8f..11e8621 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -327,7 +327,7 @@ class ContinuousGroupGenerator: return 'ContinuousGroupGenerator(\n{},\n{},\n)'.format(self.R, self.U) def apply(self, model): - """Return copy of monomials with applied infinitesimal generator. + """Return copy of model H(k) with applied infinitesimal generator. 1j * (H(k) U - U H(k)) + 1j * dH(k)/dk_i R_{ij} k_j """ R, U = self.R, self.U -- GitLab From 349454ed57dd572dfae2dc1bd0d0d342d51ca37b Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 20 May 2019 13:04:13 +0200 Subject: [PATCH 12/41] remove PerturbativeModel from docstring --- qsymm/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index d65f58a..2b44bb2 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -329,7 +329,7 @@ class Model(UserDict): return result def __matmul__(self, other): - # Multiplication by arrays and PerturbativeModel + # Multiplication by arrays and Model if isinstance(other, Model): if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") @@ -627,7 +627,7 @@ class BlochModel(Model): 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 Model. + 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 -- GitLab From de49ff675281a5513ae7bb6cc147b0c586481598 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 20 May 2019 14:10:51 +0200 Subject: [PATCH 13/41] better support for sparse and linear operator in allclose --- qsymm/linalg.py | 22 +++++++++++++++++++--- qsymm/model.py | 28 ++++++++++++++++++---------- qsymm/tests/test_model.py | 6 ++++-- 3 files changed, 41 insertions(+), 15 deletions(-) diff --git a/qsymm/linalg.py b/qsymm/linalg.py index 1ccb52f..767f4c3 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -3,6 +3,7 @@ import numpy as np import scipy.linalg as la import scipy.sparse.linalg as sla import scipy +from numbers import Number from scipy.optimize import minimize from scipy.spatial.distance import cdist from scipy.sparse.csgraph import connected_components @@ -34,13 +35,28 @@ def allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): are elementwise close. Unlike the numpy version, the relative tolerance is not elementwise, but it is relative to the largest entry in the arrays.""" - a = np.asarray(a) - b = np.asarray(b) + # If either is a list or number, recast it to ndarray, if it is + # a LinearOperator, recast it to dense matrix. This is rather + # inefficient, but LinearOperator doesn't support multiplication + # with sparse matrices. + if isinstance(a, (list, Number)): + a = np.asarray(a) + elif isinstance(a, sla.LinearOperator): + a = a @ np.eye(a.shape[-1], dtype=complex) + if isinstance(b, (list, Number)): + b = np.asarray(b) + elif isinstance(b, sla.LinearOperator): + b = b @ np.eye(a.shape[-1], dtype=complex) + # Check if empty arrays, only compare shape if a.size == 0: return a.shape == b.shape atol = atol + rtol * max(np.max(np.abs(a)), np.max(np.abs(b))) - return np.allclose(a, b, rtol=0, atol=atol, equal_nan=equal_nan) + diff = a - b + if isinstance(diff, scipy.sparse.spmatrix): + diff = diff.data + + return np.allclose(diff, 0, rtol=0, atol=atol, equal_nan=equal_nan) def mtm(a, B, c): diff --git a/qsymm/model.py b/qsymm/model.py index 2b44bb2..1d3a62f 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -524,16 +524,26 @@ class Model(UserDict): def tocsr(self): result = self.zeros_like() - result.data = {key: scipy.sparse.csr_matrix(val, dtype=complex) - for key, val in self.items()} + for key, val in self.items(): + if isinstance(val, (Number, np.ndarray, scipy.sparse.spmatrix)): + result[key] = scipy.sparse.csr_matrix(val, dtype=complex) + else: + # LinearOperator doesn't support multiplication with sparse matrix + val = scipy.sparse.csr_matrix(val @ np.eye(val.shape[-1], dtype=complex), dtype=complex) result._isarray = False return result def toarray(self): result = self.zeros_like() - result.data = {key : (val if isinstance(val, np.ndarray) or isinstance(val, Number) - else val.toarray()) - for key, val in self.items()} + for key, val in self.items(): + if isinstance(val, np.ndarray): + result[key] = val + elif isinstance(val, Number): + result[key] = np.asarray(val) + elif scipy.sparse.spmatrix: + result[key] = val.A + else: + val = val @ np.eye(val.shape[-1], dtype=complex) result._isarray = True return result @@ -591,16 +601,14 @@ class Model(UserDict): def allclose(self, other, rtol=1e-05, atol=1e-08, equal_nan=False): # Test whether two arrays are approximately equal - a = self.toarray() if other == {} or other == 0: if self.data == {}: return True else: - return all(allclose(val, 0, rtol, atol, equal_nan) for val in a.values()) + return all(allclose(val, 0, rtol, atol, equal_nan) for val in self.values()) else: - b = other.toarray() - return all(allclose(a[key], b[key], rtol, atol, equal_nan) - for key in a.keys() | b.keys()) + return all(allclose(self[key], other[key], rtol, atol, equal_nan) + for key in self.keys() | other.keys()) class BlochModel(Model): diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index bcad131..695a0cf 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -71,10 +71,11 @@ def test_sparse_algebra(): assert d.shape == (3, 3) assert d._isarray == False assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) + assert d @ d == 3 * b * d assert (d @ d) @ c == 9 * b * b * c assert (d @ d)._isarray == False - assert c.T() @ d == 3 * a * b - assert c.T() @ d @ c == 9 * a * a * b + assert c.T() @ d == 3 * b * c.T() + assert c.T() @ d @ c == 9 * a * a * b * np.eye(1) assert d.trace() == 3 * b assert d.reshape((1, 9)) @ np.ones((9,)) == 9 * b e = d @ np.ones((3,)) @@ -86,6 +87,7 @@ def test_sparse_algebra(): c = a * np.ones((3, 1)) assert d.shape == (3, 3) assert d._isarray == False + assert d @ d == 3 * b * d assert (d @ d) @ c == 9 * b * b * c assert (d @ d)._isarray == False assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) -- GitLab From ea1c12e978e3cc39aa283907628aa7de23fedfb3 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Mon, 20 May 2019 14:56:43 +0200 Subject: [PATCH 14/41] clear up isinstance and imports --- qsymm/hamiltonian_generator.py | 4 +--- qsymm/model.py | 6 +++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py index 567bde3..9f803bf 100644 --- a/qsymm/hamiltonian_generator.py +++ b/qsymm/hamiltonian_generator.py @@ -5,11 +5,9 @@ import itertools as it import scipy.linalg as la from copy import deepcopy -from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref, allclose from .model import Model, BlochModel, _commutative_momenta, e, I -from .groups import PointGroupElement, ContinuousGroupGenerator -from .groups import generate_group +from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group from . import kwant_continuum diff --git a/qsymm/model.py b/qsymm/model.py index 1d3a62f..e32e0d0 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -513,9 +513,9 @@ class Model(UserDict): vnsimplify = np.vectorize(sympy.nsimplify, otypes=[object]) result = sympy.MatAdd(*[key * sympy.Matrix(vnsimplify(val)) for key, val in self.toarray().items()]).doit() - if any([isinstance(result, matrix_type) for matrix_type in (sympy.MatrixBase, - sympy.ImmutableDenseMatrix, - sympy.ImmutableDenseNDimArray)]): + if isinstance(result, (sympy.MatrixBase, + sympy.ImmutableDenseMatrix, + sympy.ImmutableDenseNDimArray)): result = sympy.Matrix(result).reshape(*result.shape) return result -- GitLab From cf55966f4e6aebf7b3966305bf69fb526dc896bc Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 20:12:29 +0200 Subject: [PATCH 15/41] fix some inconsistencies in Model --- qsymm/model.py | 49 ++++++++++++++++++++++++++----------------------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/qsymm/model.py b/qsymm/model.py index e32e0d0..7228aa2 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -186,7 +186,7 @@ class Model(UserDict): # remove matrices == zeros monomials = {k: v for k, v in monomials.items() if not np.allclose(v, 0)} - self.data = monomials + super().__init__(monomials) restructure_dict = True # Keep track of whether this is a dense array @@ -286,12 +286,14 @@ class Model(UserDict): def __mul__(self, other): # Multiplication by numbers, sympy symbols, arrays and Model - result = self.zeros_like() if isinstance(other, Number): + result = self.zeros_like() result.data = {key: val * other for key, val in self.items()} elif isinstance(other, Basic): - result.data = {key * other: val for key, val in self.items() - if (key * other in interesting_keys or not interesting_keys)} + result = sum((type(self)({key * other: copy(val)}, interesting_keys=interesting_keys) + for key, val in self.items() + if (key * other in interesting_keys or not interesting_keys)), + self.zeros_like()) elif isinstance(other, Model): if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") @@ -308,6 +310,7 @@ class Model(UserDict): result.momenta = self.momenta.copy() else: # Otherwise try to multiply every value with other + result = self.zeros_like() result.data = {key: val * other for key, val in self.items()} result.shape = _find_shape(result.data) if result.data else (self[1] * other).shape result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) @@ -318,8 +321,10 @@ class Model(UserDict): if isinstance(other, Number): result = self.__mul__(other) elif isinstance(other, Basic): - result = self.zeros_like() - result.data = {other * key: copy(val) for key, val in self.items()} + result = sum((type(self)({other * key: copy(val)}, interesting_keys=interesting_keys) + for key, val in self.items() + if (key * other in interesting_keys or not interesting_keys)), + self.zeros_like()) else: # Otherwise try to multiply every value with other result = self.zeros_like() @@ -340,7 +345,7 @@ class Model(UserDict): # Find out the shape of the result even if it is empty if result is 0: result = self.zeros_like() - result.shape = (self[1] * other[1]).shape + result.shape = (self[1] @ other[1]).shape else: result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) result.momenta = self.momenta.copy() @@ -655,24 +660,19 @@ class BlochModel(Model): if hamiltonian is None: hamiltonian = {} if isinstance(hamiltonian, Model): - # First initialize an empty BlochModel, this is the same as init for Model - self.__init__(hamiltonian={}, - locals=locals, - momenta=hamiltonian.momenta, - interesting_keys=hamiltonian.interesting_keys) - # Initialize same as input model, so __missing__ works + # Use Model's init, only need to recast keys to BlochCoeff + super().__init__(hamiltonian=hamiltonian.data, + locals=locals, + momenta=hamiltonian.momenta, + interesting_keys=hamiltonian.interesting_keys, + symbol_normalizer=lambda key: _to_bloch_coeff(key, hamiltonian.momenta)) + # set these in case it was and empty Model self._isarray = hamiltonian._isarray self.shape = hamiltonian.shape - # Recast keys into BlochCoeffs, if some keys are different but close, - # such that BlochCoeff is the same, collect them. - for key, val in hamiltonian.items(): - self[_to_bloch_coeff(key, hamiltonian.momenta)] += val elif isinstance(hamiltonian, abc.Mapping): keys = hamiltonian.keys() symbolic = all(not isinstance(k, BlochCoeff) for k in keys) hopping = all(isinstance(k, BlochCoeff) for k in keys) - if not (symbolic or hopping): - raise ValueError('All keys must have the same type (sympy expression or BlochCoeff).') if hopping or hamiltonian == {}: # initialize as Model without any of the preprocessing super().__init__(hamiltonian, @@ -682,11 +682,14 @@ class BlochModel(Model): symbol_normalizer=lambda x: x, restructure_dict=False) elif symbolic: - # First cast it to model, then try to interpret it as BlochModel + # First cast it to model with restructuring, then try to interpret it as BlochModel self.__init__(Model(hamiltonian, locals=locals, momenta=momenta, - interesting_keys=interesting_keys)) + interesting_keys=interesting_keys, + restructure_dict=True)) + else: + raise ValueError('All keys must have the same type (sympy expression or BlochCoeff).') else: # Use Model to parse input self.__init__(Model(hamiltonian, @@ -703,7 +706,7 @@ class BlochModel(Model): assert len(momenta) == R.shape[0], (momenta, R) # do rotation on hopping vectors with transpose matrix R_T = np.array(R).astype(float).T - return BlochModel({BlochCoeff(R_T @ hop, coeff): val + return BlochModel({BlochCoeff(R_T @ hop, coeff): copy(val) for (hop, coeff), val in self.items()}, momenta=momenta) def conj(self): @@ -724,7 +727,7 @@ class BlochModel(Model): return self.tomodel(nsimplify=nsimplify).tosympy(nsimplify) def tomodel(self, nsimplify=False): - return Model({key.tosympy(self.momenta, nsimplify=nsimplify): val + return Model({key.tosympy(self.momenta, nsimplify=nsimplify): copy(val) for key, val in self.items()}, momenta=self.momenta) -- GitLab From 677bbe04004f2a71cb86cb15c3b73d7da984bf3b Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Wed, 22 May 2019 14:54:34 +0200 Subject: [PATCH 16/41] keep track of data types more rigorously in Model --- qsymm/groups.py | 1 - qsymm/model.py | 241 ++++++++++++++++++++++--------------- qsymm/tests/test_model.py | 33 +++-- qsymm/tests/test_mutual.py | 32 ++--- 4 files changed, 175 insertions(+), 132 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 11e8621..96bb4c7 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -343,7 +343,6 @@ class ContinuousGroupGenerator: return sum([sympy.diff(key, momenta[i]) * R[i, j] * momenta[j] for i, j in it.product(range(dim), repeat=2)]) result += 1j * model.transform_symbolic(trf) - assert result._isarray == model._isarray if U_nonzero: result += model @ (1j*U) + (-1j*U) @ model return result diff --git a/qsymm/model.py b/qsymm/model.py index 7228aa2..54fb3ec 100644 --- a/qsymm/model.py +++ b/qsymm/model.py @@ -66,6 +66,8 @@ class BlochCoeff(tuple): def __eq__(self, other): hop1, coeff1 = self + if other == 1: + return allclose(hop1, 0) and coeff1 == 1 hop2, coeff2 = other # test equality of hop with allclose return allclose(hop1, hop2) and coeff1 == coeff2 @@ -110,7 +112,7 @@ class Model(UserDict): __array_ufunc__ = None def __init__(self, hamiltonian=None, locals=None, momenta=(0, 1, 2), interesting_keys=None, - symbol_normalizer=None, restructure_dict=False): + symbol_normalizer=None, restructure_dict=False, shape=None, dtype=None): """ General class to efficiently store any matrix valued function. The Model represents sum(symbol * value), where symbol is a symbolic @@ -153,9 +155,21 @@ class Model(UserDict): Whether to clean input dict by splitting summands in symbols, moving numerical factors in the symbols to values, removing entries with values np.allclose to zero + 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. + dtype : class or None (default) + Type of the values in the model. Supported types are np.complex128, + and subclasses of 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 dtype is not provided, it is inferred + from the type of the values. Must be provided if hamiltonian is None + or {}. If hamiltonian is any other format, dtype is ignored an set to + np.ndarray. """ - self.momenta = _find_momenta(momenta) - if interesting_keys is not None: self.interesting_keys = {sympy.sympify(k) for k in interesting_keys} else: @@ -163,14 +177,18 @@ class Model(UserDict): if hamiltonian is None: hamiltonian = {} + if hamiltonian == {} and (shape is None or dtype is None): + raise ValueError('Must provide shape and dtype when initializing empty Model.') + if symbol_normalizer is None: + symbol_normalizer = lambda x: sympy.expand_power_exp(sympy.sympify(x)) + self.momenta = _find_momenta(momenta) if hamiltonian == {} or isinstance(hamiltonian, abc.Mapping): - if symbol_normalizer is None: - symbol_normalizer = lambda x: sympy.expand_power_exp(sympy.sympify(x)) # Initialize as dict sympifying the keys super().__init__({symbol_normalizer(k): v for k, v in hamiltonian.items() if symbol_normalizer(k) in self.interesting_keys or not self.interesting_keys}) + else: # Try to parse the input with kwant_continuum.sympify hamiltonian = kwant_continuum.sympify(hamiltonian, locals=locals) @@ -189,16 +207,30 @@ class Model(UserDict): super().__init__(monomials) restructure_dict = True - # Keep track of whether this is a dense array - self._isarray = any(isinstance(val, np.ndarray) for val in self.values()) + # Find shape and dtype + self.shape = shape + self.dtype = dtype + if self.shape is None or self.dtype is None: + val = next(iter(self.values())) + shape, dtype = _shape_and_dtype(val) + self.shape = (shape if self.shape is not None else shape) + self.dtype = (dtype if self.dtype is not None else dtype) + if shape == (): + # Recast numbers to np.complex128 + self.data = {k: np.complex128(v) for k, v in self.items()} + if not all(issubclass(type(v), dtype) for v in self.values()): + raise ValueError('All values must have the same dtype.') + if not all(v.shape == shape for v in self.values()): + raise ValueError('All values must have the same shape.') if restructure_dict: # Clean internal data by: # * splitting summands in keys # * moving numerical factors to values # * removing entries which values care np.allclose to zero - new_data = defaultdict(lambda: list()) - for key, val in self.data.items(): + old_data = {copy(key): copy(val) for key, val in self.items()} + self.data = {} + for key, val in old_data.items(): for summand in key.expand().powsimp(combine='exp').as_ordered_terms(): factors = summand.as_ordered_factors() symbols, numbers = [], [] @@ -210,32 +242,24 @@ class Model(UserDict): symbols.append(f) new_key = sympy.Mul(*symbols) new_val = complex(sympy.Mul(*numbers)) * val - new_data[new_key].append(new_val) - # translate list to single values - new_data = {k: np.sum(np.array(v), axis=0) - for k, v in new_data.items()} + self[new_key] += new_val # remove zero entries - new_data = {k: v for k, v in new_data.items() - if not allclose(v, 0)} - # overwrite internal data - self.data = new_data - - self.shape = _find_shape(self.data) + self.data = {k: v for k, v in self.items() if not allclose(v, 0)} # Defaultdict functionality def __missing__(self, key): - if self.shape is not None: - if self.shape == (): - #scalar - return np.complex128(0) - elif self._isarray: - # Return dense zero array if dense - return np.zeros(self.shape, dtype=complex) - else: - # Otherwise return a csr_matrix - return scipy.sparse.csr_matrix(self.shape, dtype=complex) - else: - return None + if self.dtype is np.complex128: + #scalar + return np.complex128(0) + elif self.dtype is np.ndarray: + # Return dense zero array if dense + return np.zeros(self.shape, dtype=complex) + elif issubclass(self.dtype, scipy.sparse.spmatrix): + # Return a zero sparse matrix of the same type + return self.dtype(self.shape, dtype=complex) + elif issubclass(self.dtype, scipy.sparse.linalg.LinearOperator): + return scipy.sparse.linalg.aslinearoperator( + scipy.sparse.csr_matrix(self.shape, dtype=complex)) def __eq__(self, other): # Call allclose with default tolerances @@ -244,26 +268,23 @@ class Model(UserDict): def __add__(self, other): # Addition of Models. It is assumed that both Models are # structured correctly, every key is in standard form. - # Define addition of 0 and {} - if (not isinstance(other, type(self)) and (other == 0 or other == {}) + # Define addition of 0 + if (not isinstance(other, type(self)) and (other == 0) or (isinstance(other, type(self)) and other.data=={})): result = self.copy() elif isinstance(other, type(self)): + if not (self.dtype is other.dtype 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") result = self.zeros_like() for key in self.keys() & other.keys(): - total = self[key] + other[key] - # If only one is sparse matrix, the result is np.matrix, recast it to np.ndarray - if isinstance(total, np.matrix): - total = total.A - result[key] = total + result[key] = self[key] + other[key] for key in self.keys() - other.keys(): result[key] = copy(self[key]) for key in other.keys() - self.keys(): result[key] = copy(other[key]) - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) else: raise NotImplementedError('Addition of {} with {} not supported'.format(type(self), type(other))) return result @@ -295,6 +316,11 @@ class Model(UserDict): if (key * other in interesting_keys or not interesting_keys)), self.zeros_like()) elif isinstance(other, Model): + if not (issubclass(self.dtype, (Number, np.ndarray)) or + issubclass(other.dtype, (Number, np.ndarray))): + raise ValueError('Elementwise multiplication only allowed for scalar ' + 'and ndarra data types. With sprse matrices use @ ' + 'for matrix multiplication.') if self.momenta != other.momenta: raise ValueError("Can only multiply Models with the same momenta") interesting_keys = self.interesting_keys | other.interesting_keys @@ -304,16 +330,15 @@ class Model(UserDict): # Find out the shape of the result even if it is empty if result is 0: result = self.zeros_like() - result.shape = (self[1] * other[1]).shape + result.shape, result.dtype = _shape_and_dtype(self[1] * other[1]) else: - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) - result.momenta = self.momenta.copy() + result.shape, result.dtype = _shape_and_dtype(next(iter(result.values()))) + result.momenta = self.momenta else: # Otherwise try to multiply every value with other result = self.zeros_like() result.data = {key: val * other for key, val in self.items()} - result.shape = _find_shape(result.data) if result.data else (self[1] * other).shape - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) + result.shape, result.dtype = _shape_and_dtype(self[1] * other) return result def __rmul__(self, other): @@ -329,8 +354,7 @@ class Model(UserDict): # Otherwise try to multiply every value with other result = self.zeros_like() result.data = {key: other * val for key, val in self.items()} - result.shape = _find_shape(result.data) if result.data else (other * self[1]).shape - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) + result.shape, result.dtype = _shape_and_dtype(other * self[1]) return result def __matmul__(self, other): @@ -345,24 +369,20 @@ class Model(UserDict): # Find out the shape of the result even if it is empty if result is 0: result = self.zeros_like() - result.shape = (self[1] @ other[1]).shape - else: - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) - result.momenta = self.momenta.copy() + result.shape, result.dtype = _shape_and_dtype(self[1] @ other[1]) + result.momenta = self.momenta else: # Otherwise try to multiply every value with other result = self.zeros_like() result.data = {key: val @ other for key, val in self.items()} - result.shape = _find_shape(result.data) if result.data else (self[1] @ other).shape - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) + result.shape, result.dtype = _shape_and_dtype(self[1] @ other) return result def __rmatmul__(self, other): # Left multiplication by arrays result = self.zeros_like() result.data = {key: other @ val for key, val in self.items()} - result.shape = _find_shape(result.data) if result.data else (other @ self[1]).shape - result._isarray = any(isinstance(val, np.ndarray) for val in result.values()) + result.shape, result.dtype = _shape_and_dtype(other @ self[1]) return result def __truediv__(self, other): @@ -384,11 +404,10 @@ class Model(UserDict): def zeros_like(self): """Return an empty model object that inherits the other properties""" - result = type(self)() - result.interesting_keys = self.interesting_keys.copy() - result.momenta = self.momenta.copy() - result.shape = self.shape - result._isarray = self._isarray + result = type(self)(shape=self.shape, + dtype=self.dtype) + result.interesting_keys=self.interesting_keys.copy() + result.momenta=self.momenta return result def transform_symbolic(self, func): @@ -397,7 +416,7 @@ class Model(UserDict): # Add possible duplicate keys that only differ in constant factors result = sum((type(self)({func(key): copy(val)}, restructure_dict=True, - momenta=self.momenta.copy()) + momenta=self.momenta) for key, val in self.items()), self.zeros_like()) return result @@ -436,21 +455,21 @@ class Model(UserDict): elif isinstance(args[0], dict): # Input is a dictionary args = ([(key, value) for key, value in args[0].items()], ) - momenta = self.momenta.copy() + momenta = self.momenta for (old, new) in args[0]: # Substitution of a momentum variable with a symbol # is a renaming of the momentum. if old in momenta and isinstance(new, sympy.Symbol): - momenta = [momentum if old is not momentum else new - for momentum in momenta] + momenta = tuple(momentum if old is not momentum else new + for momentum in momenta) # If no momenta appear in the replacement for a momentum, we consider # that momentum removed. # Replacement is not a sympy object. elif not isinstance(new, sympy.Basic): - momenta = [momentum for momentum in momenta if old is not momentum] + momenta = tuple(momentum for momentum in momenta if old is not momentum) # Replacement is a sympy object, but does not contain momenta. elif not any([momentum in new.atoms() for momentum in momenta]): - momenta = [momentum for momentum in momenta if old is not momentum] + momenta = tuple(momentum for momentum in momenta if old is not momentum) substituted = self.transform_symbolic(lambda x: x.subs(*args, **kwargs)) substituted.momenta = momenta # If there are exponentials, evaluate any numerical exponents, @@ -491,7 +510,7 @@ class Model(UserDict): def trace(self): result = self.zeros_like() result.data = {key: np.sum(val.diagonal()) for key, val in self.items()} - result.shape = () + result.shape, result.dtype = (), np.complex128 return result def value_list(self, key_list): @@ -535,7 +554,7 @@ class Model(UserDict): else: # LinearOperator doesn't support multiplication with sparse matrix val = scipy.sparse.csr_matrix(val @ np.eye(val.shape[-1], dtype=complex), dtype=complex) - result._isarray = False + result.dtype = scipy.sparse.csr_matrix return result def toarray(self): @@ -549,7 +568,7 @@ class Model(UserDict): result[key] = val.A else: val = val @ np.eye(val.shape[-1], dtype=complex) - result._isarray = True + result.dtype = np.ndarray return result def copy(self): @@ -601,7 +620,7 @@ class Model(UserDict): def reshape(self, *args, **kwargs): result = self.zeros_like() result.data = {key: val.reshape(*args, **kwargs) for key, val in self.items()} - result.shape = _find_shape(result.data) + result.shape, _ = _shape_and_dtype(self[1].reshape(*args, **kwargs)) return result def allclose(self, other, rtol=1e-05, atol=1e-08, equal_nan=False): @@ -618,7 +637,7 @@ class Model(UserDict): class BlochModel(Model): def __init__(self, hamiltonian=None, locals=None, momenta=(0, 1, 2), - interesting_keys=None): + interesting_keys=None, shape=None, dtype=None): """ Class to efficiently store matrix valued Bloch Hamiltonians. The BlochModel represents sum(BlochCoeff * value), where BlochCoeff @@ -656,6 +675,19 @@ class BlochModel(Model): 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. + dtype : 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 dtype is not provided, it is inferred from the type of the values. + If hamiltonian is any other format, dtype is ignored an set to + np.ndarray. Ignored when initialized with Model. """ if hamiltonian is None: hamiltonian = {} @@ -665,9 +697,11 @@ class BlochModel(Model): locals=locals, momenta=hamiltonian.momenta, interesting_keys=hamiltonian.interesting_keys, - symbol_normalizer=lambda key: _to_bloch_coeff(key, hamiltonian.momenta)) + symbol_normalizer=lambda key: _to_bloch_coeff(key, hamiltonian.momenta), + shape=hamiltonian.shape, + dtype=hamiltonian.dtype) # set these in case it was and empty Model - self._isarray = hamiltonian._isarray + self.dtype = hamiltonian.dtype self.shape = hamiltonian.shape elif isinstance(hamiltonian, abc.Mapping): keys = hamiltonian.keys() @@ -680,14 +714,19 @@ class BlochModel(Model): momenta=momenta, interesting_keys=interesting_keys, symbol_normalizer=lambda x: x, - restructure_dict=False) + restructure_dict=False, + shape=shape, + dtype=dtype, + ) elif symbolic: # First cast it to model with restructuring, then try to interpret it as BlochModel self.__init__(Model(hamiltonian, locals=locals, momenta=momenta, interesting_keys=interesting_keys, - restructure_dict=True)) + restructure_dict=True, + shape=shape, + dtype=dtype)) else: raise ValueError('All keys must have the same type (sympy expression or BlochCoeff).') else: @@ -695,7 +734,9 @@ class BlochModel(Model): self.__init__(Model(hamiltonian, locals=locals, momenta=momenta, - interesting_keys=interesting_keys)) + interesting_keys=interesting_keys, + shape=shape, + dtype=dtype)) def transform_symbolic(self, func): raise NotImplementedError('transform_symbolic is not implemented for BlochModel') @@ -831,28 +872,32 @@ def _to_bloch_coeff(key, momenta): return bloch_coeff def _find_momenta(momenta): - if all([type(i) is int for i in momenta]): - return [_commutative_momenta[i] for i in momenta] + if all(type(i) is int for i in momenta): + return tuple(_commutative_momenta[i] for i in momenta) + elif all(m in _commutative_momenta for m in momenta): + return tuple(momenta) else: _momenta = [kwant_continuum.sympify(k) for k in momenta] - return [kwant_continuum.make_commutative(k, k) - for k in _momenta] - -def _find_shape(data): - # Make sure every matrix has the same size - if data == {}: - return None - else: - val = next(iter(data.values())) - if isinstance(val, Number): - shape = () - if not all(isinstance(v, Number) for v in data.values()): - raise ValueError('All terms must have the same shape') - # Recast numbers to numpy complex128 - for key, val in data.items(): - data[key] = np.complex128(val) - else: - shape = val.shape - if not all(v.shape == shape for v in data.values()): - raise ValueError('All terms must have the same shape') - return shape + return tuple(kwant_continuum.make_commutative(k, k) + for k in _momenta) + +def _shape_and_dtype(val): + # Find shape and type of val + dtype = type(val) + try: + shape = val.shape + except AttributeError: + # Treat it as a scalar + shape = () + if issubclass(dtype, Number): + # Cast all numbers to np.complex128 + dtype = np.complex128 + elif issubclass(dtype, np.ndarray): + dtype = np.ndarray + elif issubclass(dtype, scipy.sparse.linalg.LinearOperator): + # Make all subclasses of LinearOperator work + dtype = scipy.sparse.linalg.LinearOperator + elif not issubclass(dtype, scipy.sparse.spmatrix): + raise ValueError('Only dtypes which are subclasses of np.ndarray, scipy.sparse.spmatrix ' + 'scipy.sparse.linalg.LinearOperator or Number are supported.') + return shape, dtype diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py index 695a0cf..477a379 100644 --- a/qsymm/tests/test_model.py +++ b/qsymm/tests/test_model.py @@ -19,7 +19,7 @@ def test_dense_algebra(): # scalar models a = Model({1: 1, 'x': 3}) assert a.shape == () - assert a._isarray == False + assert a.dtype is np.complex128 assert a * 0 == 0 * a assert 0 * a == a.zeros_like() assert 0 * a == {} @@ -39,19 +39,18 @@ def test_dense_algebra(): vec = np.ones((3,)) c = a * vec assert c.shape == (3,) - assert c._isarray == True + assert c.dtype == np.ndarray # matrix models mat = np.ones((3, 3)) d = b * mat assert d.shape == (3, 3) - assert d._isarray == True + assert d.dtype == np.ndarray assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) assert c.T() @ d == 3 * a * b assert c.T() @ d @ c == 9 * a * a * b assert d @ d == 3 * b * d - assert (d @ d)._isarray == True - # numpy broadcasting rules apply - assert d + c == (a + b) * mat + assert (d @ d).dtype == np.ndarray + # numpy elemntwise multiplication assert d * c == a * b * mat assert d.trace() == 3 * b assert d.reshape(-1) == b * np.ones((9,)) @@ -65,33 +64,33 @@ def test_sparse_algebra(): vec = csr_matrix(np.ones((3, 1))) c = a * vec assert c.shape == (3, 1) - assert c._isarray == False + assert c.dtype is csr_matrix mat = csr_matrix(np.ones((3, 3))) d = b * mat assert d.shape == (3, 3) - assert d._isarray == False + assert d.dtype is csr_matrix assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) assert d @ d == 3 * b * d assert (d @ d) @ c == 9 * b * b * c - assert (d @ d)._isarray == False + assert (d @ d).dtype is csr_matrix assert c.T() @ d == 3 * b * c.T() assert c.T() @ d @ c == 9 * a * a * b * np.eye(1) assert d.trace() == 3 * b assert d.reshape((1, 9)) @ np.ones((9,)) == 9 * b e = d @ np.ones((3,)) assert e == 3 * b * np.ones((3,)) - assert e._isarray == True + assert e.dtype is np.ndarray # Test LinearOperator d = b * LinearOperator((3, 3), matvec=lambda v: mat @ v) c = a * np.ones((3, 1)) assert d.shape == (3, 3) - assert d._isarray == False + assert d.dtype is LinearOperator assert d @ d == 3 * b * d assert (d @ d) @ c == 9 * b * b * c - assert (d @ d)._isarray == False + assert (d @ d).dtype is LinearOperator assert d @ c == 3 * Model({1: 2*vec, 'x': 6*vec, 'x**2': 3*vec, 'x**3': 9*vec}) - assert (d @ c)._isarray == True + assert (d @ c).dtype is np.ndarray assert c.T() @ (d @ c) == 9 * a * a * b @@ -172,7 +171,7 @@ def test_Model(): assert prod[keys[1]] == m1[k_x] * m2[k_y] assert (m1 * np.array([2]))[k_x] == 2*m1[k_x] - m2.momenta = [k_x] + m2.momenta = (k_x,) with raises(ValueError, message="Only multiplication of Models with identical momenta allowed."): m1*m2 @@ -212,7 +211,7 @@ def test_Model_subs(): c1 * e**(I*(4*k_x + 3*k_y)) : 2*T}, momenta=[0, 1]) u_1, u_2 = sympy.symbols('u_1 u_2', real=True) nHam = Ham.subs({k_x: u_1, k_y: 2*k_y + 1}) - assert nHam.momenta == [u_1, k_y] + assert nHam.momenta == (u_1, k_y) right_keys = [sympy.simplify(c0*e**(-2*I*k_y - I*u_1/2)), sympy.simplify(c1*e**(6*I*k_y + 4*I*u_1))] nHam_keys = [sympy.simplify(key) for key in nHam.keys()] @@ -224,7 +223,7 @@ def test_Model_subs(): nHam = Ham.subs(k_y, 1.5) - assert nHam.momenta == [k_x] + assert nHam.momenta == (k_x,) right_keys = [sympy.simplify(c0*e**(-I*k_x/2)), sympy.simplify(c1*e**(4*I*k_x))] nHam_keys = [sympy.simplify(key) for key in nHam.keys()] @@ -236,7 +235,7 @@ def test_Model_subs(): nHam = Ham.subs(k_y, sympy.sqrt(3) + u_1) - assert nHam.momenta == [k_x] + assert nHam.momenta == (k_x,) Ham = BlochModel({c0 * e**(-I*(k_x/2 + k_y )) : T, c1 * e**(I*(4*k_x + 3*k_y)) : 2*T}, momenta=[0, 1]) diff --git a/qsymm/tests/test_mutual.py b/qsymm/tests/test_mutual.py index 2b875f6..c92dbfb 100644 --- a/qsymm/tests/test_mutual.py +++ b/qsymm/tests/test_mutual.py @@ -31,10 +31,10 @@ def test_mutual_continuum(): subgroups = generate_subgroups(group) for sg, gen in subgroups.items(): families = continuum_hamiltonian(list(gen), dim, 3) - H = Model({sympify('a_' + str(i)) * k: v - for i, fam in enumerate(families) for k, v in fam.items()}, - momenta = range(dim)) - if not H == {}: + if not len(families) == 0: + H = Model({sympify('a_' + str(i)) * k: v + for i, fam in enumerate(families) for k, v in fam.items()}, + momenta = range(dim)) sg2, Ps = discrete_symmetries(H, groupnoU) # new symmetry group may bigger because of additional constraints assert sg2 >= sg, (sg2, sg) @@ -54,10 +54,10 @@ def test_mutual_continuum(): subgroups = generate_subgroups(group) for sg, gen in subgroups.items(): families = continuum_hamiltonian(list(gen), dim, 3) - H = Model({sympify('a_' + str(i)) * k: v - for i, fam in enumerate(families) for k, v in fam.items()}, - momenta = range(dim)) - if not H == {}: + if not len(families) == 0: + H = Model({sympify('a_' + str(i)) * k: v + for i, fam in enumerate(families) for k, v in fam.items()}, + momenta = range(dim)) sg2, Ps = discrete_symmetries(H, groupnoU) # new symmetry group may bigger because of additional constraints assert sg2 == sg, (sg2, sg) @@ -80,10 +80,10 @@ def test_mutual_continuum(): subgroups = generate_subgroups(group) for sg, gen in subgroups.items(): families = continuum_hamiltonian(list(gen), dim, 1) - H = Model({sympify('a_' + str(i)) * k: v - for i, fam in enumerate(families) for k, v in fam.items()}, - momenta = range(dim)) - if not H == {}: + if not len(families) == 0: + H = Model({sympify('a_' + str(i)) * k: v + for i, fam in enumerate(families) for k, v in fam.items()}, + momenta = range(dim)) sg2, Ps = discrete_symmetries(H, groupnoU) assert sg2 == sg, (sg2, sg) for g1, g2 in it.product(sg, sg2): @@ -110,10 +110,10 @@ def test_mutual_continuum(): for g in groupnoU: g.U = None families = continuum_hamiltonian(list(gens), dim, 1) - H = Model({sympify('a_' + str(i)) * k: v - for i, fam in enumerate(families) for k, v in fam.items()}, - momenta = range(dim)) - if not H == {}: + if not len(families) == 0: + H = Model({sympify('a_' + str(i)) * k: v + for i, fam in enumerate(families) for k, v in fam.items()}, + momenta = range(dim)) sg2, Ps = discrete_symmetries(H, groupnoU) assert sg2 == group, (sg2, group) for g1, g2 in it.product(group, sg2): -- GitLab From 555418e7e7813cd494649b25c543b5f9143d73b9 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Wed, 22 May 2019 20:15:33 +0200 Subject: [PATCH 17/41] make symmetry finder work with sparse matrices --- qsymm/linalg.py | 64 ++++++++++++++--------------- qsymm/symmetry_finder.py | 39 +++++++++++++----- qsymm/tests/test_symmetry_finder.py | 35 ++++++++++++++++ 3 files changed, 94 insertions(+), 44 deletions(-) diff --git a/qsymm/linalg.py b/qsymm/linalg.py index 767f4c3..8aedb41 100644 --- a/qsymm/linalg.py +++ b/qsymm/linalg.py @@ -170,12 +170,12 @@ def family_to_vectors(family, all_keys=None): return np.vstack(basis_vectors) -def nullspace(A, atol=1e-6, return_complement=False, sparse=False, k_max=-10): +def nullspace(A, atol=1e-6, return_complement=False, sparse=None, k_max=-10): """Compute an approximate basis for the nullspace of matrix A. Parameters: ----------- - A : ndarray + A : ndarray or spmatrix A should be 2-D. atol : real Absolute tolerance when deciding whether an eigen/singular @@ -184,9 +184,10 @@ def nullspace(A, atol=1e-6, return_complement=False, sparse=False, k_max=-10): return_complement : bool Whether to return the basis of the orthocomplement ot the null space as well. Treated as False if sparse is True. - sparse : bool + sparse : bool or None (default) Whether to use sparse linear algebra to find the null space. - k_max : int + The default value is inferred from the type of the input array. + k_max : int (default -10) If k_max<0, the full null-space is found by increasing the number of eigenvectors found in each step by abs(k_max). If k_max>0, k_max is the maximum number of null space vectors requested. @@ -204,13 +205,15 @@ def nullspace(A, atol=1e-6, return_complement=False, sparse=False, k_max=-10): othonormal basis of the row span of A. Only returned if return_complement=True and sparse=False. """ - + if sparse is None: + sparse = isinstance(A, scipy.sparse.spmatrix) if sparse: # Do sparse eigenvalue solving # sigma used for shift-invert, should be a small negative value sigma = -1e-1 # Make sure A is sparse matrix - A = scipy.sparse.csr_matrix(A) + if isinstance(A, np.ndarray): + A = scipy.sparse.csr_matrix(A) A.eliminate_zeros() # Treat A=0 case if np.allclose(A.data, 0, atol=atol/A.shape[1]): @@ -250,6 +253,8 @@ def nullspace(A, atol=1e-6, return_complement=False, sparse=False, k_max=-10): ns, _ = la.qr(ns, mode='economic') return ns.T else: + if isinstance(A, scipy.sparse.spmatrix): + A = A.A # Do dense SVD u, s, vh = la.svd(A, full_matrices = True) nnz = np.isclose(s, 0, atol=atol) @@ -375,14 +380,15 @@ def simult_diag(mats, tol=1e-6, checks=0): return Ps -def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False, sparse=False, k_max=-10): +def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, + conjugate=False, sparse=None, k_max=-10): """Solve for X the simultaneous matrix equatioins X HL[i] = HR[i] X for every i. It is mapped to a system of linear equations, the null space of which gives a basis for all sulutions. Parameters: ----------------- - HL : ndarray or list of ndarrays + HL : ndarray or list of arrays or sparse matrices Coefficient matrices of identical square shape, list of arrays of shape (n, n) or one array of shape (m, n, n). HR : ndarray or list of ndarrays or None @@ -394,8 +400,9 @@ def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False If True, solve X HL[i] = HR[i] X^* instead. If a list with the same length as HL and HR is provided, conjugation is applied to the equations with index corresponding to the True entries. - sparse : bool - Whether to use sparse linear algebra to find the solutions. + sparse : bool or None (default) + Whether to use sparse linear algebra to find the null space. + The default value is inferred from the type of the input array. k_max : int If k_max<0, all solutions are found by increasing the number of soulutions sought in each step by abs(k_max). @@ -407,27 +414,24 @@ def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False Returns: --------------- - ndarray of shape (l, n, n), list of linearly independent square matrices + ndarray of shape (l, n, n), or list of csr_matrix + List of linearly independent square matrices that span all solutions of the eaquations. """ - HL = np.array(HL) + if sparse is None: + sparse = isinstance(HL[0], scipy.sparse.spmatrix) if HR is None: HR = HL - else: - HR = np.array(HR) - if HL.shape != HR.shape: - raise ValueError('HL and HR must have the same shape') + if len(HL) != len(HR) or not all(l.shape == r.shape for l, r in zip(HL, HR)): + raise ValueError('HL and HR must have the same shape.') if isinstance(conjugate, bool): conjugate = [conjugate] * len(HL) if len(conjugate) != len(HL): - raise ValueError('conugate must have the same length as HL') - if len(HL.shape) == 3: - if HL.shape[1] != HL.shape[2]: - raise ValueError('HL and HR must be (a list of) square matrices') - else: - raise ValueError('HL and HR must have the shape (n, n) or (m, n, n)') + raise ValueError('Conugate must have the same length as HL.') + if not all(l.shape[0] == l.shape[1] == r.shape[0] == r.shape[1] for l, r in zip(HL, HR)): + raise ValueError('HL and HR must be a list of square matrices.') - dim = HL.shape[-1] + dim = HL[0].shape[-1] # number of basis matrices N = dim**2 - (1 if traceless else 0) if N == 0: @@ -435,14 +439,7 @@ def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False # Prepare for differences in sparse and dense algebra if sparse: - # It is worth doing sparse multiplication if the matrices are - # over 100 x 100 - if HL.shape[-1] >= 100: - HL = [scipy.sparse.csr_matrix(hL) for hL in HL] - HR = [scipy.sparse.csr_matrix(hR) for hR in HR] - basis = lambda: matrix_basis(dim, traceless=traceless, sparse=True) - else: - basis = lambda: matrix_basis(dim, traceless=traceless, sparse=False) + basis = lambda: matrix_basis(dim, traceless=traceless, sparse=True) vstack = scipy.sparse.vstack bmat = scipy.sparse.bmat # Cast it to coo format and reshape @@ -457,9 +454,10 @@ def solve_mat_eqn(HL, HR=None, hermitian=False, traceless=False, conjugate=False null_mat = [] for hL, hR, conj in zip(HL, HR, conjugate): if conj: - row = [flatten(mat.dot(hL) - hR.dot(mat.conj())) for mat in basis()] + print(next(basis)) + row = [flatten(mat @ hL - hR @ mat.conj()) for mat in basis()] else: - row = [flatten(mat.dot(hL) - hR.dot(mat)) for mat in basis()] + row = [flatten(mat @ hL - hR @ mat) for mat in basis()] null_mat.append(row) null_mat = bmat(null_mat) diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index a6e5e27..9ed8a13 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -1,5 +1,6 @@ import numpy as np import scipy.linalg as la +import scipy.sparse import itertools as it from copy import deepcopy @@ -68,6 +69,9 @@ def symmetries(model, candidates=None, continuous_rotations=False, rotation generators. The trivial conserved quantity proportional to identity is not included. """ + if not issubclass(model.dtype, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Symmetry finding is only supported for Models with ' + '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) @@ -273,7 +277,7 @@ def _reduce_hamiltonian(H, sparse=False): Parameters: ----------------- H : ndarray - 3 index array of shape (m, N, N), list of NxN Hermitian matrices + 3 index array of shape (m, N, N), or list of NxN Hermitian matrices defining the family. sparse : bool Whether to use sparse linear algebra in the calculation. @@ -294,11 +298,11 @@ def _reduce_hamiltonian(H, sparse=False): gens = solve_mat_eqn(H, hermitian=True, traceless=True, sparse=sparse) if len(gens) == 0: - return (np.array([np.eye(H.shape[-1])]),) + return (np.array([np.eye(H[0].shape[-1])]),) # Find the center gensc, _ = separate_lie_algebra(gens) if len(gensc) == 0: - Ps = [np.eye(H.shape[-1])] + Ps = [np.eye(H[0].shape[-1])] else: # Simultaneaously diagonalize central generators Ps = simult_diag(gensc) @@ -306,11 +310,12 @@ def _reduce_hamiltonian(H, sparse=False): Preduced = tuple() # Find symmetry adapted basis of restricted LA in each subspace for P in Ps: - Hr = mtm(P.T.conjugate(), H, P) + # Suppot list of sparse matrices + Hr = [P.T.conjugate() @ h @ P for h in H] gensr = solve_mat_eqn(Hr, hermitian=True, traceless=True, sparse=sparse) # Find symmetry adapted basis in subspace if len(gensr) == 0: - P2s = [np.eye(Hr.shape[-1])] + P2s = [np.eye(Hr[0].shape[-1])] else: P2s = symmetry_adapted_sun(gensr) Preduced += (np.array([np.dot(P, P2i) for P2i in P2s]),) @@ -394,11 +399,14 @@ def discrete_symmetries(model, candidates, Ps=None, generators=False, Ps : ndarray Projectors as returned by _reduce_hamiltonian. """ + if not issubclass(model.dtype, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Symmetry finding is only supported for Models with ' + 'np.ndarray or scipy.sparse.spmatrix data type.') symmetry_candidates = deepcopy(candidates) m = len(symmetry_candidates) # Reduce Hamiltonian by onsite unitaries if not Ps: - Ps = _reduce_hamiltonian(np.array(list(model.values())), sparse=sparse_linalg) + Ps = _reduce_hamiltonian(list(model.values()), sparse=sparse_linalg) # 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]) @@ -473,6 +481,9 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False): Point group operator with gr.U set to the Hilbert space action of the symmetry is found, otherwise identical to g. """ + if not issubclass(model.dtype, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Symmetry finding is only supported for Models with ' + 'np.ndarray or scipy.sparse.spmatrix data type.') if g.U is not None: raise ValueError('g.U must be None.') Rmodel = g.apply(model) @@ -486,14 +497,13 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False): matR = Rmodel[key] HL.append(matR) ev_test = ev_test and allclose(matL, matL.T.conj()) and allclose(matR, matR.T.conj()) - HR, HL = np.array(HR), np.array(HL) # Need to carry conjugation on left side through P if g.conjugate: PsL = [P.conj() for P in Ps] else: PsL = Ps - HRs = [mtm(P[0].T.conj(), HR, P[0]) for P in Ps] - HLs = [mtm(P[0].T.conj(), HL, P[0]) for P in PsL] + 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] squares_to_1 = g * g == g.identity() block_dict = _find_unitary_blocks(HLs, HRs, Ps, conjugate=g.conjugate, ev_test=ev_test, @@ -501,6 +511,7 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False): S = _construct_unitary(block_dict, Ps, conjugate=g.conjugate, squares_to_1=squares_to_1) 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: @@ -661,11 +672,14 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li symmetries : list of ContinuousGroupGenerator List of linearly independent symmetry generators. """ + if not issubclass(model.dtype, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Symmetry finding is only supported for Models with ' + 'np.ndarray or scipy.sparse.spmatrix data type.') if isinstance(model, BlochModel): # BlochModel cannot have continuous rotation symmetry return [] if not Ps: - Ps = _reduce_hamiltonian(np.array(list(model.values())), sparse=sparse_linalg) + Ps = _reduce_hamiltonian(list(model.values()), sparse=sparse_linalg) reduced_hamiltonians = _reduced_model(model, Ps) dim = len(model.momenta) if dim <= 1: @@ -730,7 +744,7 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li symmetries.append(g) # Check that it is a symmetry trf = g.apply(model) - assert trf == trf.zeros_like(), (trf, g) + assert trf.allclose(0, atol=1e-6), (trf, g) return symmetries @@ -752,6 +766,9 @@ def _reduced_model(model, Ps=None): List of reduced Hamiltonian families, each projected on the symmetry irreducible subspaces. """ + if not issubclass(model.dtype, (np.ndarray, scipy.sparse.spmatrix)): + raise ValueError('Symmetry finding is only supported for Models with ' + 'np.ndarray or scipy.sparse.spmatrix data type.') if Ps is None: Ps = _reduce_hamiltonian(np.array([H for H in model.values()])) reduced_hamiltonians = [] diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index cb5e9ee..bb41394 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -506,6 +506,24 @@ def test_continuum(): 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) + 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) + 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) + 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 + def test_bloch(): # Simple tests for Bloch models @@ -557,6 +575,23 @@ def test_bloch(): assert len(sg) == 48 assert sg == generate_group({Mx, C6, TR, PH}) + # Test sparse + H64 = Model(ham64, momenta=[0, 1]) + sg, Ps = discrete_symmetries(H64, hex_group_2D, sparse_linalg=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) + 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) + 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}) + + def test_bravais_symmetry(): # 2D -- GitLab From 66ff7e218afaf33655cfd6c6d47221598a624e6d Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Sun, 17 Feb 2019 23:12:03 +0100 Subject: [PATCH 18/41] add factory functions for fundamental discrete symmetries Closes #5. --- qsymm/__init__.py | 4 +++- qsymm/groups.py | 57 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/qsymm/__init__.py b/qsymm/__init__.py index 58d3260..92a6f88 100644 --- a/qsymm/__init__.py +++ b/qsymm/__init__.py @@ -3,7 +3,9 @@ from . import linalg from . import hamiltonian_generator from . import symmetry_finder from . import model -from .groups import PointGroupElement, ContinuousGroupGenerator +from .groups import ( + PointGroupElement, ContinuousGroupGenerator, + time_reversal_operator, particle_hole_operator, chiral_operator) from .model import Model, BlochModel from .hamiltonian_generator import continuum_hamiltonian, continuum_pairing, display_family, \ bloch_family diff --git a/qsymm/groups.py b/qsymm/groups.py index 96bb4c7..d704a04 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -294,6 +294,63 @@ def identity(dim, shape=None): return PointGroupElement(R, False, False, U) +def time_reversal_operator(realspace_dim, U=None): + """Return a time-reversal symmetry operator + + parameters + ---------- + realspace_dim : int + Realspace dimension + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + T : PointGroupElement + """ + R = ta.identity(realspace_dim, int) + return PointGroupElement(R, conjugate=True, antisymmetry=False, U=U) + + +def particle_hole_operator(realspace_dim, U=None): + """Return a particle-hole symmetry operator + + parameters + ---------- + realspace_dim : int + Realspace dimension + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + P : PointGroupElement + """ + R = ta.identity(realspace_dim, int) + return PointGroupElement(R, conjugate=True, antisymmetry=True, U=U) + + +def chiral_operator(realspace_dim, U=None): + """Return a chiral symmetry operator + + parameters + ---------- + realspace_dim : int + Realspace dimension + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + P : PointGroupElement + """ + R = ta.identity(realspace_dim, int) + return PointGroupElement(R, conjugate=False, antisymmetry=True, U=U) + + class ContinuousGroupGenerator: """A generator of a continuous group. -- GitLab From b3d2df9f2cc181e816241ed721b12e9eb4885ece Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 25 Feb 2019 16:41:18 +0100 Subject: [PATCH 19/41] use in tests --- qsymm/symmetry_finder.py | 3 ++- qsymm/tests/test_symmetry_finder.py | 36 ++++++++++++++--------------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index 9ed8a13..10e81bb 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -9,7 +9,8 @@ from .linalg import matrix_basis, nullspace, split_list, simult_diag, commutator allclose from .model import Model, BlochModel, BlochCoeff from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, \ - set_multiply, L_matrices, spin_rotation + set_multiply, L_matrices, spin_rotation, time_reversal_operator, \ + particle_hole_operator, chiral_operator from . import kwant_continuum from .kwant_linalg_lll import lll, cvp, voronoi diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index bb41394..1a15133 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -103,8 +103,8 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -173,8 +173,8 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -243,8 +243,8 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -316,8 +316,8 @@ def test_disc_finder(verbose = False): # Find the symmetry operator Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(1, 2*dim, dim), (1, 2*dim, dim)] assert len(sg) == 2 @@ -355,8 +355,8 @@ def test_disc_finder(verbose = False): assert np.allclose(H, T.dot(H.conj()).dot(T.T.conj())) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(1, 3*dim, dim), (1, 3*dim, dim), (1, 3*dim, dim)] assert len(sg) == 2 @@ -392,8 +392,8 @@ def test_disc_finder(verbose = False): assert np.allclose(H, T.dot(H.conj()).dot(T.T.conj())) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(2, 2*dim, dim)] assert len(sg) == 2 @@ -427,8 +427,8 @@ def test_disc_finder(verbose = False): h1, t_mat.dot(h1.conj()).dot(t_mat.T.conj()))) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {PointGroupElement(np.empty((0, 0)), c, a) - for c, a in [[True, False], [True, True], [False, True]]} + sg = {time_reversal_operator(0), particle_hole_operator(0), + chiral_operator(0)} sg, Ps = discrete_symmetries(Hs, sg) 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 @@ -468,8 +468,8 @@ def test_continuum(): C3 = PointGroupElement(np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]], int), False, False, None) - TR = PointGroupElement(eye, True, False, None) - PH = PointGroupElement(eye, True, True, None) + TR = time_reversal_operator(realspace_dim=3) + PH = particle_hole_operator(realspace_dim=3) cubic_gens = {I, C4, C3, TR, PH} cubic_group = generate_group(cubic_gens) assert len(cubic_group) == 192 @@ -535,8 +535,8 @@ def test_bloch(): C6 = PointGroupElement(sympy.ImmutableMatrix([[sympy.Rational(1, 2), sympy.sqrt(3)/2], [-sympy.sqrt(3)/2, sympy.Rational(1, 2)]]), False, False, None) - TR = PointGroupElement(eyesym, True, False, None) - PH = PointGroupElement(eyesym, True, True, None) + TR = time_reversal_operator(realspace_dim=2) + PH = particle_hole_operator(realspace_dim=2) gens_hex_2D ={Mx, C6, TR, PH} hex_group_2D = generate_group(gens_hex_2D) assert len(hex_group_2D) == 48 -- GitLab From 6cdeb832ff309c22babff66de3c423065ae9e912 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 01:03:18 +0200 Subject: [PATCH 20/41] better representation of PointGroupElement --- qsymm/groups.py | 143 ++++++++++++++++++++++++------------------------ 1 file changed, 73 insertions(+), 70 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index d704a04..6e9a5e4 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -107,8 +107,7 @@ class PointGroupElement: self._strict_eq = _strict_eq def __repr__(self): - return ('PointGroupElement(\n{},\n{},\n{},\n{},\n)' - .format(self.R, self.conjugate, self.antisymmetry, self.U)) + return name_PG_elements(self, full=True) def __eq__(self, other): # We do not allow mixing of PointGroupElements @@ -581,81 +580,85 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): ## Human readable group element names -def name_PG_elements(g): - """Automatically name PG operators in 3D cubic group""" - - def find_order(g): - dim = g.shape[0] - gpower = g - order = 1 - while not allclose(gpower, np.eye(dim)): - gpower = g.dot(gpower) - order += 1 - return order - - def find_axis(g): - eig = la.eig(g) - n = [n for i, n in enumerate(eig[1].T) if np.isclose(eig[0][i], 1)] - assert len(n) == 1 - n = n[0] - if not np.isclose(n[0], 0): - n = np.real(n/n[0]) - elif not np.isclose(n[1], 0): - n = np.real(n/n[1]) +def name_PG_elements(g, full=False): + """Human readable format of PointGroupElement""" + + def name_angle(theta): + frac = int(np.round(6 * theta / (np.pi), 2)) + gcd = np.gcd(frac, 6) + angle = '{}π{}'.format( + "" if gcd == frac else ("-" if gcd == -frac else str(frac//gcd)), + "" if gcd == 6 else ("/" + str(6//gcd))) + return angle + + R = np.array(g.R).astype(float) + if R.shape[0] == 1: + if R[0, 0] == 1: + rot_name = '1' else: - n = np.real(n/n[2]) - return n - - def axisname(n): - epsilon = 0.2 - return ''.join([(s if n[i] > epsilon else '') + - ('-' + s if n[i] < -epsilon else '') - for i, s in enumerate(['x', 'y', 'z'])]) - - def find_direction(g, n): - if allclose(g.dot(g), np.eye(g.shape[0])): - return '' - v = np.random.random(3) - p = np.cross(v, g.dot(v)).dot(n) - if p < 0.0001: - return '-' - else: - return '' - - name = '' - if g.conjugate and not g.antisymmetry: - name += 'T' - if g.conjugate and g.antisymmetry: - name += 'P' - - if not g.conjugate and g.antisymmetry: - name += 'TP' - g = np.array(g.R).astype(np.float_) - # if less than 3D, pad with identity - if g.shape[0] < 3: - g = la.block_diag(g, np.eye(3 - g.shape[0])) - dim = g.shape[0] - det = la.det(g) - evals = la.eigvals(g) - if np.isclose(det, 1): - if allclose(g, np.eye(dim)): - name += 'Id' + rot_name = 'M' + elif R.shape[0] == 2: + if la.det(R) == 1: + # pure rotation + theta = np.arctan2(R[1, 0], R[0, 0]) + if np.isclose(theta, 0): + rot_name = '1' + else: + rot_name = f'R({name_angle(theta)})' else: - n = find_axis(g) - name += 'C' + find_direction(g, n) + str(find_order(g)) + axisname(n) - else: - if allclose(g, -np.eye(dim)): - name += 'I' + # mirror + val, vec = la.eigh(R) + assert np.allclose(val, [-1, 1]) + rot_name = f'M({np.round(vec[0], 2)})' + elif R.shape[0] == 3: + if la.det(R) == 1: + # pure rotation + n, theta = rotation_to_angle(R) + if np.isclose(theta, 0): + rot_name = '1' + else: + rot_name = f'R({name_angle(theta)}, {np.round(n, 2)})' else: - order = find_order(-g) - n = find_axis(-g) - if order == 2: - name += 'M' + axisname(n) + # rotoinversion + n, theta = rotation_to_angle(-R) + if np.isclose(theta, 0): + # inversion + rot_name = 'I' + elif np.isclose(theta, np.pi): + # mirror + rot_name = f'M({np.round(n, 2)})' else: - name += 'S' + find_direction(g, n) + str(order) + axisname(n) + # generic rotoinversion + rot_name = f'S({name_angle(theta)}, {np.round(n, 2)})' + + if full: + name = 'U⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", + "-" if g.antisymmetry else "", + "-" if g.conjugate else "") + name += f'R = {rot_name}\n' + if g.U is not None: + name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) + else: + name = rot_name + if g.conjugate and not g.antisymmetry: + name += " T" + elif g.conjugate and g.antisymmetry: + name += " P" + elif not g.conjugate and g.antisymmetry: + name += " C" return name +def rotation_to_angle(R): + # Convert 3D rotation matrix to axis and angle + assert allclose(R, R.real) + n = np.array([1j * np.trace(L @ R) for L in L_matrices()]).real + absn = la.norm(n) + theta = np.arctan2(absn, np.trace(R).real - 1) + if not np.isclose(absn, 0): + n /= absn + return n, theta + ## Predefined representations def spin_matrices(s, include_0=False): -- GitLab From 951b9dfd62e09f3bc124704f6b691a4b215ad89f Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 10:37:52 +0200 Subject: [PATCH 21/41] some fixes in naming --- qsymm/groups.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 6e9a5e4..2cffdb7 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -581,7 +581,7 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): ## Human readable group element names def name_PG_elements(g, full=False): - """Human readable format of PointGroupElement""" + """Human readable format of PGE""" def name_angle(theta): frac = int(np.round(6 * theta / (np.pi), 2)) @@ -591,6 +591,14 @@ def name_PG_elements(g, full=False): "" if gcd == 6 else ("/" + str(6//gcd))) return angle + def round_axis(n): + # Try to find integer axis + for vec in it.product([-1, 0, 1], repeat=len(n)): + if np.isclose(vec @ n, la.norm(vec)) and not np.isclose(la.norm(vec), 0): + return vec + # otherwise round + return np.round(n, 2) + R = np.array(g.R).astype(float) if R.shape[0] == 1: if R[0, 0] == 1: @@ -608,8 +616,8 @@ def name_PG_elements(g, full=False): else: # mirror val, vec = la.eigh(R) - assert np.allclose(val, [-1, 1]) - rot_name = f'M({np.round(vec[0], 2)})' + assert allclose(val, [-1, 1]) + rot_name = f'M({round_axis(n)})' elif R.shape[0] == 3: if la.det(R) == 1: # pure rotation @@ -617,7 +625,7 @@ def name_PG_elements(g, full=False): if np.isclose(theta, 0): rot_name = '1' else: - rot_name = f'R({name_angle(theta)}, {np.round(n, 2)})' + rot_name = f'R({name_angle(theta)}, {round_axis(n)})' else: # rotoinversion n, theta = rotation_to_angle(-R) @@ -626,10 +634,10 @@ def name_PG_elements(g, full=False): rot_name = 'I' elif np.isclose(theta, np.pi): # mirror - rot_name = f'M({np.round(n, 2)})' + rot_name = f'M({round_axis(n)})' else: # generic rotoinversion - rot_name = f'S({name_angle(theta)}, {np.round(n, 2)})' + rot_name = f'S({name_angle(theta)}, {round_axis(n)})' if full: name = 'U⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", @@ -653,12 +661,21 @@ def rotation_to_angle(R): # Convert 3D rotation matrix to axis and angle assert allclose(R, R.real) n = np.array([1j * np.trace(L @ R) for L in L_matrices()]).real - absn = la.norm(n) - theta = np.arctan2(absn, np.trace(R).real - 1) - if not np.isclose(absn, 0): + # Choose direction to minimize number of minus signs + absn = la.norm(n) * np.sign(np.sum(n)) + if np.isclose(absn, 0): + # n is zero for 2-fold rotation, find eigenvector with +1 + val, vec = la.eigh(R) + assert np.isclose(val[-1], 1) + n = vec[-1] + # Choose direction to minimize number of minus signs + n /= np.sign(np.sum(n)) + else: n /= absn + theta = np.arctan2(absn, np.trace(R).real - 1) return n, theta + ## Predefined representations def spin_matrices(s, include_0=False): -- GitLab From 0ebb0d6cab2608d2cd0e79d143d26a12bc3a2868 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 11:52:32 +0200 Subject: [PATCH 22/41] use fractions for angle --- qsymm/groups.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 2cffdb7..cd0a687 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -5,6 +5,7 @@ import tinyarray as ta import scipy.linalg as la import itertools as it import functools as ft +from fractions import Fraction from collections import OrderedDict import sympy from copy import deepcopy @@ -584,11 +585,11 @@ def name_PG_elements(g, full=False): """Human readable format of PGE""" def name_angle(theta): - frac = int(np.round(6 * theta / (np.pi), 2)) - gcd = np.gcd(frac, 6) + frac = Fraction(theta / np.pi).limit_denominator(100) + num, den = frac.numerator, frac.denominator angle = '{}π{}'.format( - "" if gcd == frac else ("-" if gcd == -frac else str(frac//gcd)), - "" if gcd == 6 else ("/" + str(6//gcd))) + "" if num == 1 else ("-" if num == -1 else str(num)), + "" if den == 1 else ("/" + str(den))) return angle def round_axis(n): -- GitLab From 9a29592a7b1241bd2ef5d6132f9df4fa316747e8 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 15:07:01 +0200 Subject: [PATCH 23/41] add factory functions for other symmetries and add docstring for name_PG_elements --- qsymm/groups.py | 133 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 127 insertions(+), 6 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index cd0a687..73abb63 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -351,6 +351,89 @@ def chiral_operator(realspace_dim, U=None): return PointGroupElement(R, conjugate=False, antisymmetry=True, U=U) +def inversion(realspace_dim, U=None): + """Return an inversion operator + + parameters + ---------- + realspace_dim : int + Realspace dimension + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + P : PointGroupElement + """ + R = -ta.identity(realspace_dim, int) + return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) + + +def rotation(angle, axis=None, inversion=False, U=None): + """Return a rotation operator + + parameters + ---------- + angle : float + Rotation angle in units of 2 pi. + axis : ndarray or None (default) + Rotation axis, optional. If not provided, a 2D rotation is generated + around the axis normal to the plane. If a 3D vector is provided, + a 3D rotation is generated around this axis. Does not need to be + normalized to 1. + inversion : bool (default False) + Whether to generate a rotoinversion. By default a proper rotation + is returned. Only valid in 3D. + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + P : PointGroupElement + """ + angle = 2 * np.pi * angle + if axis is None: + # 2D + R = np.array([[np.cos(angle), -np.sin(angle)], + [np.sin(angle), np.cos(angle)]]) + elif len(axis) == 3: + # 3D + axis = np.array(axis) + R = spin_rotation(angle * axis / la.norm(axis), L_matrices(d=3, l=1)) + R *= (-1 if inversion else 1) + else: + raise ValueError('axis needs to be None or a 3D vector.') + # Try to round it to integer + R = _make_int(R) + return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) + + +def mirror(axis, U=None): + """Return a mirror operator + + parameters + ---------- + axis : ndarray + Normal of the mirror. The dimensionality of the operator is the same + as the length of axis. + U: ndarray (optional) + The unitary action on the Hilbert space. + May be None, to be able to treat symmetry candidates + + Returns + ------- + P : PointGroupElement + """ + axis = np.array(axis) + axis /= la.norm(axis) + R = np.eye(axis.shape[0]) - 2 * np.outer(axis, axis) + # Try to round it to integer + R = _make_int(R) + return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) + + class ContinuousGroupGenerator: """A generator of a continuous group. @@ -582,7 +665,43 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): ## Human readable group element names def name_PG_elements(g, full=False): - """Human readable format of PGE""" + """ + Return a human readable string representation of PointGroupElement + + Parameters + ---------- + + g : PointGroupElement + Point group element to be represented. + full : bool (default False) + Whether to return a full representation. + The default short representation only contains the real space action + and the symbol of the Altland-Zirnbauer part of the symmetry (see below). + The full representation presents the symmetry action on the Hamiltonian + and the unitary Hilbert-space action if set. + + Returns + ------- + name : string + In the short representation it is a sting rot_name + az_name. + In the long representation the first line is the action on the + Hamiltonian, the second line is rot_name and the third line + is the unitary action as a matrix, if set. + + rot_name can be (axis is not normalized): + - 1 for identity + - I for inversion (in 1D mirror is the same as inversion) + - R(angle) for 2D rotation + - R(angle, axis) for 3D rotation + - M(normal) for mirror + - S(angle, axis) for 3D + + az_name can be: + - T for time-reversal (antiunitary symmetry) + - P for particle-hole (antiunitary antisymmetry) + - C for chiral (unitary antisymmetry) + - az_name is ommited for unitary symmetries + """ def name_angle(theta): frac = Fraction(theta / np.pi).limit_denominator(100) @@ -605,7 +724,7 @@ def name_PG_elements(g, full=False): if R[0, 0] == 1: rot_name = '1' else: - rot_name = 'M' + rot_name = 'I' elif R.shape[0] == 2: if la.det(R) == 1: # pure rotation @@ -648,13 +767,15 @@ def name_PG_elements(g, full=False): if g.U is not None: name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) else: - name = rot_name if g.conjugate and not g.antisymmetry: - name += " T" + az_name = " T" elif g.conjugate and g.antisymmetry: - name += " P" + az_name = " P" elif not g.conjugate and g.antisymmetry: - name += " C" + az_name = " C" + else: + az_name = "" + name = az_name (if rot_name == '1' and az_name != "") else rot_name + az_name return name -- GitLab From 5db632b2ba51196906da41130439d2c9a5764fef Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 15:29:45 +0200 Subject: [PATCH 24/41] fix a typo and remove outdated docstring --- qsymm/groups.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 73abb63..025de54 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -53,7 +53,7 @@ class PointGroupElement: Parameters ---------- - R : sympy.ImmutableMatrix or integer array + R : sympy.ImmutableMatrix or array Real space rotation action of the operator. Square matrix with size of the number of spatial dimensions. conjugate : boolean (default False) @@ -74,13 +74,6 @@ class PointGroupElement: As U is floating point and has a phase ambiguity at least, it is ignored when comparing objects. - R must provide an exact representation, either as a - sympy.ImmutableMatrix or an integer array, as exact arithmetic is - assumed. If R is an integer array, it must be invertible over - the integers; this is always possible for crystallographic groups - in the basis of the translation vectors. - Performance is much higher if integer arrays are used. - R is the real space rotation acion. Do not include minus sign for the k-space action of antiunitary operators, such as time reversal. This minus sign will be included automatically if 'conjugate=True'. @@ -89,8 +82,6 @@ class PointGroupElement: __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_Rinv', '_strict_eq') def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False): - # Make sure that R is either an immutable sympy matrix, - # or an integer tinyarray. if isinstance(R, sympy.ImmutableMatrix): pass elif isinstance(R, ta.ndarray_int): @@ -688,13 +679,13 @@ def name_PG_elements(g, full=False): Hamiltonian, the second line is rot_name and the third line is the unitary action as a matrix, if set. - rot_name can be (axis is not normalized): + rot_name can be: - 1 for identity - I for inversion (in 1D mirror is the same as inversion) - R(angle) for 2D rotation - - R(angle, axis) for 3D rotation + - R(angle, axis) for 3D rotation (axis is not normalized) - M(normal) for mirror - - S(angle, axis) for 3D + - S(angle, axis) for 3D rotoinversion (axis is not normalized) az_name can be: - T for time-reversal (antiunitary symmetry) @@ -775,7 +766,7 @@ def name_PG_elements(g, full=False): az_name = " C" else: az_name = "" - name = az_name (if rot_name == '1' and az_name != "") else rot_name + az_name + name = az_name if (rot_name == '1' and az_name != "") else rot_name + az_name return name -- GitLab From a2c49348fd9f0efc8996e2583df0bd84fe611268 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 16:59:51 +0200 Subject: [PATCH 25/41] fix bug in naming --- qsymm/groups.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 025de54..2fbe234 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -717,7 +717,7 @@ def name_PG_elements(g, full=False): else: rot_name = 'I' elif R.shape[0] == 2: - if la.det(R) == 1: + if np.isclose(la.det(R), 1): # pure rotation theta = np.arctan2(R[1, 0], R[0, 0]) if np.isclose(theta, 0): @@ -727,10 +727,11 @@ def name_PG_elements(g, full=False): else: # mirror val, vec = la.eigh(R) - assert allclose(val, [-1, 1]) + assert allclose(val, [-1, 1]), R + n = vec[0] rot_name = f'M({round_axis(n)})' elif R.shape[0] == 3: - if la.det(R) == 1: + if np.isclose(la.det(R), 1): # pure rotation n, theta = rotation_to_angle(R) if np.isclose(theta, 0): @@ -779,7 +780,7 @@ def rotation_to_angle(R): if np.isclose(absn, 0): # n is zero for 2-fold rotation, find eigenvector with +1 val, vec = la.eigh(R) - assert np.isclose(val[-1], 1) + assert np.isclose(val[-1], 1), R n = vec[-1] # Choose direction to minimize number of minus signs n /= np.sign(np.sum(n)) -- GitLab From 5260efa9d00f56f716c970e95b38eadcd5816834 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 17:00:12 +0200 Subject: [PATCH 26/41] loosen restriction on rotation matrix types --- qsymm/groups.py | 43 ++++++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 2fbe234..530e334 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -24,8 +24,9 @@ def rmul(R1, R2): def _make_int(R): # If close to an integer array convert to integer tinyarray, else # return original array - R_int = ta.array(np.round(R), int) - if allclose(R, R_int): + R_float = (np.array(R).astype(float) if is_sympy_matrix(R) else R) + R_int = ta.array(np.round(R_float), int) + if allclose(R_float, R_int): return R_int else: return R @@ -83,11 +84,13 @@ class PointGroupElement: def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False): if isinstance(R, sympy.ImmutableMatrix): - pass + # If it is integer, recast to integer tinyarray + R = _make_int(R) elif isinstance(R, ta.ndarray_int): pass elif isinstance(R, sympy.matrices.MatrixBase): R = sympy.ImmutableMatrix(R) + R = _make_int(R) elif isinstance(R, np.ndarray): # If it is integer, recast to integer tinyarray R = _make_int(R) @@ -102,11 +105,6 @@ class PointGroupElement: return name_PG_elements(self, full=True) def __eq__(self, other): - # We do not allow mixing of PointGroupElements - # if one has a sympy spatial part R, but the other not. - if is_sympy_matrix(self.R) ^ is_sympy_matrix(other.R): - raise ValueError("Mixing of sympy with other types " - "in the spatial part R is not allowed.") # If Rs are of different type, convert it to numpy array if type(self.R) != type(other.R): Rs = np.array(self.R).astype(float) @@ -138,8 +136,8 @@ class PointGroupElement: # Sort group elements: # First by conjugate and a, then R = identity, then the rest # lexicographically - Rs = ta.array(self.R, float) - Ro = ta.array(other.R, float) + 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): @@ -162,23 +160,30 @@ class PointGroupElement: R1, c1, a1, U1 = g1.R, g1.conjugate, g1.antisymmetry, g1.U R2, c2, a2, U2 = g2.R, g2.conjugate, g2.antisymmetry, g2.U - # We do not allow mixing of PointGroupElements - # if one has a sympy spatial part R, but the other not. - if is_sympy_matrix(R1) ^ is_sympy_matrix(R2): - raise ValueError("Mixing of sympy with other types " - "in the spatial part R is not allowed.") - if (U1 is None) or (U2 is None): U = None elif c1: U = U1.dot(U2.conj()) else: U = U1.dot(U2) - # If spatial parts are sympy matrices, use cached multiplication. - if is_sympy_matrix(R1): # R1 and R2 are either both sympy or both not sympy. + + if is_sympy_matrix(R1) and is_sympy_matrix(R2): + # If spatial parts are sympy matrices, use cached multiplication. R = rmul(R1, R2) + elif not (is_sympy_matrix(R1) or is_sympy_matrix(R2)): + # If arrays, use dot + R = np.dot(R1, R2) + elif ((is_sympy_matrix(R1) or is_sympy_matrix(R2)) and + (isinstance(R1, ta.ndarray_int) or isinstance(R2, ta.ndarray_int))): + # Multiplying sympy and integer tinyarray is ok, should result in sympy + R = rmul(sympy.ImmutableMatrix(R1), sympy.ImmutableMatrix(R2)) else: - R = _make_int(np.dot(np.array(R1), np.array(R2)).astype(float)) + raise ValueError("Mixing of sympy and floating point in the spatial part R is not allowed. " + "To avoid this error, make sure that all PointGroupElements are initialized " + "with either floating point arrays or sympy matrices as rotations. " + "Integer arrays are allowed in both cases.") + # Try to round to integer + R = _make_int(R) return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq)) def __pow__(self, n): -- GitLab From 39b8ef3c0261c30b77d279e61a8c54d6568ba2e1 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Tue, 21 May 2019 17:04:06 +0200 Subject: [PATCH 27/41] modify test to the new behaviour --- qsymm/tests/test_util.py | 35 +++++++++++++++++++---------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/qsymm/tests/test_util.py b/qsymm/tests/test_util.py index 686e401..568e106 100644 --- a/qsymm/tests/test_util.py +++ b/qsymm/tests/test_util.py @@ -28,24 +28,27 @@ def test_sparse_basis(): def test_spatial_types(): - S1 = PointGroupElement(sympy.eye(2), True, True, - np.random.rand(3, 3)) + S1 = PointGroupElement(sympy.eye(2), False, False, + np.eye(3)) S2 = PointGroupElement(sympy.Matrix([[0, 1], [1, 0]]), True, False, - np.random.rand(3, 3)) + np.eye(3)) S3 = PointGroupElement(np.eye(2), False, False, - np.random.rand(3, 3)) - S4 = PointGroupElement(np.array([[0, 1.2], [1.5, 0]]), True, False, - np.random.rand(3, 3)) + 1j * np.eye(3)) + C6s = PointGroupElement(sympy.ImmutableMatrix( + [[sympy.Rational(1, 2), sympy.sqrt(3)/2], + [-sympy.sqrt(3)/2, sympy.Rational(1, 2)]] + )) + C6f = PointGroupElement(np.array( + [[1/2, np.sqrt(3)/2], + [-np.sqrt(3)/2, 1/2]] + )) - # Multiplying or comparing objects allowed if both or neither have sympy spatial parts - S = S1 * S2 + assert S2**2 == S1 assert not S1 == S2 - S = S3 * S4 - assert not S3 == S4 + assert S1 == S3 + assert C6s == C6f # Mixing sympy with other types raises an error - with raises(ValueError, message="Multiplying PointGroupElements only allowed " - "if neither or both have sympy spatial parts R."): - S = S1 * S3 - with raises(ValueError, message="Comparing PointGroupElements only allowed " - "if neither or both have sympy spatial parts R."): - S4 == S2 + with raises(ValueError, message="Multiplying PointGroupElements with sympy and " + "floating point spatial parts not allowed."): + S = C6s * C6f + -- GitLab From 80d35c3bd9c19559d5686a9af98c81e763e1381e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1niel=20Varjas?= Date: Wed, 22 May 2019 21:26:22 +0000 Subject: [PATCH 28/41] Apply suggestion to qsymm/groups.py --- qsymm/groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 530e334..b7e10b6 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -696,7 +696,7 @@ def name_PG_elements(g, full=False): - T for time-reversal (antiunitary symmetry) - P for particle-hole (antiunitary antisymmetry) - C for chiral (unitary antisymmetry) - - az_name is ommited for unitary symmetries + - missing if the symmetry is unitary """ def name_angle(theta): -- GitLab From 8b8e4d89f588e0e64c52f703620d31a3b3745d31 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 00:29:20 +0200 Subject: [PATCH 29/41] add latex display of PGE --- qsymm/groups.py | 62 ++++++++++++++++++++++++++++++++----------- symmetry_finder.ipynb | 11 ++------ 2 files changed, 48 insertions(+), 25 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index b7e10b6..888100d 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -104,6 +104,9 @@ class PointGroupElement: def __repr__(self): return name_PG_elements(self, full=True) + def _repr_latex_(self): + return name_PG_elements(self, full=False, latex=True) + def __eq__(self, other): # If Rs are of different type, convert it to numpy array if type(self.R) != type(other.R): @@ -660,7 +663,7 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): ## Human readable group element names -def name_PG_elements(g, full=False): +def name_PG_elements(g, full=False, latex=False): """ Return a human readable string representation of PointGroupElement @@ -675,6 +678,8 @@ def name_PG_elements(g, full=False): and the symbol of the Altland-Zirnbauer part of the symmetry (see below). The full representation presents the symmetry action on the Hamiltonian and the unitary Hilbert-space action if set. + latex : bool (default False) + Whether to output LateX formatted string. Returns ------- @@ -699,12 +704,22 @@ def name_PG_elements(g, full=False): - missing if the symmetry is unitary """ - def name_angle(theta): + def name_angle(theta, latex=False): frac = Fraction(theta / np.pi).limit_denominator(100) num, den = frac.numerator, frac.denominator - angle = '{}π{}'.format( - "" if num == 1 else ("-" if num == -1 else str(num)), - "" if den == 1 else ("/" + str(den))) + if latex: + if den == 1: + angle = '{}{}\pi'.format("-" if num < 0 else "", + "" if abs(num) == 1 else abs(num)) + else: + angle = '{}\\frac{{{}\pi}}{{{}}}'.format( + "-" if num < 0 else "", + "" if abs(num) == 1 else abs(num), + den) + else: + angle = '{}{}π{}'.format("-" if num < 0 else "", + "" if abs(num) == 1 else abs(num), + "" if den == 1 else ("/" + str(den))) return angle def round_axis(n): @@ -728,13 +743,19 @@ def name_PG_elements(g, full=False): if np.isclose(theta, 0): rot_name = '1' else: - rot_name = f'R({name_angle(theta)})' + if latex: + rot_name = f'R\left({name_angle(theta, latex)}\\right)' + else: + rot_name = f'R({name_angle(theta)})' else: # mirror val, vec = la.eigh(R) assert allclose(val, [-1, 1]), R n = vec[0] - rot_name = f'M({round_axis(n)})' + if latex: + rot_name = f'M\left({round_axis(n)}\\right)' + else: + rot_name = f'M({round_axis(n)})' elif R.shape[0] == 3: if np.isclose(la.det(R), 1): # pure rotation @@ -742,7 +763,10 @@ def name_PG_elements(g, full=False): if np.isclose(theta, 0): rot_name = '1' else: - rot_name = f'R({name_angle(theta)}, {round_axis(n)})' + if latex: + rot_name = f'R\left({name_angle(theta, latex)}, {round_axis(n)}\\right)' + else: + rot_name = f'R({name_angle(theta)}, {round_axis(n)})' else: # rotoinversion n, theta = rotation_to_angle(-R) @@ -751,29 +775,35 @@ def name_PG_elements(g, full=False): rot_name = 'I' elif np.isclose(theta, np.pi): # mirror - rot_name = f'M({round_axis(n)})' + if latex: + rot_name = f'M\left({round_axis(n)}\\right)' + else: + rot_name = f'M({round_axis(n)})' else: # generic rotoinversion - rot_name = f'S({name_angle(theta)}, {round_axis(n)})' + if latex: + rot_name = f'S\left({name_angle(theta, latex)}, {round_axis(n)}\\right)' + else: + rot_name = f'S({name_angle(theta)}, {round_axis(n)})' if full: - name = 'U⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", + name = '\nU⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", "-" if g.antisymmetry else "", "-" if g.conjugate else "") name += f'R = {rot_name}\n' if g.U is not None: - name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) + name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) +'\n\n' else: if g.conjugate and not g.antisymmetry: - az_name = " T" + az_name = " \mathcal{T}" if latex else " T" elif g.conjugate and g.antisymmetry: - az_name = " P" + az_name = " \mathcal{P}" if latex else " P" elif not g.conjugate and g.antisymmetry: - az_name = " C" + az_name = " \mathcal{C}" if latex else " C" else: az_name = "" name = az_name if (rot_name == '1' and az_name != "") else rot_name + az_name - return name + return '$' + name + '$' if latex else name def rotation_to_angle(R): diff --git a/symmetry_finder.ipynb b/symmetry_finder.ipynb index e65e09d..8de2af9 100644 --- a/symmetry_finder.ipynb +++ b/symmetry_finder.ipynb @@ -148,7 +148,7 @@ }, "outputs": [], "source": [ - "print([qsymm.groups.name_PG_elements(g) for g in sg])" + "[display(g) for g in sg];" ] }, { @@ -754,7 +754,7 @@ }, "outputs": [], "source": [ - "print([qsymm.groups.name_PG_elements(g) for g in sg])" + "[display(g) for g in sg];" ] }, { @@ -785,13 +785,6 @@ " %time sg, cg = qsymm.symmetries(Hkp, cubic_group)\n", " print(len(sg), len(cg))" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { -- GitLab From fceb463d5e0dff30277299e60e81eae8b7d14abe Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 14:13:08 +0200 Subject: [PATCH 30/41] rename factory functions and use them where appropriate --- qsymm/__init__.py | 2 +- qsymm/groups.py | 10 ++-- qsymm/symmetry_finder.py | 86 +++++++++++++---------------- qsymm/tests/test_symmetry_finder.py | 29 ++++------ 4 files changed, 55 insertions(+), 72 deletions(-) diff --git a/qsymm/__init__.py b/qsymm/__init__.py index 92a6f88..fd6c9ab 100644 --- a/qsymm/__init__.py +++ b/qsymm/__init__.py @@ -5,7 +5,7 @@ from . import symmetry_finder from . import model from .groups import ( PointGroupElement, ContinuousGroupGenerator, - time_reversal_operator, particle_hole_operator, chiral_operator) + time_reversal, particle_hole, chiral, inversion, rotation, mirror) from .model import Model, BlochModel from .hamiltonian_generator import continuum_hamiltonian, continuum_pairing, display_family, \ bloch_family diff --git a/qsymm/groups.py b/qsymm/groups.py index 888100d..764c94d 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -269,6 +269,7 @@ class PointGroupElement: U = None return PointGroupElement(R, False, False, U) +## Factory functions for point group elements def identity(dim, shape=None): """Return identity operator with appropriate shape. @@ -293,7 +294,7 @@ def identity(dim, shape=None): return PointGroupElement(R, False, False, U) -def time_reversal_operator(realspace_dim, U=None): +def time_reversal(realspace_dim, U=None): """Return a time-reversal symmetry operator parameters @@ -312,7 +313,7 @@ def time_reversal_operator(realspace_dim, U=None): return PointGroupElement(R, conjugate=True, antisymmetry=False, U=U) -def particle_hole_operator(realspace_dim, U=None): +def particle_hole(realspace_dim, U=None): """Return a particle-hole symmetry operator parameters @@ -331,7 +332,7 @@ def particle_hole_operator(realspace_dim, U=None): return PointGroupElement(R, conjugate=True, antisymmetry=True, U=U) -def chiral_operator(realspace_dim, U=None): +def chiral(realspace_dim, U=None): """Return a chiral symmetry operator parameters @@ -405,7 +406,7 @@ def rotation(angle, axis=None, inversion=False, U=None): else: raise ValueError('axis needs to be None or a 3D vector.') # Try to round it to integer - R = _make_int(R) + R = _make_int(R.real) return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) @@ -432,6 +433,7 @@ def mirror(axis, U=None): R = _make_int(R) return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) +## Continuous symmetry generators (conserved quantities) class ContinuousGroupGenerator: """A generator of a continuous group. diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py index 10e81bb..0916dd7 100644 --- a/qsymm/symmetry_finder.py +++ b/qsymm/symmetry_finder.py @@ -9,8 +9,8 @@ from .linalg import matrix_basis, nullspace, split_list, simult_diag, commutator allclose from .model import Model, BlochModel, BlochCoeff from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, \ - set_multiply, L_matrices, spin_rotation, time_reversal_operator, \ - particle_hole_operator, chiral_operator + set_multiply, L_matrices, spin_rotation, time_reversal, \ + particle_hole, chiral, inversion, rotation, mirror from . import kwant_continuum from .kwant_linalg_lll import lll, cvp, voronoi @@ -87,8 +87,7 @@ def symmetries(model, candidates=None, continuous_rotations=False, if candidates is None: # Make discrete onsite symmetries dim = len(model.momenta) - candidates = generate_group({PointGroupElement(np.eye(dim), True, False), - PointGroupElement(np.eye(dim), False, True)}) + candidates = generate_group({time_reversal(dim), particle_hole(dim), chiral(dim)}) if candidates: disc_sym, _ = discrete_symmetries(model, set(candidates), Ps=Ps, @@ -817,7 +816,7 @@ def bravais_point_group(periods, tr=False, ph=False, generators=False, verbose=F num_eq, sets_eq = equals(neighbors) # Everey Bravais-lattice group contains inversion - gens = {PointGroupElement(-np.eye(dim))} + gens = {inversion(dim)} if dim==1: # The only bravais lattice group in 1D only contains inversion @@ -830,10 +829,10 @@ def bravais_point_group(periods, tr=False, ph=False, generators=False, verbose=F raise NotImplementedError('Only 1, 2, and 3 dimensional translation symmetry supported.') if tr: - TR = PointGroupElement(np.eye(dim), True, False, None) + TR = time_reversal(dim) gens.add(TR) if ph: - PH = PointGroupElement(np.eye(dim), True, True, None) + PH = particle_hole(dim) gens.add(PH) assert check_bravais_symmetry(neighbors, gens) if generators: @@ -851,26 +850,26 @@ def bravais_group_2d(neighbors, num_eq, sets_eq, verbose=False): # Sqare or simple rectangular lattice name = 'simple rectangular' # Mirrors - Mx = PointGroupElement(mirror(neighbors[0])) - My = PointGroupElement(mirror(neighbors[1])) + Mx = mirror(neighbors[0]) + My = mirror(neighbors[1]) gens |= {Mx, My} if num_eq == [2]: # Square lattice, 4-fold rotation name = 'square' - C4 = PointGroupElement(np.array([[0, 1], [-1, 0]])) + C4 = rotation(1/4) gens.add(C4) elif num_eq == [1, 2] or num_eq == [3]: # Centered rectangular, mirror symmetry name = 'centered rectangular' vecs = sets_eq[-1][:2] - Mx = PointGroupElement(mirror(vecs[0] + vecs[1])) - My = PointGroupElement(mirror(vecs[0] - vecs[1])) + Mx = mirror(vecs[0] + vecs[1]) + My = mirror(vecs[0] - vecs[1]) gens.add(Mx) gens.add(My) if num_eq == [3]: name = 'hexagonal' # Hexagonal lattice, 6-fold rotation - C6 = PointGroupElement(np.array([[1/2, s3/2], [-s3/2, 1/2]])) + C6 = rotation(1/6) gens.add(C6) else: name = 'oblique' @@ -887,61 +886,61 @@ def bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False): if num_eq == [3]: # Primitive cubic, 3 4-fold axes name = 'primitive cubic' - C4s = {rotation(n, 4) for n in neighbors} + C4s = {rotation(1/4, n) for n in neighbors} gens |= C4s elif num_eq == [1, 2]: # Primitive tetragonal, find 4-fold axis name = 'primitive tetragonal' - C4 = rotation(np.cross(*sets_eq[1]), 4) + C4 = rotation(1/4, np.cross(*sets_eq[1])) gens.add(C4) - C2s = {rotation(axis, 2) for axis in sets_eq[1]} + C2s = {rotation(1/2, axis) for axis in sets_eq[1]} gens |= C2s elif num_eq == [1, 1, 1]: # Primitive orthorhombic name = 'primitive orthorhombic' - C2s = {rotation(n, 2) for n in neighbors} + C2s = {rotation(1/2, n) for n in neighbors} gens |= C2s elif num_eq == [1, 3]: # Hexagonal name = 'hexagonal' # lone one for C6 axis - C6 = rotation(sets_eq[0][0], 6) + C6 = rotation(1/6, sets_eq[0][0]) # one of the triple for C2 axis - C2 = rotation(sets_eq[1][0], 2) + C2 = rotation(1/2, sets_eq[1][0]) gens |= {C6, C2} elif num_eq == [1, 1, 2]: # Base centered orthorhombic name = 'base centered orthorhombic' vec011, vec01m1 = sets_eq[-1] - C2x = rotation(vec011 + vec01m1, 2) - C2y = rotation(vec011 - vec01m1, 2) - C2z = rotation(np.cross(vec011, vec01m1), 2) + C2x = rotation(1/2, vec011 + vec01m1) + C2y = rotation(1/2, vec011 - vec01m1) + C2z = rotation(1/2, np.cross(vec011, vec01m1)) gens |= {C2x, C2y, C2z} elif num_eq == [1, 1, 1, 1]: # Primitive Monoclinic name = 'primitive monoclinic' axis, = pick_perp(neighbors, 3) - C2 = rotation(axis, 2) + C2 = rotation(1/2, axis) gens.add(C2) elif num_eq == [3, 4]: # Body centered cubic name = 'body centered cubic' - C4s = {rotation(n, 4) for n in sets_eq[0]} + C4s = {rotation(1/4, n) for n in sets_eq[0]} gens |= C4s elif num_eq == [1, 2, 4] or num_eq == [2, 4]: # Body centered tetragonal name = 'body centered tetragonal' axes = (sets_eq[0] if len(num_eq) == 2 else sets_eq[1]) - C2s = {rotation(n, 2) for n in axes} - C4 = rotation(np.cross(*axes), 4) + C2s = {rotation(1/2, n) for n in axes} + C4 = rotation(1/4, np.cross(*axes)) gens |= C2s gens.add(C4) elif num_eq == [1, 1, 1, 4] or num_eq == [1, 1, 4]: # Body centered orthorhombic name = 'body centered orthorhombic' axes = [vec[0] for num, vec in zip(num_eq, sets_eq) if num == 1] - C2s = {rotation(n, 2) for n in axes} - C2s.add(rotation(np.cross(axes[0], axes[1]), 2)) + C2s = {rotation(1/2, n) for n in axes} + C2s.add(rotation(1/2, np.cross(axes[0], axes[1]))) gens |= C2s elif num_eq == [6]: # FCC @@ -949,9 +948,9 @@ def bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False): # pick an orthogonal pair to define C4 axes vec110 = neighbors[0] vec1m10, = pick_perp(neighbors, 1, [vec110]) - C4x = rotation(vec110 + vec1m10, 4) - C4y = rotation(vec110 - vec1m10, 4) - C4z = rotation(np.cross(vec110, vec1m10), 4) + C4x = rotation(1/4, vec110 + vec1m10) + C4y = rotation(1/4, vec110 - vec1m10) + C4z = rotation(1/4, np.cross(vec110, vec1m10)) gens |= {C4x, C4y, C4z} elif num_eq == [1, 3, 3] or num_eq == [3, 3]: # Rhombohedral with or without face toward the 3-fold axis @@ -966,11 +965,11 @@ def bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False): name = 'face centered orthorhombic' # One cubic vector has to be there vec100 = sets_eq[0][0] - C2x = rotation(vec100, 2) + C2x = rotation(1/2, vec100) # Pick the pair orthogonal to it vec011, vec01m1 = pick_perp(neighbors, 1, [vec100]) - C2y = rotation(vec011 + vec01m1, 2) - C2z = rotation(vec011 - vec01m1, 2) + C2y = rotation(1/2, vec011 + vec01m1) + C2z = rotation(1/2, vec011 - vec01m1) gens |= {C2x, C2y, C2z} elif num_eq[-1] == num_eq[-2] == 2: # Base centered monoclinic @@ -989,17 +988,6 @@ def bravais_group_3d(neighbors, num_eq, sets_eq, verbose=False): return gens -def mirror(n): - n = np.array(n) - return np.eye(len(n)) - 2 * np.outer(n, n) / (n @ n) - - -def rotation(normal, fold): - n = 2 * np.pi / fold * np.array(normal) / la.norm(normal) - Ls = L_matrices() - return PointGroupElement(np.real(spin_rotation(n, Ls))) - - def equals(vectors): # group equivalent vectors based on length and angles one = kwant_continuum.sympify('1') @@ -1045,8 +1033,8 @@ def threefold_axis(vectors, neighbors): if np.allclose([overlaps[0, 1], overlaps[0, 2], overlaps[1, 2]], 1/2): # coplanar vectors, may be 6-fold axis = np.cross(vectors[0], vectors[1]) - C3 = rotation(axis, 3) - C6 = rotation(axis, 6) + C3 = rotation(1/3, axis) + C6 = rotation(1/6, axis) if check_bravais_symmetry(neighbors, {C6}): return C6 elif check_bravais_symmetry(neighbors, {C3}): @@ -1056,7 +1044,7 @@ def threefold_axis(vectors, neighbors): for signs in it.product([1, -1], repeat=3): axis = signs @ vectors overlaps = axis @ vectors.T - C3 = rotation(axis, 3) + C3 = rotation(1/3, axis) prop, _ = prop_to_id(np.diag(np.abs(overlaps))) if prop and check_bravais_symmetry(neighbors, {C3}): return C3 @@ -1068,7 +1056,7 @@ def twofold_axis(vectors, neighbors): # Find twofold axis from vectors that leaves neighbors invariant for sign, (vec1, vec2) in it.product([1, -1], it.combinations(vectors, 2)): axis = vec1 + sign * vec2 - C2 = rotation(axis, 2) + C2 = rotation(1/2, axis) if check_bravais_symmetry(neighbors, {C2}): return C2 else: diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py index 1a15133..4c5ed4b 100644 --- a/qsymm/tests/test_symmetry_finder.py +++ b/qsymm/tests/test_symmetry_finder.py @@ -103,8 +103,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -173,8 +172,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -243,8 +241,7 @@ def test_disc_finder(verbose = False): Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) # Find symmetry operators - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) if verbose: print('sym', sg) @@ -316,8 +313,7 @@ def test_disc_finder(verbose = False): # Find the symmetry operator Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(1, 2*dim, dim), (1, 2*dim, dim)] assert len(sg) == 2 @@ -355,8 +351,7 @@ def test_disc_finder(verbose = False): assert np.allclose(H, T.dot(H.conj()).dot(T.T.conj())) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(1, 3*dim, dim), (1, 3*dim, dim), (1, 3*dim, dim)] assert len(sg) == 2 @@ -392,8 +387,7 @@ def test_disc_finder(verbose = False): assert np.allclose(H, T.dot(H.conj()).dot(T.T.conj())) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) assert [P.shape for P in Ps] == [(2, 2*dim, dim)] assert len(sg) == 2 @@ -427,8 +421,7 @@ def test_disc_finder(verbose = False): h1, t_mat.dot(h1.conj()).dot(t_mat.T.conj()))) Hs = Model({kwant_continuum.sympify('a_' + str(i)) : H for i, H in enumerate(Hs)}, momenta=[]) - sg = {time_reversal_operator(0), particle_hole_operator(0), - chiral_operator(0)} + sg = {time_reversal(0), particle_hole(0), chiral(0)} sg, Ps = discrete_symmetries(Hs, sg) 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 @@ -468,8 +461,8 @@ def test_continuum(): C3 = PointGroupElement(np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]], int), False, False, None) - TR = time_reversal_operator(realspace_dim=3) - PH = particle_hole_operator(realspace_dim=3) + TR = time_reversal(realspace_dim=3) + PH = particle_hole(realspace_dim=3) cubic_gens = {I, C4, C3, TR, PH} cubic_group = generate_group(cubic_gens) assert len(cubic_group) == 192 @@ -535,8 +528,8 @@ def test_bloch(): 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_operator(realspace_dim=2) - PH = particle_hole_operator(realspace_dim=2) + 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 -- GitLab From 790229fc807635800c9f29744fa857304f42de47 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 14:20:55 +0200 Subject: [PATCH 31/41] rename and slightly rework point group printer --- qsymm/groups.py | 86 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 56 insertions(+), 30 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 764c94d..e8c8173 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -102,10 +102,10 @@ class PointGroupElement: self._strict_eq = _strict_eq def __repr__(self): - return name_PG_elements(self, full=True) + return pretty_print_pge(self, full=True) def _repr_latex_(self): - return name_PG_elements(self, full=False, latex=True) + return pretty_print_pge(self, full=False, latex=True) def __eq__(self, other): # If Rs are of different type, convert it to numpy array @@ -665,7 +665,7 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): ## Human readable group element names -def name_PG_elements(g, full=False, latex=False): +def pretty_print_pge(g, full=False, latex=False): """ Return a human readable string representation of PointGroupElement @@ -711,10 +711,10 @@ def name_PG_elements(g, full=False, latex=False): num, den = frac.numerator, frac.denominator if latex: if den == 1: - angle = '{}{}\pi'.format("-" if num < 0 else "", + angle = r'{}{}\pi'.format("-" if num < 0 else "", "" if abs(num) == 1 else abs(num)) else: - angle = '{}\\frac{{{}\pi}}{{{}}}'.format( + angle = r'{}\frac{{{}\pi}}{{{}}}'.format( "-" if num < 0 else "", "" if abs(num) == 1 else abs(num), den) @@ -724,14 +724,6 @@ def name_PG_elements(g, full=False, latex=False): "" if den == 1 else ("/" + str(den))) return angle - def round_axis(n): - # Try to find integer axis - for vec in it.product([-1, 0, 1], repeat=len(n)): - if np.isclose(vec @ n, la.norm(vec)) and not np.isclose(la.norm(vec), 0): - return vec - # otherwise round - return np.round(n, 2) - R = np.array(g.R).astype(float) if R.shape[0] == 1: if R[0, 0] == 1: @@ -746,7 +738,7 @@ def name_PG_elements(g, full=False, latex=False): rot_name = '1' else: if latex: - rot_name = f'R\left({name_angle(theta, latex)}\\right)' + rot_name = fr'R\left({name_angle(theta, latex)}\right)' else: rot_name = f'R({name_angle(theta)})' else: @@ -755,9 +747,9 @@ def name_PG_elements(g, full=False, latex=False): assert allclose(val, [-1, 1]), R n = vec[0] if latex: - rot_name = f'M\left({round_axis(n)}\\right)' + rot_name = fr'M\left({_round_axis(n)}\right)' else: - rot_name = f'M({round_axis(n)})' + rot_name = f'M({_round_axis(n)})' elif R.shape[0] == 3: if np.isclose(la.det(R), 1): # pure rotation @@ -766,9 +758,9 @@ def name_PG_elements(g, full=False, latex=False): rot_name = '1' else: if latex: - rot_name = f'R\left({name_angle(theta, latex)}, {round_axis(n)}\\right)' + rot_name = fr'R\left({name_angle(theta, latex)}, {_round_axis(n)}\right)' else: - rot_name = f'R({name_angle(theta)}, {round_axis(n)})' + rot_name = f'R({name_angle(theta)}, {_round_axis(n)})' else: # rotoinversion n, theta = rotation_to_angle(-R) @@ -778,30 +770,39 @@ def name_PG_elements(g, full=False, latex=False): elif np.isclose(theta, np.pi): # mirror if latex: - rot_name = f'M\left({round_axis(n)}\\right)' + rot_name = fr'M\left({_round_axis(n)}\right)' else: - rot_name = f'M({round_axis(n)})' + rot_name = f'M({_round_axis(n)})' else: # generic rotoinversion if latex: - rot_name = f'S\left({name_angle(theta, latex)}, {round_axis(n)}\\right)' + rot_name = fr'S\left({name_angle(theta, latex)}, {_round_axis(n)}\right)' else: - rot_name = f'S({name_angle(theta)}, {round_axis(n)})' + rot_name = f'S({name_angle(theta)}, {_round_axis(n)})' if full: - name = '\nU⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", - "-" if g.antisymmetry else "", - "-" if g.conjugate else "") - name += f'R = {rot_name}\n' + 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 "", + "-" if g.conjugate else "") + else: + name = '\nU⋅H(k){}⋅U^-1 = {}H({}R⋅k)\n'.format("*" if g.conjugate else "", + "-" if g.antisymmetry else "", + "-" if g.conjugate else "") + name += f'R = {rot_name}' + (r'\\' if latex else '\n') if g.U is not None: - name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) +'\n\n' + if latex: + Umat = _array_to_latex(np.round(g.U, 3)) + name += 'U = {}'.format(Umat) + else: + name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) +'\n\n' else: if g.conjugate and not g.antisymmetry: - az_name = " \mathcal{T}" if latex else " T" + az_name = r" \mathcal{T}" if latex else " T" elif g.conjugate and g.antisymmetry: - az_name = " \mathcal{P}" if latex else " P" + az_name = r" \mathcal{P}" if latex else " P" elif not g.conjugate and g.antisymmetry: - az_name = " \mathcal{C}" if latex else " C" + az_name = r" \mathcal{C}" if latex else " C" else: az_name = "" name = az_name if (rot_name == '1' and az_name != "") else rot_name + az_name @@ -827,6 +828,31 @@ def rotation_to_angle(R): return n, theta +def _round_axis(n): + # Try to find integer axis + for vec in it.product([-1, 0, 1], repeat=len(n)): + if np.isclose(vec @ n, la.norm(vec)) and not np.isclose(la.norm(vec), 0): + return np.array(vec, int) + # otherwise round + return np.round(n, 2) + + +def _array_to_latex(a, precision=2): + """Returns a LaTeX bmatrix + a: numpy array + returns: LaTeX bmatrix as a string + + based on https://stackoverflow.com/questions/17129290/numpy-2d-and-1d-array-to-latex-bmatrix + """ + if len(a.shape) > 2: + raise ValueError('bmatrix can at most display two dimensions') + rv = r'\begin{bmatrix}' + for line in a: + rv += np.array2string(line, precision=precision, separator='&', suppress_small=True)[1:-1] + r'\\' + rv += r'\end{bmatrix}' + return rv + + ## Predefined representations def spin_matrices(s, include_0=False): -- GitLab From 4a1fca2749c52a56f928d829ac1d8aaa1370df12 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 15:59:14 +0200 Subject: [PATCH 32/41] define pretty printing for CGG --- qsymm/groups.py | 77 +++++++++++++++++++++++++++++++++++++ symmetry_finder.ipynb | 89 +++++++++++++++++++++++++++---------------- 2 files changed, 133 insertions(+), 33 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index e8c8173..5e16ed6 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -102,6 +102,10 @@ class PointGroupElement: self._strict_eq = _strict_eq def __repr__(self): + return ('PointGroupElement(\n{},\n{},\n{},\n{},\n)' + .format(self.R, self.conjugate, self.antisymmetry, self.U)) + + def __str__(self): return pretty_print_pge(self, full=True) def _repr_latex_(self): @@ -467,6 +471,12 @@ class ContinuousGroupGenerator: def __repr__(self): return 'ContinuousGroupGenerator(\n{},\n{},\n)'.format(self.R, self.U) + def __str__(self): + return pretty_print_cgg(self, latex=False) + + def _repr_latex_(self): + return pretty_print_cgg(self, latex=True) + def apply(self, model): """Return copy of model H(k) with applied infinitesimal generator. 1j * (H(k) U - U H(k)) + 1j * dH(k)/dk_i R_{ij} k_j @@ -809,6 +819,73 @@ def pretty_print_pge(g, full=False, latex=False): return '$' + name + '$' if latex else name +def pretty_print_cgg(g, latex=False): + """ + Return a human readable string representation of ContinuousGroupGenerator + + Parameters + ---------- + + g : ContinuousGroupGenerator + Point group element to be represented. + latex : bool (default False) + Whether to output LateX formatted string. + + Returns + ------- + name : string + The first line is the action on the Hamiltonian, the following lines + display the real space rotation as R(ϕ) for 2D rotation and + R(ϕ, axis) for 3D rotation (axis is not normalized), and the conserved + quantity as a matrix. If either of these is trivial, it is omitted. + Note that a ContinuousGroupGenerator defines a continuous group + of symmetries, so there is always a free parameter ϕ. + """ + R, L = g.R, g.U + # Find rotation axis + if R is not None and allclose(R, np.zeros(R.shape)): + R = None + if L is not None and allclose(L, np.zeros(L.shape)): + L = None + + if R is not None and R.shape[0] == 3: + n = np.array([np.trace(l @ R) for l in L_matrices()]).real + n /= la.norm(n) + n = _round_axis(n) + if latex: + rot_name = fr'R_{{\phi}} = R\left(\phi, {_round_axis(n)}\right)\\' + else: + rot_name = f'R_ϕ = R(ϕ, {_round_axis(n)})\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: + rot_name = '' + + if L is not None: + if latex: + L_name = r'L = {} \\'.format(_array_to_latex(np.round(L, 3))) + else: + L_name = 'L = {}\n'.format(str(np.round(L, 3)).replace('\n', '\n ')) + + if latex: + if R is not None: + name = r'{}H(\mathbf{{k}}){} = H(R_{{\phi}}\mathbf{{k}}) \\' + name = name.format(r'e^{-i\phi L}' if L is not None else '', + r'e^{i\phi L}' if L is not None else '') + else: + name = r'\left[ H(\mathbf{{k}}), L \right] = 0 \\' + else: + if R is not None: + 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 = r'[H(k), L] = 0' + '\n' + + name += rot_name + L_name + return '$' + name + '$' if latex else name + + def rotation_to_angle(R): # Convert 3D rotation matrix to axis and angle assert allclose(R, R.real) diff --git a/symmetry_finder.ipynb b/symmetry_finder.ipynb index 8de2af9..8f4f1cb 100644 --- a/symmetry_finder.ipynb +++ b/symmetry_finder.ipynb @@ -151,6 +151,39 @@ "[display(g) for g in sg];" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Detailed printout in LateX format." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "extensions": { + "jupyter_dashboards": { + "version": 1, + "views": { + "grid_default": { + "col": 0, + "height": 7, + "hidden": false, + "row": 14, + "width": 4 + }, + "report_default": {} + } + } + } + }, + "outputs": [], + "source": [ + "from IPython.display import display, Math\n", + "[display(Math(qsymm.groups.pretty_print_pge(g, latex=True, full=True))) for g in sg];" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -227,7 +260,7 @@ "display(H2.tosympy())\n", "sg, cg = qsymm.symmetries(H2, cubic_group)\n", "print(len(sg))\n", - "print(cg)" + "[display(g) for g in cg];" ] }, { @@ -266,7 +299,7 @@ "display(H3.tosympy())\n", "sg, cg = qsymm.symmetries(H3, cubic_group)\n", "print(len(sg))\n", - "print(cg)" + "[display(g) for g in cg];" ] }, { @@ -336,8 +369,7 @@ "outputs": [], "source": [ "sg, cg = qsymm.symmetries(H1, cubic_group)\n", - "print(len(sg))\n", - "print(cg)" + "print(len(sg))" ] }, { @@ -542,7 +574,7 @@ "sg, cg = qsymm.symmetries(H63, hex_group_2D)\n", "print(len(sg))\n", "print(set(sg) == qsymm.groups.hexagonal(ph=False))\n", - "print(cg)" + "[display(g) for g in cg];" ] }, { @@ -581,7 +613,7 @@ "sg, cg = qsymm.symmetries(H64, hex_group_2D)\n", "print(len(sg))\n", "print(set(sg) == hex_group_2D)\n", - "print(cg)" + "[display(g) for g in cg];" ] }, { @@ -637,7 +669,16 @@ "outputs": [], "source": [ "display(qsymm.sympify(ham1))\n", - "qsymm.symmetries(H1, continuous_rotations=True, prettify=True)" + "pg, cg = qsymm.symmetries(H1, continuous_rotations=True, prettify=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[display(g) for g in cg];" ] }, { @@ -729,32 +770,7 @@ "source": [ "sg, cg = qsymm.symmetries(Hkp, cubic_group)\n", "print(len(sg))\n", - "print(cg)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "extensions": { - "jupyter_dashboards": { - "version": 1, - "views": { - "grid_default": { - "col": 4, - "height": 13, - "hidden": false, - "row": 86, - "width": 4 - }, - "report_default": {} - } - } - } - }, - "outputs": [], - "source": [ - "[display(g) for g in sg];" + "[display(g) for g in cg];" ] }, { @@ -785,6 +801,13 @@ " %time sg, cg = qsymm.symmetries(Hkp, cubic_group)\n", " print(len(sg), len(cg))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { -- GitLab From a2bfe0a37473f5af1ca9a4c9c7874c582fa442b8 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 17:08:14 +0200 Subject: [PATCH 33/41] more readable __repr__ --- qsymm/groups.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 5e16ed6..6956f82 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -102,8 +102,11 @@ class PointGroupElement: self._strict_eq = _strict_eq def __repr__(self): - return ('PointGroupElement(\n{},\n{},\n{},\n{},\n)' - .format(self.R, self.conjugate, self.antisymmetry, self.U)) + return ('\nPointGroupElement(\nR = {},\nconjugate = {},\nantisymmetry = {},\nU = {})' + .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')) def __str__(self): return pretty_print_pge(self, full=True) @@ -469,7 +472,9 @@ class ContinuousGroupGenerator: self.R, self.U = R, U def __repr__(self): - return 'ContinuousGroupGenerator(\n{},\n{},\n)'.format(self.R, self.U) + return ('\nContinuousGroupGenerator(\nR = {},\nU = {})' + .format(repr(self.R).replace('\n', '\n ') if self.R is not None else 'None', + repr(self.U).replace('\n', '\n ') if self.U is not None else 'None')) def __str__(self): return pretty_print_cgg(self, latex=False) -- GitLab From 8e56cae5dd16c73a7ad27f68be32b9abaaac27ea Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 17:20:17 +0200 Subject: [PATCH 34/41] add pretty printing --- qsymm/groups.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 6956f82..712f778 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -114,6 +114,9 @@ class PointGroupElement: def _repr_latex_(self): return pretty_print_pge(self, full=False, latex=True) + def _repr_pretty_(self, pp, cycle): + pp.text(pretty_print_pge(self, full=False)) + def __eq__(self, other): # If Rs are of different type, convert it to numpy array if type(self.R) != type(other.R): @@ -482,6 +485,9 @@ class ContinuousGroupGenerator: def _repr_latex_(self): return pretty_print_cgg(self, latex=True) + def _repr_pretty_(self, pp, cycle): + pp.text(pretty_print_cgg(self)) + def apply(self, model): """Return copy of model H(k) with applied infinitesimal generator. 1j * (H(k) U - U H(k)) + 1j * dH(k)/dk_i R_{ij} k_j @@ -813,14 +819,15 @@ def pretty_print_pge(g, full=False, latex=False): name += 'U = {}'.format(str(np.round(g.U, 3)).replace('\n', '\n ')) +'\n\n' else: if g.conjugate and not g.antisymmetry: - az_name = r" \mathcal{T}" if latex else " T" + az_name = r" \mathcal{T}" if latex else "T" elif g.conjugate and g.antisymmetry: - az_name = r" \mathcal{P}" if latex else " P" + az_name = r" \mathcal{P}" if latex else "P" elif not g.conjugate and g.antisymmetry: - az_name = r" \mathcal{C}" if latex else " C" + az_name = r" \mathcal{C}" if latex else "C" else: az_name = "" - name = az_name if (rot_name == '1' and az_name != "") else rot_name + az_name + name = (az_name if (rot_name == '1' and az_name != "") + else rot_name + (" " if az_name is not "" else "") + az_name) return '$' + name + '$' if latex else name @@ -860,7 +867,7 @@ def pretty_print_cgg(g, latex=False): if latex: rot_name = fr'R_{{\phi}} = R\left(\phi, {_round_axis(n)}\right)\\' else: - rot_name = f'R_ϕ = R(ϕ, {_round_axis(n)})\n' + rot_name = f'\nR_ϕ = R(ϕ, {_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: @@ -870,7 +877,7 @@ def pretty_print_cgg(g, latex=False): if latex: L_name = r'L = {} \\'.format(_array_to_latex(np.round(L, 3))) else: - L_name = 'L = {}\n'.format(str(np.round(L, 3)).replace('\n', '\n ')) + L_name = '\nL = {}'.format(str(np.round(L, 3)).replace('\n', '\n ')) if latex: if R is not None: @@ -881,11 +888,11 @@ def pretty_print_cgg(g, latex=False): name = r'\left[ H(\mathbf{{k}}), L \right] = 0 \\' else: if R is not None: - name = r'{}H(k){} = H(R_ϕ⋅k)' + '\n' + name = '\n' + r'{}H(k){} = H(R_ϕ⋅k)' 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 = r'[H(k), L] = 0' + '\n' + name = '\n[H(k), L] = 0' name += rot_name + L_name return '$' + name + '$' if latex else name -- GitLab From c6e1e735c999e599f14005d40c83d92c099f5071 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Thu, 23 May 2019 22:19:15 +0200 Subject: [PATCH 35/41] update notebooks to use new factory fucntions and fix Model for bloch generator --- basics.ipynb | 24 ++++++++---------- bloch_generator.ipynb | 50 ++++++++++++++++++------------------ kdotp_generator.ipynb | 59 ++++++++++++------------------------------- kekule.ipynb | 20 +++++++++------ qsymm/__init__.py | 2 +- qsymm/groups.py | 12 +++------ 6 files changed, 66 insertions(+), 101 deletions(-) diff --git a/basics.ipynb b/basics.ipynb index 1bb4234..f1d605d 100644 --- a/basics.ipynb +++ b/basics.ipynb @@ -267,21 +267,17 @@ "outputs": [], "source": [ "# Identity\n", - "E = qsymm.groups.identity(3)\n", + "E = qsymm.identity(3)\n", "# Inversion\n", - "I = qsymm.PointGroupElement(-np.eye(3))\n", + "I = qsymm.inversion(3)\n", "# 4-fold rotation\n", - "C4 = qsymm.PointGroupElement(np.array([[1, 0, 0],\n", - " [0, 0, 1],\n", - " [0, -1, 0]]))\n", + "C4 = qsymm.rotation(1/4, [1, 0, 0])\n", "# 3-fold rotation\n", - "C3 = qsymm.PointGroupElement(np.array([[0, 0, 1],\n", - " [1, 0, 0],\n", - " [0, 1, 0]]))\n", + "C3 = qsymm.rotation(1/3, [1, 1, 1])\n", "# Time reversal\n", - "TR = qsymm.PointGroupElement(np.eye(3), conjugate=True)\n", + "TR = qsymm.time_reversal(3)\n", "# Particle-hole\n", - "PH = qsymm.PointGroupElement(np.eye(3), conjugate=True, antisymmetry=True)\n", + "PH = qsymm.particle_hole(3)\n", "\n", "cubic_gens = {I, C4, C3, TR, PH}" ] @@ -410,7 +406,7 @@ }, "outputs": [], "source": [ - "TR.apply(H)" + "TR.apply(H).tosympy(nsimplify=True)" ] }, { @@ -476,7 +472,7 @@ "metadata": {}, "outputs": [], "source": [ - "sz.apply(H)" + "sz.apply(H).tosympy(nsimplify=True)" ] }, { @@ -518,7 +514,7 @@ "metadata": {}, "outputs": [], "source": [ - "discrete_symm[:5]" + "discrete_symm" ] }, { @@ -592,7 +588,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.6" + "version": "3.7.1" } }, "nbformat": 4, diff --git a/bloch_generator.ipynb b/bloch_generator.ipynb index 44d73ae..a3f0843 100644 --- a/bloch_generator.ipynb +++ b/bloch_generator.ipynb @@ -45,12 +45,12 @@ "outputs": [], "source": [ "# Time reversal\n", - "TR = qsymm.PointGroupElement(sympy.eye(2), True, False, np.eye(2))\n", + "TR = qsymm.time_reversal(2, U=np.eye(2))\n", "\n", "# Chiral symmetry\n", - "C = qsymm.PointGroupElement(sympy.eye(2), False, True, np.array([[1, 0], [0, -1]]))\n", + "C = qsymm.chiral(2, U=np.array([[1, 0], [0, -1]]))\n", "\n", - "# Atom A rotates into A, B into B.\n", + "# Atom A rotates into A, B into B, use exact sympy representation\n", "sphi = 2*sympy.pi/3\n", "RC3 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", " [sympy.sin(sphi), sympy.cos(sphi)]])\n", @@ -105,7 +105,17 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Scale the bond length in terms of the graphene lattice constant, and have the function return a list of BlochModel objects." + "Scale the bond length in terms of the graphene lattice constant, and have the function return a list of BlochModel objects. For this we can use a more convenient definition of the rotation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "C3 = qsymm.rotation(1/3, U=np.eye(2))\n", + "symmetries = [C, TR, C3]" ] }, { @@ -158,10 +168,10 @@ "outputs": [], "source": [ "# Time reversal\n", - "TR = qsymm.PointGroupElement(sympy.eye(2), True, False, np.eye(3))\n", + "TR = qsymm.time_reversal(2, np.eye(3))\n", "\n", "# Mirror symmetry \n", - "Mx = qsymm.PointGroupElement(sympy.Matrix([[-1, 0], [0, 1]]), False, False, np.diag([1, -1, 1]))\n", + "Mx = qsymm.mirror([1, 0], np.diag([1, -1, 1]))\n", "\n", "# Threefold rotation on d_z^2, d_xy, d_x^2-y^2 states.\n", "C3U = np.array([[1, 0, 0],\n", @@ -175,12 +185,7 @@ "C3U2 = C3U2[mask][:, mask]\n", "assert np.allclose(C3U, C3U2)\n", "\n", - "# The spatial part of the symmetry must be provided as a symbolic matrix\n", - "# or an integer array\n", - "sphi = 2*sympy.pi/3\n", - "C3R = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", - " [sympy.sin(sphi), sympy.cos(sphi)]])\n", - "C3 = qsymm.PointGroupElement(C3R, False, False, C3U)\n", + "C3 = qsymm.rotation(1/3, U=C3U)\n", "\n", "symmetries = [TR, Mx, C3]" ] @@ -232,7 +237,7 @@ "metadata": {}, "outputs": [], "source": [ - "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs)" + "family = qsymm.bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)" ] }, { @@ -320,7 +325,7 @@ "glide = qsymm.groups.symmetry_from_permutation(np.array([[-1, 0],[0, 1]]), perm_glide, norbs, onsite_glide)\n", "\n", "# TR\n", - "time_reversal = qsymm.PointGroupElement(np.eye(2), True, False, np.eye(4))\n", + "time_reversal = qsymm.time_reversal(2, np.eye(4))\n", "\n", "gens = {glide, inversion, time_reversal}\n", "sg = qsymm.groups.generate_group(gens)\n", @@ -398,11 +403,11 @@ " perm_C4, norbs, onsite_C4)\n", "\n", "# Fermionic TR\n", - "time_reversal = qsymm.PointGroupElement(np.eye(2), True, False, \n", + "time_reversal = qsymm.time_reversal(2, \n", " np.kron(np.eye(4), qsymm.groups.spin_rotation(2*np.pi * np.array([0, 1/2, 0]), S)))\n", "\n", "# define strange PH symmetry\n", - "particle_hole = qsymm.PointGroupElement(np.eye(2), True, True, np.eye(8))\n", + "particle_hole = qsymm.particle_hole(2, np.eye(8))\n", "\n", "gens = {C4, time_reversal, particle_hole}\n", "\n", @@ -420,20 +425,13 @@ "print(len(family), [len(fam) for fam in family])\n", "qsymm.display_family(family)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "new_dev", + "display_name": "Python [conda env:kwant_conda]", "language": "python", - "name": "new_dev" + "name": "conda-env-kwant_conda-py" }, "language_info": { "codemirror_mode": { @@ -445,7 +443,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.7.1" } }, "nbformat": 4, diff --git a/kdotp_generator.ipynb b/kdotp_generator.ipynb index 6508b14..c085781 100644 --- a/kdotp_generator.ipynb +++ b/kdotp_generator.ipynb @@ -63,22 +63,16 @@ "outputs": [], "source": [ "# C3 rotational symmetry - invariant under 2pi/3\n", - "phi = 2 * np.pi / 3\n", - "C3U = qsymm.groups.spin_rotation([0, 0, phi], S)\n", - "sphi = 2*sympy.pi/3\n", - "C3R = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", - " [sympy.sin(sphi), sympy.cos(sphi)]])\n", - "C3 = qsymm.PointGroupElement(C3R, False, False, C3U)\n", + "C3U = qsymm.groups.spin_rotation([0, 0, 2 * np.pi / 3], S)\n", + "C3 = qsymm.rotation(1/3, U=C3U)\n", "\n", "# Time reversal\n", "TRU = qsymm.groups.spin_rotation([0, np.pi, 0], S)\n", - "TRR = np.eye(2)\n", - "TR = qsymm.PointGroupElement(TRR, True, False, TRU)\n", + "TR = qsymm.time_reversal(2, TRU)\n", "\n", "# Mirror symmetry in x\n", "MxU = qsymm.groups.spin_rotation([np.pi, 0, 0], S)\n", - "MxR = np.diag([-1, 1])\n", - "Mx = qsymm.PointGroupElement(MxR, False, False, MxU)\n", + "Mx = qsymm.mirror([1, 0], U=MxU)\n", "\n", "symmetries = [C3, TR, Mx]" ] @@ -176,16 +170,14 @@ " [0.0, -1.0, 0.0, 0.0],\n", " [0.0, 0.0, 1.0, 0.0],\n", " [0.0, 0.0, 0.0, -1.0]])\n", - "pR = -np.eye(2)\n", - "pS = qsymm.PointGroupElement(pR, False, False, pU)\n", + "pS = qsymm.inversion(2, U=pU)\n", "\n", "# Time reversal\n", "trU = np.array([[0.0, 0.0, -1.0, 0.0],\n", " [0.0, 0.0, 0.0, -1.0],\n", " [1.0, 0.0, 0.0, 0.0],\n", " [0.0, 1.0, 0.0, 0.0]])\n", - "trR = np.eye(2)\n", - "trS = qsymm.PointGroupElement(trR, True, False, trU)\n", + "trS = qsymm.time_reversal(2, U=trU)\n", "\n", "# Rotation\n", "phi = 2.0*np.pi/4.0 # Impose 4-fold rotational symmetry\n", @@ -193,10 +185,7 @@ " [0.0, np.exp(-1j*3*phi/2), 0.0, 0.0],\n", " [0.0, 0.0, np.exp(1j*phi/2), 0.0],\n", " [0.0, 0.0, 0.0, np.exp(1j*3*phi/2)]])\n", - "sphi = 2*sympy.pi/4\n", - "rotR = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],\n", - " [sympy.sin(sphi), sympy.cos(sphi)]])\n", - "rotS = qsymm.PointGroupElement(rotR, False, False, rotU)\n", + "rotS = qsymm.rotation(1/4, U=rotU)\n", "\n", "symmetries = [rotS, trS, pS]" ] @@ -290,20 +279,15 @@ "outputs": [], "source": [ "# Spatial inversion\n", - "pU = np.array([[1.0, 0.0, 0.0, 0.0],\n", - " [0.0, -1.0, 0.0, 0.0],\n", - " [0.0, 0.0, 1.0, 0.0],\n", - " [0.0, 0.0, 0.0, -1.0]])\n", - "pR = -np.eye(3)\n", - "pS = qsymm.PointGroupElement(pR, False, False, pU)\n", + "pU = np.diag([1, -1, 1, -1])\n", + "pS = qsymm.inversion(3, pU)\n", "\n", "# Time reversal\n", "trU = np.array([[0.0, 0.0, -1.0, 0.0],\n", " [0.0, 0.0, 0.0, -1.0],\n", " [1.0, 0.0, 0.0, 0.0],\n", " [0.0, 1.0, 0.0, 0.0]])\n", - "trR = np.eye(3)\n", - "trS = qsymm.PointGroupElement(trR, True, False, trU)\n", + "trS = qsymm.time_reversal(3, trU)\n", "\n", "# Rotation\n", "phi = 2.0*np.pi/3.0 # Impose 3-fold rotational symmetry\n", @@ -311,12 +295,7 @@ " [0.0, np.exp(1j*phi/2), 0.0, 0.0],\n", " [0.0, 0.0, np.exp(-1j*phi/2), 0.0],\n", " [0.0, 0.0, 0.0, np.exp(-1j*phi/2)]])\n", - "sphi = 2*sympy.pi/3 # Impose 3-fold rotational symmetry\n", - "rotR = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi), 0],\n", - " [sympy.sin(sphi), sympy.cos(sphi), 0],\n", - " [0, 0, 1]])\n", - "\n", - "rotS = qsymm.PointGroupElement(rotR, False, False, rotU)\n", + "rotS = qsymm.rotation(1/3, axis=[0, 0, 1], U=rotU)\n", "\n", "symmetries = [rotS, trS, pS]" ] @@ -497,25 +476,19 @@ "# C3 rotational symmetry\n", "C3U = np.kron(qsymm.groups.spin_rotation([0, 0, 2 * np.pi / 3], S), np.eye(2))\n", "sphi = 2*sympy.pi/3\n", - "C3R = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi), 0],\n", - " [sympy.sin(sphi), sympy.cos(sphi), 0],\n", - " [0, 0, 1]])\n", - "C3 = qsymm.PointGroupElement(C3R, False, False, C3U)\n", + "C3 = qsymm.rotation(1/3, axis=[0, 0, 1], U=C3U)\n", "\n", "# Time reversal\n", "TRU = np.kron(qsymm.groups.spin_rotation([0, np.pi, 0], S), np.eye(2))\n", - "TRR = np.eye(3)\n", - "TR = qsymm.PointGroupElement(TRR, True, False, TRU)\n", + "TR = qsymm.time_reversal(3, TRU)\n", "\n", "# Mirror x\n", "MxU = np.kron(qsymm.groups.spin_rotation([np.pi, 0, 0], S), np.eye(2))\n", - "MxR = np.diag([-1, 1, 1])\n", - "Mx = qsymm.PointGroupElement(MxR, False, False, MxU)\n", + "Mx = qsymm.mirror([1, 0, 0], MxU)\n", "\n", "# Inversion\n", "IU = np.kron(np.eye(2), np.diag([1, -1]))\n", - "IR = -np.eye(3)\n", - "I = qsymm.PointGroupElement(IR, False, False, IU)" + "I = qsymm.inversion(3, IU)" ] }, { @@ -669,7 +642,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.7.1" } }, "nbformat": 4, diff --git a/kekule.ipynb b/kekule.ipynb index 5550b47..86aa33c 100644 --- a/kekule.ipynb +++ b/kekule.ipynb @@ -53,8 +53,8 @@ "candidates = qsymm.groups.hexagonal()\n", "sgy, cg = qsymm.symmetries(Hkeky, candidates)\n", "print(len(sgy))\n", - "print(sorted([qsymm.groups.name_PG_elements(g) for g in sgy]))\n", - "print(cg)" + "display(sorted(sgy))\n", + "display(cg)" ] }, { @@ -84,7 +84,7 @@ "\n", "dd_sg = []\n", "for subgroup, gens in all_subgroups.items():\n", - " print(len(subgroup),' ' , end='\\r')\n", + " # print(len(subgroup))\n", " family = qsymm.continuum_hamiltonian(gens, 2, 0)\n", " family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)\n", " if family == []:\n", @@ -126,9 +126,13 @@ "source": [ "rdd_sg = []\n", "for subgroup, gens in all_subgroups.items():\n", - " print(len(subgroup),' ', end='\\r')\n", + " # print(len(subgroup),' ', end='\\r')\n", + " invs = {qsymm.inversion(2) * g for g in [qsymm.identity(2),\n", + " qsymm.time_reversal(2),\n", + " qsymm.particle_hole(2),\n", + " qsymm.chiral(2)]}\n", " # Skip subgroups that have inversion\n", - " if [g for g in subgroup if g.R == -sympy.eye(2)]:\n", + " if subgroup & invs:\n", " continue\n", " family = qsymm.continuum_hamiltonian(gens, 2, 0)\n", " family = qsymm.hamiltonian_generator.subtract_family(family, identity_terms)\n", @@ -372,10 +376,10 @@ "norbs = OrderedDict([(site, 1) for site in sites])\n", "\n", "# Time reversal\n", - "TR = qsymm.PointGroupElement(sympy.eye(2), True, False, np.eye(6))\n", + "TR = qsymm.time_reversal(2, U=np.eye(6))\n", "\n", "# Chiral symmetry\n", - "C = qsymm.PointGroupElement(sympy.eye(2), False, True, np.kron(np.eye(3), ([[1, 0], [0, -1]])))\n", + "C = qsymm.chiral(2, U=np.kron(np.eye(3), ([[1, 0], [0, -1]])))\n", "\n", "# Atom A rotates into B, B into A.\n", "sphi = 2*sympy.pi/3\n", @@ -484,7 +488,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.7.1" } }, "nbformat": 4, diff --git a/qsymm/__init__.py b/qsymm/__init__.py index fd6c9ab..a3ce9b8 100644 --- a/qsymm/__init__.py +++ b/qsymm/__init__.py @@ -4,7 +4,7 @@ from . import hamiltonian_generator from . import symmetry_finder from . import model from .groups import ( - PointGroupElement, ContinuousGroupGenerator, + PointGroupElement, ContinuousGroupGenerator, identity, time_reversal, particle_hole, chiral, inversion, rotation, mirror) from .model import Model, BlochModel from .hamiltonian_generator import continuum_hamiltonian, continuum_pairing, display_family, \ diff --git a/qsymm/groups.py b/qsymm/groups.py index 712f778..af99300 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -195,8 +195,6 @@ class PointGroupElement: "To avoid this error, make sure that all PointGroupElements are initialized " "with either floating point arrays or sympy matrices as rotations. " "Integer arrays are allowed in both cases.") - # Try to round to integer - R = _make_int(R) return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq)) def __pow__(self, n): @@ -410,14 +408,12 @@ def rotation(angle, axis=None, inversion=False, U=None): [np.sin(angle), np.cos(angle)]]) elif len(axis) == 3: # 3D - axis = np.array(axis) + axis = np.array(axis, float) R = spin_rotation(angle * axis / la.norm(axis), L_matrices(d=3, l=1)) R *= (-1 if inversion else 1) else: raise ValueError('axis needs to be None or a 3D vector.') - # Try to round it to integer - R = _make_int(R.real) - return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) + return PointGroupElement(R.real, conjugate=False, antisymmetry=False, U=U) def mirror(axis, U=None): @@ -436,11 +432,9 @@ def mirror(axis, U=None): ------- P : PointGroupElement """ - axis = np.array(axis) + axis = np.array(axis, float) axis /= la.norm(axis) R = np.eye(axis.shape[0]) - 2 * np.outer(axis, axis) - # Try to round it to integer - R = _make_int(R) return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) ## Continuous symmetry generators (conserved quantities) -- GitLab From 200e0992765459dcbe2b5942ed2fea0ceae9856f Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 24 May 2019 01:59:14 +0200 Subject: [PATCH 36/41] optimize group multiplication --- qsymm/groups.py | 53 ++++++++++++++++++++++++------------------------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index af99300..4243e8b 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -16,11 +16,27 @@ from .model import Model @ft.lru_cache(maxsize=100000) def rmul(R1, R2): - # Cached multiplication of sympy spatial parts. - # Only called if both R1 and R2 are sympy. - return R1 * R2 + # Cached multiplication of spatial parts. + if is_sympy_matrix(R1) and is_sympy_matrix(R2): + # If spatial parts are sympy matrices, use cached multiplication. + R = R1 * R2 + elif not (is_sympy_matrix(R1) or is_sympy_matrix(R2)): + # If arrays, use dot + R = ta.dot(R1, R2) + elif ((is_sympy_matrix(R1) or is_sympy_matrix(R2)) and + (isinstance(R1, ta.ndarray_int) or isinstance(R2, ta.ndarray_int))): + # Multiplying sympy and integer tinyarray is ok, should result in sympy + R = sympy.ImmutableMatrix(R1) * sympy.ImmutableMatrix(R2) + else: + raise ValueError("Mixing of sympy and floating point in the spatial part R is not allowed. " + "To avoid this error, make sure that all PointGroupElements are initialized " + "with either floating point arrays or sympy matrices as rotations. " + "Integer arrays are allowed in both cases.") + R = _make_int(R) + return R +@ft.lru_cache(maxsize=1000) def _make_int(R): # If close to an integer array convert to integer tinyarray, else # return original array @@ -88,11 +104,14 @@ class PointGroupElement: R = _make_int(R) elif isinstance(R, ta.ndarray_int): pass + elif isinstance(R, ta.ndarray_float): + R = _make_int(R) elif isinstance(R, sympy.matrices.MatrixBase): R = sympy.ImmutableMatrix(R) R = _make_int(R) elif isinstance(R, np.ndarray): # If it is integer, recast to integer tinyarray + R = ta.array(R) R = _make_int(R) else: raise ValueError('Real space rotation must be provided as a sympy matrix or an array.') @@ -125,7 +144,7 @@ class PointGroupElement: else: Rs = self.R Ro = other.R - if isinstance(Rs, np.ndarray): + if isinstance(Rs, (np.ndarray, ta.ndarray_float)): # Check equality with allclose if floating point R_eq = allclose(Rs, Ro) else: @@ -163,7 +182,7 @@ class PointGroupElement: 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, np.ndarray): + if isinstance(R, ta.ndarray_float): return hash((c, a)) else: return hash((R, c, a)) @@ -179,22 +198,7 @@ class PointGroupElement: U = U1.dot(U2.conj()) else: U = U1.dot(U2) - - if is_sympy_matrix(R1) and is_sympy_matrix(R2): - # If spatial parts are sympy matrices, use cached multiplication. - R = rmul(R1, R2) - elif not (is_sympy_matrix(R1) or is_sympy_matrix(R2)): - # If arrays, use dot - R = np.dot(R1, R2) - elif ((is_sympy_matrix(R1) or is_sympy_matrix(R2)) and - (isinstance(R1, ta.ndarray_int) or isinstance(R2, ta.ndarray_int))): - # Multiplying sympy and integer tinyarray is ok, should result in sympy - R = rmul(sympy.ImmutableMatrix(R1), sympy.ImmutableMatrix(R2)) - else: - raise ValueError("Mixing of sympy and floating point in the spatial part R is not allowed. " - "To avoid this error, make sure that all PointGroupElements are initialized " - "with either floating point arrays or sympy matrices as rotations. " - "Integer arrays are allowed in both cases.") + R = rmul(R1, R2) return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq)) def __pow__(self, n): @@ -265,12 +269,7 @@ class PointGroupElement: def identity(self): """Return identity element with the same structure as self.""" dim = self.R.shape[0] - if isinstance(self.R, sympy.matrices.MatrixBase): - R = sympy.ImmutableMatrix(sympy.eye(dim)) - elif isinstance(self.R, ta.ndarray_int): - R = ta.identity(dim, int) - else: - R = np.eye(dim, dtype=float) + R = ta.identity(dim, int) if self.U is not None: U = np.eye(self.U.shape[0]) else: -- GitLab From 1f09eabf86a47fc179440a4534b92aaf0c3665fe Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 24 May 2019 10:24:32 +0200 Subject: [PATCH 37/41] further optimization of PGE --- qsymm/groups.py | 76 +++++++++++++++++++++++++++++-------------------- 1 file changed, 45 insertions(+), 31 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 4243e8b..c8dede9 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -14,8 +14,10 @@ from .linalg import prop_to_id, _inv_int, allclose from .model import Model -@ft.lru_cache(maxsize=100000) -def rmul(R1, R2): +# Chache the operations that are potentially slow and happen a lot in +# group theory applications. Essentially store the multiplication table. +@ft.lru_cache(maxsize=10000) +def _mul(R1, R2): # Cached multiplication of spatial parts. if is_sympy_matrix(R1) and is_sympy_matrix(R2): # If spatial parts are sympy matrices, use cached multiplication. @@ -36,6 +38,33 @@ def rmul(R1, R2): return R +@ft.lru_cache(maxsize=1000) +def _inv(R): + if isinstance(R, ta.ndarray_int): + Rinv = _inv_int(R) + elif isinstance(R, ta.ndarray_float): + Rinv = la.inv(R) + elif is_sympy_matrix(R): + Rinv = R**(-1) + else: + raise ValueError('Invalid rotation matrix.') + return Rinv + + +@ft.lru_cache(maxsize=10000) +def _eq(R1, R2): + if type(R1) != type(R2): + R1 = ta.array(np.array(R1).astype(float), float) + R2 = ta.array(np.array(R2).astype(float), float) + if isinstance(R1, ta.ndarray_float): + # Check equality with allclose if floating point + R_eq = allclose(R1, R2) + else: + # If exact use exact equality + R_eq = (R1 == R2) + return R_eq + + @ft.lru_cache(maxsize=1000) def _make_int(R): # If close to an integer array convert to integer tinyarray, else @@ -88,15 +117,22 @@ class PointGroupElement: Notes ----- - As U is floating point and has a phase ambiguity at least, - it is ignored when comparing objects. + As U is floating point and has a phase ambiguity at least, hence + it is ignored when comparing objects by default. R is the real space rotation acion. Do not include minus sign for the k-space action of antiunitary operators, such as time reversal. This minus sign will be included automatically if 'conjugate=True'. + + For most uses R can be provided as a floating point array. It is + necessary to use exact sympy matrix representation if the PGE has + to act on Models with complicated momentum dependence (not polynomial), + as the function parts of models are compared exactly. If the momentum + dependence is periodic (sine, cosine and exponential), use BlochModel, + this works with floating point rotations. """ - __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_Rinv', '_strict_eq') + __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq') def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False): if isinstance(R, sympy.ImmutableMatrix): @@ -117,7 +153,6 @@ class PointGroupElement: raise ValueError('Real space rotation must be provided as a sympy matrix or an array.') self.R, self.conjugate, self.antisymmetry, self.U = R, conjugate, antisymmetry, U # Calculating sympy inverse is slow, remember it - self._Rinv = None self._strict_eq = _strict_eq def __repr__(self): @@ -137,19 +172,7 @@ class PointGroupElement: pp.text(pretty_print_pge(self, full=False)) def __eq__(self, other): - # If Rs are of different type, convert it to numpy array - if type(self.R) != type(other.R): - Rs = np.array(self.R).astype(float) - Ro = np.array(other.R).astype(float) - else: - Rs = self.R - Ro = other.R - if isinstance(Rs, (np.ndarray, ta.ndarray_float)): - # Check equality with allclose if floating point - R_eq = allclose(Rs, Ro) - else: - # If exact use exact equality - R_eq = (Rs == Ro) + R_eq = _eq(self.R, other.R) basic_eq = R_eq and ((self.conjugate, self.antisymmetry) == (other.conjugate, other.antisymmetry)) # Equality is only checked for U if _strict_eq is True and basic_eq is True. if basic_eq and (self._strict_eq is True or other._strict_eq is True): @@ -198,7 +221,7 @@ class PointGroupElement: U = U1.dot(U2.conj()) else: U = U1.dot(U2) - R = rmul(R1, R2) + R = _mul(R1, R2) return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq)) def __pow__(self, n): @@ -218,17 +241,8 @@ class PointGroupElement: else: Uinv = la.inv(U) # Check if inverse is stored, if not, calculate it - if self._Rinv is None: - if isinstance(R, sympy.matrices.MatrixBase): - self._Rinv = R**(-1) - elif isinstance(R, ta.ndarray_int): - self._Rinv = _inv_int(R) - elif isinstance(R, np.ndarray): - self._Rinv = la.inv(R) - else: # This is probably unnecessary - raise ValueError('Illegal datatype for the spatial part') - result = PointGroupElement(self._Rinv, c, a, Uinv, _strict_eq=self._strict_eq) - result._Rinv = R + Rinv = _inv(R) + result = PointGroupElement(Rinv, c, a, Uinv, _strict_eq=self._strict_eq) return result def _strictereq(self, other): -- GitLab From 1fdf88ab86620b6811518abe522d1867a45ceefc Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 24 May 2019 11:54:05 +0200 Subject: [PATCH 38/41] add error if PGE format is misused --- qsymm/groups.py | 13 ++----------- qsymm/hamiltonian_generator.py | 6 ++++++ 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index c8dede9..e6a3744 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -87,11 +87,7 @@ def _is_antisymmetric(a): def is_sympy_matrix(R): # Returns True if the input is a sympy.Matrix or sympy.ImmutableMatrix. - types = [sympy.ImmutableMatrix, sympy.matrices.MatrixBase] - if any([isinstance(R, t) for t in types]): - return True - else: - return False + return isinstance(R, (sympy.ImmutableMatrix, sympy.matrices.MatrixBase)) class PointGroupElement: @@ -263,12 +259,7 @@ class PointGroupElement: If self.U is None, U is taken as the identity. """ R, antiunitary, antisymmetry, U = self.R, self.conjugate, self.antisymmetry, self.U - if isinstance(R, sympy.matrices.MatrixBase): - R = R**(-1) - elif isinstance(R, ta.ndarray_int): - R = _inv_int(R) - else: - R = la.inv(R) + R = _inv(R) R *= (-1 if antiunitary else 1) result = model.rotate_momenta(R) if antiunitary: diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py index 9f803bf..974535e 100644 --- a/qsymm/hamiltonian_generator.py +++ b/qsymm/hamiltonian_generator.py @@ -4,6 +4,7 @@ import sympy import itertools as it import scipy.linalg as la 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, _commutative_momenta, e, I @@ -638,6 +639,11 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, conserved = [g for g in symmetries if isinstance(g, ContinuousGroupGenerator)] if not all([(g.R is None or np.allclose(g.R, np.zeros_like(g.R))) for g in conserved]): raise ValueError('Bloch Hamiltonian cannot have continuous rotation symmetry.') + if (not bloch_model) and any(isinstance(g.R, ta.ndarray_float) for g in pg): + raise ValueError('Cannot generate Bloch Hamiltonian in Model format using ' + 'floating point rotation matrices. To avoid this error, use ' + 'only PointGroupElements that are defined with exact sympy ' + 'rotation matrices or use bloch_model=True.') # Check dimensionality dim = len(hopping_vectors[0][-1]) -- GitLab From c573c3c69079a6dfbdf1c85add045f4d8d6589e5 Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Fri, 24 May 2019 15:08:16 +0200 Subject: [PATCH 39/41] allow specification of spin in factory functions. Add more docstrings --- kdotp_generator.ipynb | 40 ++---- qsymm/groups.py | 320 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 274 insertions(+), 86 deletions(-) diff --git a/kdotp_generator.ipynb b/kdotp_generator.ipynb index c085781..50ced21 100644 --- a/kdotp_generator.ipynb +++ b/kdotp_generator.ipynb @@ -26,18 +26,6 @@ "import qsymm" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "# Spin matrices\n", - "S = qsymm.groups.spin_matrices(1/2)" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -63,16 +51,13 @@ "outputs": [], "source": [ "# C3 rotational symmetry - invariant under 2pi/3\n", - "C3U = qsymm.groups.spin_rotation([0, 0, 2 * np.pi / 3], S)\n", - "C3 = qsymm.rotation(1/3, U=C3U)\n", + "C3 = qsymm.rotation(1/3, spin=1/2)\n", "\n", "# Time reversal\n", - "TRU = qsymm.groups.spin_rotation([0, np.pi, 0], S)\n", - "TR = qsymm.time_reversal(2, TRU)\n", + "TR = qsymm.time_reversal(2, spin=1/2)\n", "\n", "# Mirror symmetry in x\n", - "MxU = qsymm.groups.spin_rotation([np.pi, 0, 0], S)\n", - "Mx = qsymm.mirror([1, 0], U=MxU)\n", + "Mx = qsymm.mirror([1, 0], spin=1/2)\n", "\n", "symmetries = [C3, TR, Mx]" ] @@ -473,18 +458,16 @@ }, "outputs": [], "source": [ + "# Double spin-1/2 representation\n", + "spin = [np.kron(s, np.eye(2)) for s in qsymm.groups.spin_matrices(1/2)]\n", "# C3 rotational symmetry\n", - "C3U = np.kron(qsymm.groups.spin_rotation([0, 0, 2 * np.pi / 3], S), np.eye(2))\n", - "sphi = 2*sympy.pi/3\n", - "C3 = qsymm.rotation(1/3, axis=[0, 0, 1], U=C3U)\n", + "C3 = qsymm.rotation(1/3, axis=[0, 0, 1], spin=spin)\n", "\n", "# Time reversal\n", - "TRU = np.kron(qsymm.groups.spin_rotation([0, np.pi, 0], S), np.eye(2))\n", - "TR = qsymm.time_reversal(3, TRU)\n", + "TR = qsymm.time_reversal(3, spin=spin)\n", "\n", "# Mirror x\n", - "MxU = np.kron(qsymm.groups.spin_rotation([np.pi, 0, 0], S), np.eye(2))\n", - "Mx = qsymm.mirror([1, 0, 0], MxU)\n", + "Mx = qsymm.mirror([1, 0, 0], spin=spin)\n", "\n", "# Inversion\n", "IU = np.kron(np.eye(2), np.diag([1, -1]))\n", @@ -624,6 +607,13 @@ "qsymm.display_family(new_terms_2)\n", "print(len(new_terms_2))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/qsymm/groups.py b/qsymm/groups.py index e6a3744..83b18d6 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -6,6 +6,7 @@ import scipy.linalg as la import itertools as it import functools as ft from fractions import Fraction +from numbers import Number from collections import OrderedDict import sympy from copy import deepcopy @@ -306,7 +307,7 @@ def identity(dim, shape=None): return PointGroupElement(R, False, False, U) -def time_reversal(realspace_dim, U=None): +def time_reversal(realspace_dim, U=None, spin=None): """Return a time-reversal symmetry operator parameters @@ -315,12 +316,24 @@ def time_reversal(realspace_dim, U=None): Realspace dimension 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. + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the time reversal + operator. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of time-reversal operator is U = exp(-i π s_y). Only one of U and spin + may be provided. Returns ------- T : PointGroupElement """ + 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) @@ -382,7 +395,7 @@ def inversion(realspace_dim, U=None): return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) -def rotation(angle, axis=None, inversion=False, U=None): +def rotation(angle, axis=None, inversion=False, U=None, spin=None): """Return a rotation operator parameters @@ -400,45 +413,83 @@ def rotation(angle, axis=None, inversion=False, U=None): U: ndarray (optional) The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the + operator. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be + the rotation axis. Only one of U and spin may be provided. Returns ------- P : PointGroupElement """ + if U is not None and spin is not None: + raise ValueError('Only one of U and spin may be provided.') + angle = 2 * np.pi * angle if axis is None: # 2D R = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) + if spin is not None: + U = spin_rotation(angle * np.array([0, 0, 1]), spin) elif len(axis) == 3: # 3D - axis = np.array(axis, float) - R = spin_rotation(angle * axis / la.norm(axis), L_matrices(d=3, l=1)) + n = angle * np.array(axis, float) / la.norm(axis) + R = spin_rotation(n, L_matrices(d=3, l=1)) R *= (-1 if inversion else 1) + if spin is not None: + U = spin_rotation(n, spin) else: raise ValueError('axis needs to be None or a 3D vector.') return PointGroupElement(R.real, conjugate=False, antisymmetry=False, U=U) -def mirror(axis, U=None): +def mirror(axis, U=None, spin=None): """Return a mirror operator - parameters - ---------- + Parameters: + ----------- axis : ndarray Normal of the mirror. The dimensionality of the operator is the same as the length of axis. U: ndarray (optional) The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the + operator. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of mirror operator is U = exp(-i π n⋅s) where n is normalized to 1. In 2D the + axis is treated as x and y coordinates. Only one of U and spin may be provided. - Returns - ------- + Returns: + -------- P : PointGroupElement + + Notes: + ------ + Warning: in 2D the real space action of a mirror and and a 2-fold rotation + around an axis in the plane is identical, however the action on angular momentum + is different. Here we consider the action of the mirror, which is the same as the + action of a 2-fold rotation around the mirror axis. """ + if U is not None and spin is not None: + raise ValueError('Only one of U and spin may be provided.') + axis = np.array(axis, float) axis /= la.norm(axis) R = np.eye(axis.shape[0]) - 2 * np.outer(axis, axis) + if spin is not None: + if len(axis) == 2: + axis = np.append(axis, 0) + U = spin_rotation(np.pi * axis, spin) + return PointGroupElement(R, conjugate=False, antisymmetry=False, U=U) ## Continuous symmetry generators (conserved quantities) @@ -599,8 +650,63 @@ def generate_subgroups(group): ## Predefined point groups -def cubic(tr=True, ph=True, generators=False): - """Generate cubic point group in standard basis. +def square(tr=True, ph=True, generators=False, spin=None): + """ + Generate square point group in standard basis. + + Parameters: + ----------- + tr, ph : bool (default True) + Whether to include time-reversal and particle-hole + symmetry. + generators : bool (default false) + Only return the group generators if True. + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the + operator. If not provided, the PointGroupElements have the unitary + action set to None. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be + the rotation axis. If ph is True, spin may not be provided, as it is not + possible to deduce the unitary representation of particle-hole symmetry from + spin alone. In this case construct the particle-hole operator manually. + + Returns: + -------- + set of PointGroupElement objects with integer rotations + + Notes: + ------ + Warning: in 2D the real space action of a mirror and and a 2-fold rotation + around an axis in the plane is identical, however the action on angular momentum + is different. Here we consider the action of the mirror, which is the same as the + action of a 2-fold rotation around the mirror axis, assuming inversion acts trivially. + """ + if ph and spin is not None: + raise ValueError('If ph is True, spin may not be provided, as it is not ' + 'possible to deduce the unitary representation of particle-hole symmetry ' + 'from spin alone. In this case construct the particle-hole operator manually.') + Mx = mirror([1, 0], spin=spin) + C4 = rotation(1/4, spin=spin) + gens = {Mx, C4} + if tr: + TR = time_reversal(2, spin=spin) + gens.add(TR) + if ph: + PH = particle_hole(2) + gens.add(PH) + if generators: + return gens + else: + return generate_group(gens) + + +def cubic(tr=True, ph=True, generators=False, spin=None): + """ + Generate cubic point group in standard basis. + Parameters: ----------- tr, ph : bool (default True) @@ -608,25 +714,39 @@ def cubic(tr=True, ph=True, generators=False): symmetry. generators : bool (default false) Only return the group generators if True. + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the + operator. If not provided, the PointGroupElements have the unitary + action set to None. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of rotation operator is U = exp(-i n⋅s). If ph is True, spin may not be + provided, as it is not possible to deduce the unitary representation of + particle-hole symmetry from spin alone. In this case construct the + particle-hole operator manually. Returns: -------- - set of PointGroupElement objects with integer rotations""" - eye = ta.array(np.eye(3), int) - E = PointGroupElement(eye, False, False, None) - I = PointGroupElement(-eye, False, False, None) - C4 = PointGroupElement(ta.array([[1, 0, 0], - [0, 0, 1], - [0, -1, 0]], int), False, False, None) - C3 = PointGroupElement(ta.array([[0, 0, 1], - [1, 0, 0], - [0, 1, 0]], int), False, False, None) + set of PointGroupElement objects with integer rotations + + Notes: + ------ + We assume inversion acts trivially in spin space. + """ + if ph and spin is not None: + raise ValueError('If ph is True, spin may not be provided, as it is not ' + 'possible to deduce the unitary representation of particle-hole symmetry ' + 'from spin alone. In this case construct the particle-hole operator manually.') + I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin))) + C4 = rotation(1/4, [1, 0, 0], spin=spin) + C3 = rotation(1/3, [1, 1, 1], spin=spin) cubic_gens = {I, C4, C3} if tr: - TR = PointGroupElement(eye, True, False, None) + TR = time_reversal(3, spin=spin) cubic_gens.add(TR) if ph: - PH = PointGroupElement(eye, True, True, None) + PH = particle_hole(3) cubic_gens.add(PH) if generators: return cubic_gens @@ -634,10 +754,16 @@ def cubic(tr=True, ph=True, generators=False): return generate_group(cubic_gens) -def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): - """Generate hexagonal point group in standard basis. +def hexagonal(dim=2, tr=True, ph=True, generators=False, sympy_R=True, spin=None): + """ + Generate hexagonal point group in standard basis in 2 or 3 dimensions. + Mirror symmetries with the main coordinate axes as normals are included. + In 3D the hexagonal axis is the z axis. + Parameters: ----------- + dim : int (default 2) + Real sapce dimensionality, 2 or 3. tr, ph : bool (default True) Whether to include time-reversal and particle-hole symmetry. @@ -646,40 +772,71 @@ def hexagonal(tr=True, ph=True, generators=False, sympy_R=True): sympy_R: bool (default True) Whether the rotation matrices should be exact sympy representations. + spin : float or sequence of arrays (optional) + Spin representation to use for the unitary action of the + operator. If not provided, the PointGroupElements have the unitary + action set to None. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. The unitary action + of rotation operator is U = exp(-i n⋅s). In 2D the z axis is assumed to be + the rotation axis. If ph is True, spin may not be provided, as it is not + possible to deduce the unitary representation of particle-hole symmetry from + spin alone. In this case construct the particle-hole operator manually. Returns: -------- - set of PointGroupElement objects with Sympy rotations.""" - - if sympy_R: - eye = 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) + set of PointGroupElements + + Notes: + ------ + Warning: in 2D the real space action of a mirror and and a 2-fold rotation + around an axis in the plane is identical, however the action on angular momentum + is different. Here we consider the action of the mirror, which is the same as the + action of a 2-fold rotation around the mirror axis, assuming inversion acts trivially. + """ + if spin is not None and sympy_R: + U6 = spin_rotation(np.pi / 3 * np.array([0, 0, 1]), spin) else: - eye = np.eye(2) - Mx = PointGroupElement(np.diag([-1, 1]), False, False, None) - C6 = PointGroupElement(np.array( - [[1/2, np.sqrt(3)/2], - [-np.sqrt(3)/2, 1/2]] - ), - False, False, None) - gens_hex_2D ={Mx, C6} + U6 = None + if dim == 2: + Mx = mirror([1, 0], spin=spin) + if sympy_R: + + C6 = PointGroupElement(sympy.ImmutableMatrix( + [[sympy.Rational(1, 2), sympy.sqrt(3)/2], + [-sympy.sqrt(3)/2, sympy.Rational(1, 2)]] + ), + False, False, U6) + else: + C6 = rotation(1/6, spin=spin) + gens = {Mx, C6} + elif dim == 3: + I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin))) + C2x = rotation(1/2, [1, 0, 0], spin=spin) + if sympy_R: + C6 = PointGroupElement(sympy.ImmutableMatrix( + [[sympy.Rational(1, 2), sympy.sqrt(3)/2, 0], + [-sympy.sqrt(3)/2, sympy.Rational(1, 2), 0], + [0, 0, 1]] + ), + False, False, U6) + else: + C6 = rotation(1/6, [0, 0, 1], spin=spin) + gens = {I, C2x, C6} + else: + raise ValueError('Only 2 and 3 dimensions are supported.') + if tr: - TR = PointGroupElement(eye, True, False, None) - gens_hex_2D.add(TR) + TR = time_reversal(dim, spin=spin) + gens.add(TR) if ph: - PH = PointGroupElement(eye, True, True, None) - gens_hex_2D.add(PH) + PH = particle_hole(dim) + gens.add(PH) if generators: - return gens_hex_2D + return gens else: - return generate_group(gens_hex_2D) + return generate_group(gens) ## Human readable group element names @@ -943,9 +1100,25 @@ def _array_to_latex(a, precision=2): ## Predefined representations def spin_matrices(s, include_0=False): - """Construct spin-s matrices for any half-integer spin. - If include_0 is True, S[0] is the identity, indices 1, 2, 3 - correspond to x, y, z. Otherwise indices 0, 1, 2 are x, y, z. + """ + Construct spin-s matrices for any half-integer spin. + + Parameters: + ----------- + + s : float or int + Spin representation to use, must be integer or half-integer. + include_0 : bool (default False) + If include_0 is True, S[0] is the identity, indices 1, 2, 3 + correspond to x, y, z. Otherwise indices 0, 1, 2 are x, y, z. + + Returns: + -------- + + ndarray + Sequence of spin-s operators in the standard spin-z basis. + Array of shape (3, 2*s + 1, 2*s + 1), or if include_0 is True + (4, 2*s + 1, 2*s + 1). """ d = np.round(2*s + 1) assert np.isclose(d, int(d)) @@ -961,14 +1134,39 @@ def spin_matrices(s, include_0=False): return np.array([Sx, Sy, Sz]) -def spin_rotation(n, S, roundint=False): - """Construct the unitary spin rotation matrix for rotation generators S - with rotation specified by the vector n. It is assumed that - n and S have the same first dimension. - If roundint is True, result is converted to integer tinyarray if possible. +def spin_rotation(n, s, roundint=False): """ + Construct the unitary spin rotation matrix for rotation specified by the + vector n (in radian units) with angular momentum s, given by + U = exp(-i n⋅s). + + Parameters: + ----------- + + n : iterable + Rotation vector. Its norm is the rotation angle in radians. + s : float or sequence of arrays + Spin representation to use for the unitary action of the + operator. If float is provided, it should be integer or half-integer + specifying the spin representation in the standard basis, see spin_matrices. + Otherwise a sequence of 3 arrays of identical square size must be provided + representing 3 components of the angular momentum operator. + roundint : bool (default False) + If roundint is True, result is converted to integer tinyarray if possible. + + Returns: + -------- + U : ndarray + Unitary spin rotation matrix of the same shape as the spin matrices or + (2*s + 1, 2*s + 1). + """ + n = np.array(n) + if isinstance(s, Number): + s = spin_matrices(s) + else: + s = np.array(s) # Make matrix exponential for rotation representation - U = la.expm(-1j * np.tensordot(n, S , axes=((0), (0)))) + U = la.expm(-1j * np.tensordot(n, s, axes=((0), (0)))) if roundint: Ur = np.round(np.real(U)) assert allclose(U, Ur) -- GitLab From ca4a5fc6dc8d954aad68b6c02e50fbb6915f22ee Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Wed, 29 May 2019 13:38:19 +0200 Subject: [PATCH 40/41] fix rotation axis finding --- qsymm/groups.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index 83b18d6..c792c8a 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -921,7 +921,7 @@ def pretty_print_pge(g, full=False, latex=False): # mirror val, vec = la.eigh(R) assert allclose(val, [-1, 1]), R - n = vec[0] + n = vec[:, 0] if latex: rot_name = fr'M\left({_round_axis(n)}\right)' else: @@ -1063,7 +1063,7 @@ def rotation_to_angle(R): # n is zero for 2-fold rotation, find eigenvector with +1 val, vec = la.eigh(R) assert np.isclose(val[-1], 1), R - n = vec[-1] + n = vec[:, -1] # Choose direction to minimize number of minus signs n /= np.sign(np.sum(n)) else: -- GitLab From f020bca9b9bf610e9379a61eda9aebcf187d228b Mon Sep 17 00:00:00 2001 From: Daniel Varjas Date: Sun, 2 Jun 2019 22:58:05 +0200 Subject: [PATCH 41/41] start work on point group representations --- qsymm/groups.py | 37 ++-- qsymm/hamiltonian_generator.py | 4 +- symmetry_finder.ipynb | 388 ++++++++++++++++++++++++++++++++- 3 files changed, 410 insertions(+), 19 deletions(-) diff --git a/qsymm/groups.py b/qsymm/groups.py index c792c8a..19f0f16 100644 --- a/qsymm/groups.py +++ b/qsymm/groups.py @@ -91,6 +91,17 @@ def is_sympy_matrix(R): return isinstance(R, (sympy.ImmutableMatrix, sympy.matrices.MatrixBase)) +@ft.lru_cache(maxsize=10000) +def _U_phase_eq(U1, U2): + prop, coeff = prop_to_id(np.dot(la.inv(U1), U2)) + return (prop and np.isclose(abs(coeff), 1)) + + +@ft.lru_cache(maxsize=10000) +def _U_strict_eq(U1, U2): + return allclose(U1, U2) + + class PointGroupElement: """An element of a point group. @@ -106,11 +117,10 @@ class PointGroupElement: U : ndarray (optional) The unitary action on the Hilbert space. May be None, to be able to treat symmetry candidates - _strict_eq : boolean (default False) - Whether to test the equality of the unitary parts when comparing with + _U_eq : callable or None (default) + Function to test the equality of the unitary parts when comparing with other PointGroupElements. By default the unitary parts are ignored. - If True, PointGroupElements are considered equal, if the unitary parts - are proportional, an overall phase difference is still allowed. + Should take two tinyarrays and return a boolean. Notes ----- @@ -129,9 +139,9 @@ class PointGroupElement: this works with floating point rotations. """ - __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_strict_eq') + __slots__ = ('R', 'conjugate', 'antisymmetry', 'U', '_U_eq') - def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _strict_eq=False): + def __init__(self, R, conjugate=False, antisymmetry=False, U=None, _U_eq=None): if isinstance(R, sympy.ImmutableMatrix): # If it is integer, recast to integer tinyarray R = _make_int(R) @@ -149,8 +159,8 @@ class PointGroupElement: else: raise ValueError('Real space rotation must be provided as a sympy matrix or an array.') self.R, self.conjugate, self.antisymmetry, self.U = R, conjugate, antisymmetry, U - # Calculating sympy inverse is slow, remember it - self._strict_eq = _strict_eq + self._U_eq = _U_eq + def __repr__(self): return ('\nPointGroupElement(\nR = {},\nconjugate = {},\nantisymmetry = {},\nU = {})' @@ -171,15 +181,14 @@ class PointGroupElement: def __eq__(self, other): R_eq = _eq(self.R, other.R) basic_eq = R_eq and ((self.conjugate, self.antisymmetry) == (other.conjugate, other.antisymmetry)) - # Equality is only checked for U if _strict_eq is True and basic_eq is True. - if basic_eq and (self._strict_eq is True or other._strict_eq is True): + # Equality is only checked for U if basic_eq is True. + if basic_eq and self._U_eq is not None: if (self.U is None) and (other.U is None): U_eq = True elif (self.U is None) ^ (other.U is None): U_eq = False else: - prop, coeff = prop_to_id((self.inv() * other).U) - U_eq = (prop and np.isclose(abs(coeff), 1)) + U_eq = self._U_eq(ta.array(self.U), ta.array(other.U)) else: U_eq = True return basic_eq and U_eq @@ -219,7 +228,7 @@ class PointGroupElement: else: U = U1.dot(U2) R = _mul(R1, R2) - return PointGroupElement(R, c1^c2, a1^a2, U, _strict_eq=(self._strict_eq or g2._strict_eq)) + return PointGroupElement(R, c1^c2, a1^a2, U, _U_eq=self._U_eq) def __pow__(self, n): result = self.identity() @@ -239,7 +248,7 @@ class PointGroupElement: Uinv = la.inv(U) # Check if inverse is stored, if not, calculate it Rinv = _inv(R) - result = PointGroupElement(Rinv, c, a, Uinv, _strict_eq=self._strict_eq) + result = PointGroupElement(Rinv, c, a, Uinv, _U_eq=self._U_eq) return result def _strictereq(self, other): diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py index 974535e..16e8b30 100644 --- a/qsymm/hamiltonian_generator.py +++ b/qsymm/hamiltonian_generator.py @@ -8,7 +8,7 @@ import tinyarray as ta from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref, allclose from .model import Model, BlochModel, _commutative_momenta, e, I -from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group +from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, _U_phase_eq from . import kwant_continuum @@ -683,7 +683,7 @@ def bloch_family(hopping_vectors, symmetries, norbs, onsites=True, if symmetrize: # Make sure that group is generated while keeping track of unitary part. for g in pg: - g._strict_eq = True + g._U_eq = _U_phase_eq pg = generate_group(set(pg)) # Symmetrize every term and remove linearly dependent or zero ones family2 = [] diff --git a/symmetry_finder.ipynb b/symmetry_finder.ipynb index 8f4f1cb..501cc87 100644 --- a/symmetry_finder.ipynb +++ b/symmetry_finder.ipynb @@ -802,6 +802,388 @@ " print(len(sg), len(cg))" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import itertools as it\n", + "import scipy.linalg as la\n", + "from qsymm.linalg import prop_to_id, allclose\n", + "from qsymm.groups import PointGroupElement, set_multiply, generate_group, _U_strict_eq" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def identity_power(g, Ps=None, double_group=False):\n", + " # Append on-site U(1) symmetries to PG element to ensure that \n", + " # n'th power is exactly the identity.\n", + " # If double_group=True, if the n'th power is 2\\pi rotation,\n", + " # -identity is chosen\n", + " def find_order(r):\n", + " dim = r.shape[0]\n", + " rpower = r\n", + " order = 1\n", + " while not np.allclose(rpower, np.eye(dim)):\n", + " rpower = r.dot(rpower)\n", + " order += 1\n", + " return order\n", + " \n", + " R = np.array(g.R).astype(float)\n", + " n = find_order(R)\n", + " sign = 1\n", + " if double_group and n>1:\n", + " if np.isclose(la.det(R), 1):\n", + " # If not rotoinversion, its n'th power is 2pi rotation\n", + " sign = -1\n", + " elif find_order(-R) == n:\n", + " # here we assume R is 3d rotoinversion\n", + " # R^n is 2pi rotation\n", + " sign = -1\n", + " # if -R has a lower order (half of n) then \n", + " # R^n is a 4pi rotation or identity\n", + "\n", + " Un = np.linalg.matrix_power(g.U, n)\n", + "\n", + " if Ps is not None:\n", + " t = np.zeros_like(g.U)\n", + " for P in Ps:\n", + " # reduce Un to one block of the block\n", + " Unr = P[0].T.conj() @ Un @ P[0]\n", + " prop_id, coeff = prop_to_id(Unr)\n", + " assert prop_id\n", + " # transform U by n'th root\n", + " t += coeff ** (-1/n) * np.einsum('aij,akj->ik', P, P.conjugate())\n", + " else:\n", + " prop_id, coeff = prop_to_id(Un)\n", + " assert prop_id, (g.U, Un)\n", + " t = coeff ** (-1/n) * np.eye(Un.shape[0])\n", + "\n", + " t = t * sign ** (-1/n)\n", + " g_new = PointGroupElement(g.R, g.conjugate, g.antisymmetry, U=g.U @ t, _U_eq=_U_strict_eq)\n", + " assert allclose((g_new**n).U, sign * np.eye(g_new.U.shape[0]))\n", + " return g_new, n, sign" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[identity_power(g, double_group=True) for g in sg if g.conjugate==False]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cubic = qsymm.groups.cubic(tr=False, ph=False, generators=True)\n", + "gens = [g for g in sg if g in cubic]\n", + "gens" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def check_consistency(new_group, double_group):\n", + " for g in new_group:\n", + " _, n, sign = identity_power(g, double_group=double_group)\n", + " if n % 2 == 0 and not allclose((g**n).U, sign * np.eye(g.U.shape[0])):\n", + " # print(g, (g**n).U, sign, n)\n", + " return False\n", + " return True\n", + "\n", + "def U_double_eq(U1, U2):\n", + " if allclose(U1, U2):\n", + " return True\n", + " elif allclose(U1, -U2):\n", + " return False\n", + " else:\n", + " raise ValueError(\"Improper phase fixing detected in unitary parts of symmetry operators.\")\n", + "\n", + "def U_single_eq(U1, U2):\n", + " if allclose(U1, U2):\n", + " return True\n", + " else:\n", + " raise ValueError(\"Improper phase fixing detected in unitary parts of symmetry operators.\")\n", + "\n", + "\n", + "def fix_pg_phases(generators, double_group=False):\n", + " \"\"\"\n", + " Fix phases of unitaries such that the (double) point group generated by generators\n", + " forms a true representation.\n", + " \"\"\"\n", + " if any(g.conjugate for g in generators) or any(g.antisymmetry for g in generators):\n", + " raise ValueError('Only unitary point groups are supported. Make sure generators '\n", + " 'do not contain any elements with conjugate=True or antisymmetry=True.')\n", + " gens_orders = [identity_power(g, double_group=double_group) for g in generators]\n", + " \n", + " n_group = len(generate_group(generators))\n", + "\n", + " fixes = []\n", + " \n", + " U_eq = U_double_eq if double_group else U_single_eq\n", + "\n", + " for phases in it.product(*[range(n) for _, n, _ in gens_orders]):\n", + " new_gens = []\n", + " for (g, n, _), phase in zip(gens_orders, phases):\n", + " phase = 2 * np.pi * phase / n\n", + " g_new = PointGroupElement(g.R, U=np.exp(1j * phase) * g.U, _U_eq=U_eq)\n", + " new_gens.append(g_new)\n", + " try:\n", + " # This will raise an error if not a single/double representation\n", + " new_group = qsymm.groups.generate_group(new_gens)\n", + " except ValueError:\n", + " continue\n", + " print(len(new_group))\n", + " if ((double_group and not len(new_group) == 2 * n_group)\n", + " or (not double_group and not len(new_group) == n_group)):\n", + " continue\n", + " # Check that the phase choices are consistent\n", + " if check_consistency(new_group, double_group):\n", + " # return new_gens\n", + " fixes.append(new_gens)\n", + " # print(new_gens)\n", + " print(phases)\n", + " return fixes[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fix_pg_phases(gens, double_group=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_U_strict_eq.cache_info()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def conjugate_classes(group):\n", + " rest = set(group)\n", + " conjugate_classes = dict()\n", + " while rest:\n", + " g = rest.pop()\n", + " conjugates = {h * g * h.inv() for h in group}\n", + " conjugate_classes[g] = conjugates\n", + " rest -= conjugates\n", + " return conjugate_classes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ccs = conjugate_classes(qsymm.groups.generate_group(gens))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[(key, len(val)) for key, val in ccs.items()]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def character_table(gens, double_group=False):\n", + " \"\"\"\n", + " Find character table of group described by gens.\n", + " \"\"\"\n", + " fixed_gens = fix_pg_phases(gens, double_group)\n", + " fixed_group = generate_group(fixed_gens)\n", + " ccs = conjugate_classes(fixed_group)\n", + " charater_table = [(len(cc), rep, np.round(np.trace(rep.U), 3)) for rep, cc in ccs.items()]\n", + " return charater_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "character_table(gens, double_group=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def fix_pg_phases(generators, double_group=False):\n", + " \"\"\"\n", + " Fix phases of unitaries such that the (double) point group generated by generators\n", + " forms a true representation.\n", + " \"\"\"\n", + " if any(g.conjugate for g in generators) or any(g.antisymmetry for g in generators):\n", + " raise ValueError('Only unitary point groups are supported. Make sure generators '\n", + " 'do not contain any elements with conjugate=True or antisymmetry=True.')\n", + " gens_orders = [identity_power(g, double_group=double_group) for g in generators]\n", + "\n", + " fixes = []\n", + "\n", + " for phases in it.product(*[range(n) for _, n in gens_orders]):\n", + " new_gens = []\n", + " for (g, n), phase in zip(gens_orders, phases):\n", + " phase = 2 * np.pi * phase / n\n", + " g_new = PointGroupElement(g.R, U=np.exp(1j * phase) * g.U, _U_eq=_U_strict_eq)\n", + " new_gens.append(g_new)\n", + " new_group = qsymm.groups.generate_group(new_gens)\n", + " if len(new_group) > 96:\n", + " continue\n", + " # Check that the phase choices are consistent\n", + " # Brute force\n", + " for g, h in it.product(new_group, repeat=2):\n", + " gh = g * h\n", + " gh_new = [g for g in new_group if g == gh][0]\n", + " if not (np.allclose(gh_new.U, gh.U)\n", + " or (double_group and np.allclose(gh_new.U, -gh.U))):\n", + " # If inconsistent, break\n", + " break\n", + " # If didn't break, we found a consistent phase choice\n", + " else:\n", + " # return new_gens\n", + " fixes.append(new_gens)\n", + " # print(new_gens)\n", + " # print(phases)\n", + " return fixes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def character_table(gens, double_group=False):\n", + " \"\"\"\n", + " Find character table of group described by gens.\n", + " \"\"\"\n", + " for fixed_gens in fix_pg_phases(gens, double_group):\n", + " fixed_group = generate_group(fixed_gens)\n", + " ccs = conjugate_classes(fixed_group)\n", + " charater_table = [(len(cc), rep, np.round(np.trace(rep.U), 3)) for rep, cc in ccs.items()]\n", + " display(charater_table)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "character_table(gens, double_group=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.dot(ta.array(np.eye(2)), ta.array(np.eye(2)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import tinyarray as ta" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = ta.array(1j * np.eye(2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "type(a).mro()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from qsymm.groups import rotation_to_angle, spin_rotation\n", + "\n", + "def o3_to_su2(g):\n", + " \"\"\"\n", + " Transform O(3) spatial rotation to SU(2) x Z2 representation.\n", + " \"\"\"\n", + " R = np.array(g.R)\n", + " if np.isclose(la.det(R), 1):\n", + " inv = 1\n", + " n, theta = rotation_to_angle(R)\n", + " else:\n", + " inv = -1\n", + " n, theta = rotation_to_angle(-R)\n", + " R_new = la.block_diag(spin_rotation(theta * n, 1/2), [[inv]])\n", + " return PointGroupElement(R_new, g.conjugate, g.antisymmetry, g.U, g._U_eq)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[o3_to_su2(g) for g in gens]" + ] + }, { "cell_type": "code", "execution_count": null, @@ -831,9 +1213,9 @@ } }, "kernelspec": { - "display_name": "test_qsymm", + "display_name": "Python [conda env:kwant_conda]", "language": "python", - "name": "test_qsymm" + "name": "conda-env-kwant_conda-py" }, "language_info": { "codemirror_mode": { @@ -845,7 +1227,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.5" + "version": "3.7.1" } }, "nbformat": 4, -- GitLab