diff --git a/doc/source/reference/kwant.physics.rst b/doc/source/reference/kwant.physics.rst index 2c457f5617e405622f42602a9f01c60a945ea555..e1c6f9019ea64b05153c1f2febf7bb4139e2a571 100644 --- a/doc/source/reference/kwant.physics.rst +++ b/doc/source/reference/kwant.physics.rst @@ -11,6 +11,5 @@ Leads Bands modes selfenergy - selfenergy_from_modes PropagatingModes StabilizedModes diff --git a/kwant/builder.py b/kwant/builder.py index 7487bc7503a9bd1084037f9a77cee6fd9496e74f..b6b6c66d5ff311d8c92f9acf5c35c56c19025d2d 100644 --- a/kwant/builder.py +++ b/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() diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py index e77d0353028d3e4e9fac9754c76f8321472b570e..8fb6ad31a053883aa970ba4d853b9d2aed8a004a 100644 --- a/kwant/physics/leads.py +++ b/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): diff --git a/kwant/system.py b/kwant/system.py index cb22be259cc4a52300915c0cf825a6a6341de67f..2a69650a9439cd56802cb93a7772b3ce585145c9 100644 --- a/kwant/system.py +++ b/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))