diff --git a/example/diagram_compile.jl b/example/diagram_compile.jl index 7b2f96c..300d8c2 100644 --- a/example/diagram_compile.jl +++ b/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 diff --git a/example/diagram_compilepy.jl b/example/diagram_compilepy.jl new file mode 100644 index 0000000..d0fb756 --- /dev/null +++ b/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.") \ No newline at end of file diff --git a/example/ver4/test_PP.jl b/example/ver4/test_PP.jl index be61692..70e9671 100644 --- a/example/ver4/test_PP.jl +++ b/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] diff --git a/src/ElectronLiquid.jl b/src/ElectronLiquid.jl index cc46559..ce38f62 100644 --- a/src/ElectronLiquid.jl +++ b/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 diff --git a/src/common/counterterm.jl b/src/common/counterterm.jl index 45353ac..a5f9cdf 100644 --- a/src/common/counterterm.jl +++ b/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) diff --git a/src/common/eval.jl b/src/common/eval.jl index bab03c3..7faeb62 100644 --- a/src/common/eval.jl +++ b/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 diff --git a/src/diagram/generate.jl b/src/diagram/generate.jl index b99e634..8bcf6f1 100644 --- a/src/diagram/generate.jl +++ b/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 diff --git a/src/freeEnergy/freeEnergy.jl b/src/freeEnergy/freeEnergy.jl index d94fe39..4c0df01 100644 --- a/src/freeEnergy/freeEnergy.jl +++ b/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) diff --git a/src/freeEnergy/parquetAD.jl b/src/freeEnergy/parquetAD.jl index b5393b9..36ffa21 100644 --- a/src/freeEnergy/parquetAD.jl +++ b/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]) diff --git a/src/freeEnergy/parquetAD_Clib.jl b/src/freeEnergy/parquetAD_Clib.jl index 87ca4a0..5c72f11 100644 --- a/src/freeEnergy/parquetAD_Clib.jl +++ b/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]) diff --git a/src/sigma/sigma.jl b/src/sigma/sigma.jl index d1a8b56..b9c4d24 100644 --- a/src/sigma/sigma.jl +++ b/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, diff --git a/src/vertex3/ver3_KW.jl b/src/vertex3/ver3_KW.jl new file mode 100644 index 0000000..1fc1c75 --- /dev/null +++ b/src/vertex3/ver3_KW.jl @@ -0,0 +1,236 @@ +@inline function phase_ver3(varT, extT, nqout, nkin, β) + # println(extT) + # tq, tkin, tkout = varT[extT[1]], varT[extT[2]], varT[extT[3]] + tq, tkin, tkout = varT[extT[1]], varT[extT[2]], varT[extT[3]] + wqout, wkin = π * (2nqout) / β, π * (2nkin + 1) / β + wkout = wkin - wqout + return exp(-1im * (tkin * wkin - tq * wqout - tkout * wkout)) +end + +# @inline function phase_ver3(varT, extT, n, β) +# # println(extT) +# return phase_ver3(varT, extT, n[1], n[2], β) +# end +@inline function interactionTauNum(type::AnalyticProperty) + if type == Instant + return 1 + else + return 2 + end +end + +""" +integrand of vertex3 +""" +function integrand_ver3KW(idx, var, config) + para, kin, nkin, qout, nqout = config.userdata[1:5] + maxMomNum, extT_labels, spin_conventions, leafStat, leaf_maps, momLoopPool, root, funcGraphs! = config.userdata[6:end] + + dim, β, me, μ = para.dim, para.β, para.me, para.μ + leafval, leafType, leafOrders, leafτ_i, leafτ_o, leafMomIdx = leafStat + varK, varT = var[1], var[2] + loopNum = config.dof[idx][1] + k1 = kin[var[3][1]] + n1 = nkin[var[4][1]] + q = qout[var[5][1]] + nq = nqout[var[6][1]] + + varK.data[:, 1] .= q + varK.data[1, 2] = k1 + + FrontEnds.update(momLoopPool, varK.data[:, 1:maxMomNum]) + + for (i, lftype) in enumerate(leafType[idx]) + if lftype == 0 + continue + elseif lftype == 1 #fermionic + τ = varT[leafτ_o[idx][i]] - varT[leafτ_i[idx][i]] + kq = FrontEnds.loop(momLoopPool, leafMomIdx[idx][i]) + ϵ = dot(kq, kq) / (2me) - μ + order = leafOrders[idx][i][1] + leafval[idx][i] = Propagator.green_derive(τ, ϵ, β, order) + elseif lftype == 2 #bosonic + diagid = leaf_maps[idx][i].properties + kq = FrontEnds.loop(momLoopPool, leafMomIdx[idx][i]) + τ2, τ1 = varT[leafτ_o[idx][i]], varT[leafτ_i[idx][i]] + idorder = leafOrders[idx][i] + leafval[idx][i] = Propagator.interaction_derive(τ1, τ2, kq, para, idorder; idtype=diagid.type, tau_num=interactionTauNum(diagid.type)) + # leafval[idx][i] = Propagator.interaction_derive(τ1, τ2, kq, para, idorder; idtype=Instant, tau_num=1) + else + error("this leaftype $lftype not implemented!") + end + end + + # factor = para.NF / (2π)^(dim * (loopNum)) + factor = 1.0 / (2π)^(dim * (loopNum)) + graphfuncs! = funcGraphs![idx] + graphfuncs!(root, leafval[idx]) + wuu = zero(ComplexF64) + wud = zero(ComplexF64) + for ri in 1:length(extT_labels[idx]) + if spin_conventions[idx][ri] == UpUp + wuu += root[ri] * phase_ver3(varT, extT_labels[idx][ri], nq, n1, β) + elseif spin_conventions[idx][ri] == UpDown + wud += root[ri] * phase_ver3(varT, extT_labels[idx][ri], nq, n1, β) + end + end + + return Weight(wuu * factor, wud * factor) + + # return Weight(1.0, 1.0) +end + +function measure_ver3KW(idx, var, obs, weight, config) + KINidx = var[3][1] + NKINidx = var[4][1] + QOUTidx = var[5][1] + NQOUTidx = var[6][1] + obs[idx][1, KINidx, NKINidx, QOUTidx, NQOUTidx] += weight.d + obs[idx][2, KINidx, NKINidx, QOUTidx, NQOUTidx] += weight.e +end + +function KW(para::ParaMC, diagram; + kin=[para.kF,], #amplitude of kin + nkin=[0,], # matfreq of kin + qout=[getK(0.0, para.dim, 1),], + nqout=[0,], + neval=1e6, #number of evaluations + print=0, + alpha=3.0, #learning ratio + config=nothing, + solver=:mcmc, + integrand::Function=integrand_ver3KW, + kwargs...) + + dim, β, kF, NF = para.dim, para.β, para.kF, para.NF + partition, diagpara, FeynGraphs, extT_labels, spin_conventions = diagram + + if NoBubble in diagpara[1].filter + UEG.MCinitialize!(para, false) + else + UEG.MCinitialize!(para, true) + end + + for p in diagpara + @assert diagpara[1].filter == p.filter "filter should be the same" + end + + @assert length(diagpara) == length(FeynGraphs) == length(extT_labels) == length(spin_conventions) + + Nkin = length(kin) + Nnkin = length(nkin) + Nqout = length(qout) + Nnqout = length(nqout) + + maxMomNum = maximum([key[1] for key in partition]) + 2 + funcGraphs! = Dict{Int,Function}() + leaf_maps = Vector{Dict{Int,Graph}}() + + for (i, key) in enumerate(partition) + funcGraphs![i], leafmap = Compilers.compile(FeynGraphs[key]) + push!(leaf_maps, leafmap) + end + leafStat, loopBasis = FeynmanDiagram.leafstates(leaf_maps, maxMomNum) + + momLoopPool = FrontEnds.LoopPool(:K, dim, loopBasis) + root = zeros(Float64, maximum(length.(extT_labels))) + println("static compile has finished!") + + K = MCIntegration.FermiK(dim, kF, 0.2 * kF, 10.0 * kF, offset=2) + K.data[:, 1] .= UEG.getK(kin[1], dim, 1) + K.data[:, 2] .= UEG.getK(kin[1], dim, 1) .- qout[1] + T = MCIntegration.Continuous(0.0, β, offset=1, alpha=alpha) + T.data[1] = 0.0 + + KIN = MCIntegration.Discrete(1, Nkin, alpha=alpha) + NKIN = MCIntegration.Discrete(1, Nnkin, alpha=alpha) + QOUT = MCIntegration.Discrete(1, Nqout, alpha=alpha) + NQOUT = MCIntegration.Discrete(1, Nnqout, alpha=alpha) + + dof = [[p.innerLoopNum, p.totalTauNum - 1, 1, 1, 1, 1] for p in diagpara] # K, T, ExtKidx + obs = [zeros(ComplexF64, 2, Nkin, Nnkin, Nqout, Nnqout) for p in diagpara] + println("obs size:", size(obs[1])) + + if isnothing(config) + config = MCIntegration.Configuration(; + var=(K, T, KIN, NKIN, QOUT, NQOUT), + dof=dof, + obs=obs, + type=Weight, + # type=ComplexF64, # type of the integrand + userdata=(para, kin, nkin, qout, nqout, maxMomNum, extT_labels, + spin_conventions, leafStat, leaf_maps, momLoopPool, + root, funcGraphs!), + kwargs... + ) + end + result = integrate(integrand; measure=measure_ver3KW, config=config, solver=solver, neval=neval, print=print, kwargs...) + + if isnothing(result) == false + if print >= 0 + report(result.config) + # report(result; pick=o -> (real(o[1, 1, 1])), name="uu") + # report(result; pick=o -> (real(o[2, 1, 1])), name="ud") + end + + datadict = Dict{eltype(partition),Any}() + for k in 1:length(dof) + avg, std = result.mean[k], result.stdev[k] + r = measurement.(real(avg), real(std)) + i = measurement.(imag(avg), imag(std)) + data = Complex.(r, i) + datadict[partition[k]] = data + end + return datadict, result + else + return nothing, nothing + end + +end + +function MC_KW_angle(para; + kin=[para.kF,], nkin=[0,], + qout=[get(0.0, para.dim, 1),], nqout=[0,], + neval=1e6, filename::Union{String,Nothing}=nothing, reweight_goal=nothing, + filter=[NoHartree,], + channels=[PHr, PHEr, PPr, Alli], + partition=UEG.partition(para.order), + transferLoop=nothing, extK=nothing, optimize_level=1, + verbose=0) + kF = para.kF + diagram = Diagram.diagram_parquet_response(:vertex3, para, partition, + channels=channels, filter=filter, extK=extK, transferLoop=transferLoop, optimize_level=optimize_level) + partition = diagram[1] + println(partition) + neighbor = UEG.neighbor(partition) + + if isnothing(reweight_goal) + reweight_goal = Float64[] + for (order, sOrder, vOrder) in partition + # push!(reweight_goal, 8.0^(order + vOrder - 1)) + push!(reweight_goal, 8.0^(order - 1)) + end + push!(reweight_goal, 1.0) + end + + ver3, result = Ver3.KW(para, diagram; + kin=kin, nkin=nkin, + qout=qout, nqout=nqout, + neval=neval, print=verbose, + neighbor=neighbor, reweight_goal=reweight_goal) + + if isnothing(ver3) == false + if isnothing(filename) == false + jldopen(filename, "a+") do f + key = "$(UEG.short(para))" + if haskey(f, key) + @warn("replacing existing data for $key") + delete!(f, key) + end + f[key] = (kin, nkin, qout, nqout, ver3) + end + end + end + + return ver3, result +end \ No newline at end of file diff --git a/src/vertex3/vertex3_new.jl b/src/vertex3/vertex3_new.jl new file mode 100644 index 0000000..6cc29e9 --- /dev/null +++ b/src/vertex3/vertex3_new.jl @@ -0,0 +1,29 @@ +module Ver3 + +using JLD2, CSV +using Printf, LinearAlgebra, DataFrames +using ..StaticArrays +using ..Parameters +using ..CompositeGrids +using ..ElectronGas +using ..MCIntegration +using ..Lehmann + +using ..FeynmanDiagram +import ..FeynmanDiagram.FrontEnds: TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan +import ..FeynmanDiagram.FrontEnds: Filter, NoHartree, NoFock, DirectOnly, Wirreducible, Girreducible, NoBubble, Proper +import ..FeynmanDiagram.FrontEnds: Response, Composite, ChargeCharge, SpinSpin, UpUp, UpDown +import ..FeynmanDiagram.FrontEnds: AnalyticProperty, Instant, Dynamic +import ..FeynmanDiagram.Parquet: DiagPara, Ver4Diag +using ..Measurements +using ..Diagram +# push!(LOAD_PATH, "../common/") +using ..UEG +using ..Propagator +import ..Propagator: LeafStateADDynamic + +import ..Weight + +include("ver3_KW.jl") + +end \ No newline at end of file diff --git a/test/vertex3.jl b/test/vertex3.jl index 79cfdaf..367888f 100644 --- a/test/vertex3.jl +++ b/test/vertex3.jl @@ -1,21 +1,77 @@ + +import FeynmanDiagram.FrontEnds: Filter, NoHartree, NoFock, DirectOnly, Wirreducible, Girreducible, NoBubble, Proper + function compare(data, expect) # println(data, ", ", expect) @test isapprox(data.val, expect, atol=5 * data.err) end @testset "Vertex3" begin - ### test Yukawa interaction ########### - p = (1, 0, 0) - mass2 = 0.01 - ####################### STATIC ######################### - para = UEG.ParaMC(rs=5.0, beta=25.0, Fs=0.0, order=1, mass2=mass2, isDynamic=true) - diagram = Ver3.diagram(para, [p,]) - kin = [[para.kF, 0.0, 0.0],] - qout = [[0.0, 0.0, 0.0],] - data, result = Ver3.KW(para, diagram; neval = 1e6, kin=kin, qout=qout, nkin=[0,], nqout=[1, ]) - obs = data[p] - # for the one-loop vertex3 diagram, we expect - # \gamma_3(q=0, w->0) = - dIm\Sigma)/dw_n - expect = -0.44275 - compare(real(obs[1]), expect) -end \ No newline at end of file + + # @testset "Vertex3 init" begin + # para = UEG.ParaMC(rs=5.0, beta=25.0, Fs=0.0, order=2, mass2=0.01, isDynamic=false) + # ver3, result = Ver3.MC_KW_angle(para) + # println(ver3[(2, 0, 0)]) + # end + + @testset "Vertex3 Dynamic O(1)" begin + para = UEG.ParaMC(rs=5.0, beta=25.0, Fs=0.0, order=1, mass2=0.01, isDynamic=true) + kin = [para.kF,] + # Nth = 4 + # theta = [(i) / (Nth * π) for i in 0:Nth] # N+1 points + theta = [0.0,] + qout = [[para.kF * (1 - cos(θ)), -para.kF * sin(θ), 0.0] for θ in theta] + ver3, result = Ver3.MC_KW_angle(para; kin=kin, qout=qout, nkin=[0,], nqout=[-1,]) + obs = ver3[(1, 0, 0)][:, 1, 1, :, 1] + println(obs) + # println(ver3[(2, 0, 0)][:, 1, 1, :, 1]) + + # for the one-loop vertex3 diagram, we expect + # \gamma_3(q=0, w->0) = - dIm\Sigma)/dw_n + # using ElectronGas + # para = Parameter.rydbergUnit(1/25, 5.0, 3, Λs=0.01) + # sigma=ElectronGas.SelfEnergy.G0W0(para) + # z=SelfEnergy.zfactor(para, sigma[1]; ngrid=[-1,0])[1] + # ngrid=[-1, 0] -> -0.508 + # ngrid=[0,1] -> -0.44 + expect = -0.44275 + compare(real(obs[1]), expect) + + end + + # @testset "Vertex3 static" begin + # para = UEG.ParaMC(rs=1.0, beta=25.0, Fs=0.0, order=3, mass2=3.5, isDynamic=false) + # kin = [para.kF,] + # Nth = 4 + # theta = [(i) / (Nth * π) for i in 0:Nth] # N+1 points + # # theta = [0.0,] + # qout = [[para.kF * (1 - cos(θ)), -para.kF * sin(θ), 0.0] for θ in theta] + # transferLoop = [1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + # ver3, result = Ver3.MC_KW_angle(para; + # kin=kin, qout=qout, nkin=[0,], nqout=[0,], + # filter=[NoHartree, Proper], transferLoop=transferLoop) + # println(ver3) + # # obs = ver3[(1, 0, 0)][:, 1, 1, :, 1] + # # println(obs) + # println(ver3[(1, 0, 0)][:, 1, 1, :, 1] + ver3[(2, 0, 0)][:, 1, 1, :, 1] + ver3[(3, 0, 0)][:, 1, 1, :, 1]) + + # end + +end + +# @testset "Vertex3" begin +# ### test Yukawa interaction ########### +# p = (1, 0, 0) +# mass2 = 0.01 +# ####################### STATIC ######################### +# para = UEG.ParaMC(rs=5.0, beta=25.0, Fs=0.0, order=1, mass2=mass2, isDynamic=true) +# diagram = Ver3.diagram(para, [p,]) +# kin = [[para.kF, 0.0, 0.0],] +# qout = [[0.0, 0.0, 0.0],] +# data, result = Ver3.KW(para, diagram; neval = 1e6, kin=kin, qout=qout, nkin=[0,], nqout=[1, ]) +# obs = data[p] +# # for the one-loop vertex3 diagram, we expect +# # \gamma_3(q=0, w->0) = - dIm\Sigma)/dw_n +# expect = -0.44275 +# compare(real(obs[1]), expect) +# end \ No newline at end of file diff --git a/test/vertex4PP.jl b/test/vertex4PP.jl index e60abba..158baf7 100644 --- a/test/vertex4PP.jl +++ b/test/vertex4PP.jl @@ -29,7 +29,8 @@ end @testset "PP" begin seed = 1234 - p = (0, 0, 0) + # p = (0, 0, 0) + partition = [(0, 0, 0), (1, 0, 0)] rs = 5.0 beta = 25 mass2 = 1e-2 @@ -38,14 +39,14 @@ end UEG.MCinitialize!(para) println(para) # diagram = Ver4.diagram(para, [p,]; channel=[], filter=[]) - diagram = Diagram.diagram_parquet_response(:vertex4, para, [p,]; channels=[], filter=[]) + diagram = Diagram.diagram_parquet_response(:vertex4, para, partition; channels=[], filter=[]) ############################ generic PH one-angle average ########################### nlist = [0, 1, 2] paras = [Ver4.OneAngleAveraged(para, [para.kF, para.kF], [[0, nlist[1], -1], [0, nlist[2], -1], [0, nlist[3], -1]], :PP, 0),] data, result = Ver4.one_angle_averaged(paras, diagram; neval=neval, print=-1, seed=seed) # data, result = Ver4.one_angle_averaged_ParquetAD(paras, diagram; neval=neval, print=-1, seed=seed) - obs = data[p] + obs = data[partition[1]] # println("obs 1:", obs[:, 1, 1]) # println("obs 2:", obs[:, 2, 1]) # println("obs 3:", obs[:, 3, 1])