completed tests and added docs for SpecFunc

examples/KondoMERG.jl

@@ -1,5 +1,27 @@
-using fermions
+using fermions, CairoMakie, Measures
+include("../src/iterDiag.jl")
 
+set_theme!(merge(theme_light(), theme_latexfonts()))
+update_theme!(
+              figure_padding = 0,
+              fontsize=28,
+              ScatterLines = (
+                       linewidth = 3,
+                       markersize=10,
+                      ),
+              Lines = (
+                       linewidth = 3,
+                       markersize=20,
+                      ),
+              Scatter = (
+                       markersize=10,
+                      ),
+              Legend = (
+                        patchsize=(50,20),
+                        halign = :right,
+                        valign = :top,
+                       ),
+             )
 function RGFlow(
         energies::Vector{Float64},
         kondoJ::Float64,
@@ -30,19 +52,25 @@ dofStates = (E,D) -> (1/D)*(1 - E^2/D^2)^0.5
 energies = collect(D:-deltaD:0)
 kondoJFlow, fixedPointStep = RGFlow(energies, kondoJ, dofStates)
 
-totalSites = 11
-initStates = 5
-kondoJMatrix = zeros(totalSites, totalSites)
-kondoJMatrix[1:initStates,1:initStates] .= kondoJFlow[end]
-for i in 1:(totalSites - initStates)
-    kondoJMatrix[initStates+i, 1:initStates+i] .= kondoJFlow[end - i]
-    kondoJMatrix[1:initStates+i, initStates+i] .= kondoJFlow[end - i]
+mergResults = Dict("SF-d0" => Float64[])
+initStates = 1
+for totalSites in [3, 21]
+    kondoJMatrix = zeros(1 + 2 * totalSites, 1 + 2 * totalSites)
+    kondoJMatrix[1:1 + 2 * initStates, 1:1 + 2 * initStates] .= kondoJFlow[end]
+    for i in 1:totalSites - initStates
+        kondoJMatrix[-1 + 2 * initStates + 2 * i:1 + 2 * initStates + 2 * i, 1:1 + 2 * initStates + 2 * i] .= kondoJFlow[end - i]
+        kondoJMatrix[1:1 + 2 * initStates + 2 * i, -1 + 2 * initStates + 2 * i:1 + 2 * initStates + 2 * i] .= kondoJFlow[end - i]
+    end
+    println(kondoJMatrix)
+    println([[energies[end]]; repeat(energies[end-1:-1:end-totalSites], inner=2)])
+    hamltonian = KondoModel([[energies[end]]; repeat(energies[end-1:end-totalSites], inner=2)], kondoJMatrix)
+    hamiltonianFlow = MinceHamiltonian(hamltonian, collect((4 + 2 * initStates):2:(4 + 2 * totalSites)))
+    savePaths, resultsDict = IterDiag(hamiltonianFlow, 500;
+                         symmetries=Char['N', 'S'],
+                         correlationDefDict=Dict("SF-d0" => [("+-+-", [1,2,4,3], 1.), ("+-+-", [2,1,3,4], 1.)]),
+                         #=calculateThroughout=true,=#
+                        ) 
+    push!(mergResults["SF-d0"], resultsDict["SF-d0"])
 end
-hamltonian = KondoModel(energies[end:-1:end-totalSites+1], kondoJMatrix)
-hamiltonianFlow = MinceHamiltonian(hamltonian, collect(2 + 2 * initStates:4:2 + 2 * totalSites))
-savePaths, resultsDict = IterDiag(hamiltonianFlow, 500;
-                     symmetries=Char['N', 'S'],
-                     correlationDefDict=Dict("SF-d0" => [("+-+-", [1,2,4,3], 1.), ("+-+-", [2,1,3,4], 1.)]),
-                     calculateThroughout=true,
-                    ) 
-display(resultsDict)
+p = CairoMakie.plot(mergResults["SF-d0"])
+save("merg.pdf", p)

src/correlations.jl

@@ -279,6 +279,11 @@ end
 export ThermalAverage
 
 
+"""
+    SpecFunc(eigVals, eigVecs, probe, freqValues, standDev)
+
+Fundamental implementation of spectral function calculation given a spectrum and excitation operators(defined by probe).
+"""
 function SpecFunc(
     eigVals::Vector{Float64},
     eigVecs::Vector{Vector{Float64}},
@@ -320,16 +325,14 @@ function SpecFunc(
         specFunc = specFunc ./ areaSpecFunc
     end
     return specFunc
-
 end
 export SpecFunc
 
 
 """
     SpecFunc(eigVals, eigVecs, probe, probeDiag, freqValues, standDev)
 
-Calculate the spectral function for the excitations defined by probe
-and probeDiag.
+Extends SpecFunc above by allowing passing the spectrum in the form native to fermions.jl (dictionaries).
 
 # Examples
 ```jldoctest
@@ -437,6 +440,18 @@ end
 export SpecFunc
 
 
+"""
+    SpecFunc(eigVals, eigVecs, probe, probeDiag, freqValues, 
+             standDev, symmetries, groundStateSector)
+
+Extends SpecFunc() by allowing specific symmetry sectors from which to choose the ground state.
+
+# Examples
+```jldoctest
+julia> SpecFunc(eigenVals, eigenStates, probe, probeDag, freqValues, 1e-2, 
+                ['N']; groundStateSector=(3, 0))
+```
+"""
 function SpecFunc(
     eigVals::Vector{Float64},
     eigVecs::Vector{Dict{BitVector,Float64}},

src/iterDiag.jl

@@ -460,13 +460,9 @@ function IterDiag(
         end
     end
 
-    if calculateThroughout
-        resultsDict = Dict{String, Vector{Float64}}(name => Float64[] for name in keys(correlationDefDict))
-    else
-        resultsDict = Dict{String, Union{Nothing, Float64}}(name => nothing for name in keys(correlationDefDict))
-    end
+    resultsDict = Dict{String, Union{Nothing, Float64}}(name => nothing for name in keys(correlationDefDict))
     @assert "energyPerSite" ∉ keys(resultsDict)
-    resultsDict["energyPerSite"] = ifelse(calculateThroughout, Float64[], nothing)
+    resultsDict["energyPerSite"] = nothing
 
     pbar = Progress(length(hamltFlow); enabled=!silent)
     for (step, hamlt) in enumerate(hamltFlow)
@@ -490,11 +486,7 @@ function IterDiag(
                                            )
                      )
 
-            if calculateThroughout
-                push!(resultsDict["energyPerSite"], eigVals[1]/maximum(currentSites))
-            else
-                resultsDict["energyPerSite"] = eigVals[1]/maximum(currentSites)
-            end
+            resultsDict["energyPerSite"] = eigVals[1]/maximum(currentSites)
 
             indices = 1:length(eigVals)
             if !isnothing(corrQuantumNoReq) && !isnothing(quantumNos)
@@ -506,35 +498,14 @@ function IterDiag(
             for (name, correlationDef) in correlationDefDict
                 if !isnothing(corrOperatorDict[name])
                     correlationValue = finalState' * corrOperatorDict[name] * finalState
-                    if calculateThroughout
-                        push!(resultsDict[name], correlationValue)
-                    else
-                        resultsDict[name] = correlationValue
-                    end
+                    resultsDict[name] = correlationValue
                 end
             end
 
             next!(pbar; showvalues=[("Size", size(hamltMatrix))])
             break
         end
 
-        if calculateThroughout
-            push!(resultsDict["energyPerSite"], eigVals[1]/maximum(currentSites))
-
-            indices = 1:length(eigVals)
-            if !isnothing(corrQuantumNoReq) && !isnothing(quantumNos)
-                indices = findall(q -> corrQuantumNoReq(q, maximum(currentSites)), quantumNos)
-            end
-            currentState = rotation[:, indices[sortperm(eigVals[indices])[1]]]
-
-            # calculate correlations using the ground state of the final step
-            for (name, correlationDef) in correlationDefDict
-                if !isnothing(corrOperatorDict[name])
-                    push!(resultsDict[name], currentState' * corrOperatorDict[name] * currentState)
-                end
-            end
-        end
-
         # construct a basis and identity matrix for the new sites
         newBasis = BasisStates(length(newSitesFlow[step+1]))
 

test/iterDiag_tests.jl

@@ -31,13 +31,13 @@ function GetCorrelationOperators(totalSites)
     operators = Dict{String, Vector{Tuple{String, Vector{Int64}, Float64}}}()
 
     # 1-particle correlations
-    for _ in 1:5
+    for _ in 1:20
         indices = rand(Int, 2) .|> abs .|> n -> rem(n, totalSites) .+ 1
         operators[randstring()] = [(join([shuffle(qubitOperators)[1] for _ in indices]), indices, 1.)]
     end
 
     # 2-particle correlations
-    for _ in 1:5
+    for _ in 1:20
         indices = rand(Int, 4) .|> abs .|> n -> rem(n, totalSites) .+ 1
         operators[randstring()] = [(join([shuffle(qubitOperators)[1] for _ in indices]), indices, 1.)]
     end
@@ -48,7 +48,7 @@ end
 function GetVNEParty(totalSites)
     parties = Dict{String, Vector{Int64}}()
     for size in 1:minimum((totalSites - 1, 5))
-        for _ in 1:5
+        for _ in 1:20
             party = shuffle(1:totalSites)[1:size]
             parties[randstring()] = sort(party)
         end
@@ -61,33 +61,19 @@ function GetMIParties(totalSites)
     parties = Dict{String, NTuple{2, Vector{Int64}}}()
     for size1 in 1:minimum((totalSites - 1, 4))
         for size2 in 1:minimum((4, totalSites - size1))
-            party1 = shuffle(1:totalSites)[1:size1]
-            party2 = shuffle([i for i in 1:totalSites if i ∉ party1])[1:size2]
-            parties[randstring()] = (sort(party1), sort(party2))
+            for run in 1:10
+                party1 = shuffle(1:totalSites)[1:size1]
+                party2 = shuffle([i for i in 1:totalSites if i ∉ party1])[1:size2]
+                parties[randstring()] = (sort(party1), sort(party2))
+            end
         end
     end
     return parties
 end
 
 
-function GetSpecFuncChoices(totalSites)
-    objects = Dict{String, Vector{Tuple{String, Vector{Int64}, Float64}}}[]
-    for size in 1:2
-        for run in 1:4
-            left = shuffle(1:totalSites)[1:size]
-            right = shuffle(1:totalSites)[1:size]
-            operator = join(shuffle(['+', '-', 'n', 'h'])[1:size])
-            object = Dict("create" => [(operator, left, 1.)],
-                                         "destroy" => Dagger([(operator, right, 1.)]),                                 
-                                        )
-            push!(objects, object)
-        end
-    end
-    return objects
-end
-
 @testset "IterDiag Hubbard" begin
-    @testset for totalSites in 3:4
+    @testset for totalSites in 3:6
         hamiltonian = TestModel(totalSites)
         basis = BasisStates(totalSites)
         exactEnergies, exactEigvecs = Spectrum(hamiltonian, basis)
@@ -121,7 +107,7 @@ end
 
 @testset "IterDiag Hubbard, symmetries" begin
     chemPot = 0.01
-    @testset for totalSites in 3:4
+    @testset for totalSites in 3:6
         hamiltonian = TestModel(totalSites)
 
         basis = BasisStates(totalSites)

test/runtests.jl

@@ -2,9 +2,9 @@ using Test, LinearAlgebra, ProgressMeter, Serialization
 using fermions
 
 include("testing_helpers.jl")
-#=include("base_tests.jl")=#
-#=include("eigen_tests.jl")=#
-#=include("correlation_tests.jl")=#
-#=include("eigenstateRG_tests.jl")=#
+include("base_tests.jl")
+include("eigen_tests.jl")
+include("correlation_tests.jl")
+include("eigenstateRG_tests.jl")
 include("iterDiag_tests.jl")
-#=include("models_tests.jl")=#
+include("models_tests.jl")