test_system.py 10.4 KB
Newer Older
Joseph Weston's avatar
Joseph Weston committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# -*- coding: utf-8 -*-
# Copyright 2016 tkwant authors.
#
# This file is part of tkwant.  It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
# http://kwant-project.org/license.  A list of tkwant authors can be found in
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.

"""Test module for `tkwant.system`"""

import numpy as np
import scipy.sparse as sp
import tinyarray as ta
import kwant
import pytest

18
19
20
from .. import system, leads
from .common import (make_chain, make_simple_lead, make_complex_lead,
                     make_system_with_leads, check_boundary_hamiltonian)
Joseph Weston's avatar
Joseph Weston committed
21
22
23
24
25


def test_orbital_slices():

    def check(fsyst):
26
        slices = (system._get_orbs(fsyst, i) for i in range(len(fsyst.sites)))
Joseph Weston's avatar
Joseph Weston committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
        stop = 0
        for site, (start_orb, stop_orb) in zip(fsyst.sites, slices):
            assert start_orb == stop
            stop += site.family.norbs
            assert stop_orb == stop
            assert site.family.norbs == (stop_orb - start_orb)

    # uniform families
    lat, syst = make_chain(3, norbs=1)
    check(syst.finalized())
    lat, syst = make_chain(3, norbs=2)
    check(syst.finalized())
    # mixed families
    lat = kwant.lattice.chain(name='a', norbs=1)
    lat2 = kwant.lattice.chain(name='b', norbs=2)
    syst = kwant.Builder()
    syst[map(lat, range(3))] = 2
    syst[map(lat2, range(3))] = 2 * np.eye(2)
    check(syst.finalized())


def test_extract_matrix_elements():

50
    def onsite(i):
Joseph Weston's avatar
Joseph Weston committed
51
52
        pass

53
    def hopping(i, j):
Joseph Weston's avatar
Joseph Weston committed
54
55
        pass

56
57
58
59
60
    def td_onsite(i, time):
        pass

    def td_hopping(i, j, time):
        pass
Joseph Weston's avatar
Joseph Weston committed
61
62
63
64
65
66
67
68
69
70
71
72
73

    lat, syst = make_chain(4)
    td_sites = (lat(0), lat(1))
    td_hops = ((lat(2), lat(1)), (lat(1), lat(0)))
    # assign time-dependent parts
    for site in td_sites:
        syst[site] = td_onsite
    for hop in td_hops:
        syst[hop] = td_hopping
    syst[lat(2)] = onsite  # non time-dependent onsite function
    syst[(lat(2), lat(3))] = hopping  # non time-dependent hopping function
    fsyst = syst.finalized()

74
    got = tuple(system.extract_matrix_elements(fsyst, time_name='time'))
Joseph Weston's avatar
Joseph Weston committed
75
76
77
78
79
80
81
82
83
84
85
86
87
    assert len(got) == len(set(got))
    got = set(got)
    for site in td_sites:
        idx = fsyst.sites.index(site)
        assert (idx, idx) in got
    for hop in td_hops:
        a, b = tuple(map(fsyst.sites.index, hop))
        assert (a, b) in got


def test_extract_perturbation():

    def inner_test(fsyst, N):
88
        W = system.extract_perturbation(fsyst, time_name='time', time_start=0, params={'time': 0})
89
90
        H0 = fsyst.hamiltonian_submatrix(params={'time': 0})
        H1 = fsyst.hamiltonian_submatrix(params={'time': 1})
Joseph Weston's avatar
Joseph Weston committed
91
92
        ket = np.random.rand(N) + 1j * np.random.rand(N)
        for sparse in (True, False):
93
94
95
96
            assert np.all(W(-100, {}, sparse=sparse) == np.zeros((N, N)))
            assert np.all(W(0, {}, sparse=sparse) == np.zeros((N, N)))
            assert np.all(W(1, {}, sparse=sparse) == (H1 - H0))
            assert (type(W(1, {}, sparse=sparse)) ==
Joseph Weston's avatar
Joseph Weston committed
97
98
                    (sp.lil_matrix if sparse else np.ndarray))

