diff --git a/Project.toml b/Project.toml index 4576c55a..c95b857a 100644 --- 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 index 43490218..9a536b6f 100644 --- 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 index bfd91647..a9e0a5d3 100644 --- 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,7 +625,9 @@ _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) @@ -633,7 +635,7 @@ function nhoppings(ham::Hamiltonian) 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 index 71029640..3a0526d9 100644 --- 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 index 40d8ba3e..353e960d 100644 --- 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 LinearAlgebra: diag, norm, det +using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings, nonsites, nsites, coordination, allsitepositions @testset "basic hamiltonians" begin presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular, @@ -22,29 +22,21 @@ using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings, nonsites, nsites, h = LatticePresets.square() |> unitcell(region = RegionPresets.square(5)) |> hamiltonian(hopping(1, range = Inf)) @test nhoppings(h) == 600 - h = LatticePresets.square() |> hamiltonian(hopping(1, dn = (10,0), range = Inf)) @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 nhoppings(h) == 48 - h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 1))) @test nhoppings(h) == 12 - h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 2/√3))) @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 @@ -60,19 +52,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