Skip to content

Commit

Permalink
unitcell stress test
Browse files Browse the repository at this point in the history 
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

fix test
  • Loading branch information
@pablosanjose
pablosanjose committed Oct 16, 2020
1 parent 8758326 commit a66c2b8
Show file tree
Hide file tree
Showing 5 changed files with 101 additions and 76 deletions.
1 change: 1 addition & 0 deletions Project.toml
View file
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
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
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
102 changes: 62 additions & 40 deletions src/lattice.jl
View file
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,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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
#######################################################################
Expand Down
35 changes: 17 additions & 18 deletions test/test_hamiltonian.jl
View file
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

0 comments on commit a66c2b8

Please sign in to comment.