utilize Hugenholtz in GV

src/computational_graph/graph.jl

@@ -172,7 +172,7 @@ end
 - `c1`  first scalar multiple
 - `c2`  second scalar multiple
 """
-function linear_combination(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1)) where {F,W}
+function linear_combination(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1); properties=nothing) where {F,W}
     if length(g1.orders) > length(g2.orders)
         g2.orders = [orders(g2); zeros(Int, length(g1.orders) - length(g2.orders))]
     else
@@ -195,9 +195,9 @@ function linear_combination(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1)) wh
     end
 
     if subgraphs[1].id == subgraphs[2].id
-        g = Graph([subgraphs[1]]; subgraph_factors=[sum(subgraph_factors)], operator=Sum(), orders=orders(g1), ftype=F, wtype=W)
+        g = Graph([subgraphs[1]]; subgraph_factors=[sum(subgraph_factors)], operator=Sum(), orders=orders(g1), ftype=F, wtype=W, properties=properties)
     else
-        g = Graph(subgraphs; subgraph_factors=subgraph_factors, operator=Sum(), orders=orders(g1), ftype=F, wtype=W)
+        g = Graph(subgraphs; subgraph_factors=subgraph_factors, operator=Sum(), orders=orders(g1), ftype=F, wtype=W, properties=properties)
     end
 
     return g
@@ -222,7 +222,7 @@ where duplicate graphs in the input `graphs` are combined by summing their assoc
 # Example:
     Given graphs `g1`, `g2`, `g1` and constants `c1`, `c2`, `c3`, the function computes `(c1+c3)*g1 + c2*g2`.
 """
-function linear_combination(graphs::Vector{Graph{F,W}}, constants::AbstractVector=ones(F, length(graphs))) where {F,W}
+function linear_combination(graphs::Vector{Graph{F,W}}, constants::AbstractVector=ones(F, length(graphs)); properties=nothing) where {F,W}
     maxlen_orders = maximum(length.(orders.(graphs)))
     for g in graphs
         g.orders = [orders(g); zeros(Int, maxlen_orders - length(orders(g)))]
