From 83ea2964947980449ce8f73c28c549f7639eadcb Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Tue, 5 Jan 2021 15:34:46 +0100
Subject: [PATCH] enforce harmonics invariants

---
 src/hamiltonian.jl       | 23 ++++++++++++++---------
 src/tools.jl             |  2 --
 test/test_hamiltonian.jl |  3 +++
 3 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
index 76e887b2..6f6f9f46 100644
--- a/src/hamiltonian.jl
+++ b/src/hamiltonian.jl
@@ -32,19 +32,24 @@ struct Hamiltonian{LA<:AbstractLattice,L,M,A<:AbstractMatrix,
     lattice::LA
     harmonics::Vector{H}
     orbstruct::O
+    # Enforce sorted-dns-starting-from-zero invariant onto harmonics
+    function Hamiltonian{LA,L,M,A,H,O}(lattice, harmonics, orbstruct) where
+        {LA<:AbstractLattice,L,M,A<:AbstractMatrix, H<:HamiltonianHarmonic{L,M,A},O<:OrbitalStructure}
+        dimh = nsites(lattice)
+        all(har -> size(har.h) == (dimh, dimh), harmonics) || throw(DimensionMismatch("Harmonics don't match lattice dimensions"))
+        length(harmonics) > 0 && iszero(first(harmonics).dn) ||
+            push!(harmonics, H(zero(SVector{L,Int}), dimh, dimh))
+        sort!(harmonics, by = h -> abs.(h.dn))
+        return new(lattice, harmonics, orbstruct)
+    end
 end
 
+Hamiltonian(lat::LA, hs::Vector{H}, orb::O) where {LA<:AbstractLattice,L,M,A<:AbstractMatrix, H<:HamiltonianHarmonic{L,M,A},O<:OrbitalStructure} =
+    Hamiltonian{LA,L,M,A,H,O}(lat, hs, orb)
+
 Hamiltonian(lat, hs::Vector{H}, orbs::Tuple) where {L,M,H<:HamiltonianHarmonic{L,M}} =
     Hamiltonian(lat, hs, OrbitalStructure(orbitaltype(M), orbs, lat))
 
-function Hamiltonian(lat, hs::Vector{H}, orbs, n::Int, m::Int) where {L,M,H<:HamiltonianHarmonic{L,M}}
-    sort!(hs, by = h -> abs.(h.dn))
-    if isempty(hs) || !iszero(first(hs).dn)
-        pushfirst!(hs, H(zero(SVector{L,Int}), empty_sparse(M, n, m)))
-    end
-    return Hamiltonian(lat, hs, orbs)
-end
-
 orbitals(h::Hamiltonian) = h.orbstruct.orbitals
 
 offsets(h::Hamiltonian) = h.orbstruct.offsets
@@ -1175,7 +1180,7 @@ function hamiltonian_sparse!(builder::IJVBuilder{L,M}, lat::AbstractLattice{E,L}
     n = nsites(lat)
     HT = HamiltonianHarmonic{L,M,SparseMatrixCSC{M,Int}}
     harmonics = HT[HT(e.dn, sparse(e.i, e.j, e.v, n, n)) for e in builder.ijvs if !isempty(e)]
-    return Hamiltonian(lat, harmonics, orbs, n, n)
+    return Hamiltonian(lat, harmonics, orbs)
 end
 
 applyterms!(builder, terms...) = foreach(term -> applyterm!(builder, term), terms)
diff --git a/src/tools.jl b/src/tools.jl
index a191e575..9d0824ea 100644
--- a/src/tools.jl
+++ b/src/tools.jl
@@ -127,8 +127,6 @@ padright(v, x::T, ::Type{S}) where {E,T,S<:SVector{E,T}} =
 negative(s::SVector{L,<:Number}) where {L} = -s
 negative(s::SVector{0,<:Number}) = s
 
-empty_sparse(::Type{M}, n, m) where {M} = sparse(Int[], Int[], M[], n, m)
-
 display_as_tuple(v, prefix = "") = isempty(v) ? "()" :
     string("(", prefix, join(v, string(", ", prefix)), ")")
 
diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl
index 09481ee6..81f7df00 100644
--- a/test/test_hamiltonian.jl
+++ b/test/test_hamiltonian.jl
@@ -83,6 +83,9 @@ end
     @test nsites(h´) == nsites(h) - 1
     h = LP.square() |> hamiltonian(hopping(0)) |> unitcell(4, mincoordination = 2)
     @test nsites(h) == 0
+    # check dn=0 invariant
+    h = LP.linear() |> hamiltonian(hopping(1)) |> unitcell((1,))
+    @test length(h.harmonics) == 3 && iszero(first(h.harmonics).dn)
 end
 
 @testset "hamiltonian wrap" begin