99
100
        assert np.all(W(-100, {}, ket) == np.zeros_like(ket))
        assert np.all(W(0, {}, ket) == np.zeros_like(ket))
Joseph Weston's avatar
Joseph Weston committed
101
        # calculations not done in same order, only check if close
102
        assert np.allclose(W(1, {}, ket), np.dot(H1 - H0, ket))
Joseph Weston's avatar
Joseph Weston committed
103
        # check with wrong size
104
105
        pytest.raises(ValueError, W, 1, {}, np.zeros((N - 1),))
        pytest.raises(ValueError, W, 1, {}, np.zeros((N + 1),))
Joseph Weston's avatar
Joseph Weston committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

    def td_onsite(i, time):
        return 2 + time * kwant.digest.uniform(i.tag)

    def td_hopping(i, j, time):
        ab = ta.array(sorted((i.tag, j.tag)))
        return -1 + time * kwant.digest.uniform(ab)

    lat = kwant.lattice.chain(norbs=1)
    N = 4

    # test whole system time-dependent
    syst = kwant.Builder()
    syst[(lat(i) for i in range(N))] = td_onsite
    syst[lat.neighbors()] = td_hopping
    inner_test(syst.finalized(), N)
    # check that error is raised if norbs not given
    lat.norbs = None
    pytest.raises(RuntimeError, system.extract_perturbation,
125
                  syst.finalized(), time_name='time', time_start=0)
Joseph Weston's avatar
Joseph Weston committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    lat.norbs = 1

    # test only a partial amount time-dependent
    syst = kwant.Builder()
    syst[(lat(i) for i in range(N))] = td_onsite
    syst[(lat(i) for i in range(0, N, 2))] = 1.
    syst[lat.neighbors()] = -1
    inner_test(syst.finalized(), N)

    # test >1 orbital per site
    sigma_0 = np.eye(2)
    sigma_y = np.array([[0, -1j], [1j, 0]])

    def td_onsite2(i, time):
        U = kwant.digest.uniform(i.tag)
        return 2 + time * U * sigma_y

    def td_hopping2(i, j, time):
        U = kwant.digest.uniform(ta.array(sorted((i.tag, j.tag))))
        return -1 + time * U * sigma_0

    def td_interhopping(i, j, time):
        # *from* 2 orbital lattice *to* 1 orbital lattice
        U = kwant.digest.uniform(ta.array(sorted((i.tag, j.tag))))
        return 1j * time * U * np.array([[1, 1]])

    lat2 = kwant.lattice.chain(name='b', norbs=2)

    syst = kwant.Builder()
    syst[(lat(i) for i in range(N))] = td_onsite
    syst[(lat2(i) for i in range(N))] = td_onsite2
    syst[lat.neighbors()] = td_hopping
    syst[lat2.neighbors()] = td_hopping2
    syst[kwant.builder.HoppingKind((0,), lat, lat2)] = td_interhopping
    inner_test(syst.finalized(), 2 * N + N)
161
162
163
164
165
166
167
168
169
170
171
172
173


@pytest.mark.parametrize('lead_maker', [make_simple_lead, make_complex_lead])
def test_hamiltonian(lead_maker):
    ncells = 10
    strength = 10
    degree = 6
    # construct
    syst = make_system_with_leads(kwant.lattice.square(norbs=1),
                                  lead_maker)
    fsyst = syst.finalized()
    boundaries = [leads.SimpleBoundary(10),
                  leads.MonomialAbsorbingBoundary(ncells, strength, degree)]
