NeighborRange default for hopping range

README.md

@@ -17,7 +17,7 @@ The Quantica.jl package provides an expressive API to build arbitrary quantum sy
 # Exported API
 - `lattice`, `sublat`, `bravais`: build lattices
 - `dims`, `sitepositions`, `siteindices`: inspect lattices
-- `hopping`, `onsite`, `siteselector`, `hopselector`: build tightbinding models
+- `hopping`, `onsite`, `siteselector`, `hopselector`, `nrange`: build tightbinding models
 - `hamiltonian`: build a Hamiltonian from tightbinding model and a lattice
 - `parametric`, `@onsite!`, `@hopping!`, `parameters`: build a parametric Hamiltonian
 - `supercell`, `unitcell`, `flatten`, `wrap`, `transform!`, `combine`: build derived lattices or Hamiltonians

src/Quantica.jl

@@ -18,7 +18,7 @@ using ExprTools
 using SparseArrays: getcolptr, AbstractSparseMatrix
 
 export sublat, bravais, lattice, dims, supercell, unitcell,
-       hopping, onsite, @onsite!, @hopping!, parameters, siteselector, hopselector,
+       hopping, onsite, @onsite!, @hopping!, parameters, siteselector, hopselector, nrange,
        sitepositions, siteindices,
        ket, randomkets,
        hamiltonian, parametric, bloch, bloch!, optimize!, similarmatrix,

src/hamiltonian.jl

@@ -603,22 +603,24 @@ function applyterm!(builder::IJVBuilder{L}, term::HoppingTerm) where {L}
     for (s2, s1) in sublats(rsel)  # Each is a Pair s2 => s1
         dns = dniter(rsel)
         for dn in dns
-            foundlink = false
+            keepgoing = false
             ijv = builder[dn]
             for j in source_candidates(rsel, s2)
                 sitej = allpos[j]
                 rsource = sitej - lat.bravais.matrix * dn
                 is = targets(builder, rsel.selector.range, rsource, s1)
                 for i in is
+                    # Make sure we don't stop searching until we reach minimum range
+                    is_below_min_range((i, j), (dn, zero(dn)), rsel) && (keepgoing = true)
                     ((i, j), (dn, zero(dn))) in rsel || continue
-                    foundlink = true
+                    keepgoing = true
                     rtarget = allsitepositions(lat)[i]
                     r, dr = _rdr(rsource, rtarget)
                     v = toeltype(term(r, dr), eltype(builder), builder.orbs[s1], builder.orbs[s2])
                     push!(ijv, (i, j, v))
                 end
             end
-            foundlink && acceptcell!(dns, dn)
+            keepgoing && acceptcell!(dns, dn)
         end
     end
     return nothing

src/lattice.jl

@@ -472,6 +472,9 @@ siterange(lat::AbstractLattice, sublat) = siterange(lat.unitcell, sublat)
 
 allsitepositions(lat::AbstractLattice) = sitepositions(lat.unitcell)
 
+siteposition(i, lat::AbstractLattice) = allsitepositions(lat)[i]
+siteposition(i, dn::SVector, lat::AbstractLattice) = siteposition(i, lat) + bravais(lat) * dn
+
 offsets(lat::AbstractLattice) = offsets(lat.unitcell)
 
 sublatsites(lat::AbstractLattice) = sublatsites(lat.unitcell)

src/model.jl

@@ -1,5 +1,85 @@
 using Quantica.RegionPresets: Region
 
+#######################################################################
+# NeighborRange
+#######################################################################
+struct NeighborRange
+    n::Int
+end
+
+"""
+    nrange(n::Int)
+
+Create a `NeighborRange` that represents a hopping range to distances corresponding to the
+n-th nearest neighbors in a given lattice. Such distance is obtained by finding the n-th
+closest pairs of sites in a lattice, irrespective of their sublattice.
+
+    nrange(n::Int, lat::AbstractLattice)
+
+Obtain the actual nth-nearest-neighbot distance between sites in lattice `lat`.
+
+# See also:
+    `hopping`
+"""
+nrange(n::Int) = NeighborRange(n)
+
+function nrange(n, lat::AbstractLattice{E,L}) where {E,L}
+    sites = allsitepositions(lat)
+    T = eltype(first(sites))
+    dns = BoxIterator(zero(SVector{L,Int}))
+    br = bravais(lat)
+    # 640 is a heuristic cutoff for kdtree vs brute-force search
+    if length(sites) <= 128
+        dists = fill(T(Inf), n)
+        for dn in dns
+            iszero(dn) || ispositive(dn) || continue
+            for (i, ri) in enumerate(sites), (j, rj) in enumerate(sites)
+                j <= i && iszero(dn) && continue
+                r = ri - rj + br * dn
+                _update_dists!(dists, r'r)
+            end
+            isfinite(last(dists)) || acceptcell!(dns, dn)
+        end
+        dist = sqrt(last(dists))
+    else
+        tree = KDTree(sites)
+        dist = T(Inf)
+        for dn in dns
+            iszero(dn) || ispositive(dn) || continue
+            for r0 in sites
+                r = r0 + br * dn
+                dist = min(dist, _nrange(n, tree, r, nsites(lat)))
+            end
+            isfinite(dist) || acceptcell!(dns, dn)
+        end
+    end
+    return dist
+end
+
+function _update_dists!(dists, dist::Real)
+    len = length(dists)
+    for (n, d) in enumerate(dists)
+        isapprox(dist, d) && break
+        if dist < d
+            dists[n+1:len] .= dists[n:len-1]
+            dists[n] = dist
+            break
+        end
+    end
+    return dists
+end
+
+function _nrange(n, tree, r::AbstractVector{T}, nmax) where {T}
+    for m in n:nmax
+        _, dists = knn(tree, r, 1 + m, true)
+        popfirst!(dists)
+        unique_sorted_approx!(dists)
+        length(dists) == n && return maximum(dists)
+    end
+    return T(Inf)
+end
+
+
 #######################################################################
 # Onsite/Hopping selectors
 #######################################################################
@@ -55,7 +135,7 @@ siteselector(; region = missing, sublats = missing, indices = missing) =
     SiteSelector(region, sublats, indices)
 
 """
-    hopselector(; region = missing, sublats = missing, indices = missing, dn = missing, range = missing)
+    hopselector(; range = missing, dn = missing, sublats = missing, indices = missing, region = missing)
 
 Return a `HopSelector` object that can be used to select hops between two sites in a
 lattice. Only hops between two sites, with indices `ipair = src => dst`, at positions `r₁ =
@@ -66,6 +146,14 @@ sublattices `s₁` and `s₂` will be selected if:
 
 If any of these is `missing` it will not be used to constraint the selection.
 
+The keyword `range` admits the following possibilities
+
+    max_range                   # i.e. `norm(dr) <= max_range`
+    (min_range, max_range)      # i.e. `min_range <= norm(dr) <= max_range`
+
+Both `max_range` and `min_range` can be a `Real` or a `NeighborRange` created with
+`nrange(n)`. The latter represents the distance of `n`-th nearest neighbors.
+
 The keyword `dn` can be a `Tuple`/`Vector`/`SVector` of `Int`s, or a tuple thereof.
 
 The keyword `sublats` allows the following formats:
@@ -85,10 +173,6 @@ The keyword `indices` accepts a single `src => dest` pair or a collection thereo
     indices = [(1, 2) .=> (2, 1)]       # Broadcasted pairs, same as above
     indices = [1:10 => 20:25, 3 => 30]  # Direct product, any hopping from sites 1:10 to sites 20:25, or from 3 to 30
 
-The keyword `range` can be a number `range = max_range` or an interval `range = (min_range,
-max_range)`. In the latter case, only links with `min_range <= norm(dr) <= max_range` are
-selected.
-
 """
 hopselector(; region = missing, sublats = missing, dn = missing, range = missing, indices = missing) =
     HopSelector(region, sublats, sanitize_dn(dn), sanitize_range(range), indices)
@@ -104,8 +188,11 @@ _sanitize_dn(dn::SVector{N}) where {N} = SVector{N,Int}(dn)
 _sanitize_dn(dn::Vector) = SVector{length(dn),Int}(dn)
 
 sanitize_range(::Missing) = missing
-sanitize_range(range::Real) = ifelse(isfinite(range), float(range) + sqrt(eps(float(range))), float(range))
-sanitize_range((rmin, rmax)::Tuple{Real,Real}) = (-sanitize_range(-rmin), sanitize_range(rmax))
+sanitize_range(r) = _shift_eps(r, 1)
+sanitize_range(r::NTuple{2,Any}) = (_shift_eps(first(r), -1), _shift_eps(last(r), 1))
+
+_shift_eps(r::Real, m) = ifelse(isfinite(r), float(r) + m * sqrt(eps(float(r))), float(r))
+_shift_eps(r, m) = r
 
 # API
 
@@ -115,13 +202,18 @@ function resolve(s::SiteSelector, lat::AbstractLattice)
 end
 
 function resolve(s::HopSelector, lat::AbstractLattice)
-    s = HopSelector(s.region, resolve_sublat_pairs(s.sublats, lat), check_dn_dims(s.dns, lat), s.range, s.indices)
+    s = HopSelector(s.region, resolve_sublat_pairs(s.sublats, lat), check_dn_dims(s.dns, lat), resolve_range(s.range, lat), s.indices)
     return ResolvedSelector(s, lat)
 end
 
 resolve_sublats(::Missing, lat) = missing # must be resolved to iterate over sublats
 resolve_sublats(s, lat) = resolve_sublat_name.(s, Ref(lat))
 
+resolve_range(r::Tuple, lat) = sanitize_range(_resolve_range.(r, Ref(lat)))
+resolve_range(r, lat) = sanitize_range(_resolve_range(r, lat))
+_resolve_range(r::NeighborRange, lat) = nrange(r.n, lat)
+_resolve_range(r, lat) = r
+
 function resolve_sublat_name(name::Union{NameType,Integer}, lat)
     i = findfirst(isequal(name), lat.unitcell.names)
     return i === nothing ? 0 : i
@@ -194,6 +286,12 @@ isinrange((row, col), (dnrow, dncol), range, lat) =
 _isinrange(p, rmax::Real) = p'p <= rmax^2
 _isinrange(p, (rmin, rmax)::Tuple{Real,Real}) =  rmin^2 <= p'p <= rmax^2
 
+is_below_min_range(inds, dns, rsel::ResolvedSelector) =
+    is_below_min_range(inds, dns, rsel.selector.range, rsel.lattice)
+is_below_min_range((i, j), (dni, dnj), (rmin, rmax)::Tuple, lat) =
+    norm(siteposition(i, dni, lat) - siteposition(j, dnj, lat)) < rmin
+is_below_min_range(inds, dn, range, lat) = false
+
 # There are no sublat ranges, so supporting (:A, (:B, :C)) is not necessary
 isinsublats(s::Integer, ::Missing) = true
 isinsublats(s::Integer, sublats) = s in sublats
@@ -371,6 +469,11 @@ sublats(m::TightbindingModel) = (t -> t.selector.sublats).(terms(m))
 displayparameter(::Type{<:Function}) = "Function"
 displayparameter(::Type{T}) where {T} = "$T"
 
+displayrange(r::Real) = round(r, digits = 6)
+displayrange(::Missing) = "any"
+displayrange(nr::NeighborRange) = "NeighborRange($(nr.n))"
+displayrange(rs::Tuple) = "($(displayrange(first(rs))), $(displayrange(last(rs))))"
+
 function Base.show(io::IO, o::OnsiteTerm{F}) where {F}
     i = get(io, :indent, "")
     print(io,
@@ -385,7 +488,7 @@ function Base.show(io::IO, h::HoppingTerm{F}) where {F}
 "$(i)HoppingTerm{$(displayparameter(F))}:
 $(i)  Sublattice pairs : $(h.selector.sublats === missing ? "any" : h.selector.sublats)
 $(i)  dn cell distance : $(h.selector.dns === missing ? "any" : h.selector.dns)
-$(i)  Hopping range    : $(round(h.selector.range, digits = 6))
+$(i)  Hopping range    : $(displayrange(h.selector.range))
 $(i)  Coefficient      : $(h.coefficient)")
 end
 
@@ -480,7 +583,7 @@ _onlyonsites(s, t::HoppingTerm, args...) = (_onlyonsites(s, args...)...,)
 _onlyonsites(s) = ()
 
 """
-    hopping(t; region = missing, sublats = missing, indices = missing, dn = missing, range = 1, plusadjoint = false)
+    hopping(t; range = nrange(1), dn = missing, sublats = missing, indices = missing, region = missing, plusadjoint = false)
 
 Create an `TightbindingModel` with a single `HoppingTerm` that applies a hopping `t` to a
 `Lattice` when creating a `Hamiltonian` with `hamiltonian`.
@@ -512,9 +615,15 @@ positions `r₁ = r - dr/2` and `r₂ = r + dr`, belonging to unit cells at inte
 && dn´ in dn && norm(dr) <= range`. If any of these is `missing` it will not be used to
 constraint the selection.
 
-The keyword `range` defaults to `1` (instead of `missing`). `range` also allows a form
-`range = (min_range, max_range)`, in which case the corresponding selection condition
-becomes `min_range <= norm(dr) <= max_range`.
+The keyword `range` admits the following possibilities
+
+    max_range                   # i.e. `norm(dr) <= max_range`
+    (min_range, max_range)      # i.e. `min_range <= norm(dr) <= max_range`
+
+Both `max_range` and `min_range` can be a `Real` or a `NeighborRange` created with
+`nrange(n)`. The latter represents the distance of `n`-th nearest neighbors. Note that the
+`range` default for `hopping` (unlike for the more general `hopselector`) is `nrange(1)`,
+i.e. first-nearest-neighbors.
 
 The keyword `dn` can be a `Tuple`/`Vector`/`SVector` of `Int`s, or a tuple thereof. The
 keyword `sublats` allows the following formats:
@@ -572,9 +681,9 @@ Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
 ```
 
 # See also:
-    `onsite`
+    `onsite`, `nrange`
 """
-function hopping(t; plusadjoint = false, range = 1, kw...)
+function hopping(t; plusadjoint = false, range = nrange(1), kw...)
     hop = hopping(t, hopselector(; range = range, kw...))
     return plusadjoint ? hop + hop' : hop
 end

src/tools.jl

@@ -159,6 +159,21 @@ end
 chop(x::T, x0 = one(T)) where {T<:Real} = ifelse(abs(x) < √eps(T(x0)), zero(T), x)
 chop(x::C, x0 = one(R)) where {R<:Real,C<:Complex{R}} = chop(real(x), x0) + im*chop(imag(x), x0)
 
+function unique_sorted_approx!(v::AbstractVector{T}) where {T}
+    i = 1
+    xprev = first(v)
+    for j in 2:length(v)
+        if v[j] ≈ xprev
+            xprev = v[j]
+        else
+            i += 1
+            xprev = v[i] = v[j]
+        end
+    end
+    resize!(v, i)
+    return v
+end
+
 ############################################################################################
 
 function pushapproxruns!(runs::AbstractVector{<:UnitRange}, list::AbstractVector{T},

test/test_hamiltonian.jl

@@ -1,5 +1,5 @@
 using LinearAlgebra: diag, norm
-using Quantica: Hamiltonian, ParametricHamiltonian
+using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings
 
 @testset "basic hamiltonians" begin
     presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular,
@@ -21,31 +21,37 @@ using Quantica: Hamiltonian, ParametricHamiltonian
     # Inf range
     h = LatticePresets.square() |> unitcell(region = RegionPresets.square(5)) |>
         hamiltonian(hopping(1, range = Inf))
-    @test Quantica.nhoppings(h) == 600
+    @test nhoppings(h) == 600
 
     h = LatticePresets.square() |> hamiltonian(hopping(1, dn = (10,0), range = Inf))
-    @test Quantica.nhoppings(h) == 1
+    @test nhoppings(h) == 1
     @test isassigned(h, (10,0))
 
     h = LatticePresets.honeycomb() |> hamiltonian(onsite(1.0, sublats = :A), orbitals = (Val(1), Val(2)))
     @test Quantica.nonsites(h) == 1
 
     h = LatticePresets.square() |> unitcell(3) |> hamiltonian(hopping(1, indices = (1:8 .=> 2:9, 9=>1), range = 3, plusadjoint = true))
-    @test Quantica.nhoppings(h) == 48
+    @test nhoppings(h) == 48
 
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 1)))
-    @test Quantica.nhoppings(h) == 12
+    @test nhoppings(h) == 12
 
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 2/√3)))
-    @test Quantica.nhoppings(h) == 18
+    @test nhoppings(h) == 18
 
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (2, 1)))
     @test Quantica.nhoppings(h) == 0
