unitcell stress test

remove obsolete tests removed unitcell fastpaths for the moment shortcut remove more obsolete fastpath tests test fix more edge cases catch all outliers remove test
pablosanjose · Oct 16, 2020 · ad31ce5 · ad31ce5
1 parent 8758326
commit ad31ce5
Show file tree

Hide file tree

Showing 5 changed files with 100 additions and 67 deletions.
diff --git a/Project.toml b/Project.toml
@@ -17,6 +17,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DualNumbers = "0.6"

diff --git a/src/Quantica.jl b/src/Quantica.jl
@@ -17,6 +17,8 @@ using ExprTools
 
 using SparseArrays: getcolptr, AbstractSparseMatrix
 
+using Statistics: mean
+
 export sublat, bravais, lattice, dims, supercell, unitcell,
        hopping, onsite, @onsite!, @hopping!, parameters, siteselector, hopselector, nrange,
        sitepositions, siteindices, not,

diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
@@ -38,7 +38,7 @@ $i  Orbitals         : $(displayorbitals(ham))
 $i  Element type     : $(displayelements(ham))
 $i  Onsites          : $(nonsites(ham))
 $i  Hoppings         : $(nhoppings(ham))
-$i  Coordination     : $(nhoppings(ham) / nsites(ham))")
+$i  Coordination     : $(coordination(ham))")
     ioindent = IOContext(io, :indent => string("  "))
     issuperlattice(ham.lattice) && print(ioindent, "\n", ham.lattice.supercell)
 end
@@ -625,15 +625,17 @@ _orbitaltype(t::Type{SVector{1,Tv}}) where {Tv} = Tv
 orbitaltype(h::Hamiltonian{LA,L,M}) where {N,T,LA,L,M<:SMatrix{N,N,T}} = SVector{N,T}
 orbitaltype(h::Hamiltonian{LA,L,M}) where {LA,L,M<:Number} = M
 
-function nhoppings(ham::Hamiltonian)
+coordination(ham) = nhoppings(ham) / nsites(ham)
+
+function nhoppings(ham)
     count = 0
     for h in ham.harmonics
         count += iszero(h.dn) ? (_nnz(h.h) - _nnzdiag(h.h)) : _nnz(h.h)
     end
     return count
 end
 
-function nonsites(ham::Hamiltonian)
+function nonsites(ham)
     count = 0
     for h in ham.harmonics
         iszero(h.dn) && (count += _nnzdiag(h.h))
@@ -1253,39 +1255,41 @@ function unitcell(ham::Hamiltonian{<:Lattice}, args...; modifiers = (), kw...)
 end
 
 function unitcell(ham::Hamiltonian{LA,L}; modifiers = ()) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}}
-    lat = ham.lattice
-    sc = lat.supercell
-    isc = inv_supercell(bravais(lat), sc.matrix)
-    modifiers´ = resolve.(ensuretuple(modifiers), Ref(lat))
+    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
     mapping .= 0
-    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, lat)
+    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, slat)
     dim = nsites(sc)
     B = blocktype(ham)
     S = typeof(SparseMatrixBuilder{B}(dim, dim))
     harmonic_builders = HamiltonianHarmonic{L´,B,S}[]
-    # pinvint = pinvmultiple(sc.matrix)
-    foreach_supersite(lat) 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
-            super_dn = new_dn(target_dn, isc)
-            wrapped_dn = wrap_dn(target_dn, super_dn, sc.matrix)
-            newh = get_or_push!(harmonic_builders, super_dn, dim, newcol)
             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)
                 # 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], lat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
+                val = applymodifiers(vals[p], slat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
                 iszero(newrow) || pushtocolumn!(newh.h, newrow, val)
             end
         end
         foreach(h -> finalizecolumn!(h.h), harmonic_builders)
     end
     harmonics = [HamiltonianHarmonic(h.dn, sparse(h.h)) for h in harmonic_builders]
-    unitlat = unitcell(lat)
+    unitlat = unitcell(slat)
     orbs = ham.orbitals
     return Hamiltonian(unitlat, harmonics, orbs)
 end
@@ -1300,9 +1304,6 @@ function get_or_push!(hs::Vector{<:HamiltonianHarmonic{L,B,<:SparseMatrixBuilder
     return newh
 end
 
-inv_supercell(br, sc::SMatrix{L,L´}) where {L,L´} = inv(extended_supercell(br, sc))[SVector{L´}(1:L´), :]
-new_dn(oldn, isc) = floor.(Int, isc * oldn)
-
 wrap_dn(olddn::SVector, newdn::SVector, supercell::SMatrix) = olddn - supercell * newdn
 
 applymodifiers(val, lat, inds, dns) = val

diff --git a/src/lattice.jl b/src/lattice.jl
@@ -613,13 +613,6 @@ function supercell(lat::Lattice{E,L}, v...; seed = missing, kw...) where {E,L}
     return _superlat(lat, scmatrix, pararegion, perpselector, seed)
 end
 
-# Fast-path methods
-supercell(lat::Lattice{E,L}, factors::Vararg{Integer,L}) where {E,L} =
-    _superlat_fastpath(lat, factors)
-supercell(lat::Lattice{E,L}, factor::Integer) where {E,L} =
-    _superlat_fastpath(lat, filltuple(factor, Val(L)))
-supercell(lat::Lattice) = _superlat_fastpath(lat, ())
-
 sanitize_supercell(::Val{L}) where {L} = SMatrix{L,0,Int}()
 sanitize_supercell(::Val{L}, ::Tuple{}) where {L} = SMatrix{L,0,Int}()
 sanitize_supercell(::Val{L}, v::NTuple{L,Int}...) where {L} = toSMatrix(Int, v)
@@ -629,18 +622,33 @@ sanitize_supercell(::Val{L}, ss::Integer...) where {L} = SMatrix{L,L,Int}(Diagon
 sanitize_supercell(::Val{L}, v) where {L} =
     throw(ArgumentError("Improper supercell specification $v for an $L lattice dimensions, see `supercell`"))
 
-function supercell_regions(lat::Lattice{E,L,T}, sc::SMatrix{L,L´}) where {E,L,T,L´}
+function supercell_regions(lat::Lattice{E,L}, sc::SMatrix{L,L´}) where {E,L,L´}
+    dn_func = r_to_dn(lat, sc)
+    parainds = SVector{L´,Int}(1:L´)
+    perpinds = SVector{L-L´,Int}((1+L´):L)
+    pararegion(r) = iszero(dn_func(r)[parainds])
+    perpregion(r) = iszero(dn_func(r)[perpinds])
+    return pararegion, perpregion
+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)
+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
-    siteshift = ntuple(Val(L)) do i
-        minimum((projector * r)[i] for r in allsitepositions(lat)) - sqrt(eps(T))
-    end
-    pararegion(r) = all(i -> 0 <= (projector * r)[i] - siteshift[i] < 1, 1:L´)
-    perpregion(r) = all(i -> 0 <= (projector * r)[i] - siteshift[i] < 1, (1+L´):L)
-    return pararegion, perpregion
+    # 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])
 end
 
+supercell_center(lat::AbstractLattice{E,L,T}) where {E,L,T} =
+    mean(allsitepositions(lat)) -
+    bravais(lat) * SVector{L,T}(filltuple(1/2, Val(L))) -
+    SVector{E,T}(filltuple(sqrt(eps(T)), Val(E)))
+
 # supplements supercell with most orthogonal bravais axes
 function extended_supercell(bravais, supercell::SMatrix{L,L´}) where {L,L´}
     L == L´ && return supercell
@@ -657,7 +665,8 @@ end
 function _superlat(lat, scmatrix, pararegion, selector_perp, seed)
     br = bravais(lat)
     rsel = resolve(selector_perp, lat)
-    cells = _cell_iter(lat, pararegion, rsel, seed)
+    cells = _cell_iter(lat, scmatrix, pararegion, rsel, seed)
+    # cells = _cell_iter(lat, scmatrix)
     ns = nsites(lat)
     mask = OffsetArray(falses(ns, size(cells)...), 1:ns, cells.indices...)
     @inbounds for dn in cells
@@ -673,26 +682,51 @@ function _superlat(lat, scmatrix, pararegion, selector_perp, seed)
     return Superlattice(lat.bravais, lat.unitcell, supercell)
 end
 
-function _cell_iter(lat::Lattice{E,L}, pararegion, rsel_perp, seed) where {E,L}
+function _cell_iter(lat::Lattice{E,L}, sc::SMatrix{L,L´}, pararegion, rsel_perp, seed) where {E,L,L´}
+    # We first ensure that a full supercell bounding box is enconpassed by the iterator,
+    # the sought after cell iter is bbox only if L == L´ and all sites are inside unitcell (no outliers)
+    br = bravais(lat)
+    extsc = extended_supercell(br, sc)
+    dns = Iterators.product(ntuple(_ -> 0:1, Val(L))...)
+    bbox_min = bbox_max = zero(SVector{L,Int})
+    for dn in dns
+        isempty(dn) && break
+        bbox_min = min.(extsc * SVector(dn), bbox_min)
+        bbox_max = max.(extsc * SVector(dn), bbox_max)
+    end
+    minimum_bbox = CartesianIndices(UnitRange.(Tuple(bbox_min), Tuple(bbox_max)))
+
+    # We now iterate over a growing box of dn, ensuring that all dn are included such that
+    # r + br*dn is inside the supercell for any unitcell site at r
     seed´ = seed === missing ? zero(SVector{L,Int}) : seedcell(SVector{E}(seed), bravais(lat))
     iter = BoxIterator(seed´)
-    foundfirst = false
     counter = 0
     br = bravais(lat)
-    for dnvec in iter
+    ibr = pinverse(br)
+    for dn in iter
         found = false
         counter += 1; counter == TOOMANYITERS &&
             throw(ArgumentError("`region` seems unbounded (after $TOOMANYITERS iterations)"))
-        foundfirst || acceptcell!(iter, dnvec)
-        for i in siteindices(rsel_perp, dnvec)
-            r = siteposition(i, dnvec, lat)
-            # site i is already in perpregion through rsel. Is it also in pararegion?
-            found_in_cell = pararegion(r)
-            if found_in_cell
-                acceptcell!(iter, dnvec)
-                foundfirst = true
-                break
-            end
+        # we need to make sure we've covered at least the minimum bounding box
+        inside_minimum = CartesianIndex(Tuple(dn)) in minimum_bbox
+        inside_minimum && (acceptcell!(iter, dn); continue)
+        explored_bbox = CartesianIndices(iter)
+        # We explore all sites in the unit cell, not only `i in siteindices(rsel_perp, dn)`,
+        # because that could cause unitcell outliers to not be found
+        for (i, r) in enumerate(allsitepositions(lat))
+            r_dn = r + br * dn
+            # we now check if the dn_r of r folded onto unitcell has been explored
+            # (i.e. ensure unitcell outliers are reached)
+            dn_r = CartesianIndex(floor.(Int, Tuple(ibr * r)))
+            not_yet_found = !in(dn_r, explored_bbox)
+            # it the site has not been found this dn should be accepted. Continue to next dn
+            not_yet_found && (acceptcell!(iter, dn); break)
+            # is site i in perpregion? Otherwise go to next site
+            (i, dn) in rsel_perp || continue
+            # site i, shifted by dn, is already in perpregion through rsel. Is it also in pararegion?
+            is_in_cell = pararegion(r_dn)
+            # if it is, mark dn as accepted and continue to grow BoxIterator
+            is_in_cell && (acceptcell!(iter, dn); break)
         end
     end
     c = CartesianIndices(iter)
@@ -702,18 +736,6 @@ end
 seedcell(seed::NTuple{N,Any}, brmat) where {N} = seedcell(SVector{N}(seed), brmat)
 seedcell(seed::SVector{E}, brmat::SMatrix{E}) where {E} = round.(Int, brmat \ seed)
 
-function _superlat_fastpath(lat::Lattice{E,L}, factors) where {E,L}
-    scmatrix = sanitize_supercell(Val(L), factors...)
-    sites = 1:nsites(lat)
-    cells = _cells_fastpath(Val(L), factors)
-    mask = missing
-    supercell = Supercell(scmatrix, sites, cells, mask)
-    return Superlattice(lat.bravais, lat.unitcell, supercell)
-end
-
-_cells_fastpath(::Val, factors) = CartesianIndices((i -> 0 : i - 1).(factors))
-_cells_fastpath(::Val{L}, factors::Tuple{}) where {L} = CartesianIndices(filltuple(0:0, Val(L)))
-
 #######################################################################
 # unitcell
 #######################################################################

diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl
@@ -1,5 +1,5 @@
 using LinearAlgebra: diag, norm
-using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings, nonsites, nsites, allsitepositions
+using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings, nonsites, nsites, coordination, allsitepositions
 
 @testset "basic hamiltonians" begin
     presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular,
@@ -60,19 +60,26 @@ end
     h = LatticePresets.square() |> hamiltonian(hopping(1, range = √2)) |> unitcell(5) |> unitcell(indices = 2:2:25)
     @test nhoppings(h) == 32
     @test nsites(h) == 12
-    lat = LatticePresets.honeycomb()
-    h = lat |> hamiltonian(hopping(1)) |> unitcell
-    @test allsitepositions(h.lattice) == allsitepositions(lat)
-    h = lat |> hamiltonian(hopping(1)) |> unitcell((1,1))
-    @test allsitepositions(h.lattice) == allsitepositions(lat)
-    h = lat |> hamiltonian(hopping(1)) |> unitcell((1,-1))
-    @test allsitepositions(h.lattice) == allsitepositions(lat)
     lat = LatticePresets.honeycomb(dim = Val(3)) |> unitcell(3) |> unitcell((1,1), indices = not(1))
     h = lat |> hamiltonian(hopping(1)) |> unitcell
-    @test allsitepositions(h.lattice) == allsitepositions(lat)
     @test nsites(h) == 17
     h = unitcell(h, indices = not(1))
     @test nsites(h) == 16
+    # dim-preserving unitcell reshaping should always preserve coordination
+    lat = lattice(sublat((0.0, -0.1), (0,0), (0,1/√3), (0.0,1.6/√3)), bravais = SA[cos(pi/3) sin(pi/3); -cos(pi/3) sin(pi/3)]')
+    h = lat |> hamiltonian(hopping(1, range = 1/√3))
+    c = coordination(h)
+    iter = CartesianIndices((-3:3, -3:3))
+    for Ic´ in iter, Ic in iter
+        sc = SMatrix{2,2,Int}(Tuple(Ic)..., Tuple(Ic´)...)
+        iszero(det(sc)) && continue
+        h´ = unitcell(h, sc)
+        h´´ = unitcell(h´, 2)
+        @test coordination(h´´) ≈ coordination(h´) ≈ c
+        h´ = unitcell(h, Tuple(Ic))
+        h´´ = unitcell(h´, 2)
+        @test coordination(h´´) ≈ coordination(h´)
+    end
 end
 
 @testset "hamiltonian wrap" begin