Skip to content

Commit

Permalink
refactored and brought entanglement calculations into the direct iter…
Browse files Browse the repository at this point in the history 
…ation process
  • Loading branch information
@abhirup-m
abhirup-m committed Oct 1, 2024
1 parent 9741a98 commit 2ab0355
Show file tree
Hide file tree
Showing 2 changed files with 140 additions and 74 deletions.
29 changes: 16 additions & 13 deletions examples/Kondo1CK.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = []
Expand All @@ -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))
Expand Down
185 changes: 124 additions & 61 deletions src/iterDiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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


Expand Down Expand Up @@ -271,7 +282,10 @@ 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
)
@assert length(hamltFlow) > 1
@assert all(∈("NS"), symmetries)
Expand Down Expand Up @@ -304,7 +318,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))
Expand All @@ -321,16 +335,36 @@ function IterDiag(
create = [[c1; c2] for (c1, c2) in zip(create, corrCreate)]
basket = [[c1; c2] for (c1, c2) in zip(basket, corrBasket)]
end

display(basket)
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

display(basket)

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)]
Expand All @@ -342,9 +376,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)
Expand All @@ -354,10 +385,30 @@ 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])]
densityMatrix = gstate * gstate'

for (name, operator) in corrOperatorDict
resultsDict[name] = gstate' * operator * gstate
end

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

Expand All @@ -376,15 +427,16 @@ function IterDiag(
"newSites" => newSitesFlow[step],
"bondAntiSymmzer" => bondAntiSymmzer,
"identityEnv" => identityEnv,
"results" => resultsDict,
)
)

if !isnothing(quantumNos)
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]
Expand All @@ -407,13 +459,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)
Expand All @@ -423,7 +478,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)
Expand Down Expand Up @@ -451,38 +507,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

Expand Down Expand Up @@ -523,13 +590,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
Expand All @@ -546,29 +632,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

Expand Down

0 comments on commit 2ab0355

Please sign in to comment.