@@ -89,24 +89,18 @@ build-env:default:

   before_script:

     - source deactivate

     - source activate kwant-latest

-    - pip install qsymm

 .bleeding-edge-env: &bleeding_edge_env

   before_script:

     - source deactivate

     - conda env update -f /kwant-latest.yml

     - source activate kwant-latest

-    - pip install qsymm

 .ubuntu-env: &ubuntu_env

   image: gitlab.kwant-project.org:5005/kwant/kwant/ubuntu

-  before_script:

-    - pip3 install qsymm

 .debian-env: &debian_env

   image: gitlab.kwant-project.org:5005/kwant/kwant/debian

-  before_script:

-    - pip3 install qsymm

 ## Build Jobs

@@ -45,7 +45,7 @@ a NumPy-like Python package optimized for very small arrays,

 The following software is highly recommended though not strictly required:

  * `matplotlib <http://matplotlib.org/>`_ 1.5.1 or newer, for the module `kwant.plotter` and the tutorial,

  * `SymPy <http://sympy.org/>`_ 0.7.6 or newer, for the subpackage `kwant.continuum`.

- * `Qsymm <https://pypi.org/project/qsymm/>`_ 1.1.0 or newer, for the subpackage `kwant.qsymm`.

+ * `Qsymm <https://pypi.org/project/qsymm/>`_ 1.2.6 or newer, for the subpackage `kwant.qsymm`.

  * `MUMPS <http://graal.ens-lyon.fr/MUMPS/>`_, a sparse linear algebra library

    that will in many cases speed up Kwant several times and reduce the memory

    footprint.  (Kwant uses only the sequential, single core version

@@ -864,6 +864,10 @@ class Builder:

     amounts to creating a `Lead` object and appending it to the list of leads

     accessbile as the `~Builder.leads` attribute.

+    `conservation_law`, `time_reversal`, `particle_hole`, and `chiral`

+    affect the basis in which scattering modes derived from the builder

+    are expressed - see `~kwant.physics.DiscreteSymmetry` for details.

+

     .. warning::

         If functions are used to set values in a builder with a symmetry, then

@@ -780,6 +780,8 @@ def conductivity(hamiltonian, alpha='x', beta='x', positions=None, **kwargs):

     Parameters

     ----------

+    hamiltonian : `~kwant.system.FiniteSystem` or matrix Hamiltonian

+        If a system is passed, it should contain no leads.

     alpha, beta : str, or operators

         If ``hamiltonian`` is a kwant system, or if the ``positions``

         are provided, ``alpha`` and ``beta`` can be the directions of the

@@ -1051,10 +1053,11 @@ def _velocity(hamiltonian, params, op_type, positions):

     elif isinstance(op_type, str):

         direction = directions[op_type]

         if isinstance(hamiltonian, system.System):

-            operator = hamiltonian.hamiltonian_submatrix(params=params,

-                                                         sparse=True)

-            positions = np.array([site.pos for site in hamiltonian.sites

-                                  for iorb in range(site.family.norbs)])

+            operator, norbs, norbs = hamiltonian.hamiltonian_submatrix(

+                params=params, sparse=True, return_norb=True

+            )

+            positions = np.vstack([[hamiltonian.pos(i)] * norb

+                                   for i, norb in enumerate(norbs)])

         elif positions is not None:

             operator = coo_matrix(hamiltonian, copy=True)

         displacements = positions[operator.col] - positions[operator.row]

@@ -130,11 +130,10 @@ class PropagatingModes:

     momentum and velocity, an arbitrary orthonormal basis in the subspace of

     these modes is chosen.

-    If a conservation law is specified to block diagonalize the Hamiltonian

-    to N blocks, then `block_nmodes[i]` is the number of left or right moving

-    propagating modes in conservation law block `i`. The ordering of blocks

-    is the same as the ordering of the projectors used to block diagonalize

-    the Hamiltonian.

+    If a conservation law is specified to block diagonalize the Hamiltonian,

+    then `block_nmodes[i]` is the number of left or right moving propagating

+    modes in conservation law block `i`. The ordering of blocks is the same as

+    the ordering of the projectors used to block diagonalize the Hamiltonian.

     """

     def __init__(self, wave_functions, velocities, momenta):

         kwargs = locals()

@@ -1046,6 +1045,10 @@ def modes(h_cell, h_hop, tol=1e6, stabilization=None, *,

     Propagating modes with the same momentum are orthogonalized. All the

     propagating modes are normalized by current.

+    `projectors`, `time_reversal`, `particle_hole`, and `chiral` affect the

+    basis in which the scattering modes are expressed - see

+    `~kwant.physics.DiscreteSymmetry` for details.

+

     This function uses the most stable and efficient algorithm for calculating

     the mode decomposition that the Kwant authors are aware about. Its details

     are to be published.

@@ -37,7 +37,7 @@ _names = ['Time reversal', 'Particle-hole', 'Chiral']

 _signs = [1, -1, -1]

 class DiscreteSymmetry:

-    """A collection of discrete symmetries and conservation laws.

+    r"""A collection of discrete symmetries and conservation laws.

     Parameters

     ----------

@@ -52,22 +52,46 @@ class DiscreteSymmetry:

     Notes

     -----

-    Whenever one or more discrete symmetry is declared in conjunction with a

-    conservation law, the symmetry operators and projectors must be declared

-    in canonical form. This means that each block of the Hamiltonian is

-    transformed either to itself by a discrete symmetry or to a single

-    other block.

-

-    More formally, consider a discrete symmetry S. The symmetry projection

-    that maps from block i to block j of the Hamiltonian with projectors

-    :math:`P_i` and :math:`P_j` is :math:`S_{ji} = P_j^+ S P_i`.

-    If :math:`S_{ji}` is nonzero, a symmetry relation exists between

-    blocks i and j. Canonical form means that for each j, the block

-    :math:`S_{ji}` is nonzero at most for one i, while all other blocks vanish.

-

-    If the operators are not in canonical form, they can be made so by

-    further splitting the Hamiltonian into smaller blocks, i.e. by adding

-    more projectors.

+

+    When computing scattering modes, the representation of the

+    modes is chosen to reflect declared discrete symmetries and

+    conservation laws.

+

+    `projectors` block diagonalize the Hamiltonian, and modes are computed

+    separately in each block. The ordering of blocks is the same as of

+    `projectors`. If `conservation_law` is declared in

+    `~kwant.builder.Builder`, `projectors` is computed as the projectors

+    onto its orthogonal eigensubspaces. The projectors are stored in the

+    order of ascending eigenvalues of `conservation_law`.

+

+    Symmetrization using discrete symmetries varies depending on whether

+    a conservation law/projectors are declared. Consider the case with no

+    conservation law declared. With `time_reversal` declared, the outgoing

+    modes are chosen as the time-reversed partners of the incoming modes,

+    i.e. :math:`\psi_{out}(-k) = T \psi_{in}(k)` with k the momentum.

+    `chiral` also relates incoming and outgoing modes, such that

+    :math:`\psi_{out}(k) = C \psi_{in}(k)`. `particle_hole` gives symmetric

+    incoming and outgoing modes separately, such that

+    :math:`\psi_{in/out}(-k) = P \psi_{in/out}(k)`, except when k=-k, at

+    k = 0 or :math:`\pi`. In this case, each mode is chosen as an eigenstate

+    :math:`P \psi = \psi` if :math:`P^2=1`. If :math:`P^2=-1`, we

+    symmetrize the modes by generating pairs of orthogonal modes

+    :math:`\psi` and :math:`P\psi`. Because `chiral` and `particle_hole`

+    flip the sign of energy, they only apply at zero energy.

+

+    Discrete symmetries can be combined with a conservation law if they

+    leave each block invariant or transform it to another block. With S

+    a discrete symmetry and :math:`P_i` and :math:`P_j` projectors onto

+    blocks i and j of the Hamiltonian, :math:`S_{ji} = P_j^+ S P_i` is the

+    symmetry projection that maps from block i to block j.

+    :math:`S_{ji} = P_j^+ S P_i` must for each j be nonzero for exactly

+    one i. If S leaves block i invariant, the modes within block i are

+    symmetrized using the nonzero projection :math:`S_{ii}`, like in the

+    case without a conservation law. If S transforms between blocks i and

+    j, the modes of the block with the larger index are obtained by

+    transforming the modes of the block with the lower index. Thus, with

+    :math:`\psi_i` and :math:`\psi_j` the modes of blocks i and j, we have

+    :math:`\psi_j = S_{ji} \psi_i`.

     """

     def __init__(self, projectors=None, time_reversal=None, particle_hole=None,

                  chiral=None):

@@ -68,8 +68,8 @@ def builder_to_model(syst, momenta=None, real_space=True,

     is used in finalized kwant systems.

     """

     def term_to_model(d, par, matrix):

-        if np.allclose(matrix, 0):

-            result = BlochModel({})

+        if allclose(matrix, 0):

+            result = BlochModel({}, shape=matrix.shape, format=np.ndarray)

         else:

             result = BlochModel({BlochCoeff(d, qsymm.sympify(par)): matrix},

                                 momenta=momenta)

@@ -323,8 +323,8 @@ def model_to_builder(model, norbs, lat_vecs, atom_coords, *, coeffs=None):

     # Keep track of the hoppings and onsites by storing those

     # which have already been set.

-    hopping_dict = defaultdict(dict)

-    onsites_dict = defaultdict(dict)

+    hopping_dict = defaultdict(lambda: 0)

+    onsites_dict = defaultdict(lambda: 0)

     # Iterate over all terms in the model.

     for key, hop_mat in model.items():

@@ -80,6 +80,13 @@ def test_conductivity():

         cond = kwant.kpm.conductivity(syst, alpha=alpha, beta=beta,

                                         positions=None, **kpm_params)

+    # test when 'norbs' not provided

+    n = lat.norbs

+    lat.norbs = None

+    cond = kwant.kpm.conductivity(syst, alpha=alpha, beta=beta,

+                                  positions=None, **kpm_params)

+    lat.norbs = n

+

     # test system or hamiltonian, no positions, but velocity operators

     cond_xx = kwant.kpm.conductivity(syst, alpha='x', beta='x',

                                      positions=None, **kpm_params)

@@ -19,7 +19,7 @@ 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.model import Model, BlochModel, BlochCoeff

 from qsymm.linalg import allclose

 import kwant

@@ -41,7 +41,7 @@ def test_honeycomb():

     syst[lat.neighbors(1)] = -1

     H = builder_to_model(syst)

-    sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)

+    sg, cs = symmetries(H, hexagonal(sympy_R=False))

     assert len(sg) == 24

     assert len(cs) == 0

@@ -52,7 +52,7 @@ def test_honeycomb():

     syst[lat.neighbors(1)] = lambda site1, site2, t: t

     H = builder_to_model(syst)

-    sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)

+    sg, cs = symmetries(H, hexagonal(sympy_R=False))

     assert len(sg) == 12

     assert len(cs) == 0

@@ -93,7 +93,7 @@ def test_higher_dim():

     syst[lat(0, 0, 0), lat(1, 0, -1)] = -1

     H = builder_to_model(syst)

-    sg, cs = symmetries(H, prettify=True)

+    sg, cs = symmetries(H)

     assert len(sg) == 2

     assert len(cs) == 5

@@ -107,7 +107,7 @@ def test_higher_dim():

     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)

+    sg, cs = symmetries(H, hexagonal(sympy_R=False))

     assert len(sg) == 24

     assert len(cs) == 0

@@ -137,9 +137,9 @@ def test_graphene_to_kwant():

     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)

+    g = sympy.Symbol('g')

+    # tosympy=False to return a BlochModel

+    ham = hamiltonian_from_family(family, coeffs=[g], tosympy=False)

     syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)

     # Make the graphene Hamiltonian using kwant only

@@ -185,9 +185,9 @@ def test_graphene_to_kwant():

                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)

+    gs = list(sympy.symbols('g0:3'))

+    # tosympy=False to return a BlochModel

+    ham = hamiltonian_from_family(family, coeffs=gs, tosympy=False)

     syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)

     def onsite_A(site, c1):

@@ -277,16 +277,16 @@ def test_inverse_transform():

     # 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)

+    family = bloch_family(hopping_vectors, symmetries, norbs, bloch_model=True)

     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)

+    ham2 = builder_to_model(syst)

     # Check that it's the same as the original

-    assert fam == ham2

+    assert fam.allclose(ham2)

     # Check that the Hamiltonians are identical at random points in the Brillouin zone

     sysw = kwant.wraparound.wraparound(syst).finalized()

@@ -314,12 +314,11 @@ def test_consistency_kwant():

     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]

+    c0, c1 = sympy.symbols('c0 c1')

-    Ham = Model({c0 * e**(-I*kx): T}, momenta=['k_x'])

+    Ham = BlochModel({BlochCoeff(np.array([-1]), c0): T}, momenta=['k_x'])

     Ham += Ham.T().conj()

-    Ham += Model({c1: H}, momenta=['k_x'])

+    Ham += BlochModel({BlochCoeff(np.array([0]), c1): H}, momenta=['k_x'])

     # Two superimposed atoms, same number of orbitals on each

     norbs = OrderedDict([('A', orbs), ('B', orbs)])

@@ -367,8 +366,8 @@ def test_consistency_kwant():

                      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

+    ks = np.linspace(-np.pi, np.pi, 21)

+    for k in ks:

         assert allclose(h_model_kwant(k), h_model(coeffs[0], coeffs[1], k))

         params['k_x'] = k

         h_wrap = wsyst.hamiltonian_submatrix(params=params)

@@ -376,11 +375,11 @@ def test_consistency_kwant():

     # 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)

+    Ham1 = builder_to_model(model_syst, momenta=Ham.momenta)

     # From the pure Kwant builder

-    Ham2 = builder_to_model(kwant_syst, momenta=Ham.momenta).tomodel(nsimplify=True)

-    assert Ham == Ham1

-    assert Ham == Ham2

+    Ham2 = builder_to_model(kwant_syst, momenta=Ham.momenta)

+    assert Ham.allclose(Ham1)

+    assert Ham.allclose(Ham2)

 def test_find_builder_discrete_symmetries():

@@ -405,12 +404,15 @@ def test_find_builder_discrete_symmetries():

         bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop

         bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop

+        # We need to specify 'prettify=True' here to ensure that we do not end up with

+        # an overcomplete set of symmetries. In some badly conditioned cases sparse=True

+        # or sparse=False may affect how many symmetries are found.

         builder_symmetries_default = find_builder_symmetries(bulk, spatial_symmetries=True,

                                                              prettify=True)

         builder_symmetries_sparse = find_builder_symmetries(bulk, spatial_symmetries=True,

                                                             prettify=True, sparse=True)

         builder_symmetries_dense = find_builder_symmetries(bulk, spatial_symmetries=True,

-                                                            prettify=True, sparse=False)

+                                                           prettify=True, sparse=False)

         assert len(builder_symmetries_default) == len(builder_symmetries_sparse)

         assert len(builder_symmetries_default) == len(builder_symmetries_dense)

@@ -476,7 +478,7 @@ def test_find_cons_law():

     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)

+    builder_symmetries = find_builder_symmetries(syst, spatial_symmetries=False)

     onsites = [symm for symm in builder_symmetries if

                isinstance(symm, qsymm.ContinuousGroupGenerator) and symm.R is None]

     mham = builder_to_model(syst)

@@ -508,8 +510,7 @@ def test_basis_ordering():

         # Find the symmetries of the square

         builder_symmetries = find_builder_symmetries(square,

-                                                     spatial_symmetries=False,

-                                                     prettify=True)

+                                                     spatial_symmetries=False)

         # Finalize the square, extract Hamiltonian

         fsquare = square.finalized()

         ham = fsquare.hamiltonian_submatrix()

@@ -588,7 +588,7 @@ def main():

               'plotting': 'matplotlib >= 1.5.1',

               'continuum': 'sympy >= 0.7.6',

               # qsymm is only packaged on PyPI

-              'qsymm': 'qsymm >= 1.1.2',

+              'qsymm': 'qsymm >= 1.2.6',

           },

           classifiers=[c.strip() for c in classifiers.split('\n')])