From c3db74ea53d1b926bb4a91e1965fb81bbe487d43 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Tue, 24 May 2016 21:16:54 +0200 Subject: [PATCH 01/17] first draft of observables --- kwant/physics/__init__.py | 2 +- kwant/physics/observables.py | 518 +++++++++++++++++++++++++++++++++++ kwant/tests/test_physics.py | 197 +++++++++++++ 3 files changed, 716 insertions(+), 1 deletion(-) create mode 100644 kwant/physics/observables.py create mode 100644 kwant/tests/test_physics.py diff --git a/kwant/physics/__init__.py b/kwant/physics/__init__.py index 10edc95..3e9c32a 100644 --- a/kwant/physics/__init__.py +++ b/kwant/physics/__init__.py @@ -10,7 +10,7 @@ # Merge the public interface of all submodules. __all__ = [] -for module in ['leads', 'dispersion', 'noise']: +for module in ['leads', 'dispersion', 'noise', 'observables']: exec('from . import {0}'.format(module)) exec('from .{0} import *'.format(module)) exec('__all__.extend({0}.__all__)'.format(module)) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py new file mode 100644 index 0000000..25b3c00 --- /dev/null +++ b/kwant/physics/observables.py @@ -0,0 +1,518 @@ +# Copyright 2011-2016 Kwant authors. +# +# 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 +# http://kwant-project.org/license. A list of Kwant authors can be found in +# the file AUTHORS.rst at the top-level directory of this distribution and at +# http://kwant-project.org/authors. +"""Tools for calculating physical observables from wavefunctions. + +This module contains tools for computing observables such as charge, +current etc. inside a Kwant system. + +Local Observables +================= +In what follows we shall first define what we call "local" observables, and +then go on to look at some examples. We start by defining what we shall +call a "generalized charge" on a single site :math:`i`: + +.. math:: q_i = \psi_i^\dagger M_i \psi_i + +where :math:`\psi_i` is a wavefunction projected onto a site :math:`i` (i.e. it +is a vector of complex numbers when there is more than 1 orbital per site, or a +single complex number when there is just 1 orbital per site), and :math:`M_i` +is a square (usually Hermitian) matrix. In the simplest case :math:`M_i` is +just the identity matrix, and defines the normal charge density. If one were to +model a superconductor with electron and hole degrees of freedom attached to +each site, however, then :math:`M_i = \sigma_z` is the correct matrix to obtain +the electric charge. Analogously, if one has a spin-ful system then one could +calculate the density of spin along any quantization axis through a judicious +choice of :math:`M_i`. + +From this "generalized charge" we can define further quantities via +a continuity equation: + +.. math:: \partial_t q_i = \sum_j J_{ij} + K_i + +where :math:`J_{ij}` and :math:`K_i`, the "generalized current` and +`generalized source` are expressed as: + +.. math:: + + J_{ij} = i(\psi_j^\dagger H_{ij}^\dagger M_i \psi_i - + \psi_i^\dagger M_i H_{ij} \psi_j) + +and + +.. math:: + + K_i = \psi_i^dagger (H_{ii}^\dagger M_i - M_i H_{ii}) \psi_i. + +The "source" term goes by different names depending on the specific physical +context. In a spin-ful system the source term would correspond to a spin +torque, for example. +""" + +__all__ = ['LocalObservable', 'Charge', 'Current', 'Source', + 'charge', 'current', 'source'] + +import abc +import functools as ft +import collections +import bisect +import numpy as np + +from ..system import InfiniteSystem + + +class LocalObservable(metaclass=abc.ABCMeta): + """Abstract base class for local observables. + + Attributes + ---------- + syst : `~kwant.system.System` + The system for which this observable is defined + onsite : callable + A function that can be called with a single site (integer) and extra + arguments, and returns the representation of the observable on that + site. This should return either a scalar or a square matrix of the + same shape as that returned by the system Hamiltonian evaluated on the + same site. The extra arguments must be the same as the extra arguments + to ``syst.hamiltonian``. + where : sequence + The places where the observable should be evaluated. "Places" can be + either sites or hoppings (pairs of sites) depending on the context + check_hermiticity : bool + If True, checks that ``onsite``, as well as any relevant parts + of the Hamiltonian are hermitian. + """ + + def __call__(self, bra=None, ket=None, per_element=True, args=()): + """Return the expectation value of the observable. + + Parameters + ---------- + bra, ket : sequence of complex + Usually a `~numpy.ndarray`, must have the same length + as the number of orbitals in the system. If only one + is provided, both ``bra`` and ``ket`` are taken as equal. + per_element : bool + If True the observable is evaluated per-element (site or hopping), + otherwise sum over the contributions from all the elements. + args : tuple, optional + The arguments to pass to the system. Used to evaluate + the ``onsite`` elements and, possibly, the system Hamiltonian. + + Returns + ------- + scalar or `~numpy.ndarray` + If ``per_element`` is False, returns a scalar, else returns a sequence + with the same length as the number of sites specified by ``where``. + """ + if bra is None: + if ket is None: + raise ValueError('Either `bra` or `ket` must be provided.') + else: + ket = np.asarray(ket) + bra = ket + elif ket is None: + bra = np.asarray(bra) + ket = bra + else: + bra = np.asarray(bra) + ket = np.asarray(ket) + # TODO: vectorize this when Kwant has support for it + result = (self._act(bra, ket, s, args, inner_product=True) + for s in self.where) + if per_element: + # we use an array of `complex` to handle the case where the + # user wants non-Hermitian observables + return np.fromiter(result, dtype=np.complex, count=len(self.where)) + else: + return sum(result) + + def act(self, ket, args=()): + """Act with the observable on a wavefunction. + + Parameters + ---------- + ket : `~numpy.ndarray` + Wavefunctions defined over all the orbitals of the system. + args : tuple + The extra arguments to the Hamiltonian value functions and + the observable ``onsite`` function. + + Returns + ------- + `~numpy.ndarray` + The result of acting on the wavefunction with the observable + """ + if not self.syst.site_ranges: + raise RuntimeError('Number of orbitals not defined.\n' + 'Declare the number of orbitals using the ' + '`norbs` keyword argument when constructing ' + 'the site families (lattices).') + nsites, _, norbs = self.syst.site_ranges[-1] + out = np.zeros((norbs,), dtype=np.complex) + for s in self.where: + self._act(out, ket, s, args, inner_product=False) + return out + + @abc.abstractmethod + def _act(self, bra, ket, element, args, inner_product): + """Act with the a single element of the observable on a wavefunction. + + Parameters + ---------- + bra, ket : `~numpy.ndarray` + Wavefunctions defined over all the orbitals of the system. + element : int or pair of ints + Either a single site or a pair of sites, to be defined by the + subclass. + args : tuple + The extra arguments to the Hamiltonian value functions and + the observable ``onsite`` function. + inner_product : bool + If true, the matrix element between ``bra`` and ``ket`` should + be returned, otherwise, accumulate the action of the observable + on ``ket`` into ``bra``. + + Returns + ------- + `~numpy.ndarray` or None + Returns None if ``inner_product`` is True, otherwise returns + an array: result of acting on the wavefunction with the observable. + """ + pass + + + +################ Concrete LocalObservable subclasses + +_inf = float('inf') + +def _get_slice(syst, site): + """Return the first orbital of this site and the next""" + assert site >= 0 and site < syst.graph.num_nodes + if not syst.site_ranges: + raise RuntimeError('Number of orbitals not defined.\n' + 'Declare the number of orbitals using the ' + '`norbs` keyword argument when constructing ' + 'the site families (lattices).') + # Calculate the index of the run that contains the site. + # The `inf` is needed to avoid the fenceposting problem + run_idx = bisect.bisect(syst.site_ranges, (site, _inf)) - 1 + # calculate the slice + first_site, norbs, orb_offset = syst.site_ranges[run_idx] + orb = orb_offset + (site - first_site) * norbs + return orb, orb + norbs + + +# XXX: This has been lifted straight from `kwant.builder`; +# maybe we should put it in a common place? +def herm_conj(value): + """Calculate the hermitian conjugate of a python object. + + If the object is neither a complex number nor a matrix, the original value + is returned. + """ + if hasattr(value, 'conjugate'): + value = value.conjugate() + if hasattr(value, 'transpose'): + value = value.transpose() + return value + + +# XXX: This has been lifted from `kwant.solvers.common`; +# maybe we should put it in a common place? +def _not_herm_conj(a, b, atol=1e-300, rtol=1e-13): + tol = rtol * np.max(np.abs(a)) + atol + return np.any(np.abs(a - herm_conj(b)) > tol) + + +def _raise_non_hermitian(onsite_not_herm, ham_not_herm): + msg = [] + if onsite_not_herm: + msg.append('Observable onsite is not Hermitian. ') + if ham_not_herm: + msg.append('Hamiltonian is not Hermitian. ') + if onsite_not_herm or ham_not_herm: + msg.append('Use option `check_hermiticity=True` ' + 'if this is intentional.') + raise ValueError(''.join(msg)) + + +class Charge(LocalObservable): + """Generalized charge observable. + + Notes + ----- + Formally, this class calculates the following expression: + + .. math:: q_{a} = \psi_a^\dagger M_a \psi_a + + where :math:`\psi_a` is the wavefunction on site ``a``, + and :math:`M_a` is a square matrix. + + The above expression is sufficiently general that it can calculate a wide + class of quantities (e.g. spin projections) with the appropriate choice of + the matrix :math:`M_a`. + """ + + # TODO: move this to Cython + def _act(self, bra, ket, site, args, inner_product): + # get the orbital slice + a_s, a_e = _get_slice(self.syst, site) + M_a = self.onsite(site, *args) + + if self.check_hermiticity: + onsite_not_herm = _not_herm_conj(M_a, M_a) + _raise_non_hermitian(onsite_not_herm, False) + + r = np.dot(M_a, ket[a_s:a_e]) + if inner_product: + return np.dot(herm_conj(bra[a_s:a_e]), r) + else: + bra[a_s:a_e] += r + + +class Current(LocalObservable): + """Generalized current observable. + + Notes + ----- + Formally, this class calculates the following expression: + + .. math:: + + J_{ab} = i (\psi_b^\dagger H_{ab}^\dagger M_a \psi_a - + \psi_a^\dagger M_a H_{ab} \psi_b) + + where :math:`H_{ab}` is the Hamiltonian matrix between sites `a` and `b`, + :math:`\psi_a` is the wavefunction on site `a`, and :math:`M_a` is a + square matrix. + + The above expression is sufficiently general that it can calculate a wide + class of quantities (e.g. spin current) with the appropriate choice of + the matrix :math:`M_a`. + """ + + # TODO: move this to Cython + def _act(self, bra, ket, sites, args, inner_product): + a, b = sites + a_s, a_e = _get_slice(self.syst, a) + b_s, b_e = _get_slice(self.syst, b) + # extract the necessary matrix elements + M_a = self.onsite(a, *args) + H_ab = self.syst.hamiltonian(a, b, *args) + H_ab_dag = herm_conj(H_ab) + + if self.check_hermiticity: + onsite_not_herm = _not_herm_conj(M_a, M_a) + H_ba = self.syst.hamiltonian(b, a, *args) + ham_not_herm = _not_herm_conj(H_ab, H_ba) + _raise_non_hermitian(onsite_not_herm, ham_not_herm) + + left = ft.reduce(np.dot, (H_ab_dag, M_a, ket[a_s:a_e])) + right = ft.reduce(np.dot, (M_a, H_ab, ket[b_s:b_e])) + if inner_product: + left = np.dot(herm_conj(bra[b_s:b_e]), left) + right = np.dot(herm_conj(bra[a_s:a_e]), right) + return 1j * (left - right) + else: + bra[b_s:b_e] += 1j * left + bra[a_s:a_e] += (-1j) * right + + +class Source(LocalObservable): + """Generalized source observable. + + Notes + ----- + Formally, this class calculates the following expression: + + .. math:: + + K_a = i \psi_a^\dagger (H_{aa}^\dagger M_a \psi_a - M_a H_{aa}) \psi_a + + where :math:`H_{ab}` is the Hamiltonian matrix between sites `a` and `b`, + :math:`\psi_a` is the wavefunction on site `a`, and :math:`M_a` is a + square matrix. + + The above expression is sufficiently general that it can calculate a wide + class of quantities (e.g. spin torque) with the appropriate choice of + the matrix :math:`M_a`. + """ + + def _act(self, bra, ket, site, args, inner_product): + # get the orbital slice + a_s, a_e = _get_slice(self.syst, site) + # extract the necessary matrix elements + M_a = self.onsite(site, *args) + H_aa = self.syst.hamiltonian(site, site, *args) + H_aa_dag = herm_conj(H_aa) + K = np.dot(H_aa_dag, M_a) - np.dot(M_a, H_aa) + + if self.check_hermiticity: + onsite_not_herm = _not_herm_conj(M_a, M_a) + ham_not_herm = _not_herm_conj(H_aa, H_aa) + _raise_non_hermitian(onsite_not_herm, ham_not_herm) + + r = 1j * np.dot(K, ket[a_s:a_e]) + if inner_product: + return np.dot(herm_conj(bra[a_s:a_e]), r) + else: + bra[a_s:a_e] += r + + + +################ Constructors for use with finalized Builders + +def _normalize_where(syst, where, hoppings): + if hoppings: + if where is None: + _where = list(syst.graph) + elif callable(where): + _where = list(filter(lambda h: where(*h), syst.graph)) + else: + _where = [(syst.id_by_site[a], syst.id_by_site[b]) + for a, b in where] + + if isinstance(syst, InfiniteSystem): + if any(a > syst.cell_size or b > syst.cell_size for a, b in _where): + raise ValueError('Only intra-cell hoppings may be specified ' + 'using `where`.') + else: + if where is None: + _where = list(range(len(syst.sites))) + elif callable(where): + _where = list(filter(where, syst.sites)) + else: + _where = [syst.id_by_site[s] for s in where] + + if isinstance(syst, InfiniteSystem): + if any(w > syst.cell_size for w in _where): + raise ValueError('Only sites in the fundamental domain may be ' + 'specified using `where`.') + + return _where + + +def _normalize_onsite(syst, onsite): + _sites = syst.sites + if callable(onsite): + def _onsite(site_id, *args): + return onsite(_sites[site_id], *args) + elif isinstance(onsite, collections.Mapping): + def _onsite(site_id, *args): + return onsite[_sites[site_id].family] + else: + def _onsite(site_id, *args): + return onsite + + return _onsite + + +def charge(syst, onsite=1, where=None, check_hermiticity=True): + """Make an observable for calculating generalized charges. + + Parameters + ---------- + syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` + onsite : scalar or square matrix or dict or callable + The onsite matrix that defines the charge. If a dict is given, it + maps from site families to either scalars or square matrices. If a + function is given it must take the same arguments as the onsite + Hamiltonian functions of the system. + where : sequence of `~kwant.builder.Site`s or callable, optional + Where to evaluate the observable. If a function is provided, it + should take a single `~kwant.builder.Site` and return True or + False. If not provided, the observable will be calculated + everywhere in the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it + is not Hermitian, then an error will be raised when the + observable is evaluated. + + Returns + ------- + `~kwant.physics.observables.Charge` + An object that can be called with a wavefunction and extra + arguments (which are passed to the ``onsite`` function), + and returns the generalized charge. + """ + q = Charge() + q.syst = syst + q.onsite = _normalize_onsite(syst, onsite) + q.where = _normalize_where(syst, where, False) + q.check_hermiticity = check_hermiticity + return q + + +def current(syst, onsite=1, where=None, check_hermiticity=True): + """Make an observable for calculating generalized currents. + + Parameters + ---------- + syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` + onsite : scalar or square matrix or dict or callable, optional + The onsite matrix that defines the charge from which this current + is derived. If a dict is given, it maps from site families to + either scalars or square matrices. If a function is given it must + take the same arguments as the onsite Hamiltonian functions of the + system. + where : sequence of pairs of `~kwant.builder.Site`s, optional + The hoppings on which to evaluate the observable. If not provided, + the observable will be calculated on every hopping in the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it + is not Hermitian, then an error will be raised when the + observable is evaluated. + + Returns + ------- + `~kwant.physics.observables.Current` + An object that can be called with a wavefunction and extra + arguments (which are passed to the ``onsite`` function), + and returns the generalized current. + """ + j = Current() + j.syst = syst + j.onsite = _normalize_onsite(syst, onsite) + j.where = _normalize_where(syst, where, True) + j.check_hermiticity = check_hermiticity + return j + + +def source(syst, onsite=1, where=None, check_hermiticity=True): + """Make an observable for calculating generalized charge sources. + + Parameters + ---------- + syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` + onsite : scalar or square matrix or dict or callable + The onsite matrix that defines the charge. If a dict is given, it + maps from site families to either scalars or square matrices. If a + function is given it must take the same arguments as the onsite + Hamiltonian functions of the system. + where : sequence of `~kwant.builder.Site`s, optional + Where to evaluate the observable. If not provided, the + observable will be calculated everywhere in the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it + is not Hermitian, then an error will be raised when the + observable is evaluated. + + Returns + ------- + `~kwant.physics.observables.Charge` + An object that can be called with a wavefunction and extra + arguments (which are passed to the ``onsite`` function), + and returns the generalized charge source. + """ + s = Source() + s.syst = syst + s.onsite = _normalize_onsite(syst, onsite) + s.where = _normalize_where(syst, where, False) + s.check_hermiticity = check_hermiticity + return s diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py new file mode 100644 index 0000000..0f0e520 --- /dev/null +++ b/kwant/tests/test_physics.py @@ -0,0 +1,197 @@ +# Copyright 2011-2016 Kwant authors. +# +# 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 +# http://kwant-project.org/license. A list of Kwant authors can be found in +# the file AUTHORS.rst at the top-level directory of this distribution and at +# http://kwant-project.org/authors. + +import functools as ft +import numpy as np +import numpy.linalg as la +from nose.tools import assert_raises + +import kwant +from kwant.physics import observables as obs + + +def test_observables(): + lat = kwant.lattice.square(norbs=1) + + def random_onsite(i): + return 2 + kwant.digest.uniform(i.tag) + + def random_hopping(i, j): + return -1 + kwant.digest.uniform(i.tag + j.tag) + + syst = kwant.Builder() + syst[(lat(i, j) for i in range(3) for j in range(3))] = random_onsite + syst[lat.neighbors()] = random_hopping + fsyst = syst.finalized() + N = len(fsyst.sites) + + ## Construction tests + + # test construction failure if norbs not given + latnone = kwant.lattice.chain() + syst[latnone(0)] = 1 + for a in (obs.charge, obs.current, obs.source): + assert_raises(RuntimeError, a(syst.finalized()), np.zeros((N,))) + del syst[latnone(0)] + + # test basic construction + Q = obs.charge(fsyst) + assert Q.where == list(range(N)) + assert all(Q.onsite(i) == 1 for i in range(N)) + J = obs.current(fsyst) + assert J.where == list(fsyst.graph) + assert all(J.onsite(i) == 1 for i, j in fsyst.graph) + # test construction with dict `onsite` + for A in (J, Q): + A = obs.charge(fsyst, {lat: 1}) + assert all(A.onsite(i) == 1 for i in range(N)) + # test construction with a functional onsite + for A in (J, Q): + A = obs.charge(fsyst, lambda site: site.pos[0]) + assert all(A.onsite(i) == fsyst.sites[i].pos[0] for i in range(N)) + # test construction with `where` given by a sequence + where = [lat(2, 2), lat(1, 1)] + fwhere = [fsyst.id_by_site[s] for s in where] + A = obs.charge(fsyst, where=where) + assert A.where == fwhere + where = [(lat(2, 2), lat(1, 2)), (lat(0, 0), lat(0, 1))] + fwhere = [(fsyst.id_by_site[a], fsyst.id_by_site[b]) + for a, b in where] + A = obs.current(fsyst, where=where) + assert A.where == fwhere + # test construction with `where` given by a function + where_list = set([(1, 0), (1, 1), (1, 2)]) + def where(site): + return site.tag in where_list + A = obs.charge(fsyst, where=where) + assert all(w.tag in where_list for w in A.where) + where_list = set(kwant.HoppingKind((1, 0), lat)(syst)) + def where(a, b): + return (a, b) in where_list + A = obs.current(fsyst, where=where) + assert all(w in where_list for w in A.where) + + ## Functional tests + + def test_inner(val, A, bra, ket=None, args=()): + if not ket: + ket = bra + act_val = np.dot(bra.conj(), A.act(ket, args=args)) + inner_val = A(bra, ket, per_element=False, args=args) + assert np.isclose(act_val, val) + assert np.isclose(inner_val, val) + + def test_per_element(val, A, bra, ket=None, args=()): + if not ket: + ket = bra + assert np.allclose(A(bra, ket, args=args), val) + + + # Test closed system + ev, wfs = la.eigh(fsyst.hamiltonian_submatrix()) + + Q = obs.charge(fsyst) + J = obs.current(fsyst) + K = obs.source(fsyst) + + for wf in wfs.T: # wfs[:, i] is i'th eigenvector + assert np.isclose(np.sum(Q(wf)), 1) # eigenvectors are normalized + test_per_element(0, J, wf) + test_per_element(0, K, wf) + test_inner(1, Q, wf) + test_inner(0, J, wf) + test_inner(0, K, wf) + + # Test infinite system + lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) + lead[(lat(0, j) for j in range(3))] = 2 + lead[lat.neighbors()] = -1 + flead = lead.finalized() + J_intra = obs.current(flead, where=[(lat(0, 0), lat(0, 1)), + (lat(0, 1), lat(0, 2))]) + + # Cannot calculate current between unit cells because the wavefunction + # is only defined over a single unit cell + assert_raises(ValueError, obs.current, flead, where=[(lat(0, j), lat(1, j)) + for j in range(3)]) + prop, _ = flead.modes(energy=1.) + for wf, v in zip(prop.wave_functions.T, prop.velocities): + # unit cell is uniform, no transverse current + test_per_element(0, J_intra, wf) + + # Test open system + # Homogeneous system except for a barrier. We check that the current + # on the right of the barrier due to a mode `m` is equal to + # Σ_n |t_nm|^2. Similarly the current on the left of the barrier + # is checked against 1 - Σ_n |r_nm|^2 + + syst = kwant.Builder() + # induce some backscattering that mixes modes + syst[(lat(i, j) for i in range(3) for j in range(3))] = random_onsite + # add extra sites so we can calculate the true current + syst[(lat(i, j) for i in [-1, 3] for j in range(3))] = 2 + syst[lat.neighbors()] = -1 + lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) + lead[(lat(0, j) for j in range(3))] = 2 + lead[lat.neighbors()] = -1 + fsyst = syst.finalized() + syst.attach_lead(lead) + syst.attach_lead(lead.reversed()) + fsyst = syst.finalized() + J_right = obs.current(fsyst, where=[(lat(3, j), lat(2, j)) + for j in range(3)]) + J_left = obs.current(fsyst, where=[(lat(0, j), lat(-1, j)) + for j in range(3)]) + smatrix = kwant.smatrix(fsyst, energy=1.0) + t = smatrix.submatrix(1, 0).T # want to iterate over the columns + r = smatrix.submatrix(0, 0).T # want to iterate over the columns + wfs = kwant.wave_function(fsyst, energy=1.0)(0) + for rv, tv, wf in zip(r, t, wfs): + test_inner(np.sum(np.abs(tv)**2), J_right, wf) + test_inner(1 - np.sum(np.abs(rv)**2), J_left, wf) + + # Test with non-hermitian matrices + for A in (obs.charge, obs.current, obs.source): + assert_raises(ValueError, A(fsyst, 1j), wf) + a = A(fsyst, 1j, check_hermiticity=False) + b = a(wf, per_element=False) + test_inner(b, a, wf) + + # Test with multiple orbitals -- calculate spin current + s_x = np.array([[0, 1], [1, 0]]) + s_z = np.array([[1, 0], [0, -1]]) + + def onsite(site, B): + return 2 * (np.eye(2) + B * s_z) + + lat = kwant.lattice.chain(norbs=2) + syst = kwant.Builder() + syst[(lat(i) for i in range(2))] = onsite + syst[lat.neighbors()] = -1 * np.eye(2) + lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) + lead[lat(0)] = onsite + lead[lat.neighbors()] = -1 * np.eye(2) + syst.attach_lead(lead) + syst.attach_lead(lead.reversed()) + fsyst = syst.finalized() + args = (0.1,) + down, up = kwant.wave_function(fsyst, energy=1., args=args)(0) + + where = [(lat(1), lat(0))] + spin_current_z = obs.current(fsyst, s_z, where=where) + test_inner(1, spin_current_z, up, args=args) + test_inner(-1, spin_current_z, down, args=args) + + # calculate spin_x torque + spin_torque_x = obs.source(fsyst, s_x, where=[lat(0)]) + i = fsyst.id_by_site[lat(0)] + psi = up[2*i:2*(i+1)] + down[2*i:2*(i+1)] + H_ii = onsite(None, *args) + K = np.dot(H_ii, s_x) - np.dot(s_x, H_ii) + expect = 1j * ft.reduce(np.dot, (psi.conj(), K, psi)) + test_inner(expect, spin_torque_x, up+down, args=args) -- GitLab From 1f8bcf3b919c3cbfd4bf3596d93ebd17ca59434c Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Sun, 5 Jun 2016 13:41:50 +0200 Subject: [PATCH 02/17] clarify docs --- kwant/physics/observables.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index 25b3c00..5bb01a6 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -80,8 +80,8 @@ class LocalObservable(metaclass=abc.ABCMeta): same site. The extra arguments must be the same as the extra arguments to ``syst.hamiltonian``. where : sequence - The places where the observable should be evaluated. "Places" can be - either sites or hoppings (pairs of sites) depending on the context + The places where the observable should be evaluated. "Places" can be + either sites or hoppings (pairs of sites) depending on the context. check_hermiticity : bool If True, checks that ``onsite``, as well as any relevant parts of the Hamiltonian are hermitian. @@ -199,7 +199,7 @@ def _get_slice(syst, site): 'Declare the number of orbitals using the ' '`norbs` keyword argument when constructing ' 'the site families (lattices).') - # Calculate the index of the run that contains the site. + # Calculate the index of the range that contains the site. # The `inf` is needed to avoid the fenceposting problem run_idx = bisect.bisect(syst.site_ranges, (site, _inf)) - 1 # calculate the slice -- GitLab From 9ebde5a40258bd093e450705cca4e275a684a627 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Sun, 5 Jun 2016 13:42:37 +0200 Subject: [PATCH 03/17] fix up tests --- kwant/tests/test_physics.py | 38 +++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index 0f0e520..c607776 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -47,13 +47,13 @@ def test_observables(): assert J.where == list(fsyst.graph) assert all(J.onsite(i) == 1 for i, j in fsyst.graph) # test construction with dict `onsite` - for A in (J, Q): - A = obs.charge(fsyst, {lat: 1}) - assert all(A.onsite(i) == 1 for i in range(N)) + for A in (obs.charge, obs.current, obs.source): + B = A(fsyst, {lat: 1}) + assert all(B.onsite(i) == 1 for i in range(N)) # test construction with a functional onsite - for A in (J, Q): - A = obs.charge(fsyst, lambda site: site.pos[0]) - assert all(A.onsite(i) == fsyst.sites[i].pos[0] for i in range(N)) + for A in (obs.charge, obs.current, obs.source): + B = A(fsyst, lambda site: site.pos[0]) # x-position operator + assert all(B.onsite(i) == fsyst.sites[i].pos[0] for i in range(N)) # test construction with `where` given by a sequence where = [lat(2, 2), lat(1, 1)] fwhere = [fsyst.id_by_site[s] for s in where] @@ -101,7 +101,7 @@ def test_observables(): for wf in wfs.T: # wfs[:, i] is i'th eigenvector assert np.isclose(np.sum(Q(wf)), 1) # eigenvectors are normalized - test_per_element(0, J, wf) + test_per_element(0, J, wf) # time-reversal symmetry: no current test_per_element(0, K, wf) test_inner(1, Q, wf) test_inner(0, J, wf) @@ -121,7 +121,7 @@ def test_observables(): for j in range(3)]) prop, _ = flead.modes(energy=1.) for wf, v in zip(prop.wave_functions.T, prop.velocities): - # unit cell is uniform, no transverse current + # no transverse current test_per_element(0, J_intra, wf) # Test open system @@ -155,23 +155,25 @@ def test_observables(): test_inner(np.sum(np.abs(tv)**2), J_right, wf) test_inner(1 - np.sum(np.abs(rv)**2), J_left, wf) - # Test with non-hermitian matrices + # Test `check_hermiticity=False` for A in (obs.charge, obs.current, obs.source): assert_raises(ValueError, A(fsyst, 1j), wf) a = A(fsyst, 1j, check_hermiticity=False) b = a(wf, per_element=False) test_inner(b, a, wf) - # Test with multiple orbitals -- calculate spin current + ### Test with multiple orbitals -- calculate spin current + L = 20 s_x = np.array([[0, 1], [1, 0]]) + s_y = np.array([[0, -1j], [1j, 0]]) s_z = np.array([[1, 0], [0, -1]]) def onsite(site, B): - return 2 * (np.eye(2) + B * s_z) + return 2 * np.eye(2) + B * s_z lat = kwant.lattice.chain(norbs=2) syst = kwant.Builder() - syst[(lat(i) for i in range(2))] = onsite + syst[(lat(i) for i in range(L))] = onsite syst[lat.neighbors()] = -1 * np.eye(2) lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) lead[lat(0)] = onsite @@ -182,14 +184,14 @@ def test_observables(): args = (0.1,) down, up = kwant.wave_function(fsyst, energy=1., args=args)(0) - where = [(lat(1), lat(0))] - spin_current_z = obs.current(fsyst, s_z, where=where) - test_inner(1, spin_current_z, up, args=args) - test_inner(-1, spin_current_z, down, args=args) + x_hoppings = kwant.builder.HoppingKind((1,), lat) + spin_current_z = obs.current(fsyst, s_z, where=x_hoppings(syst)) + test_per_element(1, spin_current_z, up, args=args) + test_per_element(-1, spin_current_z, down, args=args) # calculate spin_x torque - spin_torque_x = obs.source(fsyst, s_x, where=[lat(0)]) - i = fsyst.id_by_site[lat(0)] + spin_torque_x = obs.source(fsyst, s_x, where=[lat(L//2)]) + i = fsyst.id_by_site[lat(L//2)] psi = up[2*i:2*(i+1)] + down[2*i:2*(i+1)] H_ii = onsite(None, *args) K = np.dot(H_ii, s_x) - np.dot(s_x, H_ii) -- GitLab From 4a5a4a22e82082f67a3ba49db12eae568a448bff Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 6 Jun 2016 12:12:23 +0200 Subject: [PATCH 04/17] break tests into separate functions --- kwant/tests/test_physics.py | 120 +++++++++++++++++++++--------------- 1 file changed, 72 insertions(+), 48 deletions(-) diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index c607776..15d1dad 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -15,8 +15,7 @@ import kwant from kwant.physics import observables as obs -def test_observables(): - lat = kwant.lattice.square(norbs=1) +def _random_square_system(N, norbs=1): def random_onsite(i): return 2 + kwant.digest.uniform(i.tag) @@ -24,14 +23,26 @@ def test_observables(): def random_hopping(i, j): return -1 + kwant.digest.uniform(i.tag + j.tag) + lat = kwant.lattice.square(norbs=norbs) syst = kwant.Builder() - syst[(lat(i, j) for i in range(3) for j in range(3))] = random_onsite + syst[(lat(i, j) for i in range(N) for j in range(N))] = random_onsite syst[lat.neighbors()] = random_hopping + return lat, syst + + +def _perfect_lead(N, norbs=1): + lat = kwant.lattice.square(norbs=norbs) + syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) + syst[(lat(0, j) for j in range(N))] = 2 + syst[lat.neighbors()] = -1 + return lat, syst + + +def test_observables_construction(): + lat, syst = _random_square_system(3) fsyst = syst.finalized() N = len(fsyst.sites) - ## Construction tests - # test construction failure if norbs not given latnone = kwant.lattice.chain() syst[latnone(0)] = 1 @@ -46,53 +57,61 @@ def test_observables(): J = obs.current(fsyst) assert J.where == list(fsyst.graph) assert all(J.onsite(i) == 1 for i, j in fsyst.graph) + # test construction with dict `onsite` for A in (obs.charge, obs.current, obs.source): B = A(fsyst, {lat: 1}) assert all(B.onsite(i) == 1 for i in range(N)) + # test construction with a functional onsite for A in (obs.charge, obs.current, obs.source): B = A(fsyst, lambda site: site.pos[0]) # x-position operator assert all(B.onsite(i) == fsyst.sites[i].pos[0] for i in range(N)) + # test construction with `where` given by a sequence where = [lat(2, 2), lat(1, 1)] fwhere = [fsyst.id_by_site[s] for s in where] A = obs.charge(fsyst, where=where) assert A.where == fwhere + where = [(lat(2, 2), lat(1, 2)), (lat(0, 0), lat(0, 1))] fwhere = [(fsyst.id_by_site[a], fsyst.id_by_site[b]) for a, b in where] A = obs.current(fsyst, where=where) assert A.where == fwhere + # test construction with `where` given by a function where_list = set([(1, 0), (1, 1), (1, 2)]) def where(site): return site.tag in where_list A = obs.charge(fsyst, where=where) assert all(w.tag in where_list for w in A.where) + where_list = set(kwant.HoppingKind((1, 0), lat)(syst)) def where(a, b): return (a, b) in where_list A = obs.current(fsyst, where=where) assert all(w in where_list for w in A.where) - ## Functional tests - def test_inner(val, A, bra, ket=None, args=()): - if not ket: - ket = bra - act_val = np.dot(bra.conj(), A.act(ket, args=args)) - inner_val = A(bra, ket, per_element=False, args=args) - assert np.isclose(act_val, val) - assert np.isclose(inner_val, val) +def _test_inner(val, A, bra, ket=None, args=()): + if not ket: + ket = bra + act_val = np.dot(bra.conj(), A.act(ket, args=args)) + inner_val = A(bra, ket, per_element=False, args=args) + assert np.isclose(act_val, val) + assert np.isclose(inner_val, val) + - def test_per_element(val, A, bra, ket=None, args=()): - if not ket: - ket = bra - assert np.allclose(A(bra, ket, args=args), val) +def _test_per_element(val, A, bra, ket=None, args=()): + if not ket: + ket = bra + assert np.allclose(A(bra, ket, args=args), val) - # Test closed system +def test_observables_finite(): + lat, syst = _random_square_system(3) + fsyst = syst.finalized() ev, wfs = la.eigh(fsyst.hamiltonian_submatrix()) Q = obs.charge(fsyst) @@ -101,66 +120,71 @@ def test_observables(): for wf in wfs.T: # wfs[:, i] is i'th eigenvector assert np.isclose(np.sum(Q(wf)), 1) # eigenvectors are normalized - test_per_element(0, J, wf) # time-reversal symmetry: no current - test_per_element(0, K, wf) - test_inner(1, Q, wf) - test_inner(0, J, wf) - test_inner(0, K, wf) - - # Test infinite system - lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) - lead[(lat(0, j) for j in range(3))] = 2 - lead[lat.neighbors()] = -1 + _test_per_element(0, J, wf) # time-reversal symmetry: no current + _test_per_element(0, K, wf) + _test_inner(1, Q, wf) + _test_inner(0, J, wf) + _test_inner(0, K, wf) + + +def test_observables_infinite(): + lat, lead = _perfect_lead(3) flead = lead.finalized() - J_intra = obs.current(flead, where=[(lat(0, 0), lat(0, 1)), - (lat(0, 1), lat(0, 2))]) # Cannot calculate current between unit cells because the wavefunction # is only defined over a single unit cell assert_raises(ValueError, obs.current, flead, where=[(lat(0, j), lat(1, j)) for j in range(3)]) + + transverse_kind = kwant.builder.HoppingKind((0, 1), lat) + J_intra = obs.current(flead, where=list(transverse_kind(lead))) + prop, _ = flead.modes(energy=1.) for wf, v in zip(prop.wave_functions.T, prop.velocities): # no transverse current - test_per_element(0, J_intra, wf) + _test_per_element(0, J_intra, wf) + +def test_observables_scattering(): # Test open system # Homogeneous system except for a barrier. We check that the current # on the right of the barrier due to a mode `m` is equal to # Σ_n |t_nm|^2. Similarly the current on the left of the barrier # is checked against 1 - Σ_n |r_nm|^2 - syst = kwant.Builder() - # induce some backscattering that mixes modes - syst[(lat(i, j) for i in range(3) for j in range(3))] = random_onsite - # add extra sites so we can calculate the true current - syst[(lat(i, j) for i in [-1, 3] for j in range(3))] = 2 - syst[lat.neighbors()] = -1 - lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) - lead[(lat(0, j) for j in range(3))] = 2 - lead[lat.neighbors()] = -1 - fsyst = syst.finalized() + N = 3 + lat, syst = _random_square_system(N) + # add extra sites so we can calculate the current in a region + # where there is no backscattering + syst[(lat(i, j) for i in [-1, N] for j in range(N))] = 2 + syst[((lat(-1, j), lat(0, j)) for j in range(N))] = -2 + syst[((lat(N-1, j), lat(N, j)) for j in range(N))] = -2 + + lat, lead = _perfect_lead(3) syst.attach_lead(lead) syst.attach_lead(lead.reversed()) fsyst = syst.finalized() - J_right = obs.current(fsyst, where=[(lat(3, j), lat(2, j)) + + # currents on the left and right of the disordered region + J_right = obs.current(fsyst, where=[(lat(N, j), lat(N-1, j)) for j in range(3)]) J_left = obs.current(fsyst, where=[(lat(0, j), lat(-1, j)) for j in range(3)]) + smatrix = kwant.smatrix(fsyst, energy=1.0) t = smatrix.submatrix(1, 0).T # want to iterate over the columns r = smatrix.submatrix(0, 0).T # want to iterate over the columns wfs = kwant.wave_function(fsyst, energy=1.0)(0) for rv, tv, wf in zip(r, t, wfs): - test_inner(np.sum(np.abs(tv)**2), J_right, wf) - test_inner(1 - np.sum(np.abs(rv)**2), J_left, wf) + _test_inner(np.sum(np.abs(tv)**2), J_right, wf) + _test_inner(1 - np.sum(np.abs(rv)**2), J_left, wf) # Test `check_hermiticity=False` for A in (obs.charge, obs.current, obs.source): assert_raises(ValueError, A(fsyst, 1j), wf) a = A(fsyst, 1j, check_hermiticity=False) b = a(wf, per_element=False) - test_inner(b, a, wf) + _test_inner(b, a, wf) ### Test with multiple orbitals -- calculate spin current L = 20 @@ -186,8 +210,8 @@ def test_observables(): x_hoppings = kwant.builder.HoppingKind((1,), lat) spin_current_z = obs.current(fsyst, s_z, where=x_hoppings(syst)) - test_per_element(1, spin_current_z, up, args=args) - test_per_element(-1, spin_current_z, down, args=args) + _test_per_element(1, spin_current_z, up, args=args) + _test_per_element(-1, spin_current_z, down, args=args) # calculate spin_x torque spin_torque_x = obs.source(fsyst, s_x, where=[lat(L//2)]) @@ -196,4 +220,4 @@ def test_observables(): H_ii = onsite(None, *args) K = np.dot(H_ii, s_x) - np.dot(s_x, H_ii) expect = 1j * ft.reduce(np.dot, (psi.conj(), K, psi)) - test_inner(expect, spin_torque_x, up+down, args=args) + _test_inner(expect, spin_torque_x, up+down, args=args) -- GitLab From c08327c3f32cc7adcdd53aa2b4569f8ba0d4c99b Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 6 Jun 2016 12:20:04 +0200 Subject: [PATCH 05/17] documentation fixups as per Anton's request --- kwant/physics/observables.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index 5bb01a6..856193e 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -25,9 +25,9 @@ is a square (usually Hermitian) matrix. In the simplest case :math:`M_i` is just the identity matrix, and defines the normal charge density. If one were to model a superconductor with electron and hole degrees of freedom attached to each site, however, then :math:`M_i = \sigma_z` is the correct matrix to obtain -the electric charge. Analogously, if one has a spin-ful system then one could -calculate the density of spin along any quantization axis through a judicious -choice of :math:`M_i`. +the electric charge. Analogously, if one has a spinful system then one could +calculate the density of spin along any quantization axis by defining +:math:`M_i = \sigma \cdot n` where :math:`n` is a unit vector. From this "generalized charge" we can define further quantities via a continuity equation: @@ -71,7 +71,8 @@ class LocalObservable(metaclass=abc.ABCMeta): Attributes ---------- syst : `~kwant.system.System` - The system for which this observable is defined + The system for which this observable is defined. Must have the + number of orbitals defined via the ``site_ranges`` attribute. onsite : callable A function that can be called with a single site (integer) and extra arguments, and returns the representation of the observable on that @@ -505,7 +506,7 @@ def source(syst, onsite=1, where=None, check_hermiticity=True): Returns ------- - `~kwant.physics.observables.Charge` + `~kwant.physics.observables.Source` An object that can be called with a wavefunction and extra arguments (which are passed to the ``onsite`` function), and returns the generalized charge source. -- GitLab From 5dd874a69672091651b15595b02455a655705328 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 6 Jun 2016 14:54:02 +0200 Subject: [PATCH 06/17] be more explicit about what quantities are actually being calculated --- kwant/physics/observables.py | 167 ++++++++++++++++++++++++----------- 1 file changed, 114 insertions(+), 53 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index 856193e..cd03df1 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -5,7 +5,7 @@ # http://kwant-project.org/license. A list of Kwant authors can be found in # the file AUTHORS.rst at the top-level directory of this distribution and at # http://kwant-project.org/authors. -"""Tools for calculating physical observables from wavefunctions. +r"""Tools for calculating physical observables from wavefunctions. This module contains tools for computing observables such as charge, current etc. inside a Kwant system. @@ -14,43 +14,50 @@ Local Observables ================= In what follows we shall first define what we call "local" observables, and then go on to look at some examples. We start by defining what we shall -call a "generalized charge" on a single site :math:`i`: - -.. math:: q_i = \psi_i^\dagger M_i \psi_i - -where :math:`\psi_i` is a wavefunction projected onto a site :math:`i` (i.e. it -is a vector of complex numbers when there is more than 1 orbital per site, or a -single complex number when there is just 1 orbital per site), and :math:`M_i` -is a square (usually Hermitian) matrix. In the simplest case :math:`M_i` is -just the identity matrix, and defines the normal charge density. If one were to -model a superconductor with electron and hole degrees of freedom attached to -each site, however, then :math:`M_i = \sigma_z` is the correct matrix to obtain -the electric charge. Analogously, if one has a spinful system then one could -calculate the density of spin along any quantization axis by defining -:math:`M_i = \sigma \cdot n` where :math:`n` is a unit vector. - -From this "generalized charge" we can define further quantities via -a continuity equation: - -.. math:: \partial_t q_i = \sum_j J_{ij} + K_i - -where :math:`J_{ij}` and :math:`K_i`, the "generalized current` and -`generalized source` are expressed as: +call a "generalized charge" on a single site :math:`a`, defined by the +observable :math:`Q_a`, which has matrix elements :math:`[Q_a]_{pq}`: + +.. math:: [Q_a]_{pq} = \delta_{pa} M_{a} \delta_{qb} + +where :math:`\delta_{pa}` is the Kronecker delta and :math:`M_a` is a square +(usually Hermitian) *matrix* that acts on the orbitals (degrees of freedom) +belonging to site :math:`a`. If the site has only 1 orbital, then :math:`M_a` +is a scalar. In the simplest case :math:`M_a` is just the identity matrix, and +defines the normal charge density. If one were to model a superconductor with +electron and hole degrees of freedom attached to each site, however, then +:math:`M_a = \sigma_z` is the correct matrix to obtain the electric charge. +Analogously, if one has a spinful system then one could calculate the density +of spin along any quantization axis by defining :math:`M_a = \sigma \cdot n` +where :math:`n` is a unit vector and :math:`sigma` is a vector of Pauli +matrices. + +From this "generalized charge" we can define further quantities: the +"generalized current" and "generalized source", represented by operators +:math:`J_{ab}` and :math:`K_a` respectively, with matrix elements: .. math:: - J_{ij} = i(\psi_j^\dagger H_{ij}^\dagger M_i \psi_i - - \psi_i^\dagger M_i H_{ij} \psi_j) + [J_{ab}]_{pq} = i (\delta_{pb} H^\dagger_{ab} M_{a} \delta_{qa} - + \delta_{pa} M_{a} H_{ab} \delta_{qb}) and .. math:: - K_i = \psi_i^dagger (H_{ii}^\dagger M_i - M_i H_{ii}) \psi_i. + [K_{a}]_{pq} = i (\delta_{pa} H^\dagger_{aa} M_{a} \delta_{qa} - + \delta_{pa} M_{a} H_{aa} \delta_{qa}) -The "source" term goes by different names depending on the specific physical -context. In a spin-ful system the source term would correspond to a spin -torque, for example. +where :math:`H_{ab}` is the Hamiltonian matrix element between sites :math:`a` +and :math:`b`. In the above we define :math:`[J_{ab}]_{pp} = 0` for all +:math:`p` (i.e. there are no "on-site currents"). The "source" term goes by +different names depending on the specific physical context. In a spinful system +the source term would correspond to a spin torque, for example. + +The matrix elements of these three operators are related by a continuity +equation: + +.. math:: \frac{d}{dt} \mel{\phi}{Q_a}{\psi} = \mel{\phi}{K_a}{\psi} + + \sum_{b \ne a} \mel{\phi}{J_{ab}}{\psi} """ __all__ = ['LocalObservable', 'Charge', 'Current', 'Source', @@ -246,18 +253,28 @@ def _raise_non_hermitian(onsite_not_herm, ham_not_herm): class Charge(LocalObservable): """Generalized charge observable. - Notes - ----- - Formally, this class calculates the following expression: + When this class is called in the following way:: + + Q = kwant.physics.observables.Charge(fsyst, ...) + Q(phi, psi) + + the following expression is calculated for each site ``a`` + in ``where``:: + + phi[i:j].conjugate().dot(M_a.dot(psi[i:j])) + + where ``i:j`` is a slice over all the orbitals on site ``a``, + and ``M_a`` is a square matrix associated with site ``a``. - .. math:: q_{a} = \psi_a^\dagger M_a \psi_a + When called in the following way:: - where :math:`\psi_a` is the wavefunction on site ``a``, - and :math:`M_a` is a square matrix. + phi = np.zeros_like(psi) + phi = Q.act(psi) - The above expression is sufficiently general that it can calculate a wide - class of quantities (e.g. spin projections) with the appropriate choice of - the matrix :math:`M_a`. + this is equivalent to performing the following operation + for each site ``a`` in ``where``:: + + phi[i:j] += M_a.dot(psi[i:j]) """ # TODO: move this to Cython @@ -280,22 +297,35 @@ class Charge(LocalObservable): class Current(LocalObservable): """Generalized current observable. - Notes - ----- - Formally, this class calculates the following expression: + When this class is called in the following way:: - .. math:: + J = kwant.physics.observables.Current(fsyst, ...) + J(phi, psi) - J_{ab} = i (\psi_b^\dagger H_{ab}^\dagger M_a \psi_a - - \psi_a^\dagger M_a H_{ab} \psi_b) + the following expression is calculated for each pair of sites + ``(a, b)`` in ``where``:: - where :math:`H_{ab}` is the Hamiltonian matrix between sites `a` and `b`, - :math:`\psi_a` is the wavefunction on site `a`, and :math:`M_a` is a - square matrix. + 1j * (phi[bi:bj].conjugate().dot( + H_ab.conjugate().transpose().dot(M_a.dot(psi[ai:aj]))) + - + phi[ai:aj].conjugate().dot( + M_a.dot(H_ab.dot(psi[bi:bj]))) + ) + + where ``ai:aj`` and ``bi:bj`` are slices over all the orbitals on sites + ``a`` and ``b`` respectively. ``M_a`` is a square matrix and ``H_ab`` is + the Hamiltonian matrix element between the sites. - The above expression is sufficiently general that it can calculate a wide - class of quantities (e.g. spin current) with the appropriate choice of - the matrix :math:`M_a`. + When called in the following way:: + + phi = np.zeros_like(psi) + phi = Q.act(psi) + + this is equivalent to performing the following operations + for each pair of sites ``(a, b)`` in ``where``:: + + phi[bi:bj] += 1j * H_ab.conjugate().transpose().dot(M_a.dot(psi[ai:aj])) + phi[ai:aj] += -1j * M_a.dot(H_ab.dot(psi[bi:bj])) """ # TODO: move this to Cython @@ -328,6 +358,35 @@ class Current(LocalObservable): class Source(LocalObservable): """Generalized source observable. + When this class is called in the following way:: + + K = kwant.physics.observables.Charge(fsyst, ...) + K(phi, psi) + + the following expression is calculated for each site ``a`` + in ``where``:: + + 1j * (phi[i:j].conjugate().dot( + H_aa.conjugate().transpose().dot(M_a.dot(psi[i:j]))) + - + phi[i:j].conjugate().dot( + M_a.dot(H_aa.dot(psi[i:j]))) + ) + + where ``i:j`` is a slice over all the orbitals on site ``a`` ``M_a`` is a + square matrix and ``H_aa`` is the onsite Hamiltonian matrix element for + site ``a``. + + When called in the following way:: + + phi = np.zeros_like(psi) + phi = K.act(psi) + + this is equivalent to performing the following operation + for each element in ``where``:: + + phi[i:j] += M_a.dot(psi[i:j]) + Notes ----- Formally, this class calculates the following expression: @@ -340,9 +399,11 @@ class Source(LocalObservable): :math:`\psi_a` is the wavefunction on site `a`, and :math:`M_a` is a square matrix. - The above expression is sufficiently general that it can calculate a wide - class of quantities (e.g. spin torque) with the appropriate choice of - the matrix :math:`M_a`. + then this is equivalent to performing the following operation + for each site ``a`` in ``where``:: + + phi[i:j] += 1j * (H_ab.conjugate().transpose().dot(M_a) - + M_a.dot(H_aa)).dot(psi[i:j])) """ def _act(self, bra, ket, site, args, inner_product): -- GitLab From 2b82ebab922d7cea6c25446175bdaa21763f8a3a Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 6 Jun 2016 20:09:55 +0200 Subject: [PATCH 07/17] add nontrivial test involving a random on-site gauge transformation --- kwant/tests/test_physics.py | 82 +++++++++++++++++++++++++++++++++---- 1 file changed, 73 insertions(+), 9 deletions(-) diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index 15d1dad..8b341c0 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -7,14 +7,22 @@ # http://kwant-project.org/authors. import functools as ft +from collections import deque import numpy as np import numpy.linalg as la from nose.tools import assert_raises +# needed to get round odd bug in test_mask_interpolate +from contextlib import contextmanager import kwant from kwant.physics import observables as obs +sigmax = np.array([[0, 1], [1, 0]]) +sigmay = np.array([[0, -1j], [1j, 0]]) +sigmaz = np.array([[1, 0], [0, -1]]) + + def _random_square_system(N, norbs=1): def random_onsite(i): @@ -186,15 +194,12 @@ def test_observables_scattering(): b = a(wf, per_element=False) _test_inner(b, a, wf) - ### Test with multiple orbitals -- calculate spin current - L = 20 - s_x = np.array([[0, 1], [1, 0]]) - s_y = np.array([[0, -1j], [1j, 0]]) - s_z = np.array([[1, 0], [0, -1]]) +def test_observables_spin(): def onsite(site, B): - return 2 * np.eye(2) + B * s_z + return 2 * np.eye(2) + B * sigmaz + L = 20 lat = kwant.lattice.chain(norbs=2) syst = kwant.Builder() syst[(lat(i) for i in range(L))] = onsite @@ -209,15 +214,74 @@ def test_observables_scattering(): down, up = kwant.wave_function(fsyst, energy=1., args=args)(0) x_hoppings = kwant.builder.HoppingKind((1,), lat) - spin_current_z = obs.current(fsyst, s_z, where=x_hoppings(syst)) + spin_current_z = obs.current(fsyst, sigmaz, where=x_hoppings(syst)) _test_per_element(1, spin_current_z, up, args=args) _test_per_element(-1, spin_current_z, down, args=args) # calculate spin_x torque - spin_torque_x = obs.source(fsyst, s_x, where=[lat(L//2)]) + spin_torque_x = obs.source(fsyst, sigmax, where=[lat(L//2)]) i = fsyst.id_by_site[lat(L//2)] psi = up[2*i:2*(i+1)] + down[2*i:2*(i+1)] H_ii = onsite(None, *args) - K = np.dot(H_ii, s_x) - np.dot(s_x, H_ii) + K = np.dot(H_ii, sigmax) - np.dot(sigmax, H_ii) expect = 1j * ft.reduce(np.dot, (psi.conj(), K, psi)) _test_inner(expect, spin_torque_x, up+down, args=args) + + +def test_observables_gauged(): + # Test that we get the same answer when we apply a random + # gauge (unitary) transformation to each site. We also + # adjust our definition of the current to match + + # this is to get round a bug in test_mask_interpolate (?!) that fails when + # the random number state is altered. + @contextmanager + def save_random_state(): + old_state = np.random.get_state() + yield + np.random.set_state(old_state) + + L = 20 + with save_random_state(): + Us = deque([kwant.rmt.circular(2) for i in range(L)]) + # need these to get the coupling to the leads right + Us.append(np.eye(2)) + Us.appendleft(np.eye(2)) + + H0 = 2 * np.eye(2) + 0.1 * sigmaz # onsite + V0 = -1 * np.eye(2) # hopping + + lat = kwant.lattice.chain(norbs=2) + syst = kwant.Builder() + + for i, U in enumerate(Us): + syst[lat(i)] = ft.reduce(np.dot, (U, H0, U.conjugate().transpose())) + + for a, b in kwant.builder.HoppingKind((1,), lat)(syst): + i, j = a.tag[0], b.tag[0] + syst[(a, b)] = ft.reduce(np.dot, + (Us[i], V0, Us[j].conjugate().transpose())) + + lead = kwant.Builder(kwant.TranslationalSymmetry((-1,))) + lead[lat(0)] = H0 + lead[lat.neighbors()] = V0 + syst.attach_lead(lead) + syst.attach_lead(lead.reversed()) + fsyst = syst.finalized() + down, up = kwant.wave_function(fsyst, energy=1.0)(0) + + def M_a(site): + i = site.tag[0] + return ft.reduce(np.dot, + (Us[i], sigmaz, Us[i].conjugate().transpose())) + + x_hoppings = kwant.builder.HoppingKind((1,), lat) + spin_current_gauge = obs.current(fsyst, M_a, where=x_hoppings(syst)) + print(spin_current_gauge(up)) + _test_per_element(1, spin_current_gauge, up) + _test_per_element(-1, spin_current_gauge, down) + # check the reverse is also true + minus_x_hoppings = kwant.builder.HoppingKind((-1,), lat) + spin_current_gauge = obs.current(fsyst, M_a, where=minus_x_hoppings(syst)) + _test_per_element(-1, spin_current_gauge, up) + _test_per_element(1, spin_current_gauge, down) -- GitLab From 3b66a4c3d9084537dc8dfcef986f3adebe0bc6f8 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Mon, 6 Jun 2016 21:03:52 +0200 Subject: [PATCH 08/17] add `bind` method and make `where` a tuple `bind` is useful in the case where the observable will be used repeatedly with the same arguments. --- kwant/physics/observables.py | 99 +++++++++++++++++++++++++++++++----- kwant/tests/test_physics.py | 23 +++++---- 2 files changed, 98 insertions(+), 24 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index cd03df1..6c6ce10 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -87,7 +87,7 @@ class LocalObservable(metaclass=abc.ABCMeta): same shape as that returned by the system Hamiltonian evaluated on the same site. The extra arguments must be the same as the extra arguments to ``syst.hamiltonian``. - where : sequence + where : tuple The places where the observable should be evaluated. "Places" can be either sites or hoppings (pairs of sites) depending on the context. check_hermiticity : bool @@ -95,6 +95,11 @@ class LocalObservable(metaclass=abc.ABCMeta): of the Hamiltonian are hermitian. """ + def __init__(self): + self._bound = False + self._bound_onsite = None + self._bound_hamiltonian = None + def __call__(self, bra=None, ket=None, per_element=True, args=()): """Return the expectation value of the observable. @@ -166,6 +171,16 @@ class LocalObservable(metaclass=abc.ABCMeta): self._act(out, ket, s, args, inner_product=False) return out + @abc.abstractmethod + def bind(self, args): + """Bind the given arguments to this observable. + + Returns a copy if this observable that does not need + to be passed extra arguments when calling or using the + ``act`` method. + """ + pass + @abc.abstractmethod def _act(self, bra, ket, element, args, inner_product): """Act with the a single element of the observable on a wavefunction. @@ -277,11 +292,28 @@ class Charge(LocalObservable): phi[i:j] += M_a.dot(psi[i:j]) """ + def bind(self, args): + q = self.__class__() + q._bound = True + # if/when the loop over `where` is brought into `_act`, + # this can be made into a sequence -- the lookups are the + # same speed, but lists take less space + q._bound_onsite = {site: self.onsite(site, *args) + for site in self.where} + q.syst = self.syst + q.onsite = self.onsite + q.where = self.where + q.check_hermiticity = self.check_hermiticity + return q + # TODO: move this to Cython def _act(self, bra, ket, site, args, inner_product): # get the orbital slice a_s, a_e = _get_slice(self.syst, site) - M_a = self.onsite(site, *args) + if self._bound: + M_a = self._bound_onsite[site] + else: + M_a = self.onsite(site, *args) if self.check_hermiticity: onsite_not_herm = _not_herm_conj(M_a, M_a) @@ -328,14 +360,34 @@ class Current(LocalObservable): phi[ai:aj] += -1j * M_a.dot(H_ab.dot(psi[bi:bj])) """ + def bind(self, args): + q = self.__class__() + q._bound = True + q._bound_onsite = {a: self.onsite(a, *args) + for a, b in self.where} + # if/when the loop over `where` is brought into `_act`, + # this can be made into a sequence -- the lookups are the + # same speed, but lists take less space + q._bound_hamiltonian = {(a, b): self.syst.hamiltonian(a, b, *args) + for a, b in self.where} + q.syst = self.syst + q.onsite = self.onsite + q.where = self.where + q.check_hermiticity = self.check_hermiticity + return q + # TODO: move this to Cython def _act(self, bra, ket, sites, args, inner_product): a, b = sites a_s, a_e = _get_slice(self.syst, a) b_s, b_e = _get_slice(self.syst, b) # extract the necessary matrix elements - M_a = self.onsite(a, *args) - H_ab = self.syst.hamiltonian(a, b, *args) + if self._bound: + M_a = self._bound_onsite[a] + H_ab = self._bound_hamiltonian[a, b] + else: + M_a = self.onsite(a, *args) + H_ab = self.syst.hamiltonian(a, b, *args) H_ab_dag = herm_conj(H_ab) if self.check_hermiticity: @@ -406,12 +458,33 @@ class Source(LocalObservable): M_a.dot(H_aa)).dot(psi[i:j])) """ + def bind(self, args): + q = self.__class__() + q._bound = True + # if/when the loop over `where` is brought into `_act`, + # this can be made into a sequence -- the lookups are the + # same speed, but lists take less space + q._bound_onsite = {site: self.onsite(site, *args) + for site in self.where} + q._bound_hamiltonian = {site: self.syst.hamiltonian(site, site, *args) + for site in self.where} + q.syst = self.syst + q.onsite = self.onsite + q.where = self.where + q.check_hermiticity = self.check_hermiticity + return q + + # TODO: move this to Cython def _act(self, bra, ket, site, args, inner_product): # get the orbital slice a_s, a_e = _get_slice(self.syst, site) # extract the necessary matrix elements - M_a = self.onsite(site, *args) - H_aa = self.syst.hamiltonian(site, site, *args) + if self._bound: + M_a = self._bound_onsite[site] + H_aa = self._bound_hamiltonian[site] + else: + M_a = self.onsite(site, *args) + H_aa = self.syst.hamiltonian(site, site, *args) H_aa_dag = herm_conj(H_aa) K = np.dot(H_aa_dag, M_a) - np.dot(M_a, H_aa) @@ -433,12 +506,12 @@ class Source(LocalObservable): def _normalize_where(syst, where, hoppings): if hoppings: if where is None: - _where = list(syst.graph) + _where = tuple(syst.graph) elif callable(where): - _where = list(filter(lambda h: where(*h), syst.graph)) + _where = tuple(filter(lambda h: where(*h), syst.graph)) else: - _where = [(syst.id_by_site[a], syst.id_by_site[b]) - for a, b in where] + _where = tuple((syst.id_by_site[a], syst.id_by_site[b]) + for a, b in where) if isinstance(syst, InfiniteSystem): if any(a > syst.cell_size or b > syst.cell_size for a, b in _where): @@ -446,11 +519,11 @@ def _normalize_where(syst, where, hoppings): 'using `where`.') else: if where is None: - _where = list(range(len(syst.sites))) + _where = tuple(range(len(syst.sites))) elif callable(where): - _where = list(filter(where, syst.sites)) + _where = tuple(filter(where, syst.sites)) else: - _where = [syst.id_by_site[s] for s in where] + _where = tuple(syst.id_by_site[s] for s in where) if isinstance(syst, InfiniteSystem): if any(w > syst.cell_size for w in _where): diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index 8b341c0..c6f94ac 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -60,10 +60,10 @@ def test_observables_construction(): # test basic construction Q = obs.charge(fsyst) - assert Q.where == list(range(N)) + assert Q.where == tuple(range(N)) assert all(Q.onsite(i) == 1 for i in range(N)) J = obs.current(fsyst) - assert J.where == list(fsyst.graph) + assert J.where == tuple(fsyst.graph) assert all(J.onsite(i) == 1 for i, j in fsyst.graph) # test construction with dict `onsite` @@ -78,13 +78,13 @@ def test_observables_construction(): # test construction with `where` given by a sequence where = [lat(2, 2), lat(1, 1)] - fwhere = [fsyst.id_by_site[s] for s in where] + fwhere = tuple(fsyst.id_by_site[s] for s in where) A = obs.charge(fsyst, where=where) assert A.where == fwhere where = [(lat(2, 2), lat(1, 2)), (lat(0, 0), lat(0, 1))] - fwhere = [(fsyst.id_by_site[a], fsyst.id_by_site[b]) - for a, b in where] + fwhere = tuple((fsyst.id_by_site[a], fsyst.id_by_site[b]) + for a, b in where) A = obs.current(fsyst, where=where) assert A.where == fwhere @@ -105,16 +105,18 @@ def test_observables_construction(): def _test_inner(val, A, bra, ket=None, args=()): if not ket: ket = bra - act_val = np.dot(bra.conj(), A.act(ket, args=args)) - inner_val = A(bra, ket, per_element=False, args=args) - assert np.isclose(act_val, val) - assert np.isclose(inner_val, val) + for A in (A, A.bind(args)): + act_val = np.dot(bra.conj(), A.act(ket, args=args)) + inner_val = A(bra, ket, per_element=False, args=args) + assert np.isclose(act_val, val) + assert np.isclose(inner_val, val) def _test_per_element(val, A, bra, ket=None, args=()): if not ket: ket = bra - assert np.allclose(A(bra, ket, args=args), val) + for A in (A, A.bind(args)): + assert np.allclose(A(bra, ket, args=args), val) def test_observables_finite(): @@ -277,7 +279,6 @@ def test_observables_gauged(): x_hoppings = kwant.builder.HoppingKind((1,), lat) spin_current_gauge = obs.current(fsyst, M_a, where=x_hoppings(syst)) - print(spin_current_gauge(up)) _test_per_element(1, spin_current_gauge, up) _test_per_element(-1, spin_current_gauge, down) # check the reverse is also true -- GitLab From 47279cf91fbe56e9850f0ec5701649b6fc1ab79e Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Tue, 7 Jun 2016 10:34:00 +0200 Subject: [PATCH 09/17] fix test and test description --- kwant/tests/test_physics.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index c6f94ac..f761d4f 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -156,19 +156,18 @@ def test_observables_infinite(): def test_observables_scattering(): - # Test open system - # Homogeneous system except for a barrier. We check that the current - # on the right of the barrier due to a mode `m` is equal to - # Σ_n |t_nm|^2. Similarly the current on the left of the barrier + # Disordered system with two ordered strips on the left/right. We check + # that the current on the right of the disorder due to incoming mode `m` is + # equal to Σ_n |t_nm|^2. Similarly the current on the left of the disorder # is checked against 1 - Σ_n |r_nm|^2 - N = 3 + N = 10 lat, syst = _random_square_system(N) # add extra sites so we can calculate the current in a region # where there is no backscattering syst[(lat(i, j) for i in [-1, N] for j in range(N))] = 2 - syst[((lat(-1, j), lat(0, j)) for j in range(N))] = -2 - syst[((lat(N-1, j), lat(N, j)) for j in range(N))] = -2 + syst[((lat(-1, j), lat(0, j)) for j in range(N))] = -1 + syst[((lat(N-1, j), lat(N, j)) for j in range(N))] = -1 lat, lead = _perfect_lead(3) syst.attach_lead(lead) -- GitLab From 672b9d2cfb64337ef31febd91a01b95f06121a2a Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Tue, 7 Jun 2016 10:39:33 +0200 Subject: [PATCH 10/17] check to avoid possibly common error when using `bind` --- kwant/physics/observables.py | 10 +++++++++- kwant/tests/test_physics.py | 21 ++++++++++++++------- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index 6c6ce10..e7b71dc 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -122,6 +122,10 @@ class LocalObservable(metaclass=abc.ABCMeta): If ``per_element`` is False, returns a scalar, else returns a sequence with the same length as the number of sites specified by ``where``. """ + if self._bound and args: + raise ValueError('Extra arguments are already bound to this ' + 'observable. You should call this observable ' + 'without providing `args`.') if bra is None: if ket is None: raise ValueError('Either `bra` or `ket` must be provided.') @@ -165,6 +169,10 @@ class LocalObservable(metaclass=abc.ABCMeta): 'Declare the number of orbitals using the ' '`norbs` keyword argument when constructing ' 'the site families (lattices).') + if self._bound and args: + raise ValueError('Extra arguments are already bound to this ' + 'observable. You should call this observable ' + 'without providing `args`.') nsites, _, norbs = self.syst.site_ranges[-1] out = np.zeros((norbs,), dtype=np.complex) for s in self.where: @@ -175,7 +183,7 @@ class LocalObservable(metaclass=abc.ABCMeta): def bind(self, args): """Bind the given arguments to this observable. - Returns a copy if this observable that does not need + Returns a copy of this observable that does not need to be passed extra arguments when calling or using the ``act`` method. """ diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index f761d4f..150f5b8 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -105,18 +105,25 @@ def test_observables_construction(): def _test_inner(val, A, bra, ket=None, args=()): if not ket: ket = bra - for A in (A, A.bind(args)): - act_val = np.dot(bra.conj(), A.act(ket, args=args)) - inner_val = A(bra, ket, per_element=False, args=args) - assert np.isclose(act_val, val) - assert np.isclose(inner_val, val) + act_val = np.dot(bra.conj(), A.act(ket, args=args)) + inner_val = A(bra, ket, per_element=False, args=args) + assert np.isclose(act_val, val) + assert np.isclose(inner_val, val) + # with bound args + A = A.bind(args) + act_val = np.dot(bra.conj(), A.act(ket)) + inner_val = A(bra, ket, per_element=False) + assert np.isclose(act_val, val) + assert np.isclose(inner_val, val) def _test_per_element(val, A, bra, ket=None, args=()): if not ket: ket = bra - for A in (A, A.bind(args)): - assert np.allclose(A(bra, ket, args=args), val) + assert np.allclose(A(bra, ket, args=args), val) + # with bound args + A = A.bind(args) + assert np.allclose(A(bra, ket), val) def test_observables_finite(): -- GitLab From d8188089a2e6ac91d63d30b6f7470127ae4985a9 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 17:55:16 +0200 Subject: [PATCH 11/17] fix copyright line --- kwant/tests/test_physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index 150f5b8..75d03b1 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -1,4 +1,4 @@ -# Copyright 2011-2016 Kwant authors. +# Copyright 2016 Kwant authors. # # 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 -- GitLab From 656224b89c3d0b6c56ca16b912528ca78b6e0d3c Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 17:56:46 +0200 Subject: [PATCH 12/17] do not explicitly mention low-level system API in docstring --- kwant/physics/observables.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index e7b71dc..e8284d1 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -79,7 +79,7 @@ class LocalObservable(metaclass=abc.ABCMeta): ---------- syst : `~kwant.system.System` The system for which this observable is defined. Must have the - number of orbitals defined via the ``site_ranges`` attribute. + number of orbitals defined for all site families. onsite : callable A function that can be called with a single site (integer) and extra arguments, and returns the representation of the observable on that -- GitLab From 378e2bbac31e75a31ccc95276f2b7cfca9eb0304 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 18:20:29 +0200 Subject: [PATCH 13/17] rename `per_element` to `reduce` --- kwant/physics/observables.py | 14 +++++++------- kwant/tests/test_physics.py | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index e8284d1..ebda5af 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -100,7 +100,7 @@ class LocalObservable(metaclass=abc.ABCMeta): self._bound_onsite = None self._bound_hamiltonian = None - def __call__(self, bra=None, ket=None, per_element=True, args=()): + def __call__(self, bra, ket=None, reduce=False, args=()): """Return the expectation value of the observable. Parameters @@ -109,8 +109,8 @@ class LocalObservable(metaclass=abc.ABCMeta): Usually a `~numpy.ndarray`, must have the same length as the number of orbitals in the system. If only one is provided, both ``bra`` and ``ket`` are taken as equal. - per_element : bool - If True the observable is evaluated per-element (site or hopping), + reduce : bool + If False the observable is evaluated per-element (site or hopping), otherwise sum over the contributions from all the elements. args : tuple, optional The arguments to pass to the system. Used to evaluate @@ -119,7 +119,7 @@ class LocalObservable(metaclass=abc.ABCMeta): Returns ------- scalar or `~numpy.ndarray` - If ``per_element`` is False, returns a scalar, else returns a sequence + If ``reduce`` is True, returns a scalar, else returns a sequence with the same length as the number of sites specified by ``where``. """ if self._bound and args: @@ -141,12 +141,12 @@ class LocalObservable(metaclass=abc.ABCMeta): # TODO: vectorize this when Kwant has support for it result = (self._act(bra, ket, s, args, inner_product=True) for s in self.where) - if per_element: + if reduce: + return sum(result) + else: # we use an array of `complex` to handle the case where the # user wants non-Hermitian observables return np.fromiter(result, dtype=np.complex, count=len(self.where)) - else: - return sum(result) def act(self, ket, args=()): """Act with the observable on a wavefunction. diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index 75d03b1..bb4fd6a 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -106,13 +106,13 @@ def _test_inner(val, A, bra, ket=None, args=()): if not ket: ket = bra act_val = np.dot(bra.conj(), A.act(ket, args=args)) - inner_val = A(bra, ket, per_element=False, args=args) + inner_val = A(bra, ket, reduce=True, args=args) assert np.isclose(act_val, val) assert np.isclose(inner_val, val) # with bound args A = A.bind(args) act_val = np.dot(bra.conj(), A.act(ket)) - inner_val = A(bra, ket, per_element=False) + inner_val = A(bra, ket, reduce=True) assert np.isclose(act_val, val) assert np.isclose(inner_val, val) @@ -199,7 +199,7 @@ def test_observables_scattering(): for A in (obs.charge, obs.current, obs.source): assert_raises(ValueError, A(fsyst, 1j), wf) a = A(fsyst, 1j, check_hermiticity=False) - b = a(wf, per_element=False) + b = a(wf, reduce=True) _test_inner(b, a, wf) def test_observables_spin(): -- GitLab From ce9f3ce734a8cbc04083e79ffcb4f91d8f9ecf9b Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 18:22:54 +0200 Subject: [PATCH 14/17] change `_not_herm_conj` to `_is_not_herm_conj` --- kwant/physics/observables.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index ebda5af..f7d31e3 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -256,7 +256,7 @@ def herm_conj(value): # XXX: This has been lifted from `kwant.solvers.common`; # maybe we should put it in a common place? -def _not_herm_conj(a, b, atol=1e-300, rtol=1e-13): +def _is_not_herm_conj(a, b, atol=1e-300, rtol=1e-13): tol = rtol * np.max(np.abs(a)) + atol return np.any(np.abs(a - herm_conj(b)) > tol) @@ -324,7 +324,7 @@ class Charge(LocalObservable): M_a = self.onsite(site, *args) if self.check_hermiticity: - onsite_not_herm = _not_herm_conj(M_a, M_a) + onsite_not_herm = _is_not_herm_conj(M_a, M_a) _raise_non_hermitian(onsite_not_herm, False) r = np.dot(M_a, ket[a_s:a_e]) @@ -399,9 +399,9 @@ class Current(LocalObservable): H_ab_dag = herm_conj(H_ab) if self.check_hermiticity: - onsite_not_herm = _not_herm_conj(M_a, M_a) + onsite_not_herm = _is_not_herm_conj(M_a, M_a) H_ba = self.syst.hamiltonian(b, a, *args) - ham_not_herm = _not_herm_conj(H_ab, H_ba) + ham_not_herm = _is_not_herm_conj(H_ab, H_ba) _raise_non_hermitian(onsite_not_herm, ham_not_herm) left = ft.reduce(np.dot, (H_ab_dag, M_a, ket[a_s:a_e])) @@ -497,8 +497,8 @@ class Source(LocalObservable): K = np.dot(H_aa_dag, M_a) - np.dot(M_a, H_aa) if self.check_hermiticity: - onsite_not_herm = _not_herm_conj(M_a, M_a) - ham_not_herm = _not_herm_conj(H_aa, H_aa) + onsite_not_herm = _is_not_herm_conj(M_a, M_a) + ham_not_herm = _is_not_herm_conj(H_aa, H_aa) _raise_non_hermitian(onsite_not_herm, ham_not_herm) r = 1j * np.dot(K, ket[a_s:a_e]) -- GitLab From be97b177dc606b9c5bf97a4a5b854fd9fddf9153 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 18:24:12 +0200 Subject: [PATCH 15/17] change site_ranges check to `site_ranges is None` --- kwant/physics/observables.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index f7d31e3..7793a53 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -164,7 +164,7 @@ class LocalObservable(metaclass=abc.ABCMeta): `~numpy.ndarray` The result of acting on the wavefunction with the observable """ - if not self.syst.site_ranges: + if self.syst.site_ranges is None: raise RuntimeError('Number of orbitals not defined.\n' 'Declare the number of orbitals using the ' '`norbs` keyword argument when constructing ' @@ -225,7 +225,7 @@ _inf = float('inf') def _get_slice(syst, site): """Return the first orbital of this site and the next""" assert site >= 0 and site < syst.graph.num_nodes - if not syst.site_ranges: + if syst.site_ranges is None: raise RuntimeError('Number of orbitals not defined.\n' 'Declare the number of orbitals using the ' '`norbs` keyword argument when constructing ' -- GitLab From 045a6979f72ea818111ce84a0bbe1763889097ab Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Wed, 22 Jun 2016 21:54:36 +0200 Subject: [PATCH 16/17] remove `herm_conj` and normalize everything using `ta.matrix` --- kwant/physics/observables.py | 51 ++++++++++++++---------------------- 1 file changed, 19 insertions(+), 32 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index 7793a53..bcd1f98 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -68,6 +68,7 @@ import functools as ft import collections import bisect import numpy as np +import tinyarray as ta from ..system import InfiniteSystem @@ -239,26 +240,12 @@ def _get_slice(syst, site): return orb, orb + norbs -# XXX: This has been lifted straight from `kwant.builder`; -# maybe we should put it in a common place? -def herm_conj(value): - """Calculate the hermitian conjugate of a python object. - - If the object is neither a complex number nor a matrix, the original value - is returned. - """ - if hasattr(value, 'conjugate'): - value = value.conjugate() - if hasattr(value, 'transpose'): - value = value.transpose() - return value - - # XXX: This has been lifted from `kwant.solvers.common`; # maybe we should put it in a common place? def _is_not_herm_conj(a, b, atol=1e-300, rtol=1e-13): tol = rtol * np.max(np.abs(a)) + atol - return np.any(np.abs(a - herm_conj(b)) > tol) + return np.any(np.abs(ta.matrix(a) - + ta.matrix(b).conjugate().transpose()) > tol) def _raise_non_hermitian(onsite_not_herm, ham_not_herm): @@ -306,7 +293,7 @@ class Charge(LocalObservable): # if/when the loop over `where` is brought into `_act`, # this can be made into a sequence -- the lookups are the # same speed, but lists take less space - q._bound_onsite = {site: self.onsite(site, *args) + q._bound_onsite = {site: ta.matrix(self.onsite(site, *args)) for site in self.where} q.syst = self.syst q.onsite = self.onsite @@ -321,7 +308,7 @@ class Charge(LocalObservable): if self._bound: M_a = self._bound_onsite[site] else: - M_a = self.onsite(site, *args) + M_a = ta.matrix(self.onsite(site, *args)) if self.check_hermiticity: onsite_not_herm = _is_not_herm_conj(M_a, M_a) @@ -329,7 +316,7 @@ class Charge(LocalObservable): r = np.dot(M_a, ket[a_s:a_e]) if inner_product: - return np.dot(herm_conj(bra[a_s:a_e]), r) + return np.dot(bra[a_s:a_e].conjugate(), r) else: bra[a_s:a_e] += r @@ -371,12 +358,12 @@ class Current(LocalObservable): def bind(self, args): q = self.__class__() q._bound = True - q._bound_onsite = {a: self.onsite(a, *args) + q._bound_onsite = {a: ta.matrix(self.onsite(a, *args)) for a, b in self.where} # if/when the loop over `where` is brought into `_act`, # this can be made into a sequence -- the lookups are the # same speed, but lists take less space - q._bound_hamiltonian = {(a, b): self.syst.hamiltonian(a, b, *args) + q._bound_hamiltonian = {(a, b): ta.matrix(self.syst.hamiltonian(a, b, *args)) for a, b in self.where} q.syst = self.syst q.onsite = self.onsite @@ -394,9 +381,9 @@ class Current(LocalObservable): M_a = self._bound_onsite[a] H_ab = self._bound_hamiltonian[a, b] else: - M_a = self.onsite(a, *args) - H_ab = self.syst.hamiltonian(a, b, *args) - H_ab_dag = herm_conj(H_ab) + M_a = ta.matrix(self.onsite(a, *args)) + H_ab = ta.matrix(self.syst.hamiltonian(a, b, *args)) + H_ab_dag = H_ab.conjugate().transpose() if self.check_hermiticity: onsite_not_herm = _is_not_herm_conj(M_a, M_a) @@ -407,8 +394,8 @@ class Current(LocalObservable): left = ft.reduce(np.dot, (H_ab_dag, M_a, ket[a_s:a_e])) right = ft.reduce(np.dot, (M_a, H_ab, ket[b_s:b_e])) if inner_product: - left = np.dot(herm_conj(bra[b_s:b_e]), left) - right = np.dot(herm_conj(bra[a_s:a_e]), right) + left = np.dot(bra[b_s:b_e].conjugate(), left) + right = np.dot(bra[a_s:a_e].conjugate(), right) return 1j * (left - right) else: bra[b_s:b_e] += 1j * left @@ -472,9 +459,9 @@ class Source(LocalObservable): # if/when the loop over `where` is brought into `_act`, # this can be made into a sequence -- the lookups are the # same speed, but lists take less space - q._bound_onsite = {site: self.onsite(site, *args) + q._bound_onsite = {site: ta.matrix(self.onsite(site, *args)) for site in self.where} - q._bound_hamiltonian = {site: self.syst.hamiltonian(site, site, *args) + q._bound_hamiltonian = {site: ta.matrix(self.syst.hamiltonian(site, site, *args)) for site in self.where} q.syst = self.syst q.onsite = self.onsite @@ -491,9 +478,9 @@ class Source(LocalObservable): M_a = self._bound_onsite[site] H_aa = self._bound_hamiltonian[site] else: - M_a = self.onsite(site, *args) - H_aa = self.syst.hamiltonian(site, site, *args) - H_aa_dag = herm_conj(H_aa) + M_a = ta.matrix(self.onsite(site, *args)) + H_aa = ta.matrix(self.syst.hamiltonian(site, site, *args)) + H_aa_dag = H_aa.conjugate().transpose() K = np.dot(H_aa_dag, M_a) - np.dot(M_a, H_aa) if self.check_hermiticity: @@ -503,7 +490,7 @@ class Source(LocalObservable): r = 1j * np.dot(K, ket[a_s:a_e]) if inner_product: - return np.dot(herm_conj(bra[a_s:a_e]), r) + return np.dot(bra[a_s:a_e].conjugate(), r) else: bra[a_s:a_e] += r -- GitLab From 79c59923f5e17f2d2736880fe613d7fa47d0c9b2 Mon Sep 17 00:00:00 2001 From: Joseph Weston Date: Thu, 23 Jun 2016 00:23:27 +0200 Subject: [PATCH 17/17] use constructors instead of factory functions --- kwant/physics/observables.py | 363 ++++++++++++++++++----------------- kwant/tests/test_physics.py | 55 +++--- 2 files changed, 214 insertions(+), 204 deletions(-) diff --git a/kwant/physics/observables.py b/kwant/physics/observables.py index bcd1f98..192725e 100644 --- a/kwant/physics/observables.py +++ b/kwant/physics/observables.py @@ -60,8 +60,7 @@ equation: \sum_{b \ne a} \mel{\phi}{J_{ab}}{\psi} """ -__all__ = ['LocalObservable', 'Charge', 'Current', 'Source', - 'charge', 'current', 'source'] +__all__ = ['LocalObservable', 'Charge', 'Current', 'Source'] import abc import functools as ft @@ -71,6 +70,7 @@ import numpy as np import tinyarray as ta from ..system import InfiniteSystem +from .. import builder class LocalObservable(metaclass=abc.ABCMeta): @@ -88,9 +88,8 @@ class LocalObservable(metaclass=abc.ABCMeta): same shape as that returned by the system Hamiltonian evaluated on the same site. The extra arguments must be the same as the extra arguments to ``syst.hamiltonian``. - where : tuple - The places where the observable should be evaluated. "Places" can be - either sites or hoppings (pairs of sites) depending on the context. + where : sequence of pairs of `int` + Pairs of sites where the observable should be evaluated. check_hermiticity : bool If True, checks that ``onsite``, as well as any relevant parts of the Hamiltonian are hermitian. @@ -219,7 +218,7 @@ class LocalObservable(metaclass=abc.ABCMeta): -################ Concrete LocalObservable subclasses +################ Utility functions _inf = float('inf') @@ -248,6 +247,11 @@ def _is_not_herm_conj(a, b, atol=1e-300, rtol=1e-13): ta.matrix(b).conjugate().transpose()) > tol) +def _is_finalized_builder(syst): + return isinstance(syst, (builder.FiniteSystem, + builder.InfiniteSystem)) + + def _raise_non_hermitian(onsite_not_herm, ham_not_herm): msg = [] if onsite_not_herm: @@ -260,9 +264,84 @@ def _raise_non_hermitian(onsite_not_herm, ham_not_herm): raise ValueError(''.join(msg)) +def _normalize_where(syst, where, hoppings): + if hoppings: + if where is None: + _where = tuple(syst.graph) + elif callable(where): + if _is_finalized_builder(syst): + def idx_where(hop): + a, b = hop + return where(syst.sites[a], syst.sites[b]) + _where = tuple(filter(idx_where, syst.graph)) + else: + _where = tuple(filter(lambda h: where(*h), syst.graph)) + else: + if _is_finalized_builder(syst): + _where = tuple((syst.id_by_site[a], syst.id_by_site[b]) + for a, b in where) + else: + _where = tuple(where) + + if isinstance(syst, InfiniteSystem): + if any(a > syst.cell_size or b > syst.cell_size for a, b in _where): + raise ValueError('Only intra-cell hoppings may be specified ' + 'using `where`.') + else: + if where is None: + _where = list(range(syst.graph.num_nodes)) + elif callable(where): + if _is_finalized_builder(syst): + _where = [syst.id_by_site[a] for a in filter(where, syst.sites)] + else: + _where = list(filter(where, range(syst.graph.num_nodes))) + else: + if _is_finalized_builder(syst): + _where = list(syst.id_by_site[s] for s in where) + else: + _where = list(where) + + if isinstance(syst, InfiniteSystem): + if any(w > syst.cell_size for w in _where): + raise ValueError('Only sites in the fundamental domain may be ' + 'specified using `where`.') + + # LocalObservable.where is sequence of *pairs* of sites, even if + # we are only calculating a purely onsite quantity + _where = tuple(zip(_where, _where)) + + return _where + + +def _normalize_onsite(syst, onsite): + if callable(onsite): + if _is_finalized_builder(syst): + _sites = syst.sites + def _onsite(site_id, *args): + return onsite(_sites[site_id], *args) + else: + _onsite = onsite + elif isinstance(onsite, collections.Mapping): + if not _is_finalized_builder(syst): + raise TypeError('Provide `onsite` as a value or a function for ' + 'systems that are not finalized Builders.') + _sites = syst.sites + def _onsite(site_id, *args): + return onsite[_sites[site_id].family] + else: + def _onsite(site_id, *args): + return onsite + + return _onsite + + + +################ Concrete LocalObservable subclasses + class Charge(LocalObservable): """Generalized charge observable. + When this class is called in the following way:: Q = kwant.physics.observables.Charge(fsyst, ...) @@ -287,14 +366,22 @@ class Charge(LocalObservable): phi[i:j] += M_a.dot(psi[i:j]) """ + def __init__(self, syst, onsite=1, where=None, check_hermiticity=True): + super().__init__() + self.syst = syst + self.check_hermiticity = check_hermiticity + self.onsite = _normalize_onsite(syst, onsite) + self.where = _normalize_where(syst, where, hoppings=False) + def bind(self, args): - q = self.__class__() + cls = self.__class__ + q = cls.__new__(cls) q._bound = True # if/when the loop over `where` is brought into `_act`, # this can be made into a sequence -- the lookups are the # same speed, but lists take less space - q._bound_onsite = {site: ta.matrix(self.onsite(site, *args)) - for site in self.where} + q._bound_onsite = {a: ta.matrix(self.onsite(a, *args)) + for a, _ in self.where} q.syst = self.syst q.onsite = self.onsite q.where = self.where @@ -302,13 +389,15 @@ class Charge(LocalObservable): return q # TODO: move this to Cython - def _act(self, bra, ket, site, args, inner_product): + def _act(self, bra, ket, sites, args, inner_product): # get the orbital slice - a_s, a_e = _get_slice(self.syst, site) + a, b = sites + assert a == b + a_s, a_e = _get_slice(self.syst, a) if self._bound: - M_a = self._bound_onsite[site] + M_a = self._bound_onsite[a] else: - M_a = ta.matrix(self.onsite(site, *args)) + M_a = ta.matrix(self.onsite(a, *args)) if self.check_hermiticity: onsite_not_herm = _is_not_herm_conj(M_a, M_a) @@ -324,6 +413,40 @@ class Charge(LocalObservable): class Current(LocalObservable): """Generalized current observable. + Parameters + ---------- + syst : `~kwant.system.System` + onsite : scalar or square matrix or dict or callable, optional + The onsite matrix that defines the charge from which this current + is derived. If a dict is given, it maps from site families to + either scalars or square matrices. If a function is given it must + take the same arguments as the onsite Hamiltonian functions of the + system. + where : sequence of pairs of `int` or pairs of `~kwant.builder.Site`s, + or callable, optional + Where to evaluate the observable. If `syst` is not a finalized Builder, + then this should be a sequence of pairs of integers. If a function is + provided, it should take a pair of integers or a pair of + `~kwant.builder.Site` (if `syst` is a finalized builder) and return + True or False. If not provided, the observable will be calculated over + all hoppings in the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it + is not Hermitian, then an error will be raised when the + observable is evaluated. + + Parameters + ---------- + onsite : scalar or square matrix or dict or callable + The onsite matrix that defines the charge. If a dict is given, it + maps from site families to either scalars or square matrices. If a + function is given it must take the same arguments as the onsite + Hamiltonian functions of the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it is not + Hermitian, then an error will be raised when the observable is + evaluated. + When this class is called in the following way:: J = kwant.physics.observables.Current(fsyst, ...) @@ -355,8 +478,16 @@ class Current(LocalObservable): phi[ai:aj] += -1j * M_a.dot(H_ab.dot(psi[bi:bj])) """ + def __init__(self, syst, onsite=1, where=None, check_hermiticity=True): + super().__init__() + self.syst = syst + self.check_hermiticity = check_hermiticity + self.onsite = _normalize_onsite(syst, onsite) + self.where = _normalize_where(syst, where, hoppings=True) + def bind(self, args): - q = self.__class__() + cls = self.__class__ + q = cls.__new__(cls) q._bound = True q._bound_onsite = {a: ta.matrix(self.onsite(a, *args)) for a, b in self.where} @@ -405,6 +536,25 @@ class Current(LocalObservable): class Source(LocalObservable): """Generalized source observable. + Parameters + ---------- + syst : `~kwant.system.System` + onsite : scalar or square matrix or dict or callable + The onsite matrix that defines the charge. If a dict is given, it + maps from site families to either scalars or square matrices. If a + function is given it must take the same arguments as the onsite + Hamiltonian functions of the system. + where : sequence of `int` or `~kwant.builder.Site`s, or callable, optional + Where to evaluate the observable. If `syst` is not a finalized Builder, + then this should be a sequence of integers. If a function is provided, + it should take a single `int` or `~kwant.builder.Site` (if `syst` is a + finalized builder) and return True or False. If not provided, the + observable will be calculated over all sites in the system. + check_hermiticity: bool + Check whether the provided ``onsite`` is Hermitian. If it is not + Hermitian, then an error will be raised when the observable is + evaluated. + When this class is called in the following way:: K = kwant.physics.observables.Charge(fsyst, ...) @@ -453,16 +603,24 @@ class Source(LocalObservable): M_a.dot(H_aa)).dot(psi[i:j])) """ + def __init__(self, syst, onsite=1, where=None, check_hermiticity=True): + super().__init__() + self.syst = syst + self.check_hermiticity = check_hermiticity + self.onsite = _normalize_onsite(syst, onsite) + self.where = _normalize_where(syst, where, hoppings=False) + def bind(self, args): - q = self.__class__() + cls = self.__class__ + q = cls.__new__(cls) q._bound = True # if/when the loop over `where` is brought into `_act`, # this can be made into a sequence -- the lookups are the # same speed, but lists take less space - q._bound_onsite = {site: ta.matrix(self.onsite(site, *args)) - for site in self.where} - q._bound_hamiltonian = {site: ta.matrix(self.syst.hamiltonian(site, site, *args)) - for site in self.where} + q._bound_onsite = {a: ta.matrix(self.onsite(a, *args)) + for a, _ in self.where} + q._bound_hamiltonian = {a: ta.matrix(self.syst.hamiltonian(a, a, *args)) + for a, _ in self.where} q.syst = self.syst q.onsite = self.onsite q.where = self.where @@ -470,16 +628,18 @@ class Source(LocalObservable): return q # TODO: move this to Cython - def _act(self, bra, ket, site, args, inner_product): + def _act(self, bra, ket, sites, args, inner_product): # get the orbital slice - a_s, a_e = _get_slice(self.syst, site) + a, b = sites + assert a == b + a_s, a_e = _get_slice(self.syst, a) # extract the necessary matrix elements if self._bound: - M_a = self._bound_onsite[site] - H_aa = self._bound_hamiltonian[site] + M_a = self._bound_onsite[a] + H_aa = self._bound_hamiltonian[a] else: - M_a = ta.matrix(self.onsite(site, *args)) - H_aa = ta.matrix(self.syst.hamiltonian(site, site, *args)) + M_a = ta.matrix(self.onsite(a, *args)) + H_aa = ta.matrix(self.syst.hamiltonian(a, a, *args)) H_aa_dag = H_aa.conjugate().transpose() K = np.dot(H_aa_dag, M_a) - np.dot(M_a, H_aa) @@ -493,156 +653,3 @@ class Source(LocalObservable): return np.dot(bra[a_s:a_e].conjugate(), r) else: bra[a_s:a_e] += r - - - -################ Constructors for use with finalized Builders - -def _normalize_where(syst, where, hoppings): - if hoppings: - if where is None: - _where = tuple(syst.graph) - elif callable(where): - _where = tuple(filter(lambda h: where(*h), syst.graph)) - else: - _where = tuple((syst.id_by_site[a], syst.id_by_site[b]) - for a, b in where) - - if isinstance(syst, InfiniteSystem): - if any(a > syst.cell_size or b > syst.cell_size for a, b in _where): - raise ValueError('Only intra-cell hoppings may be specified ' - 'using `where`.') - else: - if where is None: - _where = tuple(range(len(syst.sites))) - elif callable(where): - _where = tuple(filter(where, syst.sites)) - else: - _where = tuple(syst.id_by_site[s] for s in where) - - if isinstance(syst, InfiniteSystem): - if any(w > syst.cell_size for w in _where): - raise ValueError('Only sites in the fundamental domain may be ' - 'specified using `where`.') - - return _where - - -def _normalize_onsite(syst, onsite): - _sites = syst.sites - if callable(onsite): - def _onsite(site_id, *args): - return onsite(_sites[site_id], *args) - elif isinstance(onsite, collections.Mapping): - def _onsite(site_id, *args): - return onsite[_sites[site_id].family] - else: - def _onsite(site_id, *args): - return onsite - - return _onsite - - -def charge(syst, onsite=1, where=None, check_hermiticity=True): - """Make an observable for calculating generalized charges. - - Parameters - ---------- - syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` - onsite : scalar or square matrix or dict or callable - The onsite matrix that defines the charge. If a dict is given, it - maps from site families to either scalars or square matrices. If a - function is given it must take the same arguments as the onsite - Hamiltonian functions of the system. - where : sequence of `~kwant.builder.Site`s or callable, optional - Where to evaluate the observable. If a function is provided, it - should take a single `~kwant.builder.Site` and return True or - False. If not provided, the observable will be calculated - everywhere in the system. - check_hermiticity: bool - Check whether the provided ``onsite`` is Hermitian. If it - is not Hermitian, then an error will be raised when the - observable is evaluated. - - Returns - ------- - `~kwant.physics.observables.Charge` - An object that can be called with a wavefunction and extra - arguments (which are passed to the ``onsite`` function), - and returns the generalized charge. - """ - q = Charge() - q.syst = syst - q.onsite = _normalize_onsite(syst, onsite) - q.where = _normalize_where(syst, where, False) - q.check_hermiticity = check_hermiticity - return q - - -def current(syst, onsite=1, where=None, check_hermiticity=True): - """Make an observable for calculating generalized currents. - - Parameters - ---------- - syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` - onsite : scalar or square matrix or dict or callable, optional - The onsite matrix that defines the charge from which this current - is derived. If a dict is given, it maps from site families to - either scalars or square matrices. If a function is given it must - take the same arguments as the onsite Hamiltonian functions of the - system. - where : sequence of pairs of `~kwant.builder.Site`s, optional - The hoppings on which to evaluate the observable. If not provided, - the observable will be calculated on every hopping in the system. - check_hermiticity: bool - Check whether the provided ``onsite`` is Hermitian. If it - is not Hermitian, then an error will be raised when the - observable is evaluated. - - Returns - ------- - `~kwant.physics.observables.Current` - An object that can be called with a wavefunction and extra - arguments (which are passed to the ``onsite`` function), - and returns the generalized current. - """ - j = Current() - j.syst = syst - j.onsite = _normalize_onsite(syst, onsite) - j.where = _normalize_where(syst, where, True) - j.check_hermiticity = check_hermiticity - return j - - -def source(syst, onsite=1, where=None, check_hermiticity=True): - """Make an observable for calculating generalized charge sources. - - Parameters - ---------- - syst : `~kwant.builder.FiniteSystem` or `~kwant.builder.InfiniteSystem` - onsite : scalar or square matrix or dict or callable - The onsite matrix that defines the charge. If a dict is given, it - maps from site families to either scalars or square matrices. If a - function is given it must take the same arguments as the onsite - Hamiltonian functions of the system. - where : sequence of `~kwant.builder.Site`s, optional - Where to evaluate the observable. If not provided, the - observable will be calculated everywhere in the system. - check_hermiticity: bool - Check whether the provided ``onsite`` is Hermitian. If it - is not Hermitian, then an error will be raised when the - observable is evaluated. - - Returns - ------- - `~kwant.physics.observables.Source` - An object that can be called with a wavefunction and extra - arguments (which are passed to the ``onsite`` function), - and returns the generalized charge source. - """ - s = Source() - s.syst = syst - s.onsite = _normalize_onsite(syst, onsite) - s.where = _normalize_where(syst, where, False) - s.check_hermiticity = check_hermiticity - return s diff --git a/kwant/tests/test_physics.py b/kwant/tests/test_physics.py index bb4fd6a..5f68d25 100644 --- a/kwant/tests/test_physics.py +++ b/kwant/tests/test_physics.py @@ -54,52 +54,55 @@ def test_observables_construction(): # test construction failure if norbs not given latnone = kwant.lattice.chain() syst[latnone(0)] = 1 - for a in (obs.charge, obs.current, obs.source): + for a in (obs.Charge, obs.Current, obs.Source): assert_raises(RuntimeError, a(syst.finalized()), np.zeros((N,))) del syst[latnone(0)] # test basic construction - Q = obs.charge(fsyst) - assert Q.where == tuple(range(N)) + Q = obs.Charge(fsyst) + assert Q.where == tuple(zip(range(N), range(N))) assert all(Q.onsite(i) == 1 for i in range(N)) - J = obs.current(fsyst) + J = obs.Current(fsyst) assert J.where == tuple(fsyst.graph) assert all(J.onsite(i) == 1 for i, j in fsyst.graph) # test construction with dict `onsite` - for A in (obs.charge, obs.current, obs.source): + for A in (obs.Charge, obs.Current, obs.Source): B = A(fsyst, {lat: 1}) assert all(B.onsite(i) == 1 for i in range(N)) # test construction with a functional onsite - for A in (obs.charge, obs.current, obs.source): + for A in (obs.Charge, obs.Current, obs.Source): B = A(fsyst, lambda site: site.pos[0]) # x-position operator assert all(B.onsite(i) == fsyst.sites[i].pos[0] for i in range(N)) # test construction with `where` given by a sequence where = [lat(2, 2), lat(1, 1)] fwhere = tuple(fsyst.id_by_site[s] for s in where) - A = obs.charge(fsyst, where=where) + fwhere = tuple(zip(fwhere, fwhere)) + A = obs.Charge(fsyst, where=where) assert A.where == fwhere where = [(lat(2, 2), lat(1, 2)), (lat(0, 0), lat(0, 1))] fwhere = tuple((fsyst.id_by_site[a], fsyst.id_by_site[b]) for a, b in where) - A = obs.current(fsyst, where=where) + A = obs.Current(fsyst, where=where) assert A.where == fwhere # test construction with `where` given by a function - where_list = set([(1, 0), (1, 1), (1, 2)]) + tag_list = [(1, 0), (1, 1), (1, 2)] def where(site): - return site.tag in where_list - A = obs.charge(fsyst, where=where) - assert all(w.tag in where_list for w in A.where) + return site.tag in tag_list + A = obs.Charge(fsyst, where=where) + assert all(fsyst.sites[w].tag in tag_list for w, _ in A.where) where_list = set(kwant.HoppingKind((1, 0), lat)(syst)) + fwhere_list = set((fsyst.id_by_site[a], fsyst.id_by_site[b]) + for a, b in where_list) def where(a, b): return (a, b) in where_list - A = obs.current(fsyst, where=where) - assert all(w in where_list for w in A.where) + A = obs.Current(fsyst, where=where) + assert all(w in fwhere_list for w in A.where) def _test_inner(val, A, bra, ket=None, args=()): @@ -131,9 +134,9 @@ def test_observables_finite(): fsyst = syst.finalized() ev, wfs = la.eigh(fsyst.hamiltonian_submatrix()) - Q = obs.charge(fsyst) - J = obs.current(fsyst) - K = obs.source(fsyst) + Q = obs.Charge(fsyst) + J = obs.Current(fsyst) + K = obs.Source(fsyst) for wf in wfs.T: # wfs[:, i] is i'th eigenvector assert np.isclose(np.sum(Q(wf)), 1) # eigenvectors are normalized @@ -150,11 +153,11 @@ def test_observables_infinite(): # Cannot calculate current between unit cells because the wavefunction # is only defined over a single unit cell - assert_raises(ValueError, obs.current, flead, where=[(lat(0, j), lat(1, j)) + assert_raises(ValueError, obs.Current, flead, where=[(lat(0, j), lat(1, j)) for j in range(3)]) transverse_kind = kwant.builder.HoppingKind((0, 1), lat) - J_intra = obs.current(flead, where=list(transverse_kind(lead))) + J_intra = obs.Current(flead, where=list(transverse_kind(lead))) prop, _ = flead.modes(energy=1.) for wf, v in zip(prop.wave_functions.T, prop.velocities): @@ -182,9 +185,9 @@ def test_observables_scattering(): fsyst = syst.finalized() # currents on the left and right of the disordered region - J_right = obs.current(fsyst, where=[(lat(N, j), lat(N-1, j)) + J_right = obs.Current(fsyst, where=[(lat(N, j), lat(N-1, j)) for j in range(3)]) - J_left = obs.current(fsyst, where=[(lat(0, j), lat(-1, j)) + J_left = obs.Current(fsyst, where=[(lat(0, j), lat(-1, j)) for j in range(3)]) smatrix = kwant.smatrix(fsyst, energy=1.0) @@ -196,7 +199,7 @@ def test_observables_scattering(): _test_inner(1 - np.sum(np.abs(rv)**2), J_left, wf) # Test `check_hermiticity=False` - for A in (obs.charge, obs.current, obs.source): + for A in (obs.Charge, obs.Current, obs.Source): assert_raises(ValueError, A(fsyst, 1j), wf) a = A(fsyst, 1j, check_hermiticity=False) b = a(wf, reduce=True) @@ -222,12 +225,12 @@ def test_observables_spin(): down, up = kwant.wave_function(fsyst, energy=1., args=args)(0) x_hoppings = kwant.builder.HoppingKind((1,), lat) - spin_current_z = obs.current(fsyst, sigmaz, where=x_hoppings(syst)) + spin_current_z = obs.Current(fsyst, sigmaz, where=x_hoppings(syst)) _test_per_element(1, spin_current_z, up, args=args) _test_per_element(-1, spin_current_z, down, args=args) # calculate spin_x torque - spin_torque_x = obs.source(fsyst, sigmax, where=[lat(L//2)]) + spin_torque_x = obs.Source(fsyst, sigmax, where=[lat(L//2)]) i = fsyst.id_by_site[lat(L//2)] psi = up[2*i:2*(i+1)] + down[2*i:2*(i+1)] H_ii = onsite(None, *args) @@ -284,11 +287,11 @@ def test_observables_gauged(): (Us[i], sigmaz, Us[i].conjugate().transpose())) x_hoppings = kwant.builder.HoppingKind((1,), lat) - spin_current_gauge = obs.current(fsyst, M_a, where=x_hoppings(syst)) + spin_current_gauge = obs.Current(fsyst, M_a, where=x_hoppings(syst)) _test_per_element(1, spin_current_gauge, up) _test_per_element(-1, spin_current_gauge, down) # check the reverse is also true minus_x_hoppings = kwant.builder.HoppingKind((-1,), lat) - spin_current_gauge = obs.current(fsyst, M_a, where=minus_x_hoppings(syst)) + spin_current_gauge = obs.Current(fsyst, M_a, where=minus_x_hoppings(syst)) _test_per_element(-1, spin_current_gauge, up) _test_per_element(1, spin_current_gauge, down) -- GitLab