diff --git a/kwant/builder.py b/kwant/builder.py
index e8a5f34c5deb1df3db2b167611098e31d200170a..ca6729e29cffa33fed6a32a87783e051e89f9b8b 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -830,7 +830,8 @@ class Builder:
# Each edge, specified by a ``(tail, head)`` pair of nodes, holds an object
# as a value. Likewise, each tail which occurs in the graph also holds a
# value. (Nodes which only occur as heads are not required to have
- # values.)
+ # values.) Every tail node has to be in the fundamental domain of the
+ # builder's symmetry.
#
# For a given `tail` site, H[tail] is a list alternately storing
# heads and values. (The heads occupy even locations followed by the
diff --git a/kwant/tests/test_wraparound.py b/kwant/tests/test_wraparound.py
index c1233a97cb71aef6e31b240bac616081b6089f0d..18585275895f7d4c34f91c86868139ed3a1dc505 100644
--- a/kwant/tests/test_wraparound.py
+++ b/kwant/tests/test_wraparound.py
@@ -300,10 +300,9 @@ def test_fd_mismatch():
else:
def _(*args):
args = list(args)
- assert len(args) == 2
kext = args.pop(keep)
- kint = args[0]
- B = kwant.physics.Bands(syst, args=(kint,))
+ kint = args
+ B = kwant.physics.Bands(syst, args=kint)
return B(kext)
return _
@@ -312,3 +311,82 @@ def test_fd_mismatch():
# is more tricky, as we would also need to transform the k vectors
E_k = np.array([spec(0, 0) for spec in spectra]).transpose()
assert all(np.allclose(E, E[0]) for E in E_k)
+
+ # Test square lattice with oblique unit cell
+ lat = kwant.lattice.general(np.eye(2))
+ translations = kwant.lattice.TranslationalSymmetry([2, 2], [0, 2])
+ syst = kwant.Builder(symmetry=translations)
+ syst[lat.shape(lambda site: True, [0, 0])] = 1
+ syst[lat.neighbors()] = 1
+ # Check that spectra are the same at k=0.
+ spectra = [spectrum(syst, keep) for keep in (None, 0, 1)]
+ E_k = np.array([spec(0, 0) for spec in spectra]).transpose()
+ assert all(np.allclose(E, E[0]) for E in E_k)
+
+ # Test Rocksalt structure
+ # cubic lattice that contains both sublattices
+ lat = kwant.lattice.general(np.eye(3))
+ # Builder with FCC translational symmetries.
+ translations = kwant.lattice.TranslationalSymmetry([1, 1, 0], [1, 0, 1], [0, 1, 1])
+ syst = kwant.Builder(symmetry=translations)
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(0, 0, 1)] = -1
+ syst[lat.neighbors()] = 1
+ # Check that spectra are the same at k=0.
+ spectra = [spectrum(syst, keep) for keep in (None, 0, 1, 2)]
+ E_k = np.array([spec(0, 0, 0) for spec in spectra]).transpose()
+ assert all(np.allclose(E, E[0]) for E in E_k)
+ # Same with different translation vectors
+ translations = kwant.lattice.TranslationalSymmetry([1, 1, 0], [1, -1, 0], [0, 1, 1])
+ syst = kwant.Builder(symmetry=translations)
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(0, 0, 1)] = -1
+ syst[lat.neighbors()] = 1
+ # Check that spectra are the same at k=0.
+ spectra = [spectrum(syst, keep) for keep in (None, 0, 1, 2)]
+ E_k = np.array([spec(0, 0, 0) for spec in spectra]).transpose()
+ assert all(np.allclose(E, E[0]) for E in E_k)
+
+ # Test that spectrum in slab geometry is identical regardless of choice of unit
+ # cell in rocksalt structure
+ def shape(site):
+ return abs(site.tag[2]) < 4
+
+ lat = kwant.lattice.general(np.eye(3))
+ # First choice: primitive UC
+ translations = kwant.lattice.TranslationalSymmetry([1, 1, 0], [1, -1, 0], [1, 0, 1])
+ syst = kwant.Builder(symmetry=translations)
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(0, 0, 1)] = -1
+ # Set all the nearest neighbor hoppings
+ for d in [[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]]:
+ syst[(lat(0, 0, 0), lat(*d))] = 1
+
+ wrapped = kwant.wraparound.wraparound(syst, keep=2)
+ finitewrapped = kwant.Builder()
+ finitewrapped.fill(wrapped, shape, start=np.zeros(3));
+
+ sysf = finitewrapped.finalized()
+ spectrum1 = [np.linalg.eigvalsh(sysf.hamiltonian_submatrix(args=(k, 0)))
+ for k in np.linspace(-np.pi, np.pi, 5)]
+
+ # Second choice: doubled UC with third translation purely in z direction
+ translations = kwant.lattice.TranslationalSymmetry([1, 1, 0], [1, -1, 0], [0, 0, 2])
+ syst = kwant.Builder(symmetry=translations)
+ syst[lat(0, 0, 0)] = 1
+ syst[lat(0, 0, 1)] = -1
+ syst[lat(1, 0, 1)] = 1
+ syst[lat(-1, 0, 0)] = -1
+ for s in np.array([[0, 0, 0], [1, 0, 1]]):
+ for d in np.array([[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]]):
+ syst[(lat(*s), lat(*(s + d)))] = 1
+ wrapped = kwant.wraparound.wraparound(syst, keep=2)
+ finitewrapped = kwant.Builder()
+
+ finitewrapped.fill(wrapped, shape, start=np.zeros(3));
+
+ sysf = finitewrapped.finalized()
+ spectrum2 = [np.linalg.eigvalsh(sysf.hamiltonian_submatrix(args=(k, 0)))
+ for k in np.linspace(-np.pi, np.pi, 5)]
+
+ assert np.allclose(spectrum1, spectrum2)
diff --git a/kwant/wraparound.py b/kwant/wraparound.py
index 2e1cb1dd1b63b2f7cd051a9687616724b81a6759..d43d6cc8f343709d97f57cf417a8be72b80afd8a 100644
--- a/kwant/wraparound.py
+++ b/kwant/wraparound.py
@@ -269,22 +269,29 @@ def wraparound(builder, keep=None, *, coordinate_names=('x', 'y', 'z')):
mnp = -len(sym.periods) # Used by the bound functions above.
# Store lists of values, so that multiple values can be assigned to the
- # same site or hopping. We map to the FD of 'sym', as the hopping-processing
- # code assumes the sites are in the FD.
+ # same site or hopping.
for site, val in builder.site_value_pairs():
- sites[sym.to_fd(site)] = [bind_site(val) if callable(val) else val]
+ # Every 'site' is in the FD of the original symmetry.
+ # Move the sites to the FD of the remaining symmetry, this guarantees that
+ # every site in the new system is an image of an original FD site translated
+ # purely by the remaining symmetry.
+ sites[ret.symmetry.to_fd(site)] = [bind_site(val) if callable(val) else val]
for hop, val in builder.hopping_value_pairs():
- # Map hopping to FD of 'sym', as the code afterwards assumes that 'a'
- # is in this domain.
- a, b = sym.to_fd(*hop)
- b_dom = sym.which(b)
+ a, b = hop
+ # 'a' is in the FD of original symmetry.
+ # Get translation from FD of original symmetry to 'b',
+ # this is different from 'b_dom = sym.which(b)'.
+ b_dom = builder.symmetry.which(b)
+ # Throw away part that is in the remaining translation direction, so we get
+ # an element of 'sym' which is being wrapped
+ b_dom = ta.array([t for i, t in enumerate(b_dom) if i != keep])
+ # Pull back using the remainder, which is purely in the wrapped directions.
+ # This guarantees that 'b_wa' is an image of an original FD site translated
+ # purely by the remaining symmetry.
b_wa = sym.act(-b_dom, b)
- # Now map 'b_wa' to another domain of 'sym', so that the hopping
- # is compatible with 'ret.symmetry'. This is necessary when the
- # fundamental domains of the symmetries do not coincide.
- w = ret.symmetry.which(b_wa)
- b_wa = ret.symmetry.act(w, sym.to_fd(ret.symmetry.act(-w, b_wa)))
+ # Move the hopping to the FD of the remaining symmetry
+ a, b_wa = ret.symmetry.to_fd(a, b_wa)
if a == b_wa:
# The hopping gets wrapped-around into an onsite Hamiltonian.
@@ -292,13 +299,14 @@ def wraparound(builder, keep=None, *, coordinate_names=('x', 'y', 'z')):
sites[a].append(bind_hopping_as_site(b_dom, val))
else:
# The hopping remains a hopping.
- if b != b_wa or callable(val):
+ if any(b_dom) or callable(val):
# The hopping got wrapped-around or is a function.
val = bind_hopping(b_dom, val)
# Make sure that there is only one entry for each hopping
- # (pointing in one direction).
- if (b_wa, a) in hops:
+ # pointing in one direction, modulo the remaining translations.
+ b_wa_r, a_r = ret.symmetry.to_fd(b_wa, a)
+ if (b_wa_r, a_r) in hops:
assert (a, b_wa) not in hops
if callable(val):
assert not isinstance(val, HermConjOfFunc)
@@ -306,7 +314,7 @@ def wraparound(builder, keep=None, *, coordinate_names=('x', 'y', 'z')):
else:
val = herm_conj(val)
- hops[b_wa, a].append(val)
+ hops[b_wa_r, a_r].append(val)
else:
hops[a, b_wa].append(val)