generalised diagonaliser to non-hermitian matrices, added function to…

… check if state is eigenstate of operator, and added a bandstructure calculation function

src/eigen.jl

@@ -140,14 +140,17 @@ function Spectrum(
     basisStates::Vector{Dict{BitVector,Float64}};
     diagElements::Vector{Float64}=Float64[],
     tolerance::Float64=1e-14,
+    assumeHerm::Bool=true,
     )
     matrix = OperatorMatrix(basisStates, operator; tolerance=tolerance)
-    @assert maximum(abs.(matrix .- matrix')) < tolerance
+    if assumeHerm
+        @assert maximum(abs.(matrix .- matrix')) < tolerance
+    end
     if !isempty(diagElements)
         matrix += diagm(diagElements)
     end
 
-    eigenVals, eigenVecs = eigen(Hermitian(matrix))
+    eigenVals, eigenVecs = assumeHerm ? eigen(Hermitian(matrix)) : eigen(matrix)
     eigenStates = [TransformState(collect(vector), basisStates; tolerance=tolerance)
                    for vector in eachcol(eigenVecs)]
     return eigenVals, eigenStates
@@ -249,7 +252,7 @@ function Spectrum(
     classifiedBasis::Union{Dict{Tuple{Int64}, Vector{Dict{BitVector,Float64}}}, Dict{Tuple{Int64, Int64}, Vector{Dict{BitVector,Float64}}}};
     diagElements::Dict{}=Dict(),
     tolerance::Float64=1e-14,
-)
+    )
     if !isempty(diagElements)
         @assert all(x -> true, [length(E) == size(classifiedBasis[k])[1] for (k,E) in diagElements])
     end
@@ -266,3 +269,68 @@ function Spectrum(
     return classifiedEigvals, classifiedEigvecs
 end
 export Spectrum
+
+
+function BandStructure(
+        hamiltonian::Vector{Tuple{String,Vector{Int64},Float64}},
+        basisStates::Vector{Dict{BitVector,Float64}},
+        eigenOperator::Vector{Tuple{String,Vector{Int64},Float64}};
+        tolerance::Float64=1e-13
+    )
+    @assert DoesCommute(hamiltonian, eigenOperator, basisStates)
+    modifiedHamiltonian = vcat(hamiltonian, eigenOperator)
+    modifiedEigenVals, eigenStates = Spectrum(modifiedHamiltonian, basisStates)
+    diffEigenVals, _ = Spectrum(eigenOperator, basisStates; assumeHerm=true)
+    eigenVals = (modifiedEigenVals .- diffEigenVals) .|> real
+    quantumNos = trunc.(Int, length(quantumNos) .* atan.(imag.(diffEigenVals), real.(diffEigenVals)) ./ π)
+    uniqueQuantumNos = quantumNos |> unique |> sort
+    classifiedEnergies = Dict{Float64, Vector{Float64}}(Qi => Float64[] for Qi in uniqueQuantumNos)
+    for (Ei, Qi) in zip(eigenVals, quantumNos)
+        @assert Qi in keys(classifiedEnergies)
+        push!(classifiedEnergies[Qi], Ei)
+    end
+    return classifiedEnergies
+end
+export BandStructure
+
+
+function IsEigenState(
+        state::Dict{BitVector,Float64},
+        operator::Vector{Tuple{String,Vector{Int64},Float64}};
+        tolerance::Float64=1e-15,
+    )
+    @assert !isempty(state) 
+    pivotVal, pivotKey = findmax(state)
+    if pivotVal == 0
+        pivotVal, pivotKey = findmin(state)
+        @assert pivotVal ≠ 0
+    end
+    newState = ApplyOperator(operator, state)
+    for k in keys(state)
+        println((k, state[k], newState[k]))
+    end
+    if keys(newState) ≠ keys(state)
+        return false, 0.
+    end
+    for k in keys(state)
+        println((k, state[k], newState[k]))
+    end
+    if newState[pivotKey] == 0
+        if all(c -> abs(c) < tolerance, values(newState))
+            println("New state is all zero, hard to tell if eigenstate or simply the vacuum.")
+            return true, 0.
+        else
+            return false, 0.
+        end
+    else
+        eigenValue = newState[pivotKey] / state[pivotKey]
+        println(eigenValue)
+        println([newState[k] / state[k] for k in keys(state)])
+        if all(k -> abs(eigenValue - newState[k] / state[k]) < tolerance, keys(state))
+            return true, eigenValue
+        else
+            return false, 0.
+        end
+    end
+end
+export IsEigenState