+
+    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (30, 30)))
+    @test Quantica.nhoppings(h) == 12
+
+    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (10, 10.1)))
+    @test Quantica.nhoppings(h) == 48
 end
 
 @testset "hamiltonian unitcell" begin
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1/√3)) |> unitcell((1,-1), region = r -> abs(r[2])<2)
-    @test Quantica.nhoppings(h) == 22
+    @test nhoppings(h) == 22
 end
 
 @testset "hamiltonian wrap" begin
@@ -153,7 +159,6 @@ end
     @test !ishermitian(hamiltonian(lat, hopping(im, sublats = :A=>:B)))
     @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B)))
     @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (-1,0))))
-    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (1,0))))
     @test !ishermitian(hamiltonian(lat, hopping(im)))
     @test ishermitian(hamiltonian(lat, hopping(1)))
 
@@ -281,12 +286,12 @@ end
     # Issue #54. Parametric Haldane model
     sK(dr::SVector) = sK(atan(dr[2],dr[1]))
     sK(ϕ) = 2*mod(round(Int, 6*ϕ/(2π)), 2) - 1
-    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |>
+    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1)) |>
          parametric(@hopping!((t, r, dr; λ) ->  λ*im*sK(dr); sublats = :A=>:A),
                     @hopping!((t, r, dr; λ) -> -λ*im*sK(dr); sublats = :B=>:B))
     @test bloch(ph(λ=1), (π/2, -π/2)) == bloch(ph, (1, π/2, -π/2)) ≈ [4 1; 1 -4]
     # Non-numeric parameters
-    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |>
+    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1)) |>
          parametric(@hopping!((t, r, dr; λ, k) ->  λ*im*sK(dr+k); sublats = :A=>:A),
                     @hopping!((t, r, dr; λ, k) -> -λ*im*sK(dr+k); sublats = :B=>:B))
     @test bloch(ph(λ=1, k=SA[1,0]), (π/2, -π/2)) == bloch(ph, (1, SA[1,0], π/2, -π/2)) ≈ [-4 1; 1 4]
@@ -341,4 +346,17 @@ end
 
     @test Quantica.nsites(h) == 130
     @test Quantica.nsites(h3) == 64
+end
+
+@testset "nrange" begin
+    lat = LatticePresets.honeycomb(a0 = 2)
+    @test nrange(1, lat) ≈ 2/√3
+    @test nrange(2, lat) ≈ 2
+    @test nrange(3, lat) ≈ 4/√3
+    @test hamiltonian(lat, hopping(1)) |> nhoppings == 6
+    @test hamiltonian(lat, hopping(1, range = nrange(2))) |> nhoppings == 18
+    @test hamiltonian(lat, hopping(1, range = nrange(3))) |> nhoppings == 24
+    @test hamiltonian(lat, hopping(1, range = (nrange(2), nrange(3)))) |> nhoppings == 18
+    @test hamiltonian(lat, hopping(1, range = (nrange(20), nrange(20)))) |> nhoppings == 18
+    @test hamiltonian(lat, hopping(1, range = (nrange(20), nrange(21)))) |> nhoppings == 30
 end