test_solvers.py 23.3 KB
Newer Older
1
# Copyright 2011-2020 Kwant authors.
2
#
Christoph Groth's avatar
Christoph Groth committed
3
4
# This file is part of Kwant.  It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
5
# http://kwant-project.org/license.  A list of Kwant authors can be found in
Christoph Groth's avatar
Christoph Groth committed
6
# the file AUTHORS.rst at the top-level directory of this distribution and at
7
8
# http://kwant-project.org/authors.

9
import itertools
10
from math import cos, sin
11
import numpy as np
12
import pytest
Anton Akhmerov's avatar
Anton Akhmerov committed
13
from pytest import raises
14
from numpy.testing import assert_almost_equal
15

16
import kwant
17
from kwant._common import ensure_rng
18

19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
import kwant.solvers.sparse
try:
    import kwant.solvers.mumps
except ImportError:
    no_mumps = True
else:
    no_mumps = False

mumps_solver_options = [
    {},
    {'nrhs': 1},
    {'nrhs': 10},
    {'nrhs': 1, 'ordering': 'amd'},
    {'nrhs': 10, 'sparse_rhs': True},
    {'nrhs': 2, 'ordering': 'amd', 'sparse_rhs': True},
]

sparse_solver_options = [
    {},
]

40
solvers = list(itertools.chain(
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    [("mumps", opts) for opts in mumps_solver_options],
    [("sparse", opts) for opts in sparse_solver_options],
))


def solver_id(s):
    solver_name, opts = s
    args = ", ".join(f"{k}={repr(v)}" for k, v in opts.items())
    return f"{solver_name}({args})"


@pytest.fixture(scope="function", params=solvers, ids=map(solver_id, solvers))
def solver(request):
    solver_name, solver_opts = request.param
    if solver_name == "mumps" and no_mumps:
        pytest.skip("MUMPS not installed")
    solver = getattr(kwant.solvers, solver_name).Solver()
    # The sparse solver does not have "options" or "reset_options".
    if solver_name != "sparse":
        solver.options(**solver_opts)
    return solver


# These fixtures are repetitive, but doing this programatically is cryptic
# and we're not going to add any more (unless we add new methods to Solvers)

@pytest.fixture
def smatrix(solver):
    return solver.smatrix


@pytest.fixture
def greens_function(solver):
    return solver.greens_function


@pytest.fixture
def wave_function(solver):
    return solver.wave_function


@pytest.fixture
def ldos(solver):
    return solver.ldos


87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
# '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

109
n = 5
Joseph Weston's avatar
Joseph Weston committed
110
111
chain = kwant.lattice.chain(norbs=n)
sq = square = kwant.lattice.square(norbs=n)
112

113

114
class LeadWithOnlySelfEnergy:
115
116
    def __init__(self, lead):
        self.lead = lead
117
        self.parameters = frozenset()
118

119
    def selfenergy(self, energy, args=(), *, params=None):
120
        assert args == ()
121
        assert params == None
122
        return self.lead.selfenergy(energy)
123
124


125
126
127
128
129
130
131
132
def assert_modes_equal(modes1, modes2):
    assert_almost_equal(modes1.velocities, modes2.velocities)
    assert_almost_equal(modes1.momenta, modes2.momenta)
    vecs1, vecs2 = modes1.wave_functions, modes2.wave_functions
    assert_almost_equal(np.abs(np.sum(vecs1/vecs2, axis=0)),
                        vecs1.shape[0])


133
134
135
# 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
# a subblock of the full scattering matrix.
136
137
def test_output(twolead_builder, smatrix):
    fsyst = twolead_builder.finalized()
138

139
    result1 = smatrix(fsyst)
Michael Wimmer's avatar
Michael Wimmer committed
140
    s, modes1 = result1.data, result1.lead_info
