Commit e2a4343e authored by Joseph Weston's avatar Joseph Weston
Browse files

Merge branch 'feature/selfenergy-lead' into 'master'

Allow mode-space solvers to work with systems that have self-energy leads attached

Closes #368

See merge request kwant/kwant!366
parents c460cab8 65c9763d
...@@ -15,7 +15,7 @@ from numbers import Integral ...@@ -15,7 +15,7 @@ from numbers import Integral
import numpy as np import numpy as np
import scipy.sparse as sp import scipy.sparse as sp
from .._common import ensure_isinstance, deprecate_args from .._common import ensure_isinstance, deprecate_args
from .. import system from .. import system, physics
from functools import reduce from functools import reduce
# Until v0.13.0, scipy.sparse did not support making block matrices out of # Until v0.13.0, scipy.sparse did not support making block matrices out of
...@@ -139,8 +139,8 @@ class SparseSolver(metaclass=abc.ABCMeta): ...@@ -139,8 +139,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
lead_info : list of objects lead_info : list of objects
Contains one entry for each lead. If `realspace=False`, this is an Contains one entry for each lead. If `realspace=False`, this is an
instance of `~kwant.physics.PropagatingModes` with a corresponding instance of `~kwant.physics.PropagatingModes` for all leads that
format, otherwise the lead self-energy matrix. have modes, otherwise the lead self-energy matrix.
Notes Notes
----- -----
...@@ -188,55 +188,79 @@ class SparseSolver(metaclass=abc.ABCMeta): ...@@ -188,55 +188,79 @@ class SparseSolver(metaclass=abc.ABCMeta):
for leadnum, interface in enumerate(syst.lead_interfaces): for leadnum, interface in enumerate(syst.lead_interfaces):
lead = syst.leads[leadnum] lead = syst.leads[leadnum]
if not realspace: if not realspace:
prop, stab = lead.modes(energy, args, params=params) if system.is_selfenergy_lead(lead):
lead_info.append(prop) # We make equivalent checks to this in 'smatrix',
u = stab.vecs # 'greens_function' and 'wave_function' (where we
ulinv = stab.vecslmbdainv # can give more informative error messages).
nprop = stab.nmodes # Putting this check here again is just a sanity check,
svd_v = stab.sqrt_hop assert leadnum not in in_leads
sigma = np.asarray(lead.selfenergy(energy, args, params=params))
if len(u) == 0:
rhs.append(None) rhs.append(None)
continue lead_info.append(sigma)
indices.append(None)
indices.append(np.arange(lhs.shape[0], lhs.shape[0] + nprop)) # add selfenergy to the LHS
coords = np.r_[tuple(slice(offsets[i], offsets[i + 1])
u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:] for i in interface)]
u_in, ulinv_in = u[:, :nprop], ulinv[:, :nprop] if sigma.shape != 2 * coords.shape:
raise ValueError(
# Construct a matrix of 1's that translates the f"Self-energy dimension for lead {leadnum} "
# inter-cell hopping to a proper hopping "does not match the total number of orbitals "
# from the system to the lead. "of the sites for which it is defined."
iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1]) )
for i in interface)] y, x = np.meshgrid(coords, coords)
sig_sparse = splhsmat((sigma.flat, [x.flat, y.flat]),
n_lead_orbs = svd_v.shape[0] lhs.shape)
if n_lead_orbs != len(iface_orbs): lhs = lhs + sig_sparse
msg = ('Lead {0} has hopping with dimensions '
'incompatible with its interface dimension.')
raise ValueError(msg.format(leadnum))
coords = np.r_[[np.arange(len(iface_orbs))], [iface_orbs]]
transf = sp.csc_matrix((np.ones(len(iface_orbs)), coords),
shape=(iface_orbs.size, lhs.shape[0]))
v_sp = sp.csc_matrix(svd_v.T.conj()) * transf
vdaguout_sp = (transf.T *
sp.csc_matrix(np.dot(svd_v, u_out)))
lead_mat = - ulinv_out
lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]],
format=self.lhsformat)
if leadnum in in_leads and nprop > 0:
vdaguin_sp = transf.T * sp.csc_matrix(
-np.dot(svd_v, u_in))
# defer formation of the real matrix until the proper
# system size is known
rhs.append((vdaguin_sp, ulinv_in))
else: else:
rhs.append(None) prop, stab = lead.modes(energy, args, params=params)
lead_info.append(prop)
u = stab.vecs
ulinv = stab.vecslmbdainv
nprop = stab.nmodes
svd_v = stab.sqrt_hop
if len(u) == 0:
rhs.append(None)
continue
indices.append(np.arange(lhs.shape[0], lhs.shape[0] + nprop))
u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:]
u_in, ulinv_in = u[:, :nprop], ulinv[:, :nprop]
# Construct a matrix of 1's that translates the
# inter-cell hopping to a proper hopping
# from the system to the lead.
iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1])
for i in interface)]
n_lead_orbs = svd_v.shape[0]
if n_lead_orbs != len(iface_orbs):
msg = ('Lead {0} has hopping with dimensions '
'incompatible with its interface dimension.')
raise ValueError(msg.format(leadnum))
coords = np.r_[[np.arange(len(iface_orbs))], [iface_orbs]]
transf = sp.csc_matrix((np.ones(len(iface_orbs)), coords),
shape=(iface_orbs.size, lhs.shape[0]))
v_sp = sp.csc_matrix(svd_v.T.conj()) * transf
vdaguout_sp = (transf.T *
sp.csc_matrix(np.dot(svd_v, u_out)))
lead_mat = - ulinv_out
lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]],
format=self.lhsformat)
if leadnum in in_leads and nprop > 0:
vdaguin_sp = transf.T * sp.csc_matrix(
-np.dot(svd_v, u_in))
# defer formation of the real matrix until the proper
# system size is known
rhs.append((vdaguin_sp, ulinv_in))
else:
rhs.append(None)
else: else:
sigma = np.asarray(lead.selfenergy(energy, args, params=params)) sigma = np.asarray(lead.selfenergy(energy, args, params=params))
lead_info.append(sigma) lead_info.append(sigma)
...@@ -359,6 +383,15 @@ class SparseSolver(metaclass=abc.ABCMeta): ...@@ -359,6 +383,15 @@ class SparseSolver(metaclass=abc.ABCMeta):
if len(in_leads) == 0 or len(out_leads) == 0: if len(in_leads) == 0 or len(out_leads) == 0:
raise ValueError("No output is requested.") raise ValueError("No output is requested.")
for direction, leads in [("in", in_leads), ("out", out_leads)]:
for lead in leads:
if system.is_selfenergy_lead(syst.leads[lead]):
raise ValueError(
f"lead {lead} is only defined by a self-energy, "
"so we cannot calculate scattering matrix elements for it. "
f"Specify '{direction}_leads' without '{lead}'."
)
linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args, linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, args,
check_hermiticity, False, check_hermiticity, False,
params=params) params=params)
...@@ -519,7 +552,7 @@ class SparseSolver(metaclass=abc.ABCMeta): ...@@ -519,7 +552,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
"is not implemented yet.") "is not implemented yet.")
for lead in syst.leads: for lead in syst.leads:
if not hasattr(lead, 'modes') and hasattr(lead, 'selfenergy'): if system.is_selfenergy_lead(lead):
# TODO: fix this # TODO: fix this
raise NotImplementedError("ldos for leads with only " raise NotImplementedError("ldos for leads with only "
"self-energy is not implemented yet.") "self-energy is not implemented yet.")
...@@ -604,21 +637,27 @@ class WaveFunction: ...@@ -604,21 +637,27 @@ class WaveFunction:
def __init__(self, solver, sys, energy, args, check_hermiticity, params): def __init__(self, solver, sys, energy, args, check_hermiticity, params):
syst = sys # ensure consistent naming across function bodies syst = sys # ensure consistent naming across function bodies
ensure_isinstance(syst, system.System) ensure_isinstance(syst, system.System)
for lead in syst.leads:
if not hasattr(lead, 'modes'): modes_leads = [
# TODO: figure out what to do with self-energies. i for i, l in enumerate(syst.leads)
msg = ('Wave functions for leads with only self-energy' if not system.is_selfenergy_lead(l)
' are not available yet.') ]
raise NotImplementedError(msg)
linsys = solver._make_linear_sys(syst, range(len(syst.leads)), energy, linsys = solver._make_linear_sys(syst, modes_leads, energy,
args, check_hermiticity, args, check_hermiticity,
params=params)[0] params=params)[0]
self.modes_leads = modes_leads
self.solve = solver._solve_linear_sys self.solve = solver._solve_linear_sys
self.rhs = linsys.rhs self.rhs = linsys.rhs
self.factorized_h = solver._factorized(linsys.lhs) self.factorized_h = solver._factorized(linsys.lhs)
self.num_orb = linsys.num_orb self.num_orb = linsys.num_orb
def __call__(self, lead): def __call__(self, lead):
if lead not in self.modes_leads:
msg = ('Scattering wave functions for leads with only '
'self-energy are undefined.')
raise ValueError(msg)
result = self.solve(self.factorized_h, self.rhs[lead], result = self.solve(self.factorized_h, self.rhs[lead],
slice(self.num_orb)) slice(self.num_orb))
return result.transpose() return result.transpose()
...@@ -769,6 +808,8 @@ class SMatrix(BlockResult): ...@@ -769,6 +808,8 @@ class SMatrix(BlockResult):
matrix. matrix.
lead_info : list of data lead_info : list of data
a list containing `kwant.physics.PropagatingModes` for each lead. a list containing `kwant.physics.PropagatingModes` for each lead.
If a lead is a selfenergy lead, then the corresponding entry
in lead_info is a selfenergy.
out_leads, in_leads : sequence of integers out_leads, in_leads : sequence of integers
indices of the leads where current is extracted (out) or injected indices of the leads where current is extracted (out) or injected
(in). Only those are listed for which SMatrix contains the (in). Only those are listed for which SMatrix contains the
...@@ -777,7 +818,22 @@ class SMatrix(BlockResult): ...@@ -777,7 +818,22 @@ class SMatrix(BlockResult):
def __init__(self, data, lead_info, out_leads, in_leads, def __init__(self, data, lead_info, out_leads, in_leads,
current_conserving=False): current_conserving=False):
sizes = [len(i.momenta) // 2 for i in lead_info] # An equivalent condition to this is checked in 'Solver.smatrix',
# but we add it here again as a sanity check. If 'lead_info' is not
# a 'PropagatingModes' that means that the corresponding lead is a
# selfenergy lead. Scattering matrix elements to/from a selfenergy lead
# are not defined.
assert not any(
not isinstance(info, physics.PropagatingModes)
and leadnum in (out_leads + in_leads)
for leadnum, info in enumerate(lead_info)
)
# 'getattr' in order to handle self-energy leads, for which
# 'info' is just a matrix (so no 'momenta').
sizes = [
len(getattr(info, "momenta", [])) // 2 for info in lead_info
]
super().__init__(data, lead_info, out_leads, in_leads, super().__init__(data, lead_info, out_leads, in_leads,
sizes, current_conserving) sizes, current_conserving)
# in_offsets marks beginnings and ends of blocks in the scattering # in_offsets marks beginnings and ends of blocks in the scattering
......
...@@ -84,6 +84,28 @@ def ldos(solver): ...@@ -84,6 +84,28 @@ def ldos(solver):
return solver.ldos return solver.ldos
# 'scope="function"' means that mutating the returned Builder
# inside tests is won't affect other tests.
@pytest.fixture(scope="function")
def twolead_builder():
rng = ensure_rng(4)
system = kwant.Builder()
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
for b, site in [(system, chain(0)), (system, chain(1)),
(left_lead, chain(0)), (right_lead, chain(0))]:
h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
h += h.conjugate().transpose()
b[site] = h
for b, hopp in [(system, (chain(0), chain(1))),
(left_lead, (chain(0), chain(1))),
(right_lead, (chain(0), chain(1)))]:
b[hopp] = (10 * rng.random_sample((n, n)) +
1j * rng.random_sample((n, n)))
system.attach_lead(left_lead)
system.attach_lead(right_lead)
return system
n = 5 n = 5
chain = kwant.lattice.chain(norbs=n) chain = kwant.lattice.chain(norbs=n)
sq = square = kwant.lattice.square(norbs=n) sq = square = kwant.lattice.square(norbs=n)
...@@ -111,24 +133,8 @@ def assert_modes_equal(modes1, modes2): ...@@ -111,24 +133,8 @@ def assert_modes_equal(modes1, modes2):
# Test output sanity: that an error is raised if no output is requested, # Test output sanity: that an error is raised if no output is requested,
# and that solving for a subblock of a scattering matrix is the same as taking # and that solving for a subblock of a scattering matrix is the same as taking
# a subblock of the full scattering matrix. # a subblock of the full scattering matrix.
def test_output(smatrix): def test_output(twolead_builder, smatrix):
rng = ensure_rng(3) fsyst = twolead_builder.finalized()
system = kwant.Builder()
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
for b, site in [(system, chain(0)), (system, chain(1)),
(left_lead, chain(0)), (right_lead, chain(0))]:
h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
h += h.conjugate().transpose()
b[site] = h
for b, hopp in [(system, (chain(0), chain(1))),
(left_lead, (chain(0), chain(1))),
(right_lead, (chain(0), chain(1)))]:
b[hopp] = (10 * rng.random_sample((n, n)) +
1j * rng.random_sample((n, n)))
system.attach_lead(left_lead)
system.attach_lead(right_lead)
fsyst = system.finalized()
result1 = smatrix(fsyst) result1 = smatrix(fsyst)
s, modes1 = result1.data, result1.lead_info s, modes1 = result1.data, result1.lead_info
...@@ -425,25 +431,9 @@ def test_many_leads(smatrix, greens_function): ...@@ -425,25 +431,9 @@ def test_many_leads(smatrix, greens_function):
# Test equivalence between self-energy and scattering matrix representations. # Test equivalence between self-energy and scattering matrix representations.
# Also check that transmission works. # Also check that transmission works.
def test_selfenergy(greens_function, smatrix): def test_selfenergy(twolead_builder, greens_function, smatrix):
rng = ensure_rng(4) system = twolead_builder
system = kwant.Builder() fsyst = twolead_builder.finalized()
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
for b, site in [(system, chain(0)), (system, chain(1)),
(left_lead, chain(0)), (right_lead, chain(0))]:
h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
h += h.conjugate().transpose()
b[site] = h
for b, hopp in [(system, (chain(0), chain(1))),
(left_lead, (chain(0), chain(1))),
(right_lead, (chain(0), chain(1)))]:
b[hopp] = (10 * rng.random_sample((n, n)) +
1j * rng.random_sample((n, n)))
system.attach_lead(left_lead)
system.attach_lead(right_lead)
fsyst = system.finalized()
t = smatrix(fsyst, 0, (), [1], [0]).data t = smatrix(fsyst, 0, (), [1], [0]).data
eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose()) eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose())
n_eig = len(eig_should_be) n_eig = len(eig_should_be)
...@@ -471,20 +461,44 @@ def test_selfenergy(greens_function, smatrix): ...@@ -471,20 +461,44 @@ def test_selfenergy(greens_function, smatrix):
raises(ValueError, check_fsyst, fsyst.precalculate(what='modes')) raises(ValueError, check_fsyst, fsyst.precalculate(what='modes'))
def test_selfenergy_reflection(greens_function, smatrix): def test_selfenergy_lead(
rng = ensure_rng(4) twolead_builder, smatrix, greens_function, wave_function, ldos
system = kwant.Builder() ):
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) fsyst = twolead_builder.finalized()
for b, site in [(system, chain(0)), (system, chain(1)), fsyst_se = twolead_builder.finalized()
(left_lead, chain(0))]: fsyst_se.leads[0] = LeadWithOnlySelfEnergy(fsyst.leads[0])
h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
h += h.conjugate().transpose() # We cannot ask for the smatrix between selfenergy leads
b[site] = h assert pytest.raises(ValueError, smatrix, fsyst_se)
for b, hopp in [(system, (chain(0), chain(1))), assert pytest.raises(ValueError, smatrix, fsyst_se, in_leads=[0])
(left_lead, (chain(0), chain(1)))]: assert pytest.raises(ValueError, smatrix, fsyst_se, out_leads=[0])
b[hopp] = (10 * rng.random_sample((n, n)) +
1j * rng.random_sample((n, n))) # The subblocks of the smatrix between leads that have modes
system.attach_lead(left_lead) # *should* be the same.
kwargs = dict(energy=0, in_leads=[1], out_leads=[1])
assert_almost_equal(
smatrix(fsyst_se, **kwargs).data,
smatrix(fsyst, **kwargs).data,
)
psi_se = wave_function(fsyst_se, energy=0)
psi = wave_function(fsyst, energy=0)
# Scattering wavefunctions originating in the non-selfenergy
# lead should be identical.
assert_almost_equal(psi_se(1), psi(1))
# We cannot define scattering wavefunctions originating from
# a selfenergy lead.
assert pytest.raises(ValueError, psi_se, 0)
# We could in principle calculate the ldos using a mixed
# approach where some leads are specified by selfenergies
# and others by modes, but this is not done for now.
assert pytest.raises(NotImplementedError, ldos, fsyst_se, energy=0)
def test_selfenergy_reflection(twolead_builder, greens_function, smatrix):
system = twolead_builder
fsyst = system.finalized() fsyst = system.finalized()
t = smatrix(fsyst, 0, (), [0], [0]) t = smatrix(fsyst, 0, (), [0], [0])
...@@ -577,7 +591,7 @@ def test_wavefunc_ldos_consistency(wave_function, ldos): ...@@ -577,7 +591,7 @@ def test_wavefunc_ldos_consistency(wave_function, ldos):
check(fsyst) check(fsyst)
raises(ValueError, check, syst.precalculate(what='selfenergy')) raises(ValueError, check, syst.precalculate(what='selfenergy'))
syst.leads[0] = LeadWithOnlySelfEnergy(syst.leads[0]) syst.leads[0] = LeadWithOnlySelfEnergy(syst.leads[0])
raises(NotImplementedError, check, syst) raises(ValueError, check, syst)
# We need to keep testing 'args', but we don't want to see # We need to keep testing 'args', but we don't want to see
......
...@@ -955,6 +955,10 @@ def is_vectorized(syst): ...@@ -955,6 +955,10 @@ def is_vectorized(syst):
return isinstance(syst, (FiniteVectorizedSystem, InfiniteVectorizedSystem)) return isinstance(syst, (FiniteVectorizedSystem, InfiniteVectorizedSystem))
def is_selfenergy_lead(lead):
return hasattr(lead, "selfenergy") and not hasattr(lead, "modes")
def _normalize_matrix_blocks(blocks, expected_shape, *, calling_function=None): def _normalize_matrix_blocks(blocks, expected_shape, *, calling_function=None):
"""Normalize a sequence of matrices into a single 3D numpy array """Normalize a sequence of matrices into a single 3D numpy array
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment