Remove dangling bonds with mincoordination

src/hamiltonian.jl

@@ -696,6 +696,9 @@ function flatsize(h::Hamiltonian)
     return (n, n)
 end
 
+expand_supercell_mask(h::Hamiltonian{<:Superlattice}) =
+    Hamiltonian(expand_supercell_mask(h.lattice), h.harmonics, h.orbitals)
+
 Base.size(h::Hamiltonian, n) = size(first(h.harmonics).h, n)
 Base.size(h::Hamiltonian) = size(first(h.harmonics).h)
 Base.size(h::HamiltonianHarmonic, n) = size(h.h, n)
@@ -1246,51 +1249,105 @@ function supercell(ham::Hamiltonian, args...; kw...)
     return Hamiltonian(slat, ham.harmonics, ham.orbitals)
 end
 
-function unitcell(ham::Hamiltonian{<:Lattice}, args...; modifiers = (), kw...)
+function unitcell(ham::Hamiltonian{<:Lattice}, args...; modifiers = (), mincoordination = missing, kw...)
     sham = supercell(ham, args...; kw...)
-    return unitcell(sham; modifiers = modifiers)
+    return unitcell(sham; modifiers = modifiers, mincoordination = mincoordination)
 end
 
-function unitcell(ham::Hamiltonian{LA,L}; modifiers = ()) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}}
+function unitcell(ham::Hamiltonian{LA,L}; modifiers = (), mincoordination = missing) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}}
     slat = ham.lattice
     sc = slat.supercell
     supercell_dn = r_to_dn(slat, sc.matrix, SVector{L´}(1:L´))
     pos = allsitepositions(slat)
     br = bravais(slat)
     modifiers´ = resolve.(ensuretuple(modifiers), Ref(slat))
-    mapping = OffsetArray{Int}(undef, sc.sites, sc.cells.indices...) # store supersite indices newi
+
+    # Build a version of slat that has a filtered supercell mask according to mincoordination
+    slat´ = filtered_superlat!(ham, mincoordination, supercell_dn, br, sc.matrix, pos)
+    # store supersite indices newi
+    mapping = OffsetArray{Int}(undef, sc.sites, sc.cells.indices...)
     mapping .= 0
-    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, slat)
-    dim = nsites(sc)
+    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, slat´)
+
+    dim = nsites(slat´.supercell)
     B = blocktype(ham)
     S = typeof(SparseMatrixBuilder{B}(dim, dim))
     harmonic_builders = HamiltonianHarmonic{L´,B,S}[]
-    foreach_supersite(slat) do s, source_i, source_dn, newcol
+
+    foreach_supersite(slat´) do s, source_i, source_dn, newcol
         for oldh in ham.harmonics
             rows = rowvals(oldh.h)
             vals = nonzeros(oldh.h)
             target_dn = source_dn + oldh.dn
             for p in nzrange(oldh.h, source_i)
                 target_i = rows[p]
-                r = pos[target_i] + br * target_dn
-                super_dn = supercell_dn(r)
-                wrapped_dn = wrap_dn(target_dn, super_dn, sc.matrix)
+                wrapped_dn, super_dn = wrap_super_dn(target_i, target_dn, supercell_dn, br, sc.matrix, pos)
                 # check: wrapped_dn could exit bounding box along non-periodic direction
                 checkbounds(Bool, mapping, target_i, Tuple(wrapped_dn)...) || continue
                 newh = get_or_push!(harmonic_builders, super_dn, dim, newcol)
                 newrow = mapping[target_i, Tuple(wrapped_dn)...]
-                val = applymodifiers(vals[p], slat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
-                iszero(newrow) || pushtocolumn!(newh.h, newrow, val)
+                if !iszero(newrow)
+                    val = applymodifiers(vals[p], slat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
+                    pushtocolumn!(newh.h, newrow, val)
+                end
             end
         end
         foreach(h -> finalizecolumn!(h.h), harmonic_builders)
     end
     harmonics = [HamiltonianHarmonic(h.dn, sparse(h.h)) for h in harmonic_builders]
-    unitlat = unitcell(slat)
+    unitlat = unitcell(slat´)
     orbs = ham.orbitals
     return Hamiltonian(unitlat, harmonics, orbs)
 end
 
+filtered_superlat!(sham, ::Missing, args...) = sham.lattice
+filtered_superlat!(sham, mc::Int, args...) =
+    _filtered_superlat!(expand_supercell_mask(sham), mc, args...)
+
+function _filtered_superlat!(sham::Hamiltonian{LA,L}, mincoord, args...) where {LA,L}
+    slat = sham.lattice
+    sc = slat.supercell
+    mask = sc.mask
+    if mincoord > 0
+        delsites = NTuple{L+1,Int}[]
+        while true
+            foreach_supersite(sham.lattice) do _, source_i, source_dn, _
+                nn = num_neighbors_supercell(sham.harmonics, source_i, source_dn, mask, args...)
+                nn < mincoord && push!(delsites, (source_i, Tuple(source_dn)...))
+            end
+            foreach(p -> mask[p...] = false, delsites)
+            isempty(delsites) && break
+            resize!(delsites, 0)
+        end
+    end
+    sc = Supercell(sc.matrix, sc.sites, sc.cells, mask)
+    return Superlattice(slat.bravais, slat.unitcell, sc)
+end
+
+function num_neighbors_supercell(hhs, source_i, source_dn, mask, args...)
+    nn = 0
+    for hh in hhs
+        ptrs = nzrange(hh.h, source_i)
+        rows = rowvals(hh.h)
+        target_dn = source_dn + hh.dn
+        for p in nzrange(hh.h, source_i)
+            target_i = rows[p]
+            wrapped_dn, _ = wrap_super_dn(target_i, target_dn, args...)
+            isonsite = rows[p] == source_i && iszero(hh.dn)
+            isincell = isinmask(mask, rows[p], Tuple(wrapped_dn)...)
+            nn += !isonsite && isincell
+        end
+    end
+    return nn
+end
+
+function wrap_super_dn(target_i, target_dn, supercell_dn, br, smat, pos)
+    r = pos[target_i] + br * target_dn
+    super_dn = supercell_dn(r)
+    wrapped_dn = target_dn - smat * super_dn
+    return wrapped_dn, super_dn
+end
+
 function get_or_push!(hs::Vector{<:HamiltonianHarmonic{L,B,<:SparseMatrixBuilder}}, dn, dim, currentcol) where {L,B}
     for h in hs
         h.dn == dn && return h
@@ -1301,8 +1358,6 @@ function get_or_push!(hs::Vector{<:HamiltonianHarmonic{L,B,<:SparseMatrixBuilder
     return newh
 end
 
-wrap_dn(olddn::SVector, newdn::SVector, supercell::SMatrix) = olddn - supercell * newdn
-
 applymodifiers(val, lat, inds, dns) = val
 
 function applymodifiers(val, lat, inds, dns, m::UniformModifier, ms...)

src/lattice.jl

@@ -360,6 +360,10 @@ function boolean_mask(f, mask1, mask2, indranges)
     return mask
 end
 
+expand_supercell_mask(s::Supercell{L,L´,Missing}) where {L,L´} =
+    Supercell(s.matrix, s.sites, s.cells, ones(Bool, s.sites, s.cells.indices...))
+expand_supercell_mask(s::Supercell{L,L´}) where {L,L´} = s
+
 #######################################################################
 # Superlattice
 #######################################################################
@@ -414,6 +418,8 @@ function check_compatible_superlattice(s1, s2)
     return nothing
 end
 
+expand_supercell_mask(s::Superlattice) = Superlattice(s.bravais, s.unitcell, expand_supercell_mask(s.supercell))
+
 #######################################################################
 # AbstractLattice interface
 #######################################################################
@@ -469,6 +475,9 @@ ismasked(lat::Superlattice) = ismasked(lat.supercell)
 maskranges(lat::Superlattice) = (1:nsites(lat), lat.supercell.cells.indices...)
 maskranges(lat::Lattice) = (1:nsites(lat),)
 
+expand_supercell_mask(s::Superlattice) =
+    Superlattice(s.bravais, s.unitcell, expand_supercell_mask(s.supercell))
+
 """
     transform!(f::Function, lat::Lattice)
 
@@ -633,15 +642,23 @@ end
 
 # Computes δn[inds] so that r = bravais´ * δn + dr, where dr is within a supercell and
 # bravais´ = bravais * supercell, but extending supercell into a square matrix
-# Supercell center is placed at mean(allpos)
+# Supercell center is placed at mean(allpos). Enconde it as a struct for type stability
+struct PosToCell{S,V,I}
+    projector::S
+    r0::V
+    inds::I
+end
+
+(p::PosToCell)(r) = floor.(Int, (p.projector * (r - p.r0))[p.inds])
+
 function r_to_dn(lat::AbstractLattice{E,L,T}, sc::SMatrix{L}, inds = :) where {E,L,T}
     br = bravais(lat)
     extsc = extended_supercell(br, sc)
     projector = pinverse(br * extsc) # E need not be equal to L, hence pseudoinverse
     # Place mean(positions) at the center of supercell
     r0 = supercell_center(lat)
     # This results in a zero vector for all sites within the unit supercell
-    return r -> floor.(Int, (projector * (r - r0))[inds])
+    return PosToCell(projector, r0, inds)
 end
 
 supercell_center(lat::AbstractLattice{E,L,T}) where {E,L,T} =
@@ -774,20 +791,22 @@ with factors along the diagonal)
 
 Convert Superlattice `slat` into a lattice with its unit cell matching `slat`'s supercell.
 
-    unitcell(h::Hamiltonian, v...; modifiers = (), kw...)
+    unitcell(h::Hamiltonian, v...; mincoordination, modifiers = (), kw...)
 
 Transforms the `Lattice` of `h` to have a larger unitcell, while expanding the Hamiltonian
-accordingly. The modifiers (a tuple of `ElementModifier`s, either `@onsite!` or `@hopping!`
-with no free parameters) will be applied to onsite and hoppings as the hamiltonian is
-expanded. See `@onsite!` and `@hopping!` for details
+accordingly. A nonzero `mincoordination` indicates a minimum number of hopping neighbors
+required for sites to be included in the resulting unit cell. The modifiers (a tuple of
+`ElementModifier`s, either `@onsite!` or `@hopping!` with no free parameters) will be
+applied to onsite and hoppings as the hamiltonian is expanded. See `@onsite!` and
+`@hopping!` for details
 
 Note: for performance reasons, in sparse hamiltonians only the stored onsites and hoppings
 will be transformed by `ElementModifier`s, so you might want to add zero onsites or hoppings
 when building `h` to have a modifier applied to them later.
 
     lat_or_h |> unitcell(v...; kw...)
 
-Curried syntax, equivalent to `unitcell(lat_or_h, v...; kw...)
+Curried syntax, equivalent to `unitcell(lat_or_h, v...; kw...)`
 
 # Examples
 

test/test_hamiltonian.jl

@@ -72,6 +72,15 @@ end
         h´´ = unitcell(h´, 2)
         @test coordination(h´´) ≈ coordination(h´)
     end
+    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2) |> unitcell(mincoordination = 2)
+    @test nsites(h) == 6
+    h = LP.cubic() |> hamiltonian(hopping(1)) |> unitcell(4) |> unitcell(mincoordination = 4)
+    @test nsites(h) == 0
+    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(region = RP.circle(5) & !RP.circle(2)) |> unitcell(mincoordination = 2)
+    @test nsites(h) == 144
+    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(10, region = !RP.circle(2, (0,8)))
+    h´ = h |> unitcell(1, mincoordination = 2)
+    @test nsites(h´) == nsites(h) - 1
 end
 
 @testset "hamiltonian wrap" begin