update vertex4 in GV

src/FeynmanDiagram.jl

@@ -62,7 +62,7 @@ export Parquet
 include("frontend/GV.jl")
 using .GV
 export GV
-export diagdictGV, diagdict_parquet, leafstates, leafstates_diagtree
+export diagdictGV, diagdict_parquet, diagdict_parquetVer4, leafstates, leafstates_diagtree
 
 # include("strong_coupling_expansion_builder/strong_coupling_expansion.jl")
 # using .SCE

src/computational_graph/eval.jl

@@ -22,7 +22,8 @@ function eval!(g::Graph{F,W}, leafmap::Dict{Int,Int}=Dict{Int,Int}(), leaf::Vect
             if !inherit
                 if isempty(leafmap)
                     if randseed < 0
-                        node.weight = 1.0
+                        # node.weight = 1.0
+                        node.weight = -3.0
                     else
                         node.weight = rand()
                     end

src/frontend/GV.jl

@@ -17,14 +17,15 @@ import ..FrontEnds: NoBubble  # true to remove all bubble subdiagram
 import ..FrontEnds: Proper  #ver4, ver3, and polarization diagrams may require to be irreducible along the transfer momentum/frequency
 import ..FrontEnds: Response, Composite, ChargeCharge, SpinSpin, UpUp, UpDown
 import ..FrontEnds: AnalyticProperty, Instant, Dynamic
+import ..FrontEnds: TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan
 import ..FrontEnds: DiagramId, Ver4Id, Ver3Id, GreenId, SigmaId, PolarId, BareGreenId, BareInteractionId
 
 
 using AbstractTrees
 
 import ..Utility: taylorexpansion!
 
-export eachorder_diag, diagdictGV, diagdict_parquet, leafstates, leafstates_diagtree
+export eachorder_diag, diagdictGV, diagdict_parquet, diagdict_parquetVer4, leafstates, leafstates_diagtree
 
 include("GV_diagrams/readfile.jl")
 
@@ -182,30 +183,54 @@ A tuple `(diagrams, fermi_labelProd, bose_labelProd, extT_labels)` where
 """
 function eachorder_diag(type::Symbol, order::Int, GOrder::Int=0, VerOrder::Int=0;
     labelProd::Union{Nothing,LabelProduct}=nothing, spinPolarPara::Float64=0.0, tau_labels::Union{Nothing,Vector{Int}}=nothing)
-    diagtype = :polar
     if type == :spinPolar
         filename = string(@__DIR__, "/GV_diagrams/groups_spin/Polar$(order)_$(VerOrder)_$(GOrder).diag")
     elseif type == :chargePolar
         filename = string(@__DIR__, "/GV_diagrams/groups_charge/Polar$(order)_$(VerOrder)_$(GOrder).diag")
     elseif type == :sigma
-        diagtype = type
         filename = string(@__DIR__, "/GV_diagrams/groups_sigma/Sigma$(order)_$(VerOrder)_$(GOrder).diag")
     elseif type == :green
-        diagtype = type
         filename = string(@__DIR__, "/GV_diagrams/groups_green/Green$(order)_$(VerOrder)_$(GOrder).diag")
     elseif type == :freeEnergy
-        diagtype = type
         filename = string(@__DIR__, "/GV_diagrams/groups_free_energy/FreeEnergy$(order)_$(VerOrder)_$(GOrder).diag")
     end
     # println("Reading ", filename)
 
     if isnothing(labelProd)
-        return read_diagrams(filename; tau_labels=tau_labels, diagType=diagtype, spinPolarPara=spinPolarPara)
+        return read_diagrams(filename; tau_labels=tau_labels, diagType=type, spinPolarPara=spinPolarPara)
     else
-        return read_diagrams(filename; labelProd=labelProd, diagType=diagtype, spinPolarPara=spinPolarPara)
+        return read_diagrams(filename; labelProd=labelProd, diagType=type, spinPolarPara=spinPolarPara)
     end
 end
 
+function eachorder_diags(type::Symbol, order::Int)
+    innerLoopNum = order
+    if type == :spinPolar
+        filename = string(@__DIR__, "/GV_diagrams/groups_spin/Polar$(order)_0_0.diag")
+    elseif type == :chargePolar
+        filename = string(@__DIR__, "/GV_diagrams/groups_charge/Polar$(order)_0_0.diag")
+    elseif type == :sigma
+        filename = string(@__DIR__, "/GV_diagrams/groups_sigma/Sigma$(order)_0_0.diag")
+    elseif type == :green
+        innerLoopNum = order - 1
+        filename = string(@__DIR__, "/GV_diagrams/groups_green/Green$(order)_0_0.diag")
+    elseif type == :freeEnergy
+        filename = string(@__DIR__, "/GV_diagrams/groups_free_energy/FreeEnergy$(order)_0_0.diag")
+    elseif type == :vertex4
+        filename = string(@__DIR__, "/GV_diagrams/groups_vertex4/Vertex4$(order)_0_0.diag")
+    elseif type == :vertex4I
+        filename = string(@__DIR__, "/GV_diagrams/groups_vertex4/Vertex4I$(order)_0_0.diag")
+    end
+    # println("Reading ", filename)
+    para = DiagPara(type=_diagtype(type),
+        innerLoopNum=innerLoopNum,
+        hasTau=true,
+        interaction=[Interaction(ChargeCharge, Instant)],
+        filter=[NoHartree]
+    )
+    return read_diagrams(filename, para)
+end
+
 """
     function diagdict_parquet(type::Symbol, MaxOrder::Int, has_counterterm::Bool=true; MinOrder::Int=1,
         spinPolarPara::Float64=0.0, isDynamic=false, filter=[NoHartree])
@@ -392,12 +417,57 @@ function diagdict_parquet(type::Symbol, gkeys::Vector{Tuple{Int,Int,Int}}, extra
     return dict_graphs
 end
 
+function diagdict_parquetVer4(gkeys::Vector{Tuple{Int,Int,Int}};
+    spinPolarPara::Float64=0.0, isDynamic=false, filter=[NoHartree], transferLoop=nothing)
+    # spinPolarPara::Float64=0.0, isDynamic=false, channel=[PHr, PHEr, PPr], filter=[NoHartree])
+
+    spin = 2.0 / (spinPolarPara + 1)
+    dict_graphs = Dict{Tuple{Int,Int,Int},Tuple{Vector{Graph},Vector{Vector{Int}}}}()
+
+    # KinL, KoutL, KinR = zeros(16), zeros(16), zeros(16)
+    # KinL[1], KoutL[2], KinR[3] = 1.0, 1.0, 1.0
+
+    MinOrder = minimum([p[1] for p in gkeys])
+    MaxOrder = maximum([p[1] for p in gkeys])
+    for order in MinOrder:MaxOrder
+        Taylor.set_variables("x y"; orders=[MaxOrder - order, MaxOrder - order])
+        para = diagPara(Ver4Diag, isDynamic, spin, order, filter, transferLoop)
+        ver4df = Parquet.vertex4(para)
+
+        # Append fully irreducible Vertex4 diagrams
+        if 3 ≤ order ≤ 4
+            ver4I, extT_labels = eachorder_diags(:vertex4I, order)
+            responses = repeat([ChargeCharge], length(ver4I))
+            types = repeat([Dynamic], length(ver4I))
+            append!(ver4df, (response=responses, type=types, extT=extT_labels, diagram=ver4I, hash=IR.id.(ver4I)))
+        end
+        diags, extT = ver4df.diagram, ver4df.extT
+        propagator_var = Dict(BareGreenId => [true, false], BareInteractionId => [false, true]) # Specify variable dependence of fermi (first element) and bose (second element) particles.
+        taylor_vec, taylormap = taylorexpansion!(diags, propagator_var)
+
+        for t in taylor_vec
+            for (o, graph) in t.coeffs
+                key = (order, o...)
+                key ∉ gkeys && continue
+                if haskey(dict_graphs, key)
+                    push!(dict_graphs[key][1], graph)
+                else
+                    dict_graphs[key] = ([graph,], collect.(extT))
+                end
+            end
+        end
+    end
+
+    for gvec in values(dict_graphs)
+        IR.optimize!(gvec[1])
+    end
+    return dict_graphs
+end
+
 function diagPara(type, isDynamic::Bool, spin, order, filter, transferLoop=nothing)
     inter = [Interaction(ChargeCharge, isDynamic ? [Instant, Dynamic] : [Instant,]),]  #instant charge-charge interaction
     if type == VacuumDiag
         innerLoopNum = order + 1
-    elseif type == Ver4Diag
-        innerLoopNum = order - 1
     else
         innerLoopNum = order
     end
@@ -426,17 +496,19 @@ end
 
 function _diagtype(type::Symbol)
     if type == :freeEnergy
-        diagtype = VacuumDiag
+        return VacuumDiag
     elseif type == :sigma
-        diagtype = SigmaDiag
+        return SigmaDiag
     elseif type == :green
-        diagtype = GreenDiag
+        return GreenDiag
     elseif type == :chargePolar
-        diagtype = PolarDiag
+        return PolarDiag
     elseif type == :vertex3
-        diagtype = Ver3Diag
-    elseif type == :vertex4
-        diagtype = Ver4Diag
+        return Ver3Diag
+    elseif type in [:vertex4, :vertex4I]
+        return Ver4Diag
+    else
+        error("$type is not implemented")
     end
 end
 

src/frontend/GV_diagrams/diagram.py

@@ -4,7 +4,7 @@
 from numpy.linalg import matrix_rank
 import numpy as np
 import random
-
+from copy import deepcopy
 
 class diagram:
     """a feynman diagram class"""
@@ -72,7 +72,27 @@ def __init__(self, order):
 
     def GetPermu(self):
         return tuple(self.Permutation)
-
+
+    def SwapTwoVertexPairs(self, i, j, k, l):
+        # print permutation, i, j
+        Assert(i-j == k-l, "Error input two vertex pairs!")
+
+        permutation = list(self.Permutation)
+        ip, kp = permutation.index(i), permutation.index(k)
+        jp, lp = permutation.index(j), permutation.index(l)
+        permutation[ip] = k
+        permutation[kp] = i
+        permutation[jp] = l
+        permutation[lp] = j
+        permutation[i], permutation[k] = permutation[k], permutation[i]
+        permutation[j], permutation[l] = permutation[l], permutation[j]
+        # print "after", permutation
+        self.Permutation = tuple(permutation)
+
+        loopBasis = deepcopy(self.LoopBasis)
+        self.LoopBasis[:, i], self.LoopBasis[:, k] =  loopBasis[:, k], loopBasis[:, i]
+        self.LoopBasis[:, j], self.LoopBasis[:, l] =  loopBasis[:, l], loopBasis[:, j]
+
     # def GetSimpleInteractionPairs(self):
     #     if self.Type == "FreeEnergy" or self.Type == "Ver4":
     #         return [(2*i, 2*i+1) for i in range(self.Ver4Num)]
@@ -208,6 +228,22 @@ def SwapTwoVertex(permutation, i, j):
     # print "after", permutation
     return tuple(permutation)
 
+# def SwapTwoVertexPairs(permutation, i, j, k, l):
+#     # print permutation, i, j
+#     Assert(i-j == k-l, "Error input two vertex pairs!")
+
+#     permutation = list(permutation)
+#     ip, kp = permutation.index(i), permutation.index(k)
+#     jp, lp = permutation.index(j), permutation.index(l)
+#     permutation[ip] = k
+#     permutation[kp] = i
+#     permutation[jp] = l
+#     permutation[lp] = j
+#     permutation[i], permutation[k] = permutation[k], permutation[i]
+#     permutation[j], permutation[l] = permutation[l], permutation[j]
+#     # print "after", permutation
+#     return tuple(permutation)
+
 
 def Direct2Exchange(permutation, i, j):
     """change a direction interaction (i-j) to exchange interaction, or the reversed"""