make selfenergy_from_modes a method of StabilizedModes

doc/source/reference/kwant.physics.rst

@@ -11,6 +11,5 @@ Leads
    Bands
    modes
    selfenergy
-   selfenergy_from_modes
    PropagatingModes
    StabilizedModes

kwant/builder.py

@@ -479,7 +479,7 @@ class ModesLead(Lead):
 
     """
     def __init__(self, modes_func, interface):
-        self._modes_func = modes_func
+        self.modes_func = modes_func
         self.interface = tuple(interface)
 
     def finalized(self):
@@ -487,11 +487,11 @@ class ModesLead(Lead):
         return self
 
     def modes(self, energy, args=()):
-        return self._modes_func(energy, args)
+        return self.modes_func(energy, args)
 
     def selfenergy(self, energy, args=()):
-        modes = self.modes(energy, args)
-        return physics.selfenergy_from_modes(modes=modes)
+        propagating, stabilized = self.modes(energy, args)
+        return stabilized.selfenergy()
 
 
 

kwant/physics/leads.py

@@ -17,8 +17,7 @@ from .. import linalg as kla
 
 dot = np.dot
 
-__all__ = ['selfenergy', 'modes', 'PropagatingModes', 'StabilizedModes',
-           'selfenergy_from_modes']
+__all__ = ['selfenergy', 'modes', 'PropagatingModes', 'StabilizedModes']
 
 
 # Container classes
@@ -103,6 +102,25 @@ class StabilizedModes(object):
     def __getitem__(self, item):
         return (self.vecs, self.vecslmbdainv, self.nmodes, self.sqrt_hop)[item]
 
+    def selfenergy(self):
+        """
+        Compute the self-energy generated by lead modes.
+
+        Returns
+        -------
+        Sigma : numpy array, real or complex, shape (M,M)
+            The computed self-energy. Note that even if `h_cell` and `h_hop` are
+            both real, `Sigma` will typically be complex. (More precisely, if
+            there is a propagating mode, `Sigma` will definitely be complex.)
+        """
+        v = self.sqrt_hop
+        vecs = self.vecs[:, self.nmodes:]
+        vecslmbdainv = self.vecslmbdainv[:, self.nmodes:]
+        if v is not None:
+            return dot(v, dot(vecs, la.solve(vecslmbdainv, v.T.conj())))
+        else:
+            return la.solve(vecslmbdainv.T, vecs.T).T
+
 
 # Auxiliary functions that perform different parts of the calculation.
 def setup_linsys(h_cell, h_hop, tol=1e6, stabilization=None):
@@ -534,9 +552,9 @@ def modes(h_cell, h_hop, tol=1e6, stabilization=None):
     Returns
     -------
     PropagatingModes(wave_functions, momenta, velocities) : object
-        A contained for the array of the wave functions of propagating modes,
-        their momenta, and their velocities. It can be used to identify the
-        gauge in which the scattering problem is solved. See
+        Contains the array of the wave functions of propagating modes, their
+        momenta, and their velocities. It can be used to identify the gauge in
+        which the scattering problem is solved. See
         `~kwant.physics.PropagatingModes` for details.
     StabilizedModes(vecs, vecslmbdainv, nmodes, sqrt_hop) : object
         A basis of propagating and evanescent modes used by the solvers.
@@ -597,37 +615,6 @@ def modes(h_cell, h_hop, tol=1e6, stabilization=None):
     return real_space_data, StabilizedModes(vecs, vecslmbdainv, nrightmovers, v)
 
 
-def selfenergy_from_modes(lead_modes):
-    """
-    Compute the self-energy generated by lead modes.
-
-    Parameters
-    ----------
-    lead_modes : StabilizedModes(vecs, vecslmbdainv, nmodes, v) a named tuple
-        The modes in the lead, with format defined in
-        `~kwant.physics.StabilizedMods`.
-
-    Returns
-    -------
-    Sigma : numpy array, real or complex, shape (M,M)
-        The computed self-energy. Note that even if `h_cell` and `h_hop` are
-        both real, `Sigma` will typically be complex. (More precisely, if there
-        is a propagating mode, `Sigma` will definitely be complex.)
-
-    Notes
-    -----
-    For simplicity this function relies on the calculation of modes as input.
-    This may cause a small slowdown, and can be improved if necessary.
-    """
-    vecs, vecslmbdainv, nmodes, v = lead_modes
-    vecs = vecs[:, nmodes:]
-    vecslmbdainv = vecslmbdainv[:, nmodes:]
-    if v is not None:
-        return dot(v, dot(vecs, la.solve(vecslmbdainv, v.T.conj())))
-    else:
-        return la.solve(vecslmbdainv.T, vecs.T).T
-
-
 def selfenergy(h_cell, h_hop, tol=1e6):
     """
     Compute the self-energy generated by the lead.
@@ -655,7 +642,8 @@ def selfenergy(h_cell, h_hop, tol=1e6):
     For simplicity this function internally calculates the modes first.
     This may cause a small slowdown, and can be improved if necessary.
     """
-    return selfenergy_from_modes(modes(h_cell, h_hop, tol)[1])
+    propagating, stabilized = modes(h_cell, h_hop, tol)
+    return stabilized.selfenergy()
 
 
 def square_selfenergy(width, hopping, fermi_energy):

kwant/system.py

@@ -102,7 +102,9 @@ class FiniteSystem(System):
             Numbers of the leads to be precalculated. If `None`, all are
             precalculated.
         calculate_selfenergy : bool
-            Whether to calculate self-energy if modes are available.
+            Whether to calculate self-energy if modes are available.  Defaults
+            to `False`.  Disabling this saves a typically negligible amount of
+            time and memory.
 
         Returns
         -------
@@ -127,7 +129,7 @@ class FiniteSystem(System):
             try:
                 modes = lead.modes(energy, args)
                 if calculate_selfenergy:
-                    selfenergy = physics.selfenergy_from_modes(modes)
+                    selfenergy = modes[1].selfenergy()
             except AttributeError:
                 selfenergy = lead.selfenergy(energy, args)
             new_leads.append(PrecalculatedLead(modes, selfenergy))