force hermiticity in CompressedOrbitalMatrix

src/meanfield.jl

@@ -48,9 +48,9 @@ function meanfield(g::GreenFunction{T,E}, args...;
 
     lat = lattice(hamiltonian(g))
     # The sparse structure of hFock will be inherited by the evaluated mean field. Need onsite.
-    hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : Vf(dr); selector..., includeonsite = true)
+    hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : T(Vf(dr)); selector..., includeonsite = true)
     hHartree = (Uf == U && Vh === Vf) ? hFock :
-        lat |> hopping((r, dr) -> iszero(dr) ? U : Vh(dr); selector..., includeonsite = true)
+        lat |> hopping((r, dr) -> iszero(dr) ? U : T(Vh(dr)); selector..., includeonsite = true)
 
     potHartree = sum(unflat, harmonics(hHartree))
 
@@ -75,10 +75,10 @@ function meanfield(g::GreenFunction{T,E}, args...;
     return MeanField(output, potHartree, potFock, rhoHartree, rhoFock, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp)
 end
 
-sanitize_potential(x::Number) = Returns(x)
+sanitize_potential(x::Real) = Returns(x)
 sanitize_potential(x::Function) = x
 sanitize_potential(x::Nothing) = Returns(0)
-sanitize_potential(_) = argerror("Invalid potential: use a number or a function of position")
+sanitize_potential(_) = argerror("Invalid potential: use a real number or a function of position")
 
 sanitize_nambu_rotated(is_nambu_rotated, nambu, ::Type{<:SMatrix{2,2}}) = false
 sanitize_nambu_rotated(is_nambu_rotated, nambu, ::Type{<:SMatrix{4,4}}) =

src/solvers/green/schur.jl

@@ -83,8 +83,8 @@ function nearest_cell_harmonics(h)
     is_nearest ||
         argerror("Too many or too few harmonics. Perhaps try `supercell` to ensure strictly nearest-cell harmonics.")
     hm, h0, hp = h[hybrid(-1)], h[hybrid(0)], h[hybrid(1)]
-    flat(hm) == flat(hp)' ||
-        argerror("The Hamiltonian should have h[1] == h[-1]' to use the Schur solver")
+    flat(hm) ≈ flat(hp)' ||
+        argerror("The Hamiltonian should have h[1] ≈ h[-1]' to use the Schur solver")
     return hm, h0, hp
 end
 
@@ -696,7 +696,7 @@ function (s::DensityMatrixSchurSolver)(µ, kBT; params...)
     xs = sort!(ϕs ./ (2π))
     pushfirst!(xs, -0.5)
     push!(xs, 0.5)
-    xs = unique!(xs)
+    xs = [mean(view(xs, rng)) for rng in approxruns(xs)]  # elliminate approximate duplicates
     result = call!_output(s.gs)
     integrand!(x) = fermi_h!(result, s, 2π * x, µ, inv(kBT); params...)
     fx! = (y, x) -> (y .= integrand!(x))

src/specialmatrices.jl

@@ -582,24 +582,34 @@ LinearAlgebra.norm(a::AbstractOrbitalMatrix) = norm(parent(a))
 
 # decoding access
 Base.view(a::CompressedOrbitalMatrix{C}, i::AnyCellSites, j::AnyCellSites = i) where {C} =
-    checkbounds(Bool, a, i, j) ? getindex_scalar(a, only.(indexcollection(a, i, j))) : zero(decoder(a)(zero(C)))
+    checkbounds(Bool, a, i, j) ? getindex_decode(a, i, j) : zero(decoder(a)(zero(C)))
 
-@inline function getindex_scalar(a::CompressedOrbitalMatrix{C}, (i,j)::NTuple{2,Integer}) where {C}
-    dec = decoder(a)
+@inline function getindex_decode(a::CompressedOrbitalMatrix, ci, cj)
+    i, j = only.(indexcollection(a, ci, cj))
     # Only lower-triangle is stored in the hermitian case
     if ishermitian(a)
-        if i < j
-            return dec(parent(a)[j, i])'
-        elseif i == j
-            d = dec(parent(a)[i, j])
-            return 0.5 * (d + d')
-        else
-            return dec(parent(a)[i, j])
-        end
+        # we find the indices for a[sites(-ni, j), sites(0, i)],
+        # the conjugate of a[sites(ni, i), sites(0, j)],
+        ci´, cj´ = sites(-cell(ci), siteindex(cj)), sites(cell(cj), siteindex(ci))
+        i´, j´ = only.(indexcollection(a, ci´, cj´))
+        val = getindex_decode_hermitian(a, i, j)
+        val´ = getindex_decode_hermitian(a, i´, j´)'
+        return 0.5 * (val + val´)
     end
     return dec(parent(a)[i, j])
 end
 
+function getindex_decode_hermitian(a::CompressedOrbitalMatrix, i, j)
+    dec = decoder(a)
+    if i < j
+        return dec(parent(a)[j, i])'
+    elseif i == j
+        return dec(parent(a)[i, j])
+    else
+        return dec(parent(a)[i, j])
+    end
+end
+
 ## checkbounds
 
 # checkbounds axis