Ver3 new

example/diagram_compile.jl

@@ -40,7 +40,7 @@ elseif diagtype == :green || diagtype == :freeEnergy || diagtype == :chargePolar
         push!(partition, (o, sOrder, vOrder))
     end
     if diagtype == :freeEnergy
-        FeynGraphs = Diagram.diagram_freeE(para, partition, optimize_level=1)
+        FeynGraphs = Diagram.diagram_GV_freeE(para, partition, optimize_level=1)
     else
         FeynGraphs = Diagram.diagram_parquet_noresponse(diagtype, para, partition, optimize_level=1)
     end

example/diagram_compilepy.jl

@@ -0,0 +1,58 @@
+using ElectronLiquid
+using FeynmanDiagram
+
+@inline function totalMomNum(order::Int, diagtype::Symbol)
+    if diagtype == :vertex3
+        return order + 2
+    elseif diagtype == :vertex4
+        return order + 3
+    else
+        return order + 1
+    end
+end
+
+diagtype = :sigma # :sigma, :vertex3, :vertex4, :freeEnergy, :green, :chargePolar
+order = 3
+# filter = [Parquet.NoHartree, Parquet.Proper]
+filter = [Parquet.NoHartree,]
+KinL, KoutL, KinR = zeros(16), zeros(16), zeros(16)
+KinL[1], KoutL[2], KinR[3] = 1.0, 1.0, 1.0
+
+para = UEG.ParaMC(rs=1.0, beta=25, order=order, isDynamic=false)
+
+if diagtype == :chargePolar || diagtype == :sigma
+    _partition = UEG.partition(order)
+else
+    _partition = UEG.partition(order, offset=0)
+end
+
+if diagtype == :vertex4 || diagtype == :vertex3
+    partition = Vector{NTuple{3,Int}}()
+    for (o, sOrder, vOrder) in _partition
+        o == 0 && sOrder > 0 && continue
+        push!(partition, (o, sOrder, vOrder))
+    end
+
+    FeynGraphs = Diagram.diagram_parquet_response(diagtype, para, partition, optimize_level=1, filter=filter, transferLoop=KinL - KoutL)
+elseif diagtype == :green || diagtype == :freeEnergy || diagtype == :chargePolar
+    partition = Vector{NTuple{3,Int}}()
+    for (o, sOrder, vOrder) in _partition
+        o == 0 && vOrder > 0 && continue
+        push!(partition, (o, sOrder, vOrder))
+    end
+    if diagtype == :freeEnergy
+        FeynGraphs = Diagram.diagram_GV_freeE(para, partition, optimize_level=1)
+    else
+        FeynGraphs = Diagram.diagram_parquet_noresponse(diagtype, para, partition, optimize_level=1)
+    end
+else
+    partition = _partition
+    FeynGraphs = Diagram.diagram_parquet_noresponse(diagtype, para, partition, optimize_level=1)
+end
+
+# compile Python
+for p in FeynGraphs[1]
+    Compilers.compile_Python(FeynGraphs[3][p], "func_$(diagtype)_o$(p[1])$(p[2])$(p[3]).py")
+end
+
+println("Compile finished.")

example/ver4/test_PP.jl

@@ -82,7 +82,7 @@ end
     # diagram = Ver4.diagramParquet(para, partition; channel=[PHr, PHEr, PPr,], filter=[NoHartree, NoBubble])
     # data, result = Ver4.one_angle_averaged_ParquetAD(paras, diagram; neval=neval, print=-1, seed=seed)
 
-    diagram = Ver4.diagramParquet_load(para, partition; filter=[NoHartree, NoBubble])
+    diagram = Ver4.diagramParquet_load(para, partition; filter=[FeynmanDiagram.FrontEnds.NoHartree, FeynmanDiagram.FrontEnds.NoBubble])
     data, result = Ver4.one_angle_averaged_ParquetAD_Clib(paras, diagram; neval=neval, print=-1, seed=seed)
     println(data[(1, 0, 0)], data[(2, 0, 0)])
     # obs2 = data[p]

src/ElectronLiquid.jl

@@ -67,9 +67,9 @@ include("./vertex4/vertex4.jl")
 using .Ver4
 export Ver4
 
-# include("./vertex3/vertex3.jl")
-# using .Ver3
-# export Ver3
+include("./vertex3/vertex3_new.jl")
+using .Ver3
+export Ver3
 
 include("./freeEnergy/freeEnergy.jl")
 using .FreeEnergy

src/common/counterterm.jl

@@ -117,6 +117,7 @@ function chemicalpotential_renormalization(order, data, δμ; offset::Int=0)
     @assert length(δμ) + 1 >= order
     data = mergeInteraction(data)
     d = data