@@ -254,7 +254,7 @@ function linear_combination(graphs::Vector{Graph{F,W}}, constants::AbstractVecto
     if isempty(unique_graphs)
         return nothing
     end
-    g = Graph(unique_graphs; subgraph_factors=unique_factors, operator=Sum(), orders=orders(graphs[1]), ftype=F, wtype=W)
+    g = Graph(unique_graphs; subgraph_factors=unique_factors, operator=Sum(), orders=orders(graphs[1]), ftype=F, wtype=W, properties=properties)
     return g
 end
 
@@ -298,7 +298,7 @@ end
 - `c1`:  first scalar multiple (defaults to 1).
 - `c2`:  second scalar multiple (defaults to 1).
 """
-function multi_product(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1)) where {F,W}
+function multi_product(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1); properties=nothing) where {F,W}
     # @assert orders(g1) == orders(g2) "g1 and g2 have different orders."
     f1 = typeof(c1) == F ? c1 : F(c1)
     f2 = typeof(c2) == F ? c2 : F(c2)
@@ -315,14 +315,14 @@ function multi_product(g1::Graph{F,W}, g2::Graph{F,W}, c1=F(1), c2=F(1)) where {
     end
 
     if subgraphs[1].id == subgraphs[2].id
-        g = Graph([subgraphs[1]]; subgraph_factors=[prod(subgraph_factors)], operator=Power(2), orders=2 * orders(g1), ftype=F, wtype=W)
+        g = Graph([subgraphs[1]]; subgraph_factors=[prod(subgraph_factors)], operator=Power(2), orders=2 * orders(g1), ftype=F, wtype=W, properties=properties)
     else
         if length(g1.orders) > length(g2.orders)
             g2.orders = [orders(g2); zeros(Int, length(g1.orders) - length(g2.orders))]
         else
             g1.orders = [orders(g1); zeros(Int, length(g2.orders) - length(g1.orders))]
         end
-        g = Graph(subgraphs; subgraph_factors=subgraph_factors, operator=Prod(), orders=orders(g1) + orders(g2), ftype=F, wtype=W)
+        g = Graph(subgraphs; subgraph_factors=subgraph_factors, operator=Prod(), orders=orders(g1) + orders(g2), ftype=F, wtype=W, properties=properties)
     end
     return g
 end
@@ -344,7 +344,7 @@ Returns:
 # Example:
     Given graphs `g1`, `g2`, `g1` and constants `c1`, `c2`, `c3`, the function computes `(c1*c3)*(g1)^2 * c2*g2`.
 """
-function multi_product(graphs::Vector{Graph{F,W}}, constants::AbstractVector=ones(F, length(graphs))) where {F,W}
+function multi_product(graphs::Vector{Graph{F,W}}, constants::AbstractVector=ones(F, length(graphs)); properties=nothing) where {F,W}
     # @assert alleq(orders.(graphs)) "Graphs do not all have the same order."
     g1 = graphs[1]
     subgraphs = graphs
@@ -382,7 +382,7 @@ function multi_product(graphs::Vector{Graph{F,W}}, constants::AbstractVector=one
     end
 
     if length(unique_factors) == 1
-        g = Graph(unique_graphs; subgraph_factors=unique_factors, operator=Power(repeated_counts[1]), orders=g_orders, ftype=F, wtype=W)
+        g = Graph(unique_graphs; subgraph_factors=unique_factors, operator=Power(repeated_counts[1]), orders=g_orders, ftype=F, wtype=W, properties=properties)
     else
         subgraphs = Vector{Graph{F,W}}()
         for (idx, g) in enumerate(unique_graphs)
@@ -392,7 +392,7 @@ function multi_product(graphs::Vector{Graph{F,W}}, constants::AbstractVector=one
                 push!(subgraphs, Graph([g], operator=Power(repeated_counts[idx]), orders=orders(g1) * repeated_counts[idx], ftype=F, wtype=W))
             end
         end
-        g = Graph(subgraphs; subgraph_factors=unique_factors, operator=Prod(), orders=g_orders, ftype=F, wtype=W)
+        g = Graph(subgraphs; subgraph_factors=unique_factors, operator=Prod(), orders=g_orders, ftype=F, wtype=W, properties=properties)
     end
     return g
 end

src/frontend/GV_diagrams/readfile.jl

@@ -230,12 +230,13 @@ function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0
     para = (2, innerLoopNum)
 
     # Combine and merge all diagrams with the same external-tau labels
+    num_diags = length(diagrams)
     extT_labels = Vector{NTuple{4,Int}}()
     channel_vector = Vector{TwoBodyChannel}()
     DiEx = Int[]
-    sizehint!(extT_labels, length(diagrams))
-    sizehint!(channel_vector, length(diagrams))
-    sizehint!(DiEx, length(diagrams))
+    sizehint!(extT_labels, num_diags)
+    sizehint!(channel_vector, num_diags)
+    sizehint!(DiEx, num_diags)
     for prop in IR.properties.(diagrams)
         push!(extT_labels, prop.extT)
         push!(channel_vector, prop.channel)
@@ -254,15 +255,11 @@ function read_vertex4diagrams(filename::AbstractString; spinPolarPara::Float64=0
         keyEx = (extT, channel, 1)
         gId_Di = gr[keyDi][1].properties
 
-        gDi = IR.linear_combination(gr[keyDi])
-        gEx = IR.linear_combination(gr[keyEx])
-        gDi.properties = gId_Di
-        gEx.properties = gId_Di
+        gud = IR.linear_combination(gr[keyDi], properties=gId_Di)     # Direct = UpDown
+        gEx = IR.linear_combination(gr[keyEx], properties=gr[keyEx][1].properties)
 
         guuId = Ver4Id(para, UpUp, gId_Di.type, k=gId_Di.extK, t=gId_Di.extT, chan=gId_Di.channel)
-        gudId = Ver4Id(para, UpDown, gId_Di.type, k=gId_Di.extK, t=gId_Di.extT, chan=gId_Di.channel)
-        guu = Graph([gDi, gEx], properties=guuId)
-        gud = Graph([gDi], properties=gudId)
+        guu = Graph([gud, gEx], properties=guuId)
         append!(graphvec, [guu, gud])
         append!(extTs, [extT, extT])
         append!(responses, [UpUp, UpDown])
@@ -359,8 +356,19 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
             end
         end
     end
-    spinfactors_existed = Float64[]
 
+    greens = Graph{_dtype.factor,_dtype.weight}[]
+    # create all fermionic propagators
+    for (ind1, ind2) in enumerate(permutation)
+        opGType[ind1] == -2 && continue
+        diagid = BareGreenId(k=currentBasis[ind1, :], t=[tau_labels[ind1], tau_labels[ind2]])
+        push!(greens, Graph([]; properties=diagid))
+    end
+    fermi_greenProd = Graph(greens, operator=IR.Prod())
+
+    interactions_Di = Graph{_dtype.factor,_dtype.weight}[]
+    interactions_Ex = Graph{_dtype.factor,_dtype.weight}[]
+    spinfactors_existed = Float64[]
     # println("##### $permutation  $ver4Legs")
     for (iex, spinFactor) in enumerate(spinFactors)
         # create permutation and ver4Legs for each Feynman diagram from a Hugenholtz diagram
@@ -374,13 +382,6 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
         extIndex[1] = permu[1]
         extIndex[3] = permu[2]
 
-        # create all fermionic operators
-        for (ind1, ind2) in enumerate(permu)
-            opGType[ind1] == -2 && continue
-            diagid = BareGreenId(k=currentBasis[ind1, :], t=[tau_labels[ind1], tau_labels[ind2]])
-            push!(leafs, Graph([]; properties=diagid))
-        end
-
         # create all bosionic operators (relevant to interaction lines)
         for verLeg in ver4Legs_ex
             ind1, ind2 = verLeg[2] - offset, verLeg[4] - offset
@@ -391,17 +392,24 @@ function read_one_vertex4diagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int,
             push!(leafs, Graph([]; properties=diagid))
         end
 
-        para = (DiEx[iex], innerLoopNum)
-        if isDynamic
-            diagid = Ver4Id(para, ChargeCharge, Dynamic, k=extK, t=tau_labels[extIndex], chan=channel)
+        if DiEx[iex] == 0
+            push!(interactions_Di, Graph(leafs, operator=IR.Prod(), factor=spinFactor * symfactor))
         else
-            diagid = Ver4Id(para, ChargeCharge, Instant, k=extK, t=tau_labels[extIndex], chan=channel)
+            push!(interactions_Ex, Graph(leafs, operator=IR.Prod(), factor=spinFactor * symfactor))
         end
+    end
 
-        push!(graphs, Graph(leafs, operator=IR.Prod(), properties=diagid, factor=spinFactor * symfactor))
+    diagid_Di = Ver4Id((0, innerLoopNum), UpDown, isDynamic ? Dynamic : Instant, k=extK, t=tau_labels[extIndex], chan=channel)
+    diagid_Ex = Ver4Id((1, innerLoopNum), ChargeCharge, isDynamic ? Dynamic : Instant, k=extK, t=tau_labels[extIndex], chan=channel)
+    if isempty(fermi_greenProd.subgraphs)
+        g_Di = Graph(interactions_Di, operator=IR.Sum(), properties=diagid_Di)
+        g_Ex = Graph(interactions_Ex, operator=IR.Sum(), properties=diagid_Ex)
+    else
+        g_Di = IR.multi_product(fermi_greenProd, Graph(interactions_Di, operator=IR.Sum()), properties=diagid_Di)
+        g_Ex = IR.multi_product(fermi_greenProd, Graph(interactions_Ex, operator=IR.Sum()), properties=diagid_Ex)
     end
 
-    return graphs
+    return g_Di, g_Ex
 end
 
 function read_onediagram!(io::IO, GNum::Int, verNum::Int, loopNum::Int, extIndex::Vector{Int},

src/frontend/parquet/operation.jl

@@ -155,10 +155,6 @@ function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector
     # 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]
@@ -173,6 +169,9 @@ function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector
             node.id = IR.uid()
             push!(visited, node.id)
             prop = IR.properties(node)
+            if !hasproperty(prop, :extK) || !hasproperty(prop, :extT)
+                continue
+            end
             K = prop.extK
             T = prop.extT
             if prop isa Ver4Id || prop isa Ver3Id
@@ -185,7 +184,7 @@ function update_extKT!(diags::Vector{Graph}, para::DiagPara, legK::Vector{Vector
                 else
                     node.properties = FrontEnds.reconstruct(prop, :para => para)
                 end
-            else
+            elseif prop isa PropagatorId || prop isa GreenId || prop isa SigmaId || prop isa PolarId
                 original_len_K = length(K)
                 if length(K) < len_extK
                     resize!(K, len_extK)