diff --git a/src/frontend/diagram_id.jl b/src/frontend/diagram_id.jl index 7d6935d7..4f0f37ad 100644 --- a/src/frontend/diagram_id.jl +++ b/src/frontend/diagram_id.jl @@ -21,6 +21,9 @@ struct BareGreenId <: PropagatorId type::AnalyticProperty #Instant, Dynamic extK::Vector{Float64} extT::Tuple{Int,Int} #all possible extT from different interactionType + function BareGreenId(type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {T<:Real} + return new(type, mirror_symmetrize(k), t) + end function BareGreenId(type::AnalyticProperty=Dynamic; k, t) return new(type, mirror_symmetrize(k), Tuple(t)) end @@ -35,6 +38,9 @@ struct BareInteractionId <: PropagatorId # bare W-type interaction, with only on type::AnalyticProperty #Instant, Dynamic extK::Vector{Float64} extT::Tuple{Int,Int} #all possible extT from different interactionType + function BareInteractionId(response::Response, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {T<:Real} + return new(response, type, mirror_symmetrize(k), t) + end function BareInteractionId(response::Response, type::AnalyticProperty=Instant; k, t=(0, 0)) return new(response, type, mirror_symmetrize(k), Tuple(t)) end @@ -66,7 +72,7 @@ end struct GenericId{P} <: DiagramId para::P extra::Any - GenericId(para::P, extra=Nothing) where {P} = new{P}(para, extra) + GenericId(para::P, extra=nothing) where {P} = new{P}(para, extra) end Base.show(io::IO, v::GenericId) = print(io, v.extra == Nothing ? "" : "$(v.extra)") function Base.isequal(a::GenericId, b::GenericId) @@ -95,6 +101,9 @@ struct GreenId{P} <: DiagramId type::AnalyticProperty #Instant, Dynamic extK::Vector{Float64} extT::Tuple{Int,Int} #all possible extT from different interactionType + function GreenId(para::P, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real} + return new{P}(para, type, mirror_symmetrize(k), t) + end function GreenId(para::P, type::AnalyticProperty=Dynamic; k, t) where {P} return new{P}(para, type, mirror_symmetrize(k), Tuple(t)) end @@ -109,6 +118,9 @@ struct SigmaId{P} <: DiagramId type::AnalyticProperty #Instant, Dynamic extK::Vector{Float64} extT::Tuple{Int,Int} #all possible extT from different interactionType + function SigmaId(para::P, type::AnalyticProperty, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real} + return new{P}(para, type, mirror_symmetrize(k), t) + end function SigmaId(para::P, type::AnalyticProperty; k, t=(0, 0)) where {P} return new{P}(para, type, mirror_symmetrize(k), Tuple(t)) end @@ -126,7 +138,9 @@ struct PolarId{P} <: DiagramId response::Response #UpUp, UpDown, ... extK::Vector{Float64} extT::Tuple{Int,Int} #all possible extT from different interactionType - order::Vector{Int} + function PolarId(para::P, response::Response, k::Vector{T}, t::Tuple{Int,Int}) where {P,T<:Real} + return new{P}(para, response, mirror_symmetrize(k), t) + end function PolarId(para::P, response::Response; k, t=(0, 0)) where {P} return new{P}(para, response, mirror_symmetrize(k), Tuple(t)) end @@ -136,7 +150,7 @@ function Base.isequal(a::PolarId, b::PolarId) if typeof(a) != typeof(b) return false end - return a.response == b.response && a.extT == b.extT && a.order == b.order && a.extK == b.extK && a.para == b.para + return a.response == b.response && a.extT == b.extT && a.extK == b.extK && a.para == b.para end struct Ver3Id{P} <: DiagramId @@ -144,6 +158,9 @@ struct Ver3Id{P} <: DiagramId response::Response #UpUp, UpDown, ... extK::Vector{Vector{Float64}} extT::Tuple{Int,Int,Int} #all possible extT from different interactionType + function Ver3Id(para::P, response::Response, k::Vector{Vector{T}}, t::Tuple{Int,Int,Int}) where {P,T<:Real} + return new{P}(para, response, k, t) + end function Ver3Id(para::P, response::Response; k, t=(0, 0, 0)) where {P} return new{P}(para, response, k, Tuple(t)) end @@ -163,6 +180,9 @@ struct Ver4Id{P} <: DiagramId channel::TwoBodyChannel # particle-hole, particle-hole exchange, particle-particle, irreducible extK::Vector{Vector{Float64}} extT::Tuple{Int,Int,Int,Int} #all possible extT from different interactionType + function Ver4Id(para::P, response::Response, type::AnalyticProperty, chan::TwoBodyChannel, k::Vector{Vector{T}}, t::NTuple{4,Int}) where {P,T<:Real} + return new{P}(para, response, type, chan, k, t) + end function Ver4Id(para::P, response::Response, type::AnalyticProperty=Dynamic; k, t=(0, 0, 0, 0), chan::TwoBodyChannel=AnyChan) where {P} return new{P}(para, response, type, chan, k, Tuple(t)) @@ -215,7 +235,7 @@ struct BareHoppingId{P} <: PropagatorId site::Tuple{Int,Int} orbital::Tuple{Int,Int} extT::Tuple{Int,Int} - function BareHoppingId(para::P, orbital, t, r) where {P} + function BareHoppingId(para::P, r::Tuple{Int,Int}, orbital::Tuple{Int,Int}, t::Tuple{Int,Int}) where {P} return new{P}(para, r, orbital, t) end end @@ -237,6 +257,10 @@ struct BareGreenNId{P} <: PropagatorId orbital::Vector{Int} extT::Vector{Int} N::Int + function BareGreenNId(para::P, r::Int, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P} + @assert length(orbital) == length(t) == length(creation) == N + return new{P}(para, r, creation, orbital, t, N) + end function BareGreenNId(para::P; orbital=[], t=[], creation=[], r=0) where {P} @assert length(orbital) == length(t) == length(creation) return new{P}(para, r, creation, orbital, t, length(orbital)) @@ -260,6 +284,10 @@ struct GreenNId{P} <: DiagramId orbital::Vector{Int} extT::Vector{Int} N::Int + function GreenNId(para::P, r::Vector{Int}, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P} + @assert length(orbital) == length(t) == length(r) == length(creation) == N + return new{P}(para, r, creation, orbital, t, N) + end function GreenNId(para::P; orbital=[], t=[], creation=[], r=[]) where {P} @assert length(orbital) == length(t) == length(r) == length(creation) return new{P}(para, r, creation, orbital, t, length(orbital)) @@ -283,6 +311,10 @@ struct ConnectedGreenNId{P} <: DiagramId orbital::Vector{Int} extT::Vector{Int} N::Int + function ConnectedGreenNId(para::P, r::Vector{Int}, creation::Vector{Bool}, orbital::Vector{Int}, t::Vector{Int}, N::Int=length(orbital)) where {P} + @assert length(orbital) == length(t) == length(r) == length(creation) == N + return new{P}(para, r, creation, orbital, t, N) + end function ConnectedGreenNId(para::P; orbital=[], t=[], creation=[], r=[]) where {P} @assert length(orbital) == length(t) == length(r) == length(creation) return new{P}(para, r, creation, orbital, t, length(orbital)) @@ -322,4 +354,32 @@ function index(type) end end +""" + reconstruct(instance::DiagramId, updates::Pair{Symbol}...) + +Create a new instance of the same type as `instance`, with specified fields updated to new values. + +# Usage +new_instance = reconstruct(old_instance, :field1 => new_value1, :field2 => new_value2) +""" +function reconstruct(instance::DiagramId, updates::Pair{Symbol}...) + # Get the type of the instance + T = typeof(instance) + + # Extract field names and values from the instance + field_names = fieldnames(T) + field_values = [getfield(instance, fn) for fn in field_names] + + # Update fields based on the updates provided + for (field, new_value) in updates + field_idx = findfirst(==(field), field_names) + if field_idx !== nothing + field_values[field_idx] = new_value + else + throw(ArgumentError("Field $field does not exist in type $T")) + end + end + # Construct a new instance with the updated field values + return Base.typename(T).wrapper(field_values...) +end \ No newline at end of file diff --git a/src/frontend/parquet/common.jl b/src/frontend/parquet/common.jl index a7b0ae1d..586c26b0 100644 --- a/src/frontend/parquet/common.jl +++ b/src/frontend/parquet/common.jl @@ -1,10 +1,10 @@ -function build(para::DiagPara, extK=nothing, subdiagram=false) +function build(para::DiagPara, extK=nothing, subdiagram=false; channels=[PHr, PHEr, PPr, Alli]) if para.type == Ver4Diag if isnothing(extK) extK = [getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2), getK(para.totalLoopNum, 3)] end - return vertex4(para, extK, [PHr, PHEr, PPr, Alli], subdiagram) + return vertex4(para, extK, subdiagram, channels=channels) elseif para.type == SigmaDiag if isnothing(extK) extK = getK(para.totalLoopNum, 1) @@ -19,7 +19,7 @@ function build(para::DiagPara, extK=nothing, subdiagram=false) if isnothing(extK) extK = [getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2)] end - return vertex3(para, extK, subdiagram) + return vertex3(para, extK, subdiagram, channels=channels) else error("not implemented!") end diff --git a/src/frontend/parquet/ep_coupling.jl b/src/frontend/parquet/ep_coupling.jl index 050f3dc7..ac4a66bf 100644 --- a/src/frontend/parquet/ep_coupling.jl +++ b/src/frontend/parquet/ep_coupling.jl @@ -112,7 +112,7 @@ function ep_bubble!(ver4df::DataFrame, para::DiagPara, legK, chans::Vector{TwoBo LLegK, K, RLegK, Kx = legBasis(PHr, legK, LoopIdx) - Lver = vertex4(lPara, LLegK, chans, true; name=:Γf, blocks=blocks) + Lver = vertex4(lPara, LLegK, true; channels=chans, name=:Γf, blocks=blocks) isempty(Lver) && return Rver = DataFrame(response=Response[], type=AnalyticProperty[], extT=Tuple{Int,Int,Int,Int}[], diagram=Graph{Ftype,Wtype}[]) diff --git a/src/frontend/parquet/operation.jl b/src/frontend/parquet/operation.jl index e8897f6b..32149b93 100644 --- a/src/frontend/parquet/operation.jl +++ b/src/frontend/parquet/operation.jl @@ -139,62 +139,53 @@ end # mergeby(df::DataFrame; kwargs...) = mergeby(df, []; kwargs...) # mergeby(diags::Vector{Graph}; kwargs...) = mergeby(diags, []; kwargs...) +""" + update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}}) -function update_extK!(diag::Graph, extK::Vector{Vector{Float64}}) - visited = Set{Int}() - num_extK = length(extK) - len_extK = length(extK[1]) - - sumK = zeros(len_extK) - _K = zeros(len_extK) - for leaf in Leaves(diag) - if !(leaf.id in visited) - push!(visited, leaf.id) - prop = IR.properties(leaf) - K = prop.extK - - original_len_K = length(K) - if length(K) < len_extK - resize!(K, len_extK) - K[original_len_K+1:end] .= 0.0 - else - resize!(K, len_extK) - end - - _K[num_extK+1:end] .= K[num_extK+1:end] - for i in eachindex(extK) - sumK .+= K[i] * extK[i] - end - K .= sumK .+ _K - fill!(sumK, 0.0) - # K[1:end] = sum([K[i] * extK[i] for i in eachindex(extK)]) + _K - - # if prop isa BareGreenId - # new_properties = BareGreenId(prop.type, k=K, t=prop.extT) - # elseif prop isa BareInteractionId - # new_properties = BareInteractionId(prop.response, prop.type, k=K, t=prop.extT) - # else - # error("unexpected property type $prop") - # end - # IR.set_properties!(leaf, new_properties) - end - end -end + Update the external momenta (`extK`) and external times (`extT`) of all the nodes in a vector of graphs in-place. -function update_extK!(diags::Vector{Graph}, extK::Vector{Vector{Float64}}) +# Arguments +- `diags::Vector{Graph}`: A vector of `Graph` objects. +- `para::DiagPara`: parameters reconstructed in the graphs. Its `firstTauIdx` will update the `extT` of graphs. +- `legK::Vector{Vector{Float64}}`: basus of the external momenta for the legs of the diagram as [left in, left out, right in, right out]. +""" +function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}}) visited = Set{Int}() - num_extK = length(extK) - len_extK = length(extK[1]) + tauIdx = para.firstTauIdx + # len_extK = para.totalLoopNum + # num_extK = len_extK - para.innerLoopNum + # extK = [k[1:len_extK] for k in legK[1:num_extK]] + # if para.totalLoopNum - para.innerLoopNum != 3 + # println(legK) + # println(para) + # end + len_extK = length(legK[1]) + num_extK = length(legK) - 1 + extK = legK[1:end-1] sumK = zeros(len_extK) _K = zeros(len_extK) - for diag in diags - for leaf in Leaves(diag) - if !(leaf.id in visited) - push!(visited, leaf.id) - prop = IR.properties(leaf) - K = prop.extK + for graph in diags + for node in PreOrderDFS(graph) + node.id in visited && continue + node.id = IR.uid() + push!(visited, node.id) + prop = IR.properties(node) + K = prop.extK + T = prop.extT + tau_shift = tauIdx - T[1] + if prop isa Ver4Id || prop isa Ver3Id + for i in eachindex(K) + resize!(K[i], len_extK) + K[i] .= legK[i][1:len_extK] + end + if tau_shift != 0 + node.properties = FrontEnds.reconstruct(prop, :para => para, :extT => Tuple(t + tau_shift for t in T)) + else + node.properties = FrontEnds.reconstruct(prop, :para => para) + end + else original_len_K = length(K) if length(K) < len_extK resize!(K, len_extK) @@ -209,7 +200,30 @@ function update_extK!(diags::Vector{Graph}, extK::Vector{Vector{Float64}}) end K .= sumK .+ _K fill!(sumK, 0.0) + if tau_shift != 0 + node.properties = FrontEnds.reconstruct(prop, :extT => Tuple(t + tau_shift for t in T)) + end end end end +end + +""" + update_extKT(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}}) -> Vector{Graph} + + Returns a new vector of graphs with updated external momenta (`extK`) and external times (`extT`), + based on the provided graphs, parameters, and external legs' momenta. + +# Arguments +- `diags::Vector{Graph}`: A vector of `Graph` objects. +- `para::DiagPara`: parameters reconstructed in the graphs. Its `firstTauIdx` will update the `extT` of graphs. +- `legK::Vector{Vector{Float64}}`: basus of the external momenta for the legs of the diagram as [left in, left out, right in, right out]. + +# Returns +- `Vector{Graph}`: A new vector of `Graph` objects with updated `extK`, `extT`, and `para` (if existed) properties for each node. +""" +function update_extKT(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector{Float64}}) + graphs = deepcopy(diags) + update_extKT!(graphs, para, legK) + return graphs end \ No newline at end of file diff --git a/src/frontend/parquet/parquet.jl b/src/frontend/parquet/parquet.jl index 68851cb1..5260a728 100644 --- a/src/frontend/parquet/parquet.jl +++ b/src/frontend/parquet/parquet.jl @@ -9,6 +9,7 @@ Ftype, Wtype = _dtype.factor, _dtype.weight import ..Taylor: set_variables import ..Utility: taylorexpansion! +import ..FrontEnds import ..FrontEnds: TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan import ..FrontEnds: Filter, NoBubble, NoHartree, NoFock, DirectOnly, Wirreducible, Girreducible, Proper import ..FrontEnds: Response, Composite, ChargeCharge, SpinSpin, ProperChargeCharge, ProperSpinSpin, UpUp, UpDown diff --git a/src/frontend/parquet/sigma.jl b/src/frontend/parquet/sigma.jl index 2634ae56..a5d4c736 100644 --- a/src/frontend/parquet/sigma.jl +++ b/src/frontend/parquet/sigma.jl @@ -95,9 +95,9 @@ function sigma(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=false if oW == 0 # Fock-type Σ if NoHartree in paraW.filter paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=zero(K)) - ver4 = vertex4(paraW0, legK, [], true) + ver4 = vertex4(paraW0, legK, true, channels=[]) else - ver4 = vertex4(paraW, legK, [], true) + ver4 = vertex4(paraW, legK, true, channels=[]) end # println(ver4) else # composite Σ @@ -106,7 +106,7 @@ function sigma(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=false # if compact # ver4 = ep_coupling(paraW; extK=legK, subdiagram=true, name=:W, blocks=blocks) # else - ver4 = vertex4(paraW, legK, [PHr,], true; blocks=blocks, blockstoplevel=ParquetBlocks(phi=[], Γ4=[PHr, PHEr, PPr])) + ver4 = vertex4(paraW, legK, true; channels=[PHr,], blocks=blocks, blockstoplevel=ParquetBlocks(phi=[], Γ4=[PHr, PHEr, PPr])) # end end #transform extT coloum intwo extT for Σ and extT for G diff --git a/src/frontend/parquet/sigmaGV.jl b/src/frontend/parquet/sigmaGV.jl index e15826e4..29d48dc0 100644 --- a/src/frontend/parquet/sigmaGV.jl +++ b/src/frontend/parquet/sigmaGV.jl @@ -98,7 +98,7 @@ function sigmaGV(para::DiagPara, extK=getK(para.totalLoopNum, 1), subdiagram=fal if isValidG(paraG) # println(paraW) paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=zero(K)) - bareV = vertex4(paraW0, legK, [], true) + bareV = vertex4(paraW0, legK, true, channels=[]) if oW == 0 # Fock-type Σ ver4 = bareV df = transform(ver4, :extT => ByRow(x -> [(x[INL], x[OUTR]), (x[OUTL], x[INR])]) => [:extT, :GT]) diff --git a/src/frontend/parquet/to_graph.jl b/src/frontend/parquet/to_graph.jl index 7674177e..de532276 100644 --- a/src/frontend/parquet/to_graph.jl +++ b/src/frontend/parquet/to_graph.jl @@ -1,4 +1,3 @@ - """ function diagdict_parquet(type::Symbol, MaxOrder::Int, has_counterterm::Bool=true; MinOrder::Int=1, spinPolarPara::Float64=0.0, isDynamic=false, filter=[NoHartree]) @@ -20,7 +19,7 @@ The key is (order, Gorder, Vorder). The element is a Tuple (graphVector, extT_labels). """ function diagdict_parquet(type::Symbol, MaxOrder::Int, MinOrder::Int=1; - has_counterterm::Bool=true, spinPolarPara::Float64=0.0, optimize_level=0, + has_counterterm::Bool=true, spinPolarPara::Float64=0.0, optimize_level=0, channels=[PHr, PHEr, PPr, Alli], isDynamic=false, filter=[NoHartree], transferLoop=nothing, extK=nothing) diagtype = _diagtype(type) @@ -32,7 +31,7 @@ function diagdict_parquet(type::Symbol, MaxOrder::Int, MinOrder::Int=1; for order in MinOrder:MaxOrder set_variables("x y"; orders=[MaxOrder - order, MaxOrder - order]) para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop) - parquet_builder = Parquet.build(para, extK) + parquet_builder = Parquet.build(para, extK, channels=channels) diags, extT = parquet_builder.diagram, parquet_builder.extT propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true]) # Specify variable dependence of fermi (first element) and bose (second element) particles. @@ -54,7 +53,7 @@ function diagdict_parquet(type::Symbol, MaxOrder::Int, MinOrder::Int=1; set_variables("x y"; orders=[0, 0]) for order in MinOrder:MaxOrder para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop) - parquet_builder = Parquet.build(para, extK) + parquet_builder = Parquet.build(para, extK, channels=channels) diags, extT = parquet_builder.diagram, parquet_builder.extT propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true]) # Specify variable dependence of fermi (first element) and bose (second element) particles. @@ -68,17 +67,18 @@ function diagdict_parquet(type::Symbol, MaxOrder::Int, MinOrder::Int=1; dict_graphs[key] = ([graph,], extT) end end + dict_graphs[(order, 0, 0)] = (diags, collect.(extT)) end end - # for gvec in values(dict_graphs) - # IR.optimize!(gvec[1], level=optimize_level) - # end + for gvec in values(dict_graphs) + IR.optimize!(gvec[1], level=optimize_level) + end return dict_graphs end function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}; - spinPolarPara::Float64=0.0, optimize_level=0, + spinPolarPara::Float64=0.0, optimize_level=0, channels=[PHr, PHEr, PPr, Alli], isDynamic=false, filter=[NoHartree], transferLoop=nothing, extK=nothing) diagtype = _diagtype(type) @@ -91,7 +91,8 @@ function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}; for order in diag_orders set_variables("x y"; orders=[Gorder, Vorder]) para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop) - parquet_builder = Parquet.build(para, extK) + parquet_builder = Parquet.build(para, extK, channels=channels) + diags, extT = parquet_builder.diagram, parquet_builder.extT propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true]) # Specify variable dependence of fermi (first element) and bose (second element) particles. @@ -117,7 +118,7 @@ function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}; end function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}, extra_variables::Dict{String,Int}; - spinPolarPara::Float64=0.0, optimize_level=0, + spinPolarPara::Float64=0.0, optimize_level=0, channels=[PHr, PHEr, PPr, Alli], isDynamic=false, filter=[NoHartree], transferLoop=nothing, extK=nothing) diagtype = _diagtype(type) @@ -139,7 +140,7 @@ function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}, extra # set_variables("x y k"; orders=[MaxOrder - order, MaxOrder - order, 1]) set_variables("x y" * extra_varnames; orders=[Gorder, Vorder, extra_orders...]) para = diagPara(diagtype, isDynamic, spin, order, filter, transferLoop) - parquet_builder = Parquet.build(para, extK) + parquet_builder = Parquet.build(para, extK, channels=channels) diags, extT = parquet_builder.diagram, parquet_builder.extT var_dependence = Dict{Int,Vector{Bool}}() @@ -209,6 +210,10 @@ function diagPara(type, isDynamic::Bool, spin, order, filter, transferLoop=nothi innerLoopNum = order end + if Proper in filter + @assert !isnothing(transferLoop) "transferLoop must be provided if Proper is in filter" + end + if isnothing(transferLoop) return DiagPara( type=type, diff --git a/src/frontend/parquet/vertex3.jl b/src/frontend/parquet/vertex3.jl index 421c1015..8d39d637 100644 --- a/src/frontend/parquet/vertex3.jl +++ b/src/frontend/parquet/vertex3.jl @@ -1,6 +1,6 @@ """ function vertex3(para::DiagPara, _extK=[getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2)], subdiagram=false; - name=:Γ3, chan=[PHr, PHEr, PPr, Alli], resetuid=false, blocks::ParquetBlocks=ParquetBlocks()) + name=:Γ3, channels=[PHr, PHEr, PPr, Alli], resetuid=false, blocks::ParquetBlocks=ParquetBlocks()) Generate 3-vertex diagrams using Parquet Algorithm. With imaginary-time variables, all vertex3 generated has the same bosonic Tidx ``extT[1]=para.firstTauIdx`` and the incoming fermionic Tidx ``extT[2]=para.firstTauIdx+1``. @@ -10,7 +10,7 @@ With imaginary-time variables, all vertex3 generated has the same bosonic Tidx ` - `extK` : basis of external loops as a vector [bosonic leg (out), fermionic in, fermionic out], extK[1] = extK[2] - extK[3]. - `subdiagram` : a sub-vertex or not - `name` : name of the vertex -- `chan` : vector of channels of the current 4-vertex. +- `channels` : vector of channels of the current 4-vertex. - `resetuid` : restart uid count from 1 - `blocks` : building blocks of the Parquet equation. See the struct ParquetBlocks for more details. @@ -21,7 +21,7 @@ function vertex3(para::DiagPara, _extK=[getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2)], subdiagram=false; name=:Γ3, - chan=[PHr, PHEr, PPr, Alli], + channels=[PHr, PHEr, PPr, Alli], resetuid=false, blocks::ParquetBlocks=ParquetBlocks() ) @@ -74,7 +74,7 @@ function vertex3(para::DiagPara, firstLoopIdx=GoutKidx, firstTauIdx=GoutTidx) paraVer4 = reconstruct(para, type=Ver4Diag, innerLoopNum=oVer4, firstLoopIdx=Ver4Kidx, firstTauIdx=Ver4Tidx) - ver4 = vertex4(paraVer4, legK, chan, true; blocks=blocks) + ver4 = vertex4(paraVer4, legK, true; channels=channels, blocks=blocks) if isnothing(ver4) || isempty(ver4) continue end diff --git a/src/frontend/parquet/vertex4.jl b/src/frontend/parquet/vertex4.jl index acbb8d5a..432e83d4 100644 --- a/src/frontend/parquet/vertex4.jl +++ b/src/frontend/parquet/vertex4.jl @@ -1,8 +1,8 @@ """ vertex4(para::DiagPara, extK = [getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2), getK(para.totalLoopNum, 3)], - chan::AbstractVector = [PHr, PHEr, PPr, Alli], subdiagram = false; + channels::AbstractVector = [PHr, PHEr, PPr, Alli], level = 1, name = :none, resetuid = false, blocks::ParquetBlocks=ParquetBlocks(), blockstoplevel::ParquetBlocks=blocks @@ -13,8 +13,8 @@ Generate 4-vertex diagrams using Parquet Algorithm # Arguments - `para` : parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx - `extK` : basis of external loops as a vector [left in, left out, right in, right out]. -- `chan` : vector of channels of the current 4-vertex. - `subdiagram` : a sub-vertex or not +- `channels` : vector of channels of the current 4-vertex. - `name` : name of the vertex - `level` : level in the diagram tree - `resetuid` : restart uid count from 1 @@ -26,7 +26,7 @@ Generate 4-vertex diagrams using Parquet Algorithm """ function vertex4(para::DiagPara, extK=[getK(para.totalLoopNum, 1), getK(para.totalLoopNum, 2), getK(para.totalLoopNum, 3)], - chan::AbstractVector=[PHr, PHEr, PPr, Alli], subdiagram=false; + subdiagram=false; channels::AbstractVector=[PHr, PHEr, PPr, Alli], level=1, name=:none, resetuid=false, # phi_toplevel=ParquetBlocks().phi, ppi_toplevel=ParquetBlocks().ppi, Γ4_toplevel=ParquetBlocks().Γ4, blocks::ParquetBlocks=ParquetBlocks(), @@ -70,10 +70,11 @@ function vertex4(para::DiagPara, end bareVer4(ver4df, para, legK, permutation) else # loopNum>0 - for c in chan + for c in channels if c == Alli + # continue if 3 ≤ loopNum ≤ 4 - addAlli!(ver4df, para, extK) + addAlli!(ver4df, para, legK) else continue end @@ -95,9 +96,8 @@ function vertex4(para::DiagPara, end # # TODO: add envolpe diagrams end - # println(ver4df) ver4df = merge_vertex4(para, ver4df, name, legK) - # @assert all(x -> x[1] == para.firstTauIdx, ver4df.extT) "not all extT[1] are equal to the first Tau index $(para.firstTauIdx)! $ver4df" + @assert all(x -> x[1] == para.firstTauIdx, ver4df.extT) "not all extT[1] are equal to the first Tau index $(para.firstTauIdx)! $ver4df" return ver4df end @@ -115,20 +115,12 @@ function merge_vertex4(para, ver4df, name, legK) return ver4df end -function addAlli!(ver4df::DataFrame, para::DiagPara, extK::Vector{Vector{Float64}}) +function addAlli!(ver4df::DataFrame, para::DiagPara, legK::Vector{Vector{Float64}}) dict_graphs = get_ver4I() graphvec = dict_graphs[para.innerLoopNum] - - legK = [k[1:para.totalLoopNum] for k in extK[1:3]] - update_extK!(graphvec, legK) - push!(legK, legK[1] + legK[3] - legK[2]) + graphvec = update_extKT(graphvec, para, legK) for ver4diag in graphvec Id = ver4diag.properties - for node in PostOrderDFS(ver4diag) - if !isleaf(node) - node.properties = Ver4Id(para, node.properties.response, node.properties.type, k=legK, t=Id.extT, chan=Alli) - end - end push!(ver4df, (response=Id.response, type=Id.type, extT=Id.extT, diagram=ver4diag)) end end @@ -155,6 +147,9 @@ function bubble!(ver4df::DataFrame, para::DiagPara, legK, chan::TwoBodyChannel, LfirstTauIdx, G0firstTauIdx, RfirstTauIdx, GxfirstTauIdx = idx @assert maxTau == maxVer4TauIdx(para) "Partition $partition with tauNum configuration $idx. maxTau = $maxTau, yet $(maxTauIdx(para)) is expected!" + # println(partition) + # println(idx, " ", maxTau) + lPara = reconstruct(para, type=Ver4Diag, innerLoopNum=oL, firstLoopIdx=LfirstLoopIdx, firstTauIdx=LfirstTauIdx) rPara = reconstruct(para, type=Ver4Diag, innerLoopNum=oR, firstLoopIdx=RfirstLoopIdx, firstTauIdx=RfirstTauIdx) gxPara = reconstruct(para, type=GreenDiag, innerLoopNum=oGx, firstLoopIdx=GxfirstLoopIdx, firstTauIdx=GxfirstTauIdx) @@ -175,9 +170,9 @@ function bubble!(ver4df::DataFrame, para::DiagPara, legK, chan::TwoBodyChannel, LLegK, K, RLegK, Kx = legBasis(chan, legK, LoopIdx) # println(K, ", ", Kx) - Lver = vertex4(lPara, LLegK, Γi, true; level=level + 1, name=:Γi, blocks=blocks) + Lver = vertex4(lPara, LLegK, true; channels=Γi, level=level + 1, name=:Γi, blocks=blocks) isempty(Lver) && return - Rver = vertex4(rPara, RLegK, Γf, true; level=level + 1, name=:Γf, blocks=blocks) + Rver = vertex4(rPara, RLegK, true; channels=Γf, level=level + 1, name=:Γf, blocks=blocks) isempty(Rver) && return ver8 = Dict{Any,Any}() diff --git a/test/front_end.jl b/test/front_end.jl index 977af7b9..d370ad66 100644 --- a/test/front_end.jl +++ b/test/front_end.jl @@ -478,7 +478,7 @@ end varT = [rand() for i in 1:para.totalTauNum] #################### DiagTree #################################### - diags = Parquet.vertex4(para, legK, chan) + diags = Parquet.vertex4(para, legK, channels=chan) diags = mergeby(diags, :response) # DiagTreeNew.plot_tree(diags[1]) # DiagTreeNew.plot_tree(diags[2])