Merge branch 'fix/greens-function-reflection' into 'master'

Fix reflection calculation for kwant.greens_function

Closes #398

See merge request !377

kwant/solvers/common.py

@@ -989,17 +989,20 @@ class GreensFunction(BlockResult):
         result = np.trace(attdagainv).real
         if lead_out == lead_in:
             # For reflection we have to be more careful
-            gamma = 1j * (self.lead_info[lead_in] -
-                          self.lead_info[lead_in].conj().T)
+            sigma = self.lead_info[lead_in]
+            gamma = 1j * (sigma - sigma.conj().T)
             gf = self.submatrix(lead_out, lead_in)
 
             # The number of channels is given by the number of
             # nonzero eigenvalues of Gamma
             # rationale behind the threshold from
             # Golub; van Loan, chapter 5.5.8
-            eps = np.finfo(gamma.dtype).eps * 1000
-            N = np.sum(np.linalg.eigvalsh(gamma) >
-                       eps * np.linalg.norm(gamma, np.inf))
+            # We use ‖Σ‖, not ‖Γ‖, for the tolerance as ‖Γ‖~0 when there
+            # are no open modes.
+            eps = (
+                1e6 * np.finfo(gamma.dtype).eps * np.linalg.norm(sigma, np.inf)
+            )
+            N = np.sum(np.linalg.eigvalsh(gamma) > eps)
 
             result += 2 * np.trace(np.dot(gamma, gf)).imag + N
 

kwant/solvers/tests/test_solvers.py

@@ -6,7 +6,7 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from itertools import chain
+import itertools
 from math import cos, sin
 import numpy as np
 import pytest
@@ -37,7 +37,7 @@ sparse_solver_options = [
     {},
 ]
 
-solvers = list(chain(
+solvers = list(itertools.chain(
     [("mumps", opts) for opts in mumps_solver_options],
     [("sparse", opts) for opts in sparse_solver_options],
 ))
@@ -272,6 +272,29 @@ def test_two_equal_leads(smatrix):
     raises(ValueError, check_fsyst, fsyst.precalculate(what='selfenergy'))
 
 
+# Regression test against
+# https://gitlab.kwant-project.org/kwant/kwant/-/issues/398
+# that the reflection calculation in 'greens_function' does not
+# mis-count the number of open modes when there are none.
+def test_reflection_no_open_modes(greens_function):
+    # Build system
+    syst = kwant.Builder()
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+    syst[(square(i, j) for i in range(3) for j in range(3))] = 4
+    lead[(square(0, j) for j in range(3))] = 4
+    syst[square.neighbors()] = -1
+    lead[square.neighbors()] = -1
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+    syst = syst.finalized()
+
+    # Sanity check; no open modes at 0 energy
+    _, m = syst.leads[0].modes(energy=0)
+    assert m.nmodes == 0
+
+    assert np.isclose(greens_function(syst).transmission(0, 0), 0)
+
+
 # Test a more complicated graph with non-singular hopping.
 def test_graph_system(smatrix):
     rng = ensure_rng(11)