merge fix for setup_linsys from branch stable_py2

AUTHORS.rst

@@ -11,10 +11,11 @@ The principal developers of Kwant are
 
 The authors can be reached at authors@kwant-project.org.
 
-Other contributors to Kwant include
+Contributors to Kwant include
 
 * Mathieu Istas (INAC/CEA Grenoble)
 * Daniel Jaschke (INAC/CEA Grenoble)
+* Bas Nijholt (TU Delft)
 * Michał Nowak (TU Delft)
 * Viacheslav Ostroukh (Leiden University)
 * Adrien Sorgniard (INAC/CEA Grenoble)

kwant/physics/leads.py

@@ -198,22 +198,22 @@ def setup_linsys(h_cell, h_hop, tol=1e6, stabilization=None):
     if (n_nonsing == n and stabilization is None):
         # The hopping matrix is well-conditioned and can be safely inverted.
         # Hence the regular transfer matrix may be used.
-        hop_inv = la.inv(h_hop)
+        h_hop_sqrt = sqrt(np.linalg.norm(h_hop))
+        A = h_hop / h_hop_sqrt
+        B = h_hop_sqrt
+        B_H_inv = 1.0 / B     # just a real scalar here
+        A_inv = la.inv(A)
 
-        A = np.zeros((2*n, 2*n), dtype=np.common_type(h_cell, h_hop))
-        A[:n, :n] = dot(hop_inv, -h_cell)
-        A[:n, n:] = -hop_inv
-        A[n:, :n] = h_hop.T.conj()
-
-        # The function that can extract the full wave function psi from the
-        # projected one. Here it is almost trivial, but used for simplifying
-        # the logic.
+        lhs = np.zeros((2*n, 2*n), dtype=np.common_type(h_cell, h_hop))
+        lhs[:n, :n] = -dot(A_inv, h_cell) * B_H_inv
+        lhs[:n, n:] = -A_inv * B
+        lhs[n:, :n] = A.T.conj() * B_H_inv
 
         def extract_wf(psi, lmbdainv):
-            return np.copy(psi[:n])
+            return B_H_inv * np.copy(psi[:n])
 
-        matrices = (A, None)
-        v_out = None
+        matrices = (lhs, None)
+        v_out = h_hop_sqrt * np.eye(n)
     else:
         if stabilization is None:
             stabilization = [None, False]

kwant/physics/tests/test_leads.py

@@ -224,7 +224,7 @@ def test_modes():
     v = 2 * t * np.sin(k)
     prop, stab = leads.modes(np.array([[h]]), np.array([[t]]))
     assert stab.nmodes == 1
-    assert stab.sqrt_hop is None
+    assert stab.sqrt_hop[0] == np.sqrt(np.linalg.norm(t))
     np.testing.assert_almost_equal(prop.velocities, [-v, v])
     np.testing.assert_almost_equal(prop.momenta, [k, -k])
     # Test for normalization by current.
@@ -255,7 +255,7 @@ def test_algorithm_equivalence():
     u, v = u * np.sqrt(s), vh.T.conj() * np.sqrt(s)
     prop_vecs = []
     evan_vecs = []
-    algos = [None] + list(product(*(2 * [(True, False)])))
+    algos = [None, (True, True), (True, False), (False, True), (False, False)]
     for algo in algos:
         result = leads.modes(h, t, stabilization=algo)[1]
 
@@ -266,7 +266,9 @@ def test_algorithm_equivalence():
             vecs = np.dot(v, vecs)
             np.testing.assert_almost_equal(result.sqrt_hop, v)
         else:
-            vecslmbdainv = np.dot(v.T.conj(), vecslmbdainv)
+            vecslmbdainv = (np.dot(v.T.conj(), vecslmbdainv) /
+                            np.sqrt(np.linalg.norm(t)))
+            vecs = vecs * np.sqrt(np.linalg.norm(t))
         full_vecs = np.r_[vecslmbdainv, vecs]
 
         prop_vecs.append(full_vecs[:, : 2 * result.nmodes])
@@ -340,3 +342,18 @@ def test_zero_hopping():
     assert all(np.alltrue(getattr(actual[0], attr) ==
                           getattr(expected[0], attr)) for attr
                    in ('wave_functions', 'velocities', 'momenta'))
+
+
+def make_clean_lead(W, E, t):
+    syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
+    lat = kwant.lattice.square()
+    syst[(lat(0, j) for j in range(W))] = E
+    syst[lat.neighbors()] = -t
+    return syst.finalized()
+
+
+def test_momenta():
+    """Test whether the two systems have the same momenta,
+    these should not change when the Hamiltonian is scaled."""
+    momenta = [make_clean_lead(10, s, s).modes()[0].momenta for s in [1, 1e20]]
+    assert_almost_equal(*momenta)