+    println(keys(d))
     # println("size: ", size(d[(1, 0)]))
     # z = Vector{eltype(values(d))}(undef, order)
     sample = collect(values(d))[1]
@@ -328,8 +329,8 @@ function _inverse(z::AbstractVector{T}) where {T}
         zi[5] = -z[1] .^ 5 + 4z[1] .^ 3 .* z[2] - 3z[1] .* z[2] .^ 2 - 3z[1] .^ 2 .* z[3] + 2z[2] .* z[3] + 2z[1] .* z[4] - z[5]
     end
     if order >= 6
-        zi[6] = z[1] .^ 6 - 5z[1] .^ 4 .* z[2] + 6z[1] .^ 2 .* z[2] .^ 2 + 4z[3] .* z[1] .^ 3
-        -3z[1] .^ 2 .* z[4] - 6z[1] .* z[2] .* z[3] - z[2] .^ 3 + z[3] .^ 2 + 2z[2] .* z[4] + 2z[1] .* z[5] - z[6]
+        zi[6] = z[1] .^ 6 - 5z[1] .^ 4 .* z[2] + 6z[1] .^ 2 .* z[2] .^ 2 + 4z[3] .* z[1] .^ 3 -
+                3z[1] .^ 2 .* z[4] - 6z[1] .* z[2] .* z[3] - z[2] .^ 3 + z[3] .^ 2 + 2z[2] .* z[4] + 2z[1] .* z[5] - z[6]
     end
     if order >= 7
         error("order must be <= 5")
@@ -533,6 +534,7 @@ function getSigma(df::DataFrame, paraid::Dict, order::Int)
 
     mu = Dict()
     for P in _partition
+        # P[3] >0 && continue
         v = filter(r -> r["partition"] == compactPartition(P), df)[1, "μ"]
         err = filter(r -> r["partition"] == compactPartition(P), df)[1, "μ.err"]
         mu[P] = measurement(v, err)
@@ -543,6 +545,7 @@ function getSigma(df::DataFrame, paraid::Dict, order::Int)
 
     sw = Dict()
     for P in _partition
+        # P[3] >0 && continue
         v = filter(r -> r["partition"] == compactPartition(P), df)[1, "Σw"]
         err = filter(r -> r["partition"] == compactPartition(P), df)[1, "Σw.err"]
         sw[P] = measurement(v, err)

src/common/eval.jl

@@ -33,6 +33,10 @@ end
 function interaction_derive(τ1, τ2, K, p::ParaMC, idorder; idtype=Instant, tau_num=1, isLayered=false)
     dim, e0, ϵ0, mass2 = p.dim, p.e0, p.ϵ0, p.mass2
     qd = sqrt(dot(K, K))
+    if p.isDynamic == true
+        # for dynamic case Instant also has 2 tau
+        tau_num = 2
+    end
     if idorder[2] == 0
         if idtype == Instant
             if tau_num == 1

src/diagram/generate.jl

