diff --git a/doc/source/pre/whatsnew/1.5.rst b/doc/source/pre/whatsnew/1.5.rst
index 9340cec96666ac0b9f02b7fa92ab1471ea63dc78..a509ccc0a16dc525e8330671790459ae5eaa2040 100644
--- a/doc/source/pre/whatsnew/1.5.rst
+++ b/doc/source/pre/whatsnew/1.5.rst
@@ -3,6 +3,30 @@ What's new in Kwant 1.5
This article explains the user-visible changes in Kwant 1.5.0.
+Automatic addition of Peierls phase terms to Builders
+-----------------------------------------------------
+Kwant 1.4 introduced `kwant.physics.magnetic_gauge` for computing Peierls
+phases for arbitrary geometries and for systems with leads. Using this
+functionality requires that the system value functions are equipped to
+take the required Peierls phase parameters, which is not possible when
+you are not in direct control of the value functions (e.g. discretized
+systems). In Kwant 1.5 we have added the missing piece of functionality,
+`kwant.builder.add_peierls_phase`, which properly adds the Peierls phase
+parameters to a Builder's value functions::
+
+ syst = make_system()
+ lead = make_lead()
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ fsyst, phase = kwant.builder.add_peierls_phase(syst)
+
+ def B_syst(pos):
+ return np.exp(-np.sum(pos * pos))
+
+ kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
+
+
Improved tutorial building
--------------------------
Improving or adding to Kwant's tutorial is now much simpler. Now
diff --git a/doc/source/reference/kwant.builder.rst b/doc/source/reference/kwant.builder.rst
index af02a1bc27e3bbeb3adc4f82e0594c7604b7ab07..327afca9e9f44c797378b918baa0799ef30af45f 100644
--- a/doc/source/reference/kwant.builder.rst
+++ b/doc/source/reference/kwant.builder.rst
@@ -26,3 +26,10 @@ Abstract base classes
SiteFamily
Symmetry
Lead
+
+Functions
+---------
+.. autosummary::
+ :toctree: generated/
+
+ add_peierls_phase
diff --git a/kwant/_common.py b/kwant/_common.py
index e4623b8d09c0d8a307fcab8dff2b6571fb4112c3..2492e6daa1272ddc24c4fe6b66290612901035ec 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -13,6 +13,7 @@ import inspect
import warnings
import importlib
import functools
+import collections
from contextlib import contextmanager
__all__ = ['KwantDeprecationWarning', 'UserCodeError']
@@ -191,3 +192,27 @@ class lazy_import:
package = sys.modules[self.__package]
setattr(package, self.__module, mod)
return getattr(mod, name)
+
+
+def _hashable(obj):
+ return isinstance(obj, collections.abc.Hashable)
+
+
+def memoize(f):
+ """Decorator to memoize a function that works even with unhashable args.
+
+ This decorator will even work with functions whose args are not hashable.
+ The cache key is made up by the hashable arguments and the ids of the
+ non-hashable args. It is up to the user to make sure that non-hashable
+ args do not change during the lifetime of the decorator.
+
+ This decorator will keep reevaluating functions that return None.
+ """
+ def lookup(*args):
+ key = tuple(arg if _hashable(arg) else id(arg) for arg in args)
+ result = cache.get(key)
+ if result is None:
+ cache[key] = result = f(*args)
+ return result
+ cache = {}
+ return lookup
diff --git a/kwant/builder.py b/kwant/builder.py
index 2ca81a946faf5c9c50201c67a34398a69d46d41c..89ec7c24c68393a9916d6d314925e3783168d7b7 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -13,6 +13,7 @@ import collections
import copy
from functools import total_ordering, wraps, update_wrapper
from itertools import islice, chain
+import textwrap
import inspect
import tinyarray as ta
import numpy as np
@@ -20,13 +21,14 @@ from scipy import sparse
from . import system, graph, KwantDeprecationWarning, UserCodeError
from .linalg import lll
from .operator import Density
-from .physics import DiscreteSymmetry
+from .physics import DiscreteSymmetry, magnetic_gauge
from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
- interleave, deprecate_args)
+ interleave, deprecate_args, memoize)
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
- 'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead']
+ 'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead',
+ 'add_peierls_phase']
################ Sites and site families
@@ -1803,6 +1805,169 @@ class Builder:
inter_cell_hopping = cell_hamiltonian = precalculated = \
_require_system
+
+def _add_peierls_phase(syst, peierls_parameter):
+
+ @memoize
+ def wrap_hopping(hop):
+ def f(*args):
+ a, b, *args, phi = args
+ return hop(a, b, *args) * phi(a, b)
+
+ params = ('_site0', '_site1')
+ if callable(hop):
+ params += get_parameters(hop)[2:] # cut off both site parameters
+ params += (peierls_parameter,)
+
+ f.__signature__ = inspect.Signature(
+ inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
+ for name in params)
+
+ return f
+
+ const_hoppings = {}
+
+ def const_hopping(a, b, phi):
+ return const_hoppings[a, b] * phi(a, b)
+
+ # fix the signature
+ params = ('_site0', '_site1', peierls_parameter)
+ const_hopping.__signature__ = inspect.Signature(
+ inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
+ for name in params)
+
+ ret = Builder(
+ syst.symmetry,
+ # No time reversal, as magnetic fields break it
+ conservation_law=syst.conservation_law,
+ particle_hole=syst.particle_hole,
+ chiral=syst.chiral,
+ )
+ for i, lead in enumerate(syst.leads):
+ ret.leads.append(
+ BuilderLead(
+ _add_peierls_phase(
+ lead.builder,
+ peierls_parameter + '_lead{}'.format(i),
+ ),
+ lead.interface,
+ lead.padding,
+ )
+ )
+
+ for site, val in syst.site_value_pairs():
+ ret[site] = val
+
+ for hop, val in syst.hopping_value_pairs():
+ if callable(val):
+ ret[hop] = wrap_hopping(val)
+ else:
+ const_hoppings[hop] = val
+ ret[hop] = const_hopping
+
+ return ret
+
+
+def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
+ """Add a Peierls phase parameter to a Builder.
+
+ Parameters
+ ----------
+ syst : kwant.Builder
+ The system to which to add a Peierls phase.
+ peierls_parameter : str
+ The name of the Peierls phase parameter to value functions.
+ fix_gauge : bool (default: True)
+ If True, fix the gauge of the system and return it also
+
+ Returns
+ -------
+ syst : kwant.builder.Finitesystem or kwant.builder.InfiniteSystem
+ A system where all hoppings are given by value functions
+ that have an additional parameter 'phi' that should be
+ the Peierls phase returned by `kwant.physics.magnetic_gauge`.
+ Any leads have similar parameters named 'phi_lead0',
+ 'phi_lead1' etc.
+ gauge : callable
+ Only returned if 'fix_gauge' is True. Called with magnetic
+ field(s), and returns a parameter dictionary that can be
+ passed to the system as 'params' (see the example below).
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import kwant
+ >>>
+ >>> syst = make_system()
+ >>> lead = make_lead()
+ >>> syst.attach_lead(lead)
+ >>> syst.attach_lead(lead.reversed())
+ >>>
+ >>> fsyst, phase = kwant.builder.add_peierls_phase(syst)
+ >>>
+ >>> def B_syst(pos):
+ >>> return np.exp(-np.sum(pos * pos))
+ >>>
+ >>> kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
+ """
+
+ def wrap_gauge(gauge):
+ phase_names = [peierls_parameter]
+ phase_names += [peierls_parameter + '_lead{}'.format(i)
+ for i in range(len(syst.leads))]
+
+ doc = textwrap.dedent("""\
+ Return the Peierls phase for a magnetic field configuration.
+
+ Parameters
+ ----------
+ syst_field : scalar, vector or callable
+ The magnetic field to apply to the scattering region.
+ If callable, takes a position and returns the
+ magnetic field at that position. Can be a scalar if
+ the system is 1D or 2D, otherwise must be a vector.
+ Magnetic field is expressed in units :math:`φ₀ / l²`,
+ where :math:`φ₀` is the magnetic flux quantum and
+ :math:`l` is the unit of length.
+ *lead_fields : scalar, vector or callable
+ The magnetic fields to apply to each of the leads, in
+ the same format as 'syst_field'. In addition, if a callable
+ is provided, then the magnetic field must have the symmetry
+ of the associated lead.
+ tol : float, default: 1E-8
+ The tolerance to which to calculate the flux through each
+ hopping loop in the system.
+ average : bool, default: False
+ If True, estimate the magnetic flux through each hopping loop
+ in the system by evaluating the magnetic field at a single
+ position inside the loop and multiplying it by the area of the
+ loop. If False, then ``scipy.integrate.quad`` is used to integrate
+ the magnetic field. This parameter is only used when 'syst_field'
+ or 'lead_fields' are callable.
+
+ Returns
+ -------
+ phases : dict of callables
+ To be passed directly as the parameters of a system wrapped with
+ 'kwant.builder.add_peierls_phase'.
+ """)
+
+ @wraps(gauge)
+ def f(*args, **kwargs):
+ phases = gauge(*args, **kwargs)
+ return dict(zip(phase_names, phases))
+
+ f.__doc__ = doc
+
+ return f
+
+ ret = _add_peierls_phase(syst, peierls_parameter).finalized()
+
+ if fix_gauge:
+ return ret, wrap_gauge(magnetic_gauge(ret))
+ else:
+ return ret
+
################ Finalized systems
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 84e3ef024721e9812d4b7e502a09eb440853dd83..3885940a6147bbaed2720949c0814459286abe09 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1416,3 +1416,28 @@ def test_finalization_preserves_padding():
syst = syst.finalized()
# The order is guaranteed because the paddings are sorted.
assert [syst.sites[i] for i in syst.lead_paddings[0]] == padding[:-1]
+
+
+def test_add_peierls_phase():
+
+ lat = kwant.lattice.square()
+ syst = kwant.Builder()
+ syst[(lat(i, j) for i in range(5) for j in range(5))] = 4
+ syst[lat.neighbors()] = lambda a, b, t: -t
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+ lead[(lat(0, j) for j in range(5))] = 4
+ lead[lat.neighbors()] = -1
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ syst, phase = builder.add_peierls_phase(syst)
+
+ assert isinstance(syst, builder.FiniteSystem)
+
+ params = phase(1, 0, 0)
+
+ assert all(p in params for p in ('phi', 'phi_lead0', 'phi_lead1'))
+
+ kwant.smatrix(syst, energy=0.1, params=dict(t=-1, **params))
diff --git a/kwant/wraparound.py b/kwant/wraparound.py
index c3c5a38cccb9f3a8562bb894df52d7c791d3a955..25217826e4c653325072c7d325f6eae90a92d2fb 100644
--- a/kwant/wraparound.py
+++ b/kwant/wraparound.py
@@ -20,36 +20,12 @@ from . import builder, system, plotter
from .linalg import lll
from .builder import herm_conj, HermConjOfFunc
from .lattice import TranslationalSymmetry
-from ._common import get_parameters
+from ._common import get_parameters, memoize
__all__ = ['wraparound', 'plot_2d_bands']
-def _hashable(obj):
- return isinstance(obj, collections.abc.Hashable)
-
-
-def _memoize(f):
- """Decorator to memoize a function that works even with unhashable args.
-
- This decorator will even work with functions whose args are not hashable.
- The cache key is made up by the hashable arguments and the ids of the
- non-hashable args. It is up to the user to make sure that non-hashable
- args do not change during the lifetime of the decorator.
-
- This decorator will keep reevaluating functions that return None.
- """
- def lookup(*args):
- key = tuple(arg if _hashable(arg) else id(arg) for arg in args)
- result = cache.get(key)
- if result is None:
- cache[key] = result = f(*args)
- return result
- cache = {}
- return lookup
-
-
def _set_signature(func, params):
"""Set the signature of 'func'.
@@ -103,7 +79,7 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
format. It will be deprecated in the 2.0 release of Kwant.
"""
- @_memoize
+ @memoize
def bind_site(val):
def f(*args):
a, *args = args
@@ -113,7 +89,7 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
_set_signature(f, get_parameters(val) + momenta)
return f
- @_memoize
+ @memoize
def bind_hopping_as_site(elem, val):
def f(*args):
a, *args = args
@@ -128,7 +104,7 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
_set_signature(f, params + momenta)
return f
- @_memoize
+ @memoize
def bind_hopping(elem, val):
def f(*args):
a, b, *args = args
@@ -142,7 +118,7 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
_set_signature(f, params + momenta)
return f
- @_memoize
+ @memoize
def bind_sum(num_sites, *vals):
"""Construct joint signature for all 'vals'."""