diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 3335af8adfbf3348ec230800996633f9e973be4c..5f8a4262d3e021803cc13f9bef5e607fb24f4a4c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -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
diff --git a/INSTALL.rst b/INSTALL.rst
index 88209df56c007829dfb910216b58811cbd8485ec..a32658ebbcdcea4c46d78d6f0b8287974c0cd71a 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -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
diff --git a/kwant/builder.py b/kwant/builder.py
index b31edc6beeefb0a308ed530f5d79efbf3babcd49..390ca0ba3c54d30d080d67c5002b43a49b74f4bf 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -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
diff --git a/kwant/kpm.py b/kwant/kpm.py
index a897c8f9f00b949fcaf694b7f5866b191c625b0c..5da09b097e394bc57c07c7de97fb958a3665fc3b 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -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]
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index fc799edc25705cd7c2d483a59a3bec7fce70cfcd..384cc8d0d57b44bfc611985d0093d3b28e6e9a4b 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -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.
diff --git a/kwant/physics/symmetry.py b/kwant/physics/symmetry.py
index 9420f52d4967b66c22cfb971a7b948a10a377581..6c648a95d1687acd2b4c6318277e635ae6e6171a 100644
--- a/kwant/physics/symmetry.py
+++ b/kwant/physics/symmetry.py
@@ -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):
diff --git a/kwant/qsymm.py b/kwant/qsymm.py
index bb752e19314a18cbc7ec53e53eb163df5bf17f80..fef2fd6438bd7028ad7dd22f3eaecf65dc2f300d 100644
--- a/kwant/qsymm.py
+++ b/kwant/qsymm.py
@@ -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():
diff --git a/kwant/tests/test_kpm.py b/kwant/tests/test_kpm.py
index 7ae187339e51163549637cf44306643b7f816bf4..524cc57d78a97434ed6605f36b72d5833c5fd410 100644
--- a/kwant/tests/test_kpm.py
+++ b/kwant/tests/test_kpm.py
@@ -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)
diff --git a/kwant/tests/test_qsymm.py b/kwant/tests/test_qsymm.py
index bf7e45d31286b30e5dcbd12247db96bbc279b6b0..7714c1c397d49171c8c6699860752ba333ed444c 100644
--- a/kwant/tests/test_qsymm.py
+++ b/kwant/tests/test_qsymm.py
@@ -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()
diff --git a/setup.py b/setup.py
index 5d4c30873c1a420bbf9221ec2c4f239991e3044d..cef8c9f545c81c4eb49868ff125cde816237c2ec 100755
--- a/setup.py
+++ b/setup.py
@@ -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')])