From 4c1b6d75710b7bc22cb0727ceaf50e655352160d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?T=C3=B3mas=20Rosdahl?= <torosdahl@gmail.com>
Date: Wed, 6 Feb 2019 18:06:53 +0100
Subject: [PATCH] add kwant.qsymm module
Co-authored-by: Daniel Varjas <dvarjas@gmail.com>
---
kwant/qsymm.py | 459 +++++++++++++++++++++++++++++++++++
kwant/tests/test_qsymm.py | 487 ++++++++++++++++++++++++++++++++++++++
2 files changed, 946 insertions(+)
create mode 100644 kwant/qsymm.py
create mode 100644 kwant/tests/test_qsymm.py
diff --git a/kwant/qsymm.py b/kwant/qsymm.py
new file mode 100644
index 00000000..291f7e2b
--- /dev/null
+++ b/kwant/qsymm.py
@@ -0,0 +1,459 @@
+# Copyright 2011-2018 Kwant authors.
+#
+# This file is part of Kwant. It is subject to the license terms in the file
+# LICENSE.rst found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+
+__all__ = ['builder_to_model', 'model_to_builder', 'find_builder_symmetries']
+
+import itertools as it
+from collections import OrderedDict, defaultdict
+
+import numpy as np
+import tinyarray as ta
+import scipy.linalg as la
+
+try:
+ import sympy
+ import qsymm
+ from qsymm.model import Model, BlochModel, BlochCoeff
+ from qsymm.groups import PointGroupElement, ContinuousGroupGenerator
+ from qsymm.symmetry_finder import bravais_point_group
+ from qsymm.linalg import allclose
+ from qsymm.hamiltonian_generator import hamiltonian_from_family
+except ImportError as error:
+ msg = ("'kwant.qsymm' is not available because one or more of its "
+ "dependencies is not installed.")
+ raise ImportError(msg) from error
+
+from kwant import lattice, builder
+from kwant._common import get_parameters
+
+
+def builder_to_model(syst, momenta=None, real_space=True,
+ params=None):
+ """Make a qsymm.BlochModel out of a Builder.
+
+ Parameters
+ ----------
+ syst: kwant.Builder
+ May have translational symmetries.
+ momenta: list of strings or None
+ Names of momentum variables. If None, 'k_x', 'k_y', ... is used.
+ real_space: bool (default True)
+ If False, use the unit cell convention for Bloch basis, the
+ exponential has the difference in the unit cell coordinates and
+ k is expressed in the reciprocal lattice basis. This is consistent
+ with kwant.wraparound.
+ If True, the difference in the real space coordinates is used
+ and k is given in an absolute basis.
+ Only the default choice guarantees that qsymm is able to find
+ nonsymmorphic symmetries.
+ params: dict, optional
+ Dictionary of parameter names and their values; used when
+ evaluating the Hamiltonian matrix elements.
+
+ Returns
+ -------
+ qsymm.BlochModel
+ Model representing the tight-binding Hamiltonian.
+
+ Notes
+ -----
+ The sites in the the builder are in lexicographical order, i.e. ordered
+ first by their family and then by their tag. This is the same ordering that
+ is used in finalized kwant systems.
+ """
+ def term_to_model(d, par, matrix):
+ if np.allclose(matrix, 0):
+ result = BlochModel({})
+ else:
+ result = BlochModel({BlochCoeff(d, qsymm.sympify(par)): matrix},
+ momenta=momenta)
+ return result
+
+ def hopping_to_model(hop, value, proj, params):
+ site1, site2 = hop
+ if real_space:
+ d = proj @ np.array(site2.pos - site1.pos)
+ else:
+ # site in the FD
+ d = np.array(syst.symmetry.which(site2))
+
+ slice1, slice2 = slices[to_fd(site1)], slices[to_fd(site2)]
+ if callable(value):
+ return sum(term_to_model(d, par, set_block(slice1, slice2, val))
+ for par, val in function_to_terms(hop, value, params))
+ else:
+ matrix = set_block(slice1, slice2, value)
+ return term_to_model(d, '1', matrix)
+
+ def onsite_to_model(site, value, params):
+ d = np.zeros((dim, ))
+ slice1 = slices[to_fd(site)]
+ if callable(value):
+ return sum(term_to_model(d, par, set_block(slice1, slice1, val))
+ for par, val in function_to_terms(site, value, params))
+ else:
+ return term_to_model(d, '1', set_block(slice1, slice1, value))
+
+ def function_to_terms(site_or_hop, value, fixed_params):
+ assert callable(value)
+ parameters = get_parameters(value)
+ # remove site or site1, site2 parameters
+ if isinstance(site_or_hop, builder.Site):
+ parameters = parameters[1:]
+ site_or_hop = (site_or_hop,)
+ else:
+ parameters = parameters[2:]
+ free_parameters = (par for par in parameters
+ if par not in fixed_params.keys())
+ # first set all free parameters to 0
+ args = ((fixed_params[par] if par in fixed_params.keys() else 0)
+ for par in parameters)
+ h_0 = value(*site_or_hop, *args)
+ # set one of the free parameters to 1 at a time, the rest 0
+ terms = []
+ for p in free_parameters:
+ args = ((fixed_params[par] if par in fixed_params.keys() else
+ (1 if par == p else 0)) for par in parameters)
+ terms.append((p, value(*site_or_hop, *args) - h_0))
+ return terms + [('1', h_0)]
+
+ def orbital_slices(syst):
+ orbital_slices = {}
+ start_orb = 0
+
+ for site in sorted(syst.sites()):
+ n = site.family.norbs
+ if n is None:
+ raise ValueError('norbs must be provided for every lattice.')
+ orbital_slices[site] = slice(start_orb, start_orb + n)
+ start_orb += n
+ return orbital_slices, start_orb
+
+ def set_block(slice1, slice2, val):
+ matrix = np.zeros((N, N), dtype=complex)
+ matrix[slice1, slice2] = val
+ return matrix
+
+ if params is None:
+ params = dict()
+
+ periods = np.array(syst.symmetry.periods)
+ dim = len(periods)
+ to_fd = syst.symmetry.to_fd
+ if momenta is None:
+ momenta = ['k_x', 'k_y', 'k_z'][:dim]
+ # If the system is higher dimensional than the number of translation
+ # vectors, we need to project onto the subspace spanned by the
+ # translation vectors.
+ if dim == 0:
+ proj = np.empty((0, len(list(syst.sites())[0].pos)))
+ elif dim < len(list(syst.sites())[0].pos):
+ proj, r = la.qr(np.array(periods).T, mode='economic')
+ sign = np.diag(np.diag(np.sign(r)))
+ proj = sign @ proj.T
+ else:
+ proj = np.eye(dim)
+
+ slices, N = orbital_slices(syst)
+
+ one_way_hoppings = [hopping_to_model(hop, value, proj, params)
+ for hop, value in syst.hopping_value_pairs()]
+ other_way_hoppings = [term.T().conj() for term in one_way_hoppings]
+ hoppings = one_way_hoppings + other_way_hoppings
+
+ onsites = [onsite_to_model(site, value, params)
+ for site, value in syst.site_value_pairs()]
+
+ result = sum(onsites) + sum(hoppings)
+
+ return result
+
+
+def model_to_builder(model, norbs, lat_vecs, atom_coords, *, coeffs=None):
+ """Make a Builder out of qsymm.Models or qsymm.BlochModels.
+
+ Parameters
+ ----------
+ model: qsymm.model.Model, qsymm.model.BlochModel, or an iterable thereof
+ The Hamiltonian (or terms of the Hamiltonian) to convert to a
+ Builder.
+ norbs: OrderedDict or sequence of pairs
+ Maps sites to the number of orbitals per site in a unit cell.
+ lat_vecs: list of arrays
+ Lattice vectors of the underlying tight binding lattice.
+ atom_coords: list of arrays
+ Positions of the sites (or atoms) within a unit cell.
+ The ordering of the atoms is the same as in norbs.
+ coeffs: list of sympy.Symbol, default None.
+ Constant prefactors for the individual terms in model, if model
+ is a list of multiple objects. If model is a single Model or BlochModel
+ object, this argument is ignored. By default assigns the coefficient
+ c_n to element model[n].
+
+ Returns
+ -------
+ syst: kwant.Builder
+ The unfinalized Kwant system representing the qsymm Model(s).
+
+ Notes
+ -----
+ Onsite terms that are not provided in the input model are set
+ to zero by default.
+ The input Model(s) representing the tight binding Hamiltonian in
+ Bloch form should follow the convention where the difference in the real
+ space atomic positions appear in the Bloch factors.
+ """
+
+ def make_int(R):
+ # If close to an integer array convert to integer tinyarray, else
+ # return None
+ R_int = ta.array(np.round(R), int)
+ if qsymm.linalg.allclose(R, R_int):
+ return R_int
+ else:
+ return None
+
+ def term_onsite(onsites_dict, hopping_dict, hop_mat, atoms,
+ sublattices, coords_dict):
+ """Find the Kwant onsites and hoppings in a qsymm.BlochModel term
+ that has no lattice translation in the Bloch factor.
+ """
+ for atom1, atom2 in it.product(atoms, atoms):
+ # Subblock within the same sublattice is onsite
+ hop = hop_mat[ranges[atom1], ranges[atom2]]
+ if sublattices[atom1] == sublattices[atom2]:
+ onsites_dict[atom1] += Model({coeff: hop}, momenta=momenta)
+ # Blocks between sublattices are hoppings between sublattices
+ # at the same position.
+ # Only include nonzero hoppings
+ elif not allclose(hop, 0):
+ if not allclose(np.array(coords_dict[atom1]),
+ np.array(coords_dict[atom2])):
+ raise ValueError(
+ "Position of sites not compatible with qsymm model.")
+ lat_basis = np.array(zer)
+ hop = Model({coeff: hop}, momenta=momenta)
+ hop_dir = builder.HoppingKind(-lat_basis, sublattices[atom1],
+ sublattices[atom2])
+ hopping_dict[hop_dir] += hop
+ return onsites_dict, hopping_dict
+
+ def term_hopping(hopping_dict, hop_mat, atoms,
+ sublattices, coords_dict):
+ """Find Kwant hoppings in a qsymm.BlochModel term that has a lattice
+ translation in the Bloch factor.
+ """
+ # Iterate over combinations of atoms, set hoppings between each
+ for atom1, atom2 in it.product(atoms, atoms):
+ # Take the block from atom1 to atom2
+ hop = hop_mat[ranges[atom1], ranges[atom2]]
+ # Only include nonzero hoppings
+ if allclose(hop, 0):
+ continue
+ # Adjust hopping vector to Bloch form basis
+ r_lattice = (
+ r_vec
+ + np.array(coords_dict[atom1])
+ - np.array(coords_dict[atom2])
+ )
+ # Bring vector to basis of lattice vectors
+ lat_basis = np.linalg.solve(np.vstack(lat_vecs).T, r_lattice)
+ lat_basis = make_int(lat_basis)
+ # Should only have hoppings that are integer multiples of
+ # lattice vectors
+ if lat_basis is not None:
+ hop_dir = builder.HoppingKind(-lat_basis,
+ sublattices[atom1],
+ sublattices[atom2])
+ # Set the hopping as the matrix times the hopping amplitude
+ hopping_dict[hop_dir] += Model({coeff: hop}, momenta=momenta)
+ else:
+ raise RuntimeError('A nonzero hopping not matching a '
+ 'lattice vector was found.')
+ return hopping_dict
+
+ # Disambiguate single model instances from iterables thereof. Because
+ # Model is itself iterable (subclasses dict) this is a bit cumbersome.
+ if isinstance(model, Model):
+ # BlochModel can't yet handle getting a Blochmodel as input
+ if not isinstance(model, BlochModel):
+ model = BlochModel(model)
+ else:
+ model = BlochModel(hamiltonian_from_family(
+ model, coeffs=coeffs, nsimplify=False, tosympy=False))
+
+
+ # 'momentum' and 'zer' are used in the closures defined above, so don't
+ # move these declarations down.
+ momenta = model.momenta
+ if len(momenta) != len(lat_vecs):
+ raise ValueError("Dimension of the lattice and number of "
+ "momenta do not match.")
+ zer = [0] * len(momenta)
+
+
+ # Subblocks of the Hamiltonian for different atoms.
+ N = 0
+ if not any([isinstance(norbs, OrderedDict), isinstance(norbs, list),
+ isinstance(norbs, tuple)]):
+ raise ValueError('norbs must be OrderedDict, tuple, or list.')
+ else:
+ norbs = OrderedDict(norbs)
+ ranges = dict()
+ for a, n in norbs.items():
+ ranges[a] = slice(N, N + n)
+ N += n
+
+ # Extract atoms and number of orbitals per atom,
+ # store the position of each atom
+ atoms, orbs = zip(*norbs.items())
+ coords_dict = dict(zip(atoms, atom_coords))
+
+ # Make the kwant lattice
+ lat = lattice.general(lat_vecs, atom_coords, norbs=orbs)
+ # Store sublattices by name
+ sublattices = dict(zip(atoms, lat.sublattices))
+
+ # Keep track of the hoppings and onsites by storing those
+ # which have already been set.
+ hopping_dict = defaultdict(dict)
+ onsites_dict = defaultdict(dict)
+
+ # Iterate over all terms in the model.
+ for key, hop_mat in model.items():
+ # Determine whether this term is an onsite or a hopping, extract
+ # overall symbolic coefficient if any, extract the exponential
+ # part describing the hopping if present.
+ r_vec, coeff = key
+ # Onsite term; modifies onsites_dict and hopping_dict in-place
+ if allclose(r_vec, 0):
+ term_onsite(
+ onsites_dict, hopping_dict, hop_mat,
+ atoms, sublattices, coords_dict)
+ # Hopping term; modifies hopping_dict in-place
+ else:
+ term_hopping(hopping_dict, hop_mat, atoms,
+ sublattices, coords_dict)
+
+ # If some onsite terms are not set, we set them to zero.
+ for atom in atoms:
+ if atom not in onsites_dict:
+ onsites_dict[atom] = Model(
+ {sympy.numbers.One(): np.zeros((norbs[atom], norbs[atom]))},
+ momenta=momenta)
+
+ # Make the Kwant system, and set all onsites and hoppings.
+
+ sym = lattice.TranslationalSymmetry(*lat_vecs)
+ syst = builder.Builder(sym)
+
+ # Iterate over all onsites and set them
+ for atom, onsite in onsites_dict.items():
+ syst[sublattices[atom](*zer)] = onsite.lambdify(onsite=True)
+
+ # Finally, iterate over all the hoppings and set them
+ for direction, hopping in hopping_dict.items():
+ syst[direction] = hopping.lambdify(hopping=True)
+
+ return syst
+
+
+# This may be useful in the future, so we'll keep it as internal for now,
+# and can make it part of the API in the future if we wish.
+def _get_builder_symmetries(builder):
+ """Extract the declared symmetries of a Kwant builder.
+
+ Parameters
+ ----------
+ builder: kwant.Builder
+
+ Returns
+ -------
+ builder_symmetries: dict
+ Dictionary of the discrete symmetries that the builder has.
+ The symmetries can be particle-hole, time-reversal or chiral,
+ which are returned as qsymm.PointGroupElements, or
+ a conservation law, which is returned as a
+ qsymm.ContinuousGroupGenerators.
+ """
+
+ dim = len(np.array(builder.symmetry.periods))
+ symmetry_names = ['time_reversal', 'particle_hole', 'chiral',
+ 'conservation_law']
+ builder_symmetries = {name: getattr(builder, name)
+ for name in symmetry_names
+ if getattr(builder, name) is not None}
+ for name, symmetry in builder_symmetries.items():
+ if name == 'time_reversal':
+ builder_symmetries[name] = PointGroupElement(np.eye(dim),
+ True, False, symmetry)
+ elif name == 'particle_hole':
+ builder_symmetries[name] = PointGroupElement(np.eye(dim),
+ True, True, symmetry)
+ elif name == 'chiral':
+ builder_symmetries[name] = PointGroupElement(np.eye(dim),
+ False, True, symmetry)
+ elif name == 'conservation_law':
+ builder_symmetries[name] = ContinuousGroupGenerator(R=None,
+ U=symmetry)
+ else:
+ raise ValueError("Invalid symmetry name.")
+ return builder_symmetries
+
+
+def find_builder_symmetries(builder, momenta=None, params=None,
+ spatial_symmetries=True, prettify=True):
+ """Finds the symmetries of a Kwant system using qsymm.
+
+ Parameters
+ ----------
+ builder: kwant.Builder
+ momenta: list of strings or None
+ Names of momentum variables, if None 'k_x', 'k_y', ... is used.
+ params: dict, optional
+ Dictionary of parameter names and their values; used when
+ evaluating the Hamiltonian matrix elements.
+ spatial_symmetries: bool (default True)
+ If True, search for all symmetries.
+ If False, only searches for the symmetries that are declarable in
+ kwant.Builder objects, i.e. time-reversal symmetry, particle-hole
+ symmetry, chiral symmetry, or conservation laws. This can save
+ computation time.
+ prettify : bool (default True)
+ Whether to carry out sparsification of the continuous symmetry
+ generators, in general an arbitrary linear combination of the
+ symmetry generators is returned.
+
+ Returns
+ -------
+ symmetries: list of qsymm.PointGroupElements and/or
+ qsymm.ContinuousGroupElement
+ The symmetries of the Kwant system.
+ """
+
+ if params is None:
+ params = dict()
+
+ ham = builder_to_model(builder, momenta=momenta,
+ real_space=False, params=params)
+
+ dim = len(np.array(builder.symmetry.periods))
+
+ if spatial_symmetries:
+ candidates = bravais_point_group(builder.symmetry.periods, tr=True,
+ ph=True, generators=False,
+ verbose=False)
+ else:
+ candidates = [
+ qsymm.PointGroupElement(np.eye(dim), True, False, None), # T
+ qsymm.PointGroupElement(np.eye(dim), True, True, None), # P
+ qsymm.PointGroupElement(np.eye(dim), False, True, None)] # C
+ sg, cg = qsymm.symmetries(ham, candidates, prettify=prettify,
+ continuous_rotations=False)
+ return list(sg) + list(cg)
diff --git a/kwant/tests/test_qsymm.py b/kwant/tests/test_qsymm.py
new file mode 100644
index 00000000..7f0d1825
--- /dev/null
+++ b/kwant/tests/test_qsymm.py
@@ -0,0 +1,487 @@
+# Copyright 2011-2018 Kwant authors.
+#
+# This file is part of Kwant. It is subject to the license terms in the file
+# LICENSE.rst found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from collections import OrderedDict
+
+import numpy as np
+import sympy
+
+import qsymm
+from qsymm.symmetry_finder import symmetries
+from qsymm.hamiltonian_generator import bloch_family, hamiltonian_from_family
+from qsymm.groups import (hexagonal, PointGroupElement, spin_matrices,
+ spin_rotation, ContinuousGroupGenerator)
+from qsymm.model import Model, e, I, _commutative_momenta
+from qsymm.linalg import allclose
+
+import kwant
+from kwant._common import ensure_rng
+from kwant.lattice import TranslationalSymmetry
+from kwant.qsymm import (builder_to_model, model_to_builder,
+ _get_builder_symmetries, find_builder_symmetries)
+
+
+def test_honeycomb():
+ lat = kwant.lattice.honeycomb(norbs=1)
+
+ # Test simple honeycomb model with constant terms
+ # Add discrete symmetries to the kwant builder as well, to check that they are
+ # returned as well.
+ syst = kwant.Builder(symmetry=TranslationalSymmetry(*lat.prim_vecs))
+ syst[lat.a(0, 0)] = 1
+ syst[lat.b(0, 0)] = 1
+ syst[lat.neighbors(1)] = -1
+
+ H = builder_to_model(syst)
+ sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+ assert len(sg) == 24
+ assert len(cs) == 0
+
+ # Test simple honeycomb model with value functions
+ syst = kwant.Builder(symmetry=TranslationalSymmetry(*lat.prim_vecs))
+ syst[lat.a(0, 0)] = lambda site, ma: ma
+ syst[lat.b(0, 0)] = lambda site, mb: mb
+ syst[lat.neighbors(1)] = lambda site1, site2, t: t
+
+ H = builder_to_model(syst)
+ sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+ assert len(sg) == 12
+ assert len(cs) == 0
+
+
+def test_get_builder_symmetries():
+ syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0), (0, 1)),
+ particle_hole=np.eye(2),
+ conservation_law=2*np.eye(2))
+ builder_symmetries = _get_builder_symmetries(syst)
+
+ assert len(builder_symmetries) == 2
+ P = builder_symmetries['particle_hole']
+ assert isinstance(P, PointGroupElement)
+ assert allclose(P.U, np.eye(2))
+ assert P.conjugate and P.antisymmetry
+ assert allclose(P.R, np.eye(2))
+
+ cons = builder_symmetries['conservation_law']
+ assert isinstance(cons, ContinuousGroupGenerator)
+ assert allclose(cons.U, 2*np.eye(2))
+ assert cons.R is None
+
+ syst = kwant.Builder()
+ builder_symmetries = _get_builder_symmetries(syst)
+ assert len(builder_symmetries) == 0
+
+
+def test_higher_dim():
+ # Test 0D finite system
+ lat = kwant.lattice.cubic(norbs=1)
+ syst = kwant.Builder()
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(1, 1, 0)] = 1
+ syst[lat(0, 1, 1)] = 1
+ syst[lat(1, 0, -1)] = 1
+ syst[lat(0, 0, 0), lat(1, 1, 0)] = -1
+ syst[lat(0, 0, 0), lat(0, 1, 1)] = -1
+ syst[lat(0, 0, 0), lat(1, 0, -1)] = -1
+
+ H = builder_to_model(syst)
+ sg, cs = symmetries(H, prettify=True)
+ assert len(sg) == 2
+ assert len(cs) == 5
+
+ # Test triangular lattice system embedded in 3D
+ sym = TranslationalSymmetry([1, 1, 0], [0, 1, 1])
+ lat = kwant.lattice.cubic(norbs=1)
+ syst = kwant.Builder(symmetry=sym)
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(0, 0, 0), lat(1, 1, 0)] = -1
+ syst[lat(0, 0, 0), lat(0, 1, 1)] = -1
+ syst[lat(0, 0, 0), lat(1, 0, -1)] = -1
+
+ H = builder_to_model(syst)
+ sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+ assert len(sg) == 24
+ assert len(cs) == 0
+
+
+def test_graphene_to_kwant():
+
+ norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each
+ hopping_vectors = [('A', 'B', [1, 0])] # Hopping between neighbouring A and B atoms
+ # Atomic coordinates within the unit cell
+ atom_coords = [(0, 0), (1, 0)]
+ # We set the interatom distance to 1, so the lattice vectors have length sqrt(3)
+ lat_vecs = [(3/2, np.sqrt(3)/2), (3/2, -np.sqrt(3)/2)]
+
+ # Time reversal
+ TR = PointGroupElement(sympy.eye(2), True, False, np.eye(2))
+ # Chiral symmetry
+ C = PointGroupElement(sympy.eye(2), False, True, np.array([[1, 0], [0, -1]]))
+ # Atom A rotates into A, B into B.
+ sphi = 2*sympy.pi/3
+ RC3 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],
+ [sympy.sin(sphi), sympy.cos(sphi)]])
+ C3 = PointGroupElement(RC3, False, False, np.eye(2))
+
+ # Generate graphene Hamiltonian in Kwant from qsymm
+ symmetries = [C, TR, C3]
+ # Generate using a family
+ family = bloch_family(hopping_vectors, symmetries, norbs)
+ syst_from_family = model_to_builder(family, norbs, lat_vecs, atom_coords, coeffs=None)
+ # Generate using a single Model object
+ g = sympy.Symbol('g', real=True)
+ ham = hamiltonian_from_family(family, coeffs=[g])
+ ham = Model(hamiltonian=ham, momenta=family[0].momenta)
+ syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)
+
+ # Make the graphene Hamiltonian using kwant only
+ atoms, orbs = zip(*[(atom, norb) for atom, norb in
+ norbs.items()])
+ # Make the kwant lattice
+ lat = kwant.lattice.general(lat_vecs,
+ atom_coords,
+ norbs=orbs)
+ # Store sublattices by name
+ sublattices = {atom: sublat for atom, sublat in
+ zip(atoms, lat.sublattices)}
+
+ sym = kwant.TranslationalSymmetry(*lat_vecs)
+ bulk = kwant.Builder(sym)
+
+ bulk[ [sublattices['A'](0, 0), sublattices['B'](0, 0)] ] = 0
+
+ def hop(site1, site2, c0):
+ return c0
+
+ bulk[lat.neighbors()] = hop
+
+ fsyst_family = kwant.wraparound.wraparound(syst_from_family).finalized()
+ fsyst_model = kwant.wraparound.wraparound(syst_from_model).finalized()
+ fsyst_kwant = kwant.wraparound.wraparound(bulk).finalized()
+
+ # Check that the energies are identical at random points in the Brillouin zone
+ coeff = 0.5 + np.random.rand()
+ for _ in range(20):
+ kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+ params = dict(c0=coeff, k_x=kx, k_y=ky)
+ hamiltonian1 = fsyst_kwant.hamiltonian_submatrix(params=params, sparse=False)
+ hamiltonian2 = fsyst_family.hamiltonian_submatrix(params=params, sparse=False)
+ assert allclose(hamiltonian1, hamiltonian2)
+ params = dict(g=coeff, k_x=kx, k_y=ky)
+ hamiltonian3 = fsyst_model.hamiltonian_submatrix(params=params, sparse=False)
+ assert allclose(hamiltonian2, hamiltonian3)
+
+ # Include random onsites as well
+ one = sympy.numbers.One()
+ onsites = [Model({one: np.array([[1, 0], [0, 0]])}, momenta=family[0].momenta),
+ Model({one: np.array([[0, 0], [0, 1]])}, momenta=family[0].momenta)]
+ family = family + onsites
+ syst_from_family = model_to_builder(family, norbs, lat_vecs, atom_coords, coeffs=None)
+ gs = list(sympy.symbols('g0:%d'%3, real=True))
+ ham = hamiltonian_from_family(family, coeffs=gs)
+ ham = Model(hamiltonian=ham, momenta=family[0].momenta)
+ syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)
+
+ def onsite_A(site, c1):
+ return c1
+
+ def onsite_B(site, c2):
+ return c2
+
+ bulk[[sublattices['A'](0, 0)]] = onsite_A
+ bulk[[sublattices['B'](0, 0)]] = onsite_B
+
+ fsyst_family = kwant.wraparound.wraparound(syst_from_family).finalized()
+ fsyst_model = kwant.wraparound.wraparound(syst_from_model).finalized()
+ fsyst_kwant = kwant.wraparound.wraparound(bulk).finalized()
+
+ # Check equivalence of the Hamiltonian at random points in the BZ
+ coeffs = 0.5 + np.random.rand(3)
+ for _ in range(20):
+ kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+ params = dict(c0=coeffs[0], c1=coeffs[1], c2=coeffs[2], k_x=kx, k_y=ky)
+ hamiltonian1 = fsyst_kwant.hamiltonian_submatrix(params=params, sparse=False)
+ hamiltonian2 = fsyst_family.hamiltonian_submatrix(params=params, sparse=False)
+ assert allclose(hamiltonian1, hamiltonian2)
+ params = dict(g0=coeffs[0], g1=coeffs[1], g2=coeffs[2], k_x=kx, k_y=ky)
+ hamiltonian3 = fsyst_model.hamiltonian_submatrix(params=params, sparse=False)
+ assert allclose(hamiltonian2, hamiltonian3)
+
+
+def test_wraparound_convention():
+ # Test that it matches exactly kwant.wraparound convention
+ # Make the graphene Hamiltonian using kwant only
+ norbs = OrderedDict([('A', 1), ('B', 1)]) # A and B atom per unit cell, one orbital each
+ atoms, orbs = zip(*[(atom, norb) for atom, norb in
+ norbs.items()])
+ # Atomic coordinates within the unit cell
+ atom_coords = [(0, 0), (1, 0)]
+ # We set the interatom distance to 1, so the lattice vectors have length sqrt(3)
+ lat_vecs = [(3/2, np.sqrt(3)/2), (3/2, -np.sqrt(3)/2)]
+ # Make the kwant lattice
+ lat = kwant.lattice.general(lat_vecs,
+ atom_coords,
+ norbs=orbs)
+ # Store sublattices by name
+ sublattices = {atom: sublat for atom, sublat in
+ zip(atoms, lat.sublattices)}
+
+ sym = kwant.TranslationalSymmetry(*lat_vecs)
+ bulk = kwant.Builder(sym)
+
+ bulk[ [sublattices['A'](0, 0), sublattices['B'](0, 0)] ] = 0
+
+ def hop(site1, site2, c0):
+ return c0
+
+ bulk[lat.neighbors()] = hop
+
+ wrapped = kwant.wraparound.wraparound(bulk).finalized()
+ ham2 = builder_to_model(bulk, real_space=False)
+ # Check that the Hamiltonians are identical at random points in the Brillouin zone
+ H1 = wrapped.hamiltonian_submatrix
+ H2 = ham2.lambdify()
+ coeffs = 0.5 + np.random.rand(1)
+ for _ in range(20):
+ kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+ params = dict(c0=coeffs[0], k_x=kx, k_y=ky)
+ h1, h2 = H1(params=params), H2(**params)
+ assert allclose(h1, h2), (h1, h2)
+
+
+
+def test_inverse_transform():
+ # Define family on square lattice
+ s = spin_matrices(1/2)
+ # Time reversal
+ TR = PointGroupElement(np.eye(2), True, False,
+ spin_rotation(2 * np.pi * np.array([0, 1/2, 0]), s))
+ # Mirror symmetry
+ Mx = PointGroupElement(np.array([[-1, 0], [0, 1]]), False, False,
+ spin_rotation(2 * np.pi * np.array([1/2, 0, 0]), s))
+ # Fourfold
+ C4 = PointGroupElement(np.array([[0, 1], [-1, 0]]), False, False,
+ spin_rotation(2 * np.pi * np.array([0, 0, 1/4]), s))
+ symmetries = [TR, Mx, C4]
+
+ # One site per unit cell
+ norbs = OrderedDict([('A', 2)])
+ # Hopping to a neighbouring atom one primitive lattice vector away
+ hopping_vectors = [('A', 'A', [1, 0])]
+ # Make family
+ family = bloch_family(hopping_vectors, symmetries, norbs)
+ fam = hamiltonian_from_family(family, tosympy=False)
+ # Atomic coordinates within the unit cell
+ atom_coords = [(0, 0)]
+ lat_vecs = [(1, 0), (0, 1)]
+ syst = model_to_builder(fam, norbs, lat_vecs, atom_coords)
+ # Convert it back
+ ham2 = builder_to_model(syst).tomodel(nsimplify=True)
+ # Check that it's the same as the original
+ assert fam == ham2
+
+ # Check that the Hamiltonians are identical at random points in the Brillouin zone
+ sysw = kwant.wraparound.wraparound(syst).finalized()
+ H1 = sysw.hamiltonian_submatrix
+ H2 = ham2.lambdify()
+ H3 = fam.lambdify()
+ coeffs = 0.5 + np.random.rand(3)
+ for _ in range(20):
+ kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+ params = dict(c0=coeffs[0], c1=coeffs[1], c2=coeffs[2], k_x=kx, k_y=ky)
+ assert allclose(H1(params=params), H2(**params))
+ assert allclose(H1(params=params), H3(**params))
+
+
+def test_consistency_kwant():
+ """Make a random 1D Model, convert it to a builder, and compare
+ the Bloch representation of the Model with that which Kwant uses
+ in wraparound and in Bands. Then, convert the builder back to a Model
+ and compare with the original Model.
+ For comparison, we also make the system using Kwant only.
+ """
+ orbs = 4
+ T = np.random.rand(2*orbs, 2*orbs) + 1j*np.random.rand(2*orbs, 2*orbs)
+ H = np.random.rand(2*orbs, 2*orbs) + 1j*np.random.rand(2*orbs, 2*orbs)
+ H += H.T.conj()
+
+ # Make the 1D Model manually using only qsymm features.
+ c0, c1 = sympy.symbols('c0 c1', real=True)
+ kx = _commutative_momenta[0]
+
+ Ham = Model({c0 * e**(-I*kx): T}, momenta=[0])
+ Ham += Ham.T().conj()
+ Ham += Model({c1: H}, momenta=[0])
+
+ # Two superimposed atoms, same number of orbitals on each
+ norbs = OrderedDict([('A', orbs), ('B', orbs)])
+ atom_coords = [(0.3, ), (0.3, )]
+ lat_vecs = [(1, )] # Lattice vector
+
+ # Make a Kwant builder out of the qsymm Model
+ model_syst = model_to_builder(Ham, norbs, lat_vecs, atom_coords)
+ fmodel_syst = model_syst.finalized()
+
+ # Make the same system manually using only Kwant features.
+ lat = kwant.lattice.general(np.array([[1.]]),
+ [(0., )],
+ norbs=2*orbs)
+ kwant_syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
+
+ def onsite(site, c1):
+ return c1*H
+
+ def hopping(site1, site2, c0):
+ return c0*T
+
+ sublat = lat.sublattices[0]
+ kwant_syst[sublat(0,)] = onsite
+ hopp = kwant.builder.HoppingKind((1, ), sublat)
+ kwant_syst[hopp] = hopping
+ fkwant_syst = kwant_syst.finalized()
+
+ # Make sure we are consistent with bands calculations in kwant
+ # The Bloch Hamiltonian used in Kwant for the bands computation
+ # is h(k) = exp(-i*k)*hop + onsite + exp(i*k)*hop.T.conj.
+ # We also check that all is consistent with wraparound
+ coeffs = (0.7, 1.2)
+ params = dict(c0 = coeffs[0], c1 = coeffs[1])
+ kwant_hop = fkwant_syst.inter_cell_hopping(params=params)
+ kwant_onsite = fkwant_syst.cell_hamiltonian(params=params)
+ model_kwant_hop = fmodel_syst.inter_cell_hopping(params=params)
+ model_kwant_onsite = fmodel_syst.cell_hamiltonian(params=params)
+
+ assert allclose(model_kwant_hop, coeffs[0]*T)
+ assert allclose(model_kwant_hop, kwant_hop)
+ assert allclose(model_kwant_onsite, kwant_onsite)
+
+ h_model_kwant = (lambda k: np.exp(-1j*k)*model_kwant_hop + model_kwant_onsite +
+ np.exp(1j*k)*model_kwant_hop.T.conj()) # As in kwant.Bands
+ h_model = Ham.lambdify()
+ wsyst = kwant.wraparound.wraparound(model_syst).finalized()
+ for _ in range(20):
+ k = (np.random.rand() - 0.5)*2*np.pi
+ assert allclose(h_model_kwant(k), h_model(coeffs[0], coeffs[1], k))
+ params['k_x'] = k
+ h_wrap = wsyst.hamiltonian_submatrix(params=params)
+ assert allclose(h_model(coeffs[0], coeffs[1], k), h_wrap)
+
+ # Get the model back from the builder
+ # From the Kwant builder based on original Model
+ Ham1 = builder_to_model(model_syst, momenta=Ham.momenta).tomodel(nsimplify=True)
+ # From the pure Kwant builder
+ Ham2 = builder_to_model(kwant_syst, momenta=Ham.momenta).tomodel(nsimplify=True)
+ assert Ham == Ham1
+ assert Ham == Ham2
+
+
+def test_find_builder_discrete_symmetries():
+ symm_class = ['AI', 'D', 'AIII', 'BDI']
+ class_dict = {'AI': ['time_reversal'],
+ 'D': ['particle_hole'],
+ 'AIII': ['chiral'],
+ 'BDI': ['time_reversal', 'particle_hole', 'chiral']}
+ sym_dict = {'time_reversal': qsymm.PointGroupElement(np.eye(2), True, False, None),
+ 'particle_hole': qsymm.PointGroupElement(np.eye(2), True, True, None),
+ 'chiral': qsymm.PointGroupElement(np.eye(2), False, True, None)}
+ n = 4
+ rng = 11
+ for sym in symm_class:
+ # Random Hamiltonian in the symmetry class
+ h_ons = kwant.rmt.gaussian(n, sym, rng=rng)
+ h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
+ # Make a Kwant builder in the symmetry class and find its symmetries
+ lat = kwant.lattice.square(norbs=n)
+ bulk = kwant.Builder(TranslationalSymmetry([1, 0], [0, 1]))
+ bulk[lat(0, 0)] = h_ons
+ bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop
+ bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop
+ builder_symmetries = find_builder_symmetries(bulk, spatial_symmetries=True, prettify=True)
+
+ # Equality of symmetries ignores unitary part
+ fourfold_rotation = qsymm.PointGroupElement(np.array([[0, 1],[1, 0]]), False, False, None)
+ assert fourfold_rotation in builder_symmetries
+ class_symmetries = class_dict[sym]
+ for class_symmetry in class_symmetries:
+ assert sym_dict[class_symmetry] in builder_symmetries
+
+
+def random_onsite_hop(n, rng=0):
+ rng = ensure_rng(rng)
+ onsite = rng.randn(n, n) + 1j * rng.randn(n, n)
+ onsite = onsite + onsite.T.conj()
+ hop = rng.rand(n, n) + 1j * rng.rand(n, n)
+ return onsite, hop
+
+
+def test_find_cons_law():
+ sy = np.array([[0, -1j], [1j, 0]])
+
+ n = 3
+ lat = kwant.lattice.chain(norbs=2*n)
+ syst = kwant.Builder()
+ rng = 1337
+ ons, hop = random_onsite_hop(n, rng=rng)
+
+ syst[lat(0)] = np.kron(sy, ons)
+ syst[lat(1)] = np.kron(sy, ons)
+ syst[lat(1), lat(0)] = np.kron(sy, hop)
+
+ builder_symmetries = find_builder_symmetries(syst, spatial_symmetries=False, prettify=True)
+ onsites = [symm for symm in builder_symmetries if
+ isinstance(symm, qsymm.ContinuousGroupGenerator) and symm.R is None]
+ mham = builder_to_model(syst)
+ assert not any([len(symm.apply(mham)) for symm in onsites])
+
+def test_basis_ordering():
+ symm_class = ['AI', 'D', 'AIII', 'BDI']
+ n = 2
+ rng = 12
+ for sym in symm_class:
+ # Make a small finite system in the symmetry class, finalize it and check
+ # that the basis is consistent.
+ h_ons = kwant.rmt.gaussian(n, sym, rng=rng)
+ h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
+ lat = kwant.lattice.square(norbs=n)
+ bulk = kwant.Builder(TranslationalSymmetry([1, 0], [0, 1]))
+ bulk[lat(0, 0)] = h_ons
+ bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop
+ bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop
+
+ def rect(site):
+ x, y = site.pos
+ return (0 <= x < 2) and (0 <= y < 3)
+
+ square = kwant.Builder()
+ square.fill(bulk, lambda site: rect(site), (0, 0),
+ max_sites=float('inf'))
+
+ # Find the symmetries of the square
+ builder_symmetries = find_builder_symmetries(square,
+ spatial_symmetries=False,
+ prettify=True)
+ # Finalize the square, extract Hamiltonian
+ fsquare = square.finalized()
+ ham = fsquare.hamiltonian_submatrix()
+
+ # Check manually that the found symmetries are in the same basis as the
+ # finalized system
+ for symmetry in builder_symmetries:
+ U = symmetry.U
+ if isinstance(symmetry, qsymm.ContinuousGroupGenerator):
+ assert symmetry.R is None
+ assert allclose(U.dot(ham), ham.dot(U))
+ else:
+ if symmetry.conjugate:
+ left = U.dot(ham.conj())
+ else:
+ left = U.dot(ham)
+ if symmetry.antisymmetry:
+ assert allclose(left, -ham.dot(U))
+ else:
+ assert allclose(left, ham.dot(U))
--
GitLab