Commit d29bc4be authored by Bas Nijholt's avatar Bas Nijholt
Browse files

fix setup_linsys algo without stabilization and add a test

parent a82d9e88
......@@ -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]
......
......@@ -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)
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