174
    Hext = system.hamiltonian_with_boundaries(fsyst, boundaries, params={})
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221

    # test shapes
    norbs_central = len(fsyst.sites)
    norbs_leads = [lead.cell_size * ncells for lead in fsyst.leads]
    assert Hext.hamiltonian.shape[0] == Hext.hamiltonian.shape[1]
    assert Hext.hamiltonian.shape[0] == norbs_central + sum(norbs_leads)

    # test absorbing regions and coupling

    def absorbing(i):
        n = degree
        return -1j * np.eye(norbs) * ((n + 1) * strength *
                                      i**n / ncells**(n + 1))

    loop = zip(fsyst.leads, fsyst.lead_interfaces, Hext.boundary_slices,
               boundaries)
    for lead, lead_interface, slc, bdy in loop:
        norbs = lead.cell_size
        norbs_iface = lead.graph.num_nodes - lead.cell_size
        V = lead.inter_cell_hopping()
        V_dag = V.conjugate().transpose()
        if isinstance(bdy, leads.MonomialAbsorbingBoundary):
            def onsite(i):
                return lead.cell_hamiltonian() + absorbing(i)
        else:
            def onsite(i):
                return lead.cell_hamiltonian()
        check_boundary_hamiltonian(Hext.hamiltonian[slc, slc],
                                   norbs, norbs_iface, onsite, lambda i: V)

        # test coupling -- uses fact that norbs=1
        for i, iface_site in enumerate(lead_interface):
            bdy_iface = slice(slc.start, slc.start + norbs)
            iface_slice = slice(iface_site, iface_site + 1)
            assert np.allclose(Hext.hamiltonian[bdy_iface, iface_slice].todense(),
                               V[:, i:(i + 1)])
            assert np.allclose(Hext.hamiltonian[iface_slice, bdy_iface].todense(),
                               V_dag[i:(i + 1), :])

    # test is_valid and should_continue
    psi_test = np.empty(Hext.hamiltonian.shape[0], dtype=complex)
    assert Hext.solution_is_valid(psi_test)  # should always return True
    vmax = np.linalg.norm(fsyst.leads[0].inter_cell_hopping(), ord=2)
    tmax = ncells / (2 * vmax)
    eps = 1E-8
    assert Hext.time_is_valid(tmax - eps)
    assert not Hext.time_is_valid(tmax + eps)
222
223


224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
@pytest.mark.parametrize('lead_maker', [make_simple_lead])
def test_hamiltonian_with_boundaries_exceptions(lead_maker):

    # construct
    syst = make_system_with_leads(kwant.lattice.square(norbs=1), lead_maker)
    fsyst = syst.finalized()

    # too less boundary conditions
    boundaries = [leads.SimpleBoundary(10)]
    with pytest.raises(ValueError) as exc:
        system.hamiltonian_with_boundaries(fsyst, boundaries, params={})
    assert 'Number of leads= 2 does not match the number of boundaries provided= 1' in str(exc.value)

    # too many boundary conditions
    boundaries = [leads.SimpleBoundary(10)] * 3
    with pytest.raises(ValueError) as exc:
        system.hamiltonian_with_boundaries(fsyst, boundaries, params={})
    assert 'Number of leads= 2 does not match the number of boundaries provided= 3' in str(exc.value)


244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def test_is_time_dependent_function():
    def time_dependent(zeit):
        pass
    not_a_function = 'time'
    assert system._is_time_dependent_function(time_dependent, 'zeit')
    assert system._is_time_dependent_function(time_dependent, 'time') is False
    assert system._is_time_dependent_function(not_a_function, 'time') is False


def test_add_time_to_params():
    params_example = {'some_key': 2}
    time_name = 'my_time'
    time_value = 42

    # check against reference result
    tparams = system.add_time_to_params(params_example, time_name=time_name,
                                        time=time_value)
    assert isinstance(tparams, dict)
    assert tparams == {'some_key': 2, 'my_time': 42}

    # check for possible inputs
    for params in [None, {}, params_example]:
        tparams = system.add_time_to_params(params, time_name=time_name,
                                            time=time_value)
        assert isinstance(tparams, dict)
        assert time_name in tparams
        assert tparams[time_name] == time_value

        # check that having a dict with already the time input gives an error
        with pytest.raises(KeyError) as exc:
            system.add_time_to_params(tparams, time_name=time_name,
                                      time=time_value)
        error_string = str(exc.value).lstrip('\"').rstrip('\"')  # remove etra quotes
        assert error_string == "'params' must not contain my_time"

    # check that the (mutable) input dict is really copied and no reference created
    params_copy = params_example.copy()
    tparams = system.add_time_to_params(params_copy, time_name=time_name,
                                        time=time_value)
    tparams = {}
    assert params_copy == params_example