Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unitcell stress test: coordination should not change with full-rank supercells #117

Merged
merged 2 commits into from
Oct 16, 2020
+108 −78
Merged

Unitcell stress test: coordination should not change with full-rank supercells #117

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a66c2b8
unitcell stress test
pablosanjose Oct 16, 2020
01c6489
need to normalize in extended_supercell
pablosanjose Oct 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
2 changes: 2 additions & 0 deletions src/Quantica.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
37 changes: 19 additions & 18 deletions src/hamiltonian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
111 changes: 69 additions & 42 deletions src/lattice.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -629,35 +622,56 @@ 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
bravais_new = bravais * supercell
bravais_new_norm = normalize_cols(bravais * supercell)
bravais_norm = normalize_cols(bravais)
# νnorm are the L projections of old bravais on new bravais axis subspace
ν = bravais' * bravais_new # L×L´
ν = bravais_norm' * bravais_new_norm # L×L´
νnorm = SVector(ntuple(row -> norm(ν[row,:]), Val(L)))
νorder = sortperm(νnorm)
ext_axes = hcat(ntuple(i -> unitvector(SVector{L,Int}, νorder[i]), Val(L-L´))...)
ext_supercell = hcat(supercell, ext_axes)
return ext_supercell
end

normalize_cols(s::SMatrix{L,0}) where {L} = s
normalize_cols(s::SMatrix{0,L}) where {L} = s
normalize_cols(s) = mapslices(v -> v/norm(v), s, dims = 1)

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
Expand All @@ -673,26 +687,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)
Expand All @@ -702,18 +741,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
#######################################################################
Expand Down
35 changes: 17 additions & 18 deletions test/test_hamiltonian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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
Expand All @@ -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
Expand Down