@@ -124,7 +124,7 @@ function diagdict_parquet(diagtype::Union{DiagramType,Symbol}, _partition::Vecto
     return dict_graphs
 end
 
-function diagram_freeE(paramc::ParaMC, _partition::Vector{T},
+function diagram_GV_freeE(paramc::ParaMC, _partition::Vector{T},
     leaf_dep_funcs::Vector{Function}=[pr -> pr isa FrontEnds.BareGreenId, pr -> pr isa FrontEnds.BareInteractionId];
     filter=[NoHartree], optimize_level=0
 ) where {T}
@@ -166,6 +166,50 @@ function diagram_freeE(paramc::ParaMC, _partition::Vector{T},
     return (partition, diagpara, dict_graphs)
 end
 
+function diagram_GV_noresponse(diagtype::Symbol, paramc::ParaMC, _partition::Vector{T},
+    leaf_dep_funcs::Vector{Function}=[pr -> pr isa FrontEnds.BareGreenId, pr -> pr isa FrontEnds.BareInteractionId];
+    filter=[NoHartree], optimize_level=0
+) where {T}
+    deriv_num = length(leaf_dep_funcs)
+    @assert length(_partition[1]) == deriv_num + 1 "partition should have $deriv_num+1 entries"
+
+    diagpara = Vector{DiagPara}()
+    MaxOrder = maximum([p[1] for p in _partition])
+    MinOrder = minimum([p[1] for p in _partition])
+    dict_graphs = Dict{NTuple{deriv_num + 1,Int},Vector{Graph}}()
+    spinPolarPara = 2.0 / paramc.spin - 1
+
+    for order in MinOrder:MaxOrder
+        partition_ind = findall(p -> p[1] == order, _partition)
+        partition_order = [_partition[i] for i in partition_ind]
+
+        diagrams = GV.diagsGV(diagtype, order, spinPolarPara=spinPolarPara, filter=filter)
+        optimize!(diagrams, level=optimize_level)
+        optimize!(diagrams, level=optimize_level)
+
+        renormalization_orders = Int[]
+        for i in 1:deriv_num
+            push!(renormalization_orders, maximum([p[i+1] for p in partition_order]))
+        end
+        # renormalization_orders = [Max_GD_o, Max_ID_o, extra_deriv_orders...]
+        dict_graph_order = taylorAD(diagrams, renormalization_orders, leaf_dep_funcs)
+        for key in keys(dict_graph_order)
+            p = (order, key...)
+            if p in _partition
+                dict_graphs[p] = dict_graph_order[key]
+            end
+        end
+    end
+
+    extT_labels = Vector{Vector{Int}}[]
+    partition = sort(collect(keys(dict_graphs)))
+    for p in partition
+        push!(diagpara, diagPara(_diagtype(diagtype), paramc.isDynamic, p[1], paramc.spin, filter))
+        push!(extT_labels, [collect(g.properties.extT) for g in dict_graphs[p]])
+    end
+    return (partition, diagpara, dict_graphs, extT_labels)
+end
+
 function _diagtype(type::Symbol)
     if type == :freeEnergy
         return Parquet.VacuumDiag

src/freeEnergy/freeEnergy.jl

@@ -27,7 +27,7 @@ function MC(para; neval=1e6, reweight_goal=nothing,
     isLayered2D=false, # whether to use the screened Coulomb interaction in 2D or not 
     filter=[NoHartree], optimize_level=1, verbose=-1
 )
-    diagram = Diagram.diagram_freeE(para, partition, filter=filter, optimize_level=optimize_level)
+    diagram = Diagram.diagram_GV_freeE(para, partition, filter=filter, optimize_level=optimize_level)
 
     partition = diagram[1]
     neighbor = UEG.neighbor(partition)

src/freeEnergy/parquetAD.jl

@@ -5,7 +5,6 @@ function integrand_parquetAD(idx, vars, config)
     leafval, leafType, leafOrders, leafτ_i, leafτ_o, leafMomIdx = leafStat
     dim, β, me, μ = para.dim, para.β, para.me, para.μ
     loopNum = config.dof[idx][1]
-    # is_zero_order = partition[idx][1] == 0 ? true : false
     is_zero_order = partition[idx] == (0, 0, 0) ? true : false
 
     FrontEnds.update(momLoopPool, varK.data[:, 1:maxMomNum])

src/freeEnergy/parquetAD_Clib.jl

@@ -6,7 +6,6 @@ function integrand_parquetAD_Clib(idx, vars, config)
 
     dim, β, me, μ = para.dim, para.β, para.me, para.μ
     loopNum = config.dof[idx][1]
-    # is_zero_order = partition[idx][1] == 0 ? true : false
     is_zero_order = partition[idx] == (0, 0, 0) ? true : false
 
     FrontEnds.update(momLoopPool, varK.data[:, 1:maxMomNum])

src/sigma/sigma.jl

@@ -91,6 +91,7 @@ function MC(para; kgrid=[para.kF,], ngrid=[0], neval=1e6, reweight_goal=nothing,
     # spinPolarPara::Float64=0.0, # spin-polarization parameter (n_up - n_down) / (n_up + n_down) ∈ [0,1]
     filename::Union{String,Nothing}=nothing, partition=UEG.partition(para.order),
     isLayered2D=false, # whether to use the screened Coulomb interaction in 2D or not 
+    diag_generator::Symbol=:parquet,
     filter=[NoHartree], extK=nothing, optimize_level=1, verbose=-1
 )
     kF = para.kF
@@ -112,7 +113,13 @@ function MC(para; kgrid=[para.kF,], ngrid=[0], neval=1e6, reweight_goal=nothing,
         push!(reweight_goal, 4.0)
     end
 
-    diagram = Diagram.diagram_parquet_noresponse(:sigma, para, partition, filter=filter, extK=extK, optimize_level=optimize_level)
+    if diag_generator == :Parquet
+        diagram = Diagram.diagram_parquet_noresponse(:sigma, para, partition, filter=filter, extK=extK, optimize_level=optimize_level)
+    elseif diag_generator == :GV
+        diagram = Diagram.diagram_GV_noresponse(:sigma, para, partition, filter=filter, optimize_level=optimize_level)
+    else
+        error("Unknown diag_generator: $diag_generator")
+    end
     sigma, result = Sigma.ParquetAD(para, diagram;
         isLayered2D=isLayered2D, print=verbose,
         neighbor=neighbor, reweight_goal=reweight_goal,