diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index bde55660ce941f41247f46aba3d26c1008dbb121..189abcd80bc9aa4769207139ff1885c1bcaa0190 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -75,6 +75,13 @@ if bdiag_broken: # skip coverage
return b_mat
+def lstsq(a, b):
+ """Least squares version that works also with 0-shaped matrices."""
+ if a.shape[1] == 0:
+ return np.empty((0, 0), dtype=np.common_type(a, b))
+ return np.linalg.lstsq(a, b)[0]
+
+
def nonzero_symm_projection(matrix):
"""Check whether a discrete symmetry relation between two blocks of the
Hamiltonian vanishes or not.
@@ -755,9 +762,10 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
if (np.abs(lmbdainv[indx].real - 1) < eps).all():
lmbdainv[indx] = 1
else:
- # Momenta are the negative arguments of the translation eigenvalues,
- # as computed below using np.angle. np.angle of -1 is pi, so this
- # assigns k = -pi to modes with translation eigenvalue -1.
+ # Momenta are the negative arguments of the translation
+ # eigenvalues, as computed below using np.angle. np.angle of -1
+ # is pi, so this assigns k = -pi to modes with translation
+ # eigenvalue -1.
lmbdainv[indx] = -1
# Original wave functions
@@ -814,6 +822,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
if chiral is not None and time_reversal is None:
out_orig = full_psi[:, indx[len(indx)//2:]]
out = chiral.dot(full_psi[:, indx[:len(indx)//2]])
+ # No least squares below because the modes should be orthogonal.
rot = out_orig.T.conj().dot(out)
full_psi[:, indx[len(indx)//2:]] = out
psi[:, indx[len(indx)//2:]] = psi[:, indx[len(indx)//2:]].dot(rot)
@@ -859,7 +868,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
# reverse the order of the product at the end.
wf_neg_k = particle_hole.dot(
(full_psi[:, :N][:, positive_k]).conj())[:, ::-1]
- rot = orig_neg_k.T.conj().dot(wf_neg_k)
+ rot = lstsq(orig_neg_k, wf_neg_k)
full_psi[:, :N][:, positive_k[::-1]] = wf_neg_k
psi[:, :N][:, positive_k[::-1]] = \
psi[:, :N][:, positive_k[::-1]].dot(rot)
@@ -874,7 +883,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
# Reverse order at the end to match momenta of opposite sign.
wf_neg_k = particle_hole.dot(
full_psi[:, N:][:, positive_k].conj())[:, ::-1]
- rot = orig_neg_k.T.conj().dot(wf_neg_k)
+ rot = lstsq(orig_neg_k, wf_neg_k)
full_psi[:, N:][:, positive_k[::-1]] = wf_neg_k
psi[:, N:][:, positive_k[::-1]] = \
psi[:, N:][:, positive_k[::-1]].dot(rot)
@@ -886,7 +895,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
# of propagating modes, not either left or right movers.
out_orig = full_psi[:, nmodes//2:]
out = time_reversal.dot(full_psi[:, :nmodes//2].conj())
- rot = out_orig.T.conj().dot(out)
+ rot = lstsq(out_orig, out)
full_psi[:, nmodes//2:] = out
psi[:, nmodes//2:] = psi[:, nmodes//2:].dot(rot)
@@ -954,29 +963,38 @@ def transform_modes(modes_data, unitary=None, time_reversal=None,
"""
wave_functions, momenta, velocities, vecs, vecslmbdainv, v = modes_data
+
+ # Copy to not overwrite modes from previous blocks
+ wave_functions = wave_functions.copy()
+ momenta = momenta.copy()
+ velocities = velocities.copy()
+ vecs = vecs.copy()
+ vecslmbdainv = vecslmbdainv.copy()
+ v = v.copy()
+
nmodes = wave_functions.shape[1] // 2
if unitary is not None:
perm = np.arange(2*nmodes)
conj = False
- flip_velocity = False
+ flip_energy = False
elif time_reversal is not None:
unitary = time_reversal
perm = np.arange(2*nmodes)[::-1]
conj = True
- flip_velocity = True
+ flip_energy = False
elif particle_hole is not None:
unitary = particle_hole
perm = ((-1-np.arange(2*nmodes)) % nmodes +
nmodes * (np.arange(2*nmodes) // nmodes))
conj = True
- flip_velocity = False
+ flip_energy = True
elif chiral is not None:
unitary = chiral
perm = (np.arange(2*nmodes) % nmodes +
nmodes * (np.arange(2*nmodes) < nmodes))
conj = False
- flip_velocity = True
+ flip_energy = True
else: # skip coverage
raise ValueError("No relation between blocks was provided.")
@@ -987,17 +1005,18 @@ def transform_modes(modes_data, unitary=None, time_reversal=None,
wave_functions = wave_functions.conj()
v = v.conj()
+ if flip_energy:
+ vecslmbdainv *= -1
+
+ if flip_energy != conj:
+ velocities *= -1
+
wave_functions = unitary.dot(wave_functions)[:, perm]
v = unitary.dot(v)
- # Copy to not overwrite modes from previous blocks
- vecs = vecs.copy()
vecs[:, :2*nmodes] = vecs[:, perm]
- vecslmbdainv = vecslmbdainv.copy()
vecslmbdainv[:, :2*nmodes] = vecslmbdainv[:, perm]
velocities = velocities[perm]
momenta = momenta[perm]
- if flip_velocity:
- velocities *= -1
return (wave_functions, momenta, velocities, vecs, vecslmbdainv, v)
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index c4bb4b2c554df1c1aeff64613def79aee93b7b1f..c1a1d151509c7ff13d72e0ba2353c87b095d726e 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -974,10 +974,12 @@ def test_block_relations_cons_PHS():
offsets[1:]]).T]
# Need both incident and outgoing stabilized modes to make the
# comparison between blocks
- prop_modes = [(vecs[:, :nmodes], vecs[:, nmodes:2*nmodes]),
+ prop_modes = [(vecs[:, :nmodes], vecs[:, nmodes:2*nmodes], 1),
(vecslmbdainv[:, :nmodes],
- vecslmbdainv[:, nmodes:2*nmodes])]
- for (in_modes, out_modes) in prop_modes:
+ vecslmbdainv[:, nmodes:2*nmodes], -1)]
+ # P is antiunitary, such that vecslmbdainv changes sign when
+ # used between blocks to construct modes.
+ for (in_modes, out_modes, vecs_sign) in prop_modes:
# Coordinates of blocks 2 and 3
rows2, cols2 = block_rows[2], block_cols[2]
cols3 = block_cols[3]
@@ -999,7 +1001,8 @@ def test_block_relations_cons_PHS():
# they are specified. Block 2 is computed before block 3, so
# block 3 is obtained by particle-hole transforming the modes
# of block 2. Check that it is so.
- assert_almost_equal(P_mat.dot(modes2.conj())[:, perm], modes3)
+ assert_almost_equal(P_mat.dot(modes2.conj())[:, perm], vecs_sign*modes3)
+
def check_symm_ham(h_cell, h_hop, sym_op, trans_sign, sym):
"""Check that the symmetry operator and Hamiltonian are properly defined"""
@@ -1128,10 +1131,12 @@ def test_blocks_symm_complex_projectors():
# Need both incident and outgoing stabilized modes to make the
# comparison between blocks
- prop_modes = [(vecs[:, :nmodes], vecs[:, nmodes:2*nmodes]),
+ prop_modes = [(vecs[:, :nmodes], vecs[:, nmodes:2*nmodes], 1),
(vecslmbdainv[:, :nmodes],
- vecslmbdainv[:, nmodes:2*nmodes])]
- for (in_modes, out_modes) in prop_modes:
+ vecslmbdainv[:, nmodes:2*nmodes], trans_sign)]
+ # Symmetries that flip the sign of energy change the sign of
+ # vecslmbdainv when used to construct modes between blocks.
+ for (in_modes, out_modes, vecs_sign) in prop_modes:
rows0, cols0 = block_rows[0], block_cols[0]
rows1, cols1 = block_rows[1], block_cols[1]
# Make sure the blocks are not empty
@@ -1147,9 +1152,9 @@ def test_blocks_symm_complex_projectors():
# block 1 is obtained by symmetry transforming the modes of
# block 0. Check that it is so.
if sym in ['AI', 'AII', 'C', 'D']:
- assert_almost_equal(S_t.dot(modes0.conj())[:, perm], modes1)
+ assert_almost_equal(S_t.dot(modes0.conj())[:, perm], vecs_sign*modes1)
elif sym in ['AIII']:
- assert_almost_equal(S_t.dot(modes0)[:, perm], modes1)
+ assert_almost_equal(S_t.dot(modes0)[:, perm], vecs_sign*modes1)
# If first block is empty, so is the second one.
else:
assert not in_modes[rows1, cols1].size