diff --git a/examples/Kondo1CK.jl b/examples/Kondo1CK.jl index fa1ca22..246fd33 100644 --- a/examples/Kondo1CK.jl +++ b/examples/Kondo1CK.jl @@ -43,31 +43,34 @@ p1 = plot(thickness_scaling=1.6, leftmargin=-6mm, bottommargin=-3mm, label="appr p2 = plot(thickness_scaling=1.6, leftmargin=-6mm, bottommargin=-3mm, label="approx.", xlabel="sites", ylabel="energy per site") spinFlipCorrd2 = Tuple{String, Vector{Int64}, Float64}[("+-+-", [1, 2, 8, 7], 1.0), ("+-+-", [2, 1, 7, 8], 1.0)] spinFlipCorrd0 = Tuple{String, Vector{Int64}, Float64}[("+-+-", [1, 2, 4, 3], 1.0), ("+-+-", [2, 1, 3, 4], 1.0)] -savePaths, resultsDict = IterDiag(hamFlow, maxSize; +vneReq = Dict( + "SEE_Imp" => [1, 2], + "SEE_d2" => [1, 2, 7, 8], + "SEE_0" => [3, 4] + ) +@time savePaths, resultsDict = IterDiag(hamFlow, maxSize; symmetries=Char['N', 'S'], # symmetries=Char['S'], # symmetries=Char['N'], magzReq=(m, N) -> -1 ≤ m ≤ 3, occReq=(x, N) -> div(N, 2) - 3 ≤ x ≤ div(N, 2) + 3, correlationDefDict=Dict("SF0" => spinFlipCorrd0, "SF2" => spinFlipCorrd2), + vneSets=vneReq, + corrMagzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0), + corrOccReq=(x,N)->x==div(N, 2), ) -vneReq = Dict( - "SEE_Imp" => [1, 2], - "SEE_d2" => [1, 2, 7, 8], - "SEE_0" => [3, 4] - ) mutInfoReq = Dict( "MI_d0" => ([1,2],[3,4]), "MI_d1" => ([1,2],[5,6]), ) -@time vneResults = IterVNE(savePaths, vneReq; occReq=(x,N)->x==div(N, 2), magzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0)) +#=@time IterCorrelation(savePaths, spinFlipCorrd2, Dict(); occReq=(x,N)->x==div(N, 2), magzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0))=# +#=@time vneResults = IterVNE(savePaths, vneReq; occReq=(x,N)->x==div(N, 2), magzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0))=# # @time mutInfoResults = IterMutInfo(savePaths, mutInfoReq; occReq=(x,N)->x==div(N, 2), magzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0)) -display([deserialize(path)["results"]["SF2"] for path in savePaths[1:end-1] if !isnothing(deserialize(path)["results"]["SF2"])]) -display([deserialize(path)["eigVals"][1] / (initSites + i) for (i, path) in enumerate(savePaths[1:end-1])]) -display(vneResults) +# @time display(IterCorrelation(savePaths, spinFlipCorrd2, vneReq; occReq=(x,N)->x==div(N, 2), magzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0))) +display(resultsDict) # display(mutInfoResults) -scatter!(p1, [deserialize(path)["results"]["SF2"] for path in savePaths[1:end-1] if !isnothing(deserialize(path)["results"]["SF2"])]) -scatter!(p2, initSites:totalSites, [deserialize(path)["eigVals"][1] / (initSites + i) for (i, path) in enumerate(savePaths[1:end-1])]) +# scatter!(p1, [deserialize(path)["results"]["SF2"] for path in savePaths[1:end-1] if !isnothing(deserialize(path)["results"]["SF2"])]) +# scatter!(p2, initSites:totalSites, [deserialize(path)["eigVals"][1] / (initSites + i) for (i, path) in enumerate(savePaths[1:end-1])]) corrExact = [] energyExact = [] @@ -77,7 +80,7 @@ mutInfoExact = Dict(k => [] for k in keys(mutInfoReq)) basis = BasisStates(2 * (1 + num); totOccReq=[1 + num], magzReq=[ifelse(isodd(i), 0, 1)], localCriteria=x->x[1]+x[2]==1) fullHam = vcat(hamFlow[1:i]...) eigenVals, eigenStates = Spectrum(fullHam, basis) - push!(energyExact, minimum(eigenVals)/(1 + num)) + push!(energyExact, minimum(eigenVals)/(2*(1 + num))) for (k, m) in vneReq try push!(vneExact[k], VonNEntropy(eigenStates[1], m)) diff --git a/examples/KondoKSpace.jl b/examples/KondoKSpace.jl index 2790a95..e6bae92 100644 --- a/examples/KondoKSpace.jl +++ b/examples/KondoKSpace.jl @@ -6,7 +6,7 @@ totalSites = 7 initSites = 1 kondoJ = 1. bathInt = -1.5 -maxSize = 500 +maxSize = 50 dispersion = [ifelse(i % 2 == 1, -1., 1.) for i in 1:totalSites] function getHamFlow(initSites::Int64, totalSites::Int64, dispersion::Vector{Float64}, kondoJ::Float64, bathInt::Float64) @@ -58,13 +58,21 @@ p1 = plot(thickness_scaling=1.6, leftmargin=-6mm, bottommargin=-3mm, label="appr p2 = plot(thickness_scaling=1.6, leftmargin=-6mm, bottommargin=-3mm, label="approx.", xlabel="sites", ylabel="energy per site") spinFlipCorrd2 = Tuple{String, Vector{Int64}, Float64}[("+-+-", [1, 2, 8, 7], 1.0), ("+-+-", [2, 1, 7, 8], 1.0)] spinFlipCorrd0 = Tuple{String, Vector{Int64}, Float64}[("+-+-", [1, 2, 4, 3], 1.0), ("+-+-", [2, 1, 3, 4], 1.0)] + +vneReq = Dict( + "SEE_Imp" => [1, 2], + "SEE_1" => [5, 6], + ) savePaths, resultsDict = IterDiag(hamFlow, maxSize; symmetries=Char['N', 'S'], # symmetries=Char['S'], # symmetries=Char['N'], - #=magzReq=(m, N) -> -1 ≤ m ≤ 1,=# - #=occReq=(x, N) -> div(N, 2) - 3 ≤ x ≤ div(N, 2) + 3,=# + magzReq=(m, N) -> -1 ≤ m ≤ 3, + occReq=(x, N) -> div(N, 2) - 3 ≤ x ≤ div(N, 2) + 3, correlationDefDict=Dict("SF0" => spinFlipCorrd0, "SF2" => spinFlipCorrd2), + vneSets=vneReq, + corrMagzReq=(m,N)->m == ifelse(isodd(div(N, 2)), 1, 0), + corrOccReq=(x,N)->x==div(N, 2), ) display(vneResults) display(resultsDict["energyPerSite"]) diff --git a/src/iterDiag.jl b/src/iterDiag.jl index cdfe691..9c3de44 100644 --- a/src/iterDiag.jl +++ b/src/iterDiag.jl @@ -76,7 +76,9 @@ end function CombineRequirements( occReq::Union{Nothing,Function}, - magzReq::Union{Nothing,Function} + magzReq::Union{Nothing,Function}, + corrOccReq::Union{Nothing,Function}, + corrMagzReq::Union{Nothing,Function}, ) quantumNoReq = nothing if !isnothing(occReq) && isnothing(magzReq) @@ -86,7 +88,16 @@ function CombineRequirements( elseif !isnothing(occReq) && !isnothing(magzReq) quantumNoReq = (q, N) -> occReq(q[1], N) && magzReq(q[2], N) end - return quantumNoReq + + corrQuantumNoReq = nothing + if !isnothing(corrOccReq) && isnothing(corrMagzReq) + corrQuantumNoReq = (q, N) -> corrOccReq(q[1], N) + elseif isnothing(corrOccReq) && !isnothing(corrMagzReq) + corrQuantumNoReq = (q, N) -> corrMagzReq(q[1], N) + elseif !isnothing(corrOccReq) && !isnothing(corrMagzReq) + corrQuantumNoReq = (q, N) -> corrOccReq(q[1], N) && corrMagzReq(q[2], N) + end + return quantumNoReq, corrQuantumNoReq end @@ -271,8 +282,17 @@ function IterDiag( degenTol::Float64=1e-10, dataDir::String="data-iterdiag", correlationDefDict::Dict{String, Vector{Tuple{String, Vector{Int64}, Float64}}}=Dict{String, Tuple{String, Vector{Int64}, Float64}[]}(), + vneSets::Dict{String, Vector{Int64}}=Dict{String, Vector{Int64}}(), silent::Bool=false, + corrOccReq::Union{Nothing,Function}=nothing,#(o,N)->ifelse(div(N,2)-2 ≤ o ≤ div(N,2)+2, true, false), ## close to half-filling + corrMagzReq::Union{Nothing,Function}=nothing,#(m,N)->ifelse(m == N, true, false), ## maximally polarised states ) + + println(vneSets) + for hamlt in hamltFlow + println(hamlt) + end + println("-----") @assert length(hamltFlow) > 1 @assert all(∈("NS"), symmetries) if !isnothing(occReq) @@ -304,7 +324,7 @@ function IterDiag( end end - quantumNoReq = CombineRequirements(occReq, magzReq) + quantumNoReq, corrQuantumNoReq = CombineRequirements(occReq, magzReq, corrOccReq, corrMagzReq) currentSites = collect(1:maximum(maximum.([opMembers for (_, opMembers, _) in hamltFlow[1]]))) initBasis = BasisStates(maximum(currentSites)) @@ -321,16 +341,34 @@ function IterDiag( create = [[c1; c2] for (c1, c2) in zip(create, corrCreate)] basket = [[c1; c2] for (c1, c2) in zip(basket, corrBasket)] end + + requiredOperators = Dict() + for (name, parties) in vneSets + requiredOperators[name] = OperatorsRDM(parties) + vneCreate, vneBasket, newSitesFlow = UpdateRequirements(requiredOperators[name], newSitesFlow) + + for i in reverse(eachindex(vneCreate)) + if isempty(vneCreate[i]) && !isempty(vneCreate[i-1]) + vneBasket[i:end] .= repeat([vneCreate[i-1]], length(vneBasket) - i + 1) + break + end + end + + create = [[c1; c2] for (c1, c2) in zip(create, vneCreate)] + basket = [[c1; c2] for (c1, c2) in zip(basket, vneBasket)] + end + + operators = CreateProductOperator(create[1], operators, newSitesFlow[1]) quantumNos = InitiateQuantumNos(currentSites, symmetries, initBasis) corrOperatorDict = Dict{String, Union{Nothing, Matrix{Float64}}}(name => nothing for name in keys(correlationDefDict)) - resultsDict = Dict{String, Union{Nothing,Float64}}(name => nothing for name in keys(correlationDefDict)) + resultsDict = Dict{String, Union{Nothing,Float64}}() @assert "energyPerSite" ∉ keys(resultsDict) resultsDict["energyPerSite"] = nothing - pbar = Progress(length(hamltFlow); enabled=!silent)#, showvalues=[("Size", size(hamltMatrix))]) + pbar = Progress(length(hamltFlow); enabled=!silent) for (step, hamlt) in enumerate(hamltFlow) for (type, members, strength) in hamlt hamltMatrix += strength * operators[(type, members)] @@ -342,9 +380,6 @@ function IterDiag( if isnothing(corrOperatorDict[name]) && all(∈(keys(operators)), [(type, members) for (type, members, _) in correlationDef]) corrOperatorDict[name] = sum([coupling * operators[(type, members)] for (type, members, coupling) in correlationDef]) end - if !isnothing(corrOperatorDict[name]) - resultsDict[name] = rotation[:, 1]' * corrOperatorDict[name] * rotation[:, 1] - end end if step == length(hamltFlow) @@ -354,13 +389,44 @@ function IterDiag( "currentSites" => currentSites, "newSites" => newSitesFlow[step], "bondAntiSymmzer" => bondAntiSymmzer, - "results" => resultsDict, ) ) next!(pbar) + + indices = 1:size(rotation)[2] + if !isnothing(corrQuantumNoReq) + indices = findall(q -> corrQuantumNoReq(q, maximum(currentSites)), quantumNos) + @assert !isempty(indices) + end + gstate = rotation[:,argmin(eigVals[indices])] + + for (name, operator) in corrOperatorDict + resultsDict[name] = gstate' * operator * gstate + end + + densityMatrix = gstate * gstate' + for (name, parties) in collect(vneSets) + reducedDM = sum([operators[operator][indices, :]' * densityMatrix[indices,indices] * operators[operator][indices, :] for operator in requiredOperators[name]]) + reducedDM = reducedDM ./ tr(reducedDM) + eigVals = eigvals(Hermitian(reducedDM)) + filter!(isposdef, eigVals) + resultsDict[name] = -sum(log.(eigVals) .* eigVals) + end + break end + densityMatrix = rotation[:, 1] * rotation[:, 1]' + println(densityMatrix) + println(maximum(abs.(densityMatrix))) + for (name, parties) in collect(vneSets) + reducedDM = sum([operators[operator]' * densityMatrix * operators[operator] for operator in requiredOperators[name]]) + reducedDM = reducedDM ./ tr(reducedDM) + rdmeigVals = eigvals(Hermitian(reducedDM)) + filter!(isposdef, rdmeigVals) + resultsDict[name] = -sum(log.(rdmeigVals) .* rdmeigVals) + end + if length(eigVals) > maxSize rotation, eigVals, quantumNos = TruncateSpectrum(quantumNoReq, rotation, eigVals, maxSize, degenTol, currentSites, quantumNos) end @@ -376,7 +442,6 @@ function IterDiag( "newSites" => newSitesFlow[step], "bondAntiSymmzer" => bondAntiSymmzer, "identityEnv" => identityEnv, - "results" => resultsDict, ) ) @@ -384,7 +449,9 @@ function IterDiag( quantumNos = UpdateQuantumNos(newBasis, symmetries, newSitesFlow[step+1], quantumNos) end - hamltMatrix, operators, bondAntiSymmzer, corrOperator = UpdateOldOperators(eigVals, identityEnv, basket[step+1], operators, rotation, bondAntiSymmzer, corrOperatorDict) + hamltMatrix, operators, bondAntiSymmzer, corrOperator = UpdateOldOperators(eigVals, identityEnv, basket[step+1], operators, + rotation, bondAntiSymmzer, corrOperatorDict, + ) # define the qbit operators for the new sites for site in newSitesFlow[step+1] @@ -407,13 +474,16 @@ export IterDiag function IterCorrelation( savePaths::Vector{String}, - corrDef::Vector{Tuple{String, Vector{Int64}, Float64}}; + corrDef::Vector{Tuple{String, Vector{Int64}, Float64}}, + vneSets::Dict{String, Vector{Int64}}; occReq::Union{Nothing, Function}=nothing, magzReq::Union{Nothing, Function}=nothing, ) + + requiredOperators = Dict{String, Vector{Tuple{String, Vector{Int64}}}}(name => OperatorsRDM(parties) for (name, parties) in vneSets) newSitesFlow = [deserialize(savePath)["newSites"] for savePath in savePaths[1:end-1]] - create, basket, newSitesFlow = UpdateRequirements(corrDef, newSitesFlow) + create, basket, newSitesFlow = UpdateRequirements(vcat([(x,y) for (x,y,z) in values(corrDef)], vcat(values(requiredOperators)...)), newSitesFlow) quantumNoReq = nothing if !isnothing(occReq) && isnothing(magzReq) @@ -423,7 +493,8 @@ function IterCorrelation( elseif !isnothing(occReq) && !isnothing(magzReq) quantumNoReq = (q, N) -> occReq(q[1], N) && magzReq(q[2], N) end - corrVals = Float64[] + corrVals = nothing + vneResults = Dict{String, Vector{Float64}}(name => Float64[] for name in keys(vneSets)) metadata = deserialize(savePaths[end]) symmetries = metadata["symmetries"] if !isnothing(occReq) @@ -451,38 +522,49 @@ function IterCorrelation( eigVals = f["eigVals"] quantumNos = f["quantumNos"] currentSites = f["currentSites"] - if !isnothing(corrOperator) - if isnothing(quantumNoReq) - gstate = basis[:,argmin(eigVals)] - else + if step == length(savePaths) - 1 + indices = 1:size(basis)[2] + if !isnothing(quantumNoReq) indices = findall(q -> quantumNoReq(q, maximum(currentSites)), quantumNos) - gstate = basis[:, indices][:,argmin(eigVals[indices])] - end - push!(corrVals, gstate' * corrOperator * gstate) - if step < length(savePaths) - 1 - corrOperator = kron(basis' * corrOperator * basis, f["identityEnv"]) - end - else - for (k, v) in collect(operators) - operators[k] = kron(basis' * v * basis, f["identityEnv"]) + @assert !isempty(indices) unique(quantumNos) end + gstate = basis[:,argmin(eigVals[indices])] - newBasis = BasisStates(length(newSitesFlow[step+1])) + corrVals = gstate' * corrOperator * gstate - # define the qbit operators for the new sites - for site in newSitesFlow[step+1] - operators[("+", [site])] = kron(basis' * f["bondAntiSymmzer"] * basis, OperatorMatrix(newBasis, [("+", [site - length(currentSites)], 1.0)])) - operators = CreateDNH(operators, site) - end + densityMatrix = gstate * gstate' - operators = CreateProductOperator(create[step+1], operators, newSitesFlow[step+1]) - if all(∈(keys(operators)), [(type, members) for (type, members, _) in corrDef]) - corrOperator = sum([coupling * operators[(type, members)] for (type, members, coupling) in corrDef]) + for (name, parties) in collect(vneSets) + reducedDM = sum([operators[operator][indices,:]' * densityMatrix[indices,indices] * operators[operator][indices,:] for operator in requiredOperators[name]])#[indices,indices] + reducedDM = reducedDM ./ tr(reducedDM) + eigVals = eigvals(Hermitian(reducedDM)) + filter!(isposdef, eigVals) + push!(vneResults[name], -sum(log.(eigVals) .* eigVals)) end + break + end + + for (k, v) in collect(operators) + operators[k] = kron(basis' * v * basis, f["identityEnv"]) + end + + newBasis = BasisStates(length(newSitesFlow[step+1])) + # define the qbit operators for the new sites + for site in newSitesFlow[step+1] + operators[("+", [site])] = kron(basis' * f["bondAntiSymmzer"] * basis, OperatorMatrix(newBasis, [("+", [site - length(currentSites)], 1.0)])) + operators = CreateDNH(operators, site) + end + + operators = CreateProductOperator(create[step+1], operators, newSitesFlow[step+1]) + if all(∈(keys(operators)), [(type, members) for (type, members, _) in corrDef]) + corrOperator = sum([coupling * operators[(type, members)] for (type, members, coupling) in corrDef]) end + end - return corrVals + + return corrVals, vneResults + end export IterCorrelation @@ -523,13 +605,32 @@ function IterVNE( operators = CreateProductOperator(create[1], operators, newSitesFlow[1]) - for (step, savePath) in enumerate(savePaths[1:end-2]) + for (step, savePath) in enumerate(savePaths[1:end-1]) f = deserialize(savePath) basis = f["basis"] eigVals = f["eigVals"] quantumNos = f["quantumNos"] currentSites = f["currentSites"] + if step == length(savePaths)-1 + indices = 1:size(basis)[2] + if !isnothing(quantumNoReq) + indices = findall(q -> quantumNoReq(q, maximum(currentSites)), quantumNos) + @assert !isempty(indices) unique(quantumNos) + end + gstate = basis[:,argmin(eigVals[indices])] + densityMatrix = gstate * gstate' + + for (name, parties) in collect(vneSets) + reducedDM = sum([operators[operator][indices,:]' * densityMatrix[indices,indices] * operators[operator][indices,:] for operator in requiredOperators[name]])#[indices,indices] + reducedDM = reducedDM ./ tr(reducedDM) + eigVals = eigvals(Hermitian(reducedDM)) + filter!(isposdef, eigVals) + push!(vneResults[name], -sum(log.(eigVals) .* eigVals)) + end + break + end + @time for (k, v) in collect(operators) operators[k] = kron(basis' * v * basis, f["identityEnv"]) end @@ -546,29 +647,6 @@ function IterVNE( end - @time begin - f = deserialize(savePaths[end-1]) - basis = f["basis"] - eigVals = f["eigVals"] - quantumNos = f["quantumNos"] - currentSites = f["currentSites"] - - indices = 1:size(basis)[2] - if !isnothing(quantumNoReq) - indices = findall(q -> quantumNoReq(q, maximum(currentSites)), quantumNos) - @assert !isempty(indices) unique(quantumNos) - end - gstate = basis[:,argmin(eigVals[indices])] - densityMatrix = gstate * gstate' - - for (name, parties) in collect(vneSets) - reducedDM = sum([operators[operator][indices,:]' * densityMatrix[indices,indices] * operators[operator][indices,:] for operator in requiredOperators[name]])#[indices,indices] - reducedDM = reducedDM ./ tr(reducedDM) - eigVals = eigvals(Hermitian(reducedDM)) - filter!(isposdef, eigVals) - push!(vneResults[name], -sum(log.(eigVals) .* eigVals)) - end - end return vneResults end