141
    assert s.shape == 2 * (sum(len(i.momenta) for i in modes1) // 2,)
142
    s1 = result1.submatrix(1, 0)
143
    result2 = smatrix(fsyst, 0, (), [1], [0])
Michael Wimmer's avatar
Michael Wimmer committed
144
    s2, modes2 = result2.data, result2.lead_info
145
146
    assert s2.shape == (len(modes2[1].momenta) // 2,
                        len(modes2[0].momenta) // 2)
147
    assert_almost_equal(abs(s1), abs(s2))
148
    assert_almost_equal(np.dot(s.T.conj(), s),
149
                        np.identity(s.shape[0]))
Anton Akhmerov's avatar
Anton Akhmerov committed
150
    raises(ValueError, smatrix, fsyst, out_leads=[])
151
152
153
    modes = smatrix(fsyst).lead_info
    h = fsyst.leads[0].cell_hamiltonian()
    t = fsyst.leads[0].inter_cell_hopping()
154
    modes1 = kwant.physics.modes(h, t)[0]
155
156
    h = fsyst.leads[1].cell_hamiltonian()
    t = fsyst.leads[1].inter_cell_hopping()
157
    modes2 = kwant.physics.modes(h, t)[0]
158
159
    assert_modes_equal(modes1, modes[0])
    assert_modes_equal(modes2, modes[1])
160
161
162


# Test that a system with one lead has unitary scattering matrix.
163
def test_one_lead(smatrix):
164
    rng = ensure_rng(3)
165
    system = kwant.Builder()
166
    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
167
168
    for b, site in [(system, chain(0)), (system, chain(1)),
                    (system, chain(2)), (lead, chain(0))]:
169
        h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
170
171
172
173
174
        h += h.conjugate().transpose()
        b[site] = h
    for b, hopp in [(system, (chain(0), chain(1))),
                    (system, (chain(1), chain(2))),
                    (lead, (chain(0), chain(1)))]:
175
176
        b[hopp] = (10 * rng.random_sample((n, n)) +
                   1j * rng.random_sample((n, n)))
177
    system.attach_lead(lead)
178
    fsyst = system.finalized()
179

180
181
    for syst in (fsyst, fsyst.precalculate(), fsyst.precalculate(what='all')):
        s = smatrix(syst).data
182
183
        assert_almost_equal(np.dot(s.conjugate().transpose(), s),
                            np.identity(s.shape[0]))
184

Anton Akhmerov's avatar
Anton Akhmerov committed
185
    raises(ValueError, smatrix, fsyst.precalculate(what='selfenergy'))
186
187
188

# Test that a system with one lead with no propagating modes has a
# 0x0 S-matrix.
189
def test_smatrix_shape(smatrix):
Joseph Weston's avatar
Joseph Weston committed
190
191
    chain = kwant.lattice.chain(norbs=1)

192
    system = kwant.Builder()
193
194
    lead0 = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
    lead1 = kwant.Builder(kwant.TranslationalSymmetry((1,)))
195
196
197
198
199
200
201
202
203
204
205
206
207
    for b, site in [(system, chain(0)), (system, chain(1)),
                    (system, chain(2))]:
        b[site] = 2
    lead0[chain(0)] = lambda site: lead0_val
    lead1[chain(0)] = lambda site: lead1_val

    for b, hopp in [(system, (chain(0), chain(1))),
                    (system, (chain(1), chain(2))),
                    (lead0, (chain(0), chain(1))),
                    (lead1, (chain(0), chain(1)))]:
        b[hopp] = -1
    system.attach_lead(lead0)
    system.attach_lead(lead1)
208
    fsyst = system.finalized()
209
210
211

    lead0_val = 4
    lead1_val = 4
212
    s = smatrix(fsyst, 1.0, (), [1], [0]).data
213
214
215
216
    assert s.shape == (0, 0)

    lead0_val = 2
    lead1_val = 2
217
    s = smatrix(fsyst, 1.0, (), [1], [0]).data
218
219
220
221
    assert s.shape == (1, 1)

    lead0_val = 4
    lead1_val = 2
222
    s = smatrix(fsyst, 1.0, (), [1], [0]).data
223
224
225
226
    assert s.shape == (1, 0)

    lead0_val = 2
    lead1_val = 4
227
    s = smatrix(fsyst, 1.0, (), [1], [0]).data
228
229
    assert s.shape == (0, 1)

230

231
232
# Test that a translationally invariant system with two leads has only
# transmission and that transmission does not mix modes.
233
def test_two_equal_leads(smatrix):
234
235
    def check_fsyst(fsyst):
        sol = smatrix(fsyst)
Michael Wimmer's avatar
Michael Wimmer committed
236
        s, leads = sol.data, sol.lead_info
237
238
        assert_almost_equal(np.dot(s.conjugate().transpose(), s),
                            np.identity(s.shape[0]))
239
240
        n_modes = len(leads[0].momenta) // 2
        assert len(leads[1].momenta) // 2 == n_modes
241
242
243
244
245
246
        assert_almost_equal(s[: n_modes, : n_modes], 0)
        t_elements = np.sort(abs(np.asarray(s[n_modes :, : n_modes])),
                             axis=None)
        t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
        assert_almost_equal(t_elements, t_el_should_be)
        assert_almost_equal(sol.transmission(1,0), n_modes)
247
    rng = ensure_rng(11)
248
    system = kwant.Builder()
249
    lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
250
    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
251
252
    h += h.conjugate().transpose()
    h *= 0.8
253
    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
254
255
256
257
    lead[chain(0)] = system[chain(0)] = h
    lead[chain(0), chain(1)] = t
    system.attach_lead(lead)
    system.attach_lead(lead.reversed())
258
259
260
    fsyst = system.finalized()
    for syst in (fsyst, fsyst.precalculate(), fsyst.precalculate(what='all')):
        check_fsyst(syst)
Anton Akhmerov's avatar
Anton Akhmerov committed
261
    raises(ValueError, check_fsyst, fsyst.precalculate(what='selfenergy'))
262
263
264
265
266
267
268

    # Test the same, but with a larger scattering region.
    system = kwant.Builder()
    system[[chain(0), chain(1)]] = h
    system[chain(0), chain(1)] = t
    system.attach_lead(lead)
    system.attach_lead(lead.reversed())
269
270
271
    fsyst = system.finalized()
    for syst in (fsyst, fsyst.precalculate(), fsyst.precalculate(what='all')):
        check_fsyst(syst)
Anton Akhmerov's avatar
Anton Akhmerov committed
272
    raises(ValueError, check_fsyst, fsyst.precalculate(what='selfenergy'))
273
274


275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
# 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)


298
# Test a more complicated graph with non-singular hopping.
299
def test_graph_system(smatrix):
300
    rng = ensure_rng(11)
301
    system = kwant.Builder()
302
    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
303
    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
304
305
    h += h.conjugate().transpose()
    h *= 0.8
306
307
    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
    t1 = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
308
309
    lead[sq(0, 0)] = system[[sq(0, 0), sq(1, 0)]] = h
    lead[sq(0, 1)] = system[[sq(0, 1), sq(1, 1)]] = 4 * h
310
    for builder in [system, lead]:
311
312
313
        builder[sq(0, 0), sq(1, 0)] = t
        builder[sq(0, 1), sq(1, 0)] = t1
        builder[sq(0, 1), sq(1, 1)] = 1.1j * t1
314
315
    system.attach_lead(lead)
    system.attach_lead(lead.reversed())
316
    fsyst = system.finalized()
317

318
    result = smatrix(fsyst)
Michael Wimmer's avatar
Michael Wimmer committed
319
    s, leads = result.data, result.lead_info
320
321
    assert_almost_equal(np.dot(s.conjugate().transpose(), s),
                        np.identity(s.shape[0]))
322
    n_modes = len(leads[0].momenta) // 2
323
    assert len(leads[1].momenta) // 2 == n_modes
324
    assert_almost_equal(s[: n_modes, : n_modes], 0)
325
    t_elements = np.sort(abs(np.asarray(s[n_modes:, :n_modes])),
326
327
328
329
330
331
                         axis=None)
    t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
    assert_almost_equal(t_elements, t_el_should_be)


# Test a system with singular hopping.
332
def test_singular_graph_system(smatrix):
333
    rng = ensure_rng(11)
334
335

    system = kwant.Builder()
336
    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
337
    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
338
339
    h += h.conjugate().transpose()
    h *= 0.8
340
341
    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
    t1 = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
342
343
    lead[sq(0, 0)] = system[[sq(0, 0), sq(1, 0)]] = h
    lead[sq(0, 1)] = system[[sq(0, 1), sq(1, 1)]] = 4 * h
344
    for builder in [system, lead]:
345
346
        builder[sq(0, 0), sq(1, 0)] = t
        builder[sq(0, 1), sq(1, 0)] = t1
347
348
    system.attach_lead(lead)
    system.attach_lead(lead.reversed())
349
    fsyst = system.finalized()
350

351
    result = smatrix(fsyst)
Michael Wimmer's avatar
Michael Wimmer committed
352
    s, leads = result.data, result.lead_info
353
354
    assert_almost_equal(np.dot(s.conjugate().transpose(), s),
                        np.identity(s.shape[0]))
355
356
    n_modes = len(leads[0].momenta) // 2
    assert len(leads[1].momenta) // 2 == n_modes
357
358
359
360
361
362
363
    assert_almost_equal(s[: n_modes, : n_modes], 0)
    t_elements = np.sort(abs(np.asarray(s[n_modes :, : n_modes])),
                         axis=None)
    t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
    assert_almost_equal(t_elements, t_el_should_be)


364
# This test features inside the cell Hamiltonian a hopping matrix with more
Michael Wimmer's avatar
Michael Wimmer committed
365
# zero eigenvalues than the lead hopping matrix. Older version of the
366
# sparse solver failed here.
367
def test_tricky_singular_hopping(smatrix):
Joseph Weston's avatar
Joseph Weston committed
368
369
    sq = kwant.lattice.square(norbs=1)

370
    system = kwant.Builder()
371
    lead = kwant.Builder(kwant.TranslationalSymmetry((4, 0)))
372

373
    interface = []
Joseph Weston's avatar
Joseph Weston committed
374
    for i in range(n):
375
        site = sq(-1, i)
376
        interface.append(site)
377
        system[site] = 0
Joseph Weston's avatar
Joseph Weston committed
378
        for j in range(4):
379
            lead[sq(j, i)] = 0
Joseph Weston's avatar
Joseph Weston committed
380
    for i in range(n-1):
381
        system[sq(-1, i), sq(-1, i+1)] = -1
Joseph Weston's avatar
Joseph Weston committed
382
        for j in range(4):
383
            lead[sq(j, i), sq(j, i+1)] = -1
Joseph Weston's avatar
Joseph Weston committed
384
385
    for i in range(n):
        for j in range(4):
386
387
            lead[sq(j, i), sq(j+1, i)] = -1
    del lead[sq(1, 0), sq(2, 0)]
388

389
    system.leads.append(kwant.builder.BuilderLead(lead, interface))
390
    fsyst = system.finalized()
391

392
    s = smatrix(fsyst, -1.3).data
393
394
395
396
    assert_almost_equal(np.dot(s.conjugate().transpose(), s),
                        np.identity(s.shape[0]))


397
398
# Test the consistency of transmission and conductance_matrix for a four-lead
# system without time-reversal symmetry.
399
400
def test_many_leads(smatrix, greens_function):
    factories = (greens_function, smatrix)
Joseph Weston's avatar
Joseph Weston committed
401
    sq = kwant.lattice.square(norbs=1)
402
403
404
405
406
407
408
409
410
411
    E=2.1
    B=0.01

    def phase(a, b):
        ap = a.pos
        bp = b.pos
        phase = -B * (0.5 * (ap[1] + bp[1]) * (bp[0] - ap[0]))
        return -complex(cos(phase), sin(phase))

    # Build a square system with four leads and a hole in the center.
412
413
    syst = kwant.Builder()
    syst[(sq(x, y) for x in range(-4, 4) for y in range(-4, 4)
414
         if x**2 + y**2 >= 2)] = 3
415
    syst[sq.neighbors()] = phase
Joseph Weston's avatar
Joseph Weston committed
416
    for r in [range(-4, -1), range(4)]:
417
418
419
        lead = kwant.Builder(kwant.TranslationalSymmetry([-1, 0]))
        lead[(sq(0, y) for y in r)] = 4
        lead[sq.neighbors()] = phase
420
421
422
        syst.attach_lead(lead)
        syst.attach_lead(lead.reversed())
    syst = syst.finalized()
423

Joseph Weston's avatar
Joseph Weston committed
424
    r4 = list(range(4))
425
    br = greens_function(syst, E, out_leads=r4, in_leads=r4)
426
427
428
429
430
431
432
433
434
435
436
437
438
439
    trans = np.array([[br._transmission(i, j) for j in r4] for i in r4])
    cmat = br.conductance_matrix()
    assert_almost_equal(cmat.sum(axis=0), [0] * 4)
    assert_almost_equal(cmat.sum(axis=1), [0] * 4)
    for i in r4:
        for j in r4:
            assert_almost_equal(
                (br.num_propagating(i) if i == j else 0) - cmat[i, j],
                trans[i, j])

    for out_leads, in_leads in [(r4, r4), ((1, 2, 3), r4), (r4, (0, 2, 3)),
                                ((1, 3), (1, 2, 3)), ((0, 2), (0, 1, 2)),
                                ((0, 1,), (1, 2)), ((3,), (3,))]:
        for f in factories:
440
            br = f(syst, E, out_leads=out_leads, in_leads=in_leads)
441
442
443
444
445
446
447
448
449
            if len(out_leads) == 3:
                out_leads = r4
            if len(in_leads) == 3:
                in_leads = r4
            for i in r4:
                for j in r4:
                    if i in out_leads and j in in_leads:
                        assert_almost_equal(br.transmission(i, j), trans[i, j])
                    else:
Anton Akhmerov's avatar
Anton Akhmerov committed
450
                        raises(ValueError, br.transmission, i, j)
451
452
453
454
            if len(out_leads) == len(in_leads) == 4:
                assert_almost_equal(br.conductance_matrix(), cmat)


455
456
# Test equivalence between self-energy and scattering matrix representations.
# Also check that transmission works.
457
458
459
def test_selfenergy(twolead_builder, greens_function, smatrix):
    system = twolead_builder
    fsyst = twolead_builder.finalized()
460
    t = smatrix(fsyst, 0, (), [1], [0]).data
461
462
463
    eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose())
    n_eig = len(eig_should_be)

464
465
    def check_fsyst(fsyst):
        sol = greens_function(fsyst, 0, (), [1], [0])
466
467
468
469
470
471
472
473
        ttdagnew = sol._a_ttdagger_a_inv(1, 0)
        eig_are = np.linalg.eigvals(ttdagnew)
        t_should_be = np.sum(eig_are)
        assert_almost_equal(eig_are.imag, 0)
        assert_almost_equal(np.sort(eig_are.real)[-n_eig:],
                            np.sort(eig_should_be.real))
        assert_almost_equal(t_should_be, sol.transmission(1, 0))

474
475
    fsyst.leads[1] = LeadWithOnlySelfEnergy(fsyst.leads[1])
    check_fsyst(fsyst)
476

477
478
    fsyst.leads[0] = LeadWithOnlySelfEnergy(fsyst.leads[0])
    check_fsyst(fsyst)
479

480
481
482
483
    fsyst = system.finalized()
    for syst in (fsyst, fsyst.precalculate(what='selfenergy'),
                fsyst.precalculate(what='all')):
        check_fsyst(syst)
Anton Akhmerov's avatar
Anton Akhmerov committed
484
    raises(ValueError, check_fsyst, fsyst.precalculate(what='modes'))
485

486

487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
def test_selfenergy_lead(
    twolead_builder, smatrix, greens_function, wave_function, ldos
):
    fsyst = twolead_builder.finalized()
    fsyst_se = twolead_builder.finalized()
    fsyst_se.leads[0] = LeadWithOnlySelfEnergy(fsyst.leads[0])

    # We cannot ask for the smatrix between selfenergy leads
    assert pytest.raises(ValueError, smatrix, fsyst_se)
    assert pytest.raises(ValueError, smatrix, fsyst_se, in_leads=[0])
    assert pytest.raises(ValueError, smatrix, fsyst_se, out_leads=[0])

    # The subblocks of the smatrix between leads that have modes
    # *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)


523
524
def test_selfenergy_reflection(twolead_builder, greens_function, smatrix):
    system = twolead_builder
525
    fsyst = system.finalized()
526

527
    t = smatrix(fsyst, 0, (), [0], [0])
528

529
530
    fsyst.leads[0] = LeadWithOnlySelfEnergy(fsyst.leads[0])
    sol = greens_function(fsyst, 0, (), [0], [0])
531
532
    assert_almost_equal(sol.transmission(0,0), t.transmission(0,0))

533
534
535
536
    fsyst = system.finalized()
    for syst in (fsyst.precalculate(what='selfenergy'),
                fsyst.precalculate(what='all')):
        sol = greens_function(syst, 0, (), [0], [0])
537
        assert_almost_equal(sol.transmission(0,0), t.transmission(0,0))
Anton Akhmerov's avatar
Anton Akhmerov committed
538
    raises(ValueError, greens_function, fsyst.precalculate(what='modes'),
539
540
                  0, (), [0], [0])

541

542
def test_very_singular_leads(smatrix):
543
    syst = kwant.Builder()
544
    chain = kwant.lattice.chain(norbs=2)
545
546
    left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
    right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
547
    syst[chain(0)] = left_lead[chain(0)] = right_lead[chain(0)] = np.identity(2)
548
549
    left_lead[chain(0), chain(1)] = np.zeros((2, 2))
    right_lead[chain(0), chain(1)] = np.identity(2)
550
551
552
553
    syst.attach_lead(left_lead)
    syst.attach_lead(right_lead)
    fsyst = syst.finalized()
    leads = smatrix(fsyst).lead_info
554
    assert [len(i.momenta) for i in leads] == [0, 4]
555

556

557
def test_ldos(ldos):
558
    syst = kwant.Builder()
559
    chain = kwant.lattice.chain(norbs=1)
560
    lead = kwant.Builder(kwant.TranslationalSymmetry(chain.vec((1,))))
561
562
563
564
565
566
567
568
569
    syst[chain(0)] = syst[chain(1)] = lead[chain(0)] = 0
    syst[chain(0), chain(1)] = lead[chain(0), chain(1)] = 1
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())
    fsyst = syst.finalized()

    for finsyst in (fsyst, fsyst.precalculate(what='modes'),
                   fsyst.precalculate(what='all')):
        assert_almost_equal(ldos(finsyst, 0),
570
                            np.array([1, 1]) / (2 * np.pi))
Anton Akhmerov's avatar
Anton Akhmerov committed
571
    raises(ValueError, ldos, fsyst.precalculate(what='selfenergy'), 0)
572
    fsyst.leads[0] = LeadWithOnlySelfEnergy(fsyst.leads[0])
Anton Akhmerov's avatar
Anton Akhmerov committed
573
    raises(NotImplementedError, ldos, fsyst, 0)
574
575


576
def test_wavefunc_ldos_consistency(wave_function, ldos):
577
578
579
    L = 2
    W = 3

580
    rng = ensure_rng(31)
581
    syst = kwant.Builder()
582
583
    left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
    top_lead = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
584
    for b, sites in [(syst, [square(x, y)
585
586
587
588
                               for x in range(L) for y in range(W)]),
                     (left_lead, [square(0, y) for y in range(W)]),
                     (top_lead, [square(x, 0) for x in range(L)])]:
        for site in sites:
589
            h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
590
591
            h += h.conjugate().transpose()
            b[site] = h
592
        for hopping_kind in square.neighbors():
593
            for hop in hopping_kind(b):
594
595
                b[hop] = (10 * rng.random_sample((n, n)) +
                          1j * rng.random_sample((n, n)))
596
597
598
    syst.attach_lead(left_lead)
    syst.attach_lead(top_lead)
    syst = syst.finalized()
599

600
    def check(syst):
601
        for energy in [0, 1000]:
602
            wf = wave_function(syst, energy)
603
            ldos2 = np.zeros(wf.num_orb, float)
604
            for lead in range(len(syst.leads)):
605
606
607
608
609
                temp = abs(wf(lead))
                temp **= 2
                ldos2 += temp.sum(axis=0)
            ldos2 *= (0.5 / np.pi)

610
            assert_almost_equal(ldos2, ldos(syst, energy))
611

612
613
614
    for fsyst in (syst, syst.precalculate(what='modes'),
                 syst.precalculate(what='all')):
        check(fsyst)
Anton Akhmerov's avatar
Anton Akhmerov committed
615
    raises(ValueError, check, syst.precalculate(what='selfenergy'))
616
    syst.leads[0] = LeadWithOnlySelfEnergy(syst.leads[0])
617
    raises(ValueError, check, syst)
618
619


620
621
622
# We need to keep testing 'args', but we don't want to see
# all the deprecation warnings in the test logs
@pytest.mark.filterwarnings("ignore:.*'args' parameter")
623
624
625
626
627
628
629
630
def test_arg_passing(wave_function, ldos, smatrix):

    def onsite(site, a, b):
        return site.pos[0] + site.pos[1] + a + b

    def hopping(site1, site2, a, b):
        return b - a

Joseph Weston's avatar
Joseph Weston committed
631
    square = kwant.lattice.square(norbs=1)
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
    W = 3
    L = 4

    syst = kwant.Builder()
    syst[(square(i, j) for i in range(L) for j in range(W))] = onsite
    syst[square.neighbors()] = hopping

    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
    lead[(square(0, j) for j in range(W))] = onsite
    lead[square.neighbors()] = hopping

    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    fsyst = syst.finalized()

    # compare results to those when we pass `args` only
    args = (1, 3)
    params = dict(a=1, b=3)
    np.testing.assert_array_equal(
        wave_function(fsyst, args=args)(0),
        wave_function(fsyst, params=params)(0))
    np.testing.assert_array_equal(
        ldos(fsyst, args=args),
        ldos(fsyst, params=params))
    np.testing.assert_array_equal(
        smatrix(fsyst, args=args).data,
        smatrix(fsyst, params=params).data)