Dev/keldysh dlr new imaginary time Green's function type with DLR representation using Lehmann.jl

Project.toml

@@ -12,6 +12,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Keldysh = "50f2bc7e-5fd7-11e9-13c7-85fb88b4f34f"
 KeldyshED = "675b6b9c-7c2f-11e9-3bf3-dfd4d61640f7"
+Lehmann = "95bf888a-8996-4655-9f35-1c0506bdfefe"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"

docs/src/api.md

@@ -88,6 +88,19 @@ SplineInterpolatedGF(::GFType; ::kd.TimeGridPoint) where {
         T <: Number, scalar, GFType <: kd.AbstractTimeGF{T, scalar}}
 ```
 
+## [`QInchworm.keldysh_dlr`](@id api:QInchworm.keldysh_dlr)
+
+```@meta
+CurrentModule = QInchworm.keldysh_dlr
+```
+```@docs
+keldysh_dlr
+```
+```@docs
+DLRImaginaryTimeGrid
+DLRImaginaryTimeGF
+```
+
 ## [`QInchworm.utility`](@id api:QInchworm.utility)
 
 ```@meta

docs/src/index.md

@@ -43,6 +43,7 @@ two-point correlation functions for it.
 - [`QInchworm.exact_atomic_ppgf`](@ref QInchworm.exact_atomic_ppgf)
 - [`QInchworm.sector_block_matrix`](@ref QInchworm.sector_block_matrix)
 - [`QInchworm.spline_gf`](@ref QInchworm.spline_gf)
+- [`QInchworm.keldysh_dlr`](@ref QInchworm.keldysh_dlr)
 - [`QInchworm.utility`](@ref QInchworm.utility)
 - [`QInchworm.mpi`](@ref QInchworm.mpi)
 - [`QInchworm.scrambled_sobol`](@ref QInchworm.scrambled_sobol)

src/QInchworm.jl

@@ -27,6 +27,7 @@ include("utility.jl")
 include("sector_block_matrix.jl")
 include("mpi.jl")
 include("diagrammatics.jl")
+include("keldysh_dlr.jl")
 include("spline_gf.jl")
 include("ppgf.jl")
 include("exact_atomic_ppgf.jl")

src/keldysh_dlr.jl

@@ -0,0 +1,248 @@
+# QInchworm.jl
+#
+# Copyright (C) 2021-2024 I. Krivenko, H. U. R. Strand and J. Kleinhenz
+#
+# QInchworm.jl is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# QInchworm.jl is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# QInchworm.jl. If not, see <http://www.gnu.org/licenses/>.
+#
+# Authors: Hugo U. R. Strand, Igor Krivenko
+
+"""
+Extension of Keldysh.jl defining imaginary time Green's functions
+represented using the Discrete Lehmann Representation as implemented
+in Lehmann.jl.
+
+# Exports
+$(EXPORTS)
+"""
+module keldysh_dlr
+
+using DocStringExtensions
+
+using Lehmann; le = Lehmann
+
+using Keldysh; kd = Keldysh
+using Keldysh: AbstractTimeGrid
+using Keldysh: ImaginaryContour, TimeGridPoint, PeriodicStorage, GFSignEnum, BranchPoint
+
+import QInchworm.utility: ph_conj
+
+export DLRImaginaryTimeGrid, DLRImaginaryTimeGF
+export ph_conj
+
+#
+# DLRImaginaryTimeGrid
+#
+
+"""
+    $(TYPEDEF)
+
+Wrapper around Lehmann.jl describing a Discrete Lehmann Representation imaginary time grid
+conforming to the interface of Keldysh.jl TimeGrids.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+struct DLRImaginaryTimeGrid <: AbstractTimeGrid
+    contour::ImaginaryContour
+    points::Vector{TimeGridPoint}
+    branch_bounds::NTuple{1, Pair{TimeGridPoint, TimeGridPoint}}
+    ntau::Int
+    dlr::le.DLRGrid
+
+    function DLRImaginaryTimeGrid(c::ImaginaryContour, dlr::le.DLRGrid)
+        points::Vector{TimeGridPoint} = []
+        τ_branch = c.branches[1]
+        @assert τ_branch.domain == kd.imaginary_branch
+        for (idx, τ) in enumerate(dlr.τ)
+            point = TimeGridPoint(1, idx, τ_branch(τ/dlr.β))
+            push!(points, point)
+        end
+        τ_0 = TimeGridPoint(-1, -1, τ_branch(0.))
+        τ_β = TimeGridPoint(-1, -1, τ_branch(1.))
+
+        branch_bounds = ( Pair(τ_0, τ_β), )
+        ntau = length(dlr.τ)
+        return new(c, points, branch_bounds, ntau, dlr)
+    end
+end
+
+Base.:isequal(A::T, B::T) where T <: DLRImaginaryTimeGrid = all(A.points .== B.points)
+
+"""
+    $(TYPEDEF)
+
+Wrapper around Lehmann.jl describing a Discrete Lehmann Representation imaginary time Green's
+function conforming to the interface of Keldysh.jl AbstractTimeGF.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+struct DLRImaginaryTimeGF{T, scalar} <: AbstractTimeGF{T, scalar}
+    grid::DLRImaginaryTimeGrid
+    mat::PeriodicStorage{T,scalar}
+    ξ::GFSignEnum
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Make a [`DLRImaginaryTimeGF`](@ref) from a [`DLRImaginaryTimeGrid`](@ref)
+following the API of Keldysh.ImaginarTimeGF.
+
+"""
+function DLRImaginaryTimeGF(::Type{T},
+                            grid::DLRImaginaryTimeGrid,
+                            norb=1,
+                            ξ::GFSignEnum=fermionic,
+                            scalar=false) where T <: Number
+    ntau = grid.ntau
+    mat = PeriodicStorage(T, ntau, norb, scalar)
+    return DLRImaginaryTimeGF(grid, mat, ξ)
+end
+
+DLRImaginaryTimeGF(grid::DLRImaginaryTimeGrid,
+                   norb=1,
+                   ξ::GFSignEnum=fermionic,
+                   scalar=false) = DLRImaginaryTimeGF(ComplexF64, grid, norb, ξ, scalar)
+
+norbitals(G::DLRImaginaryTimeGF) = G.mat.norb
+
+#
+# Arithmetic operations
+#
+
+function Base.:+(A::T, B::T) where T <: DLRImaginaryTimeGF
+    @assert A.grid == B.grid
+    return T(A.grid, A.mat + B.mat, A.ξ)
+end
+
+function Base.:-(A::T, B::T) where T <: DLRImaginaryTimeGF
+    @assert A.grid == B.grid
+    return T(A.grid, A.mat - B.mat, A.ξ)
+end
+
+Base.:-(A::T) where T <: DLRImaginaryTimeGF = T(A.grid, -A.mat, A.ξ)
+
+Base.:*(G::T, α::Number) where T <: DLRImaginaryTimeGF = T(G.grid, α * G.mat, G.ξ)
+Base.:*(α::Number, G::T) where T <: DLRImaginaryTimeGF = G * α
+
+function Base.:isapprox(A::T, B::T) where T <: DLRImaginaryTimeGF
+    return (A.mat.data ≈ B.mat.data) && (A.grid == B.grid)
+end
+
+#
+# Matrix valued Gf interpolator interface
+#
+
+function (G::DLRImaginaryTimeGF{T, false})(t1::BranchPoint, t2::BranchPoint) where T
+  norb = norbitals(G)
+  x = zeros(T, norb, norb)
+  return interpolate!(x, G, t1, t2)
+end
+
+function interpolate!(x, G::DLRImaginaryTimeGF{T, false},
+                      t1::BranchPoint, t2::BranchPoint) where T
+    dlr = G.grid.dlr
+    τ = dlr.β*(t1.ref - t2.ref)
+    @assert 0. <= τ <= dlr.β
+    x[:] = le.dlr2tau(dlr, g.mat.data, [τ], axis=3)
+    return x
+end
+
+#
+# Scalar valued Gf interpolator interface
+#
+
+function (G::DLRImaginaryTimeGF{T, true})(t1::BranchPoint, t2::BranchPoint) where T
+    return interpolate(G, t1, t2)
+end
+
+function interpolate(G::DLRImaginaryTimeGF{T, true},
+                     t1::BranchPoint, t2::BranchPoint)::T where T
+    dlr = G.grid.dlr
+    τ = dlr.β*(t1.ref - t2.ref)
+    @assert 0. <= τ <= dlr.β
+    return le.dlr2tau(dlr, G.mat.data, [τ], axis=3)[1, 1, 1]
+end
+
+#
+# Green's function creation from function
+#
+
+"""
+    $(TYPEDSIGNATURES)
+
+Make a [`DLRImaginaryTimeGF`](@ref) from a function
+
+"""
+DLRImaginaryTimeGF(f::Function,
+                   grid::DLRImaginaryTimeGrid,
+                   norb=1,
+                   ξ::GFSignEnum=fermionic,
+                   scalar=false) = DLRImaginaryTimeGF(f, ComplexF64, grid, norb, ξ, scalar)
+
+function DLRImaginaryTimeGF(f::Function,
+                            ::Type{T},
+                            grid::DLRImaginaryTimeGrid,
+                            norb=1,
+                            ξ::GFSignEnum=fermionic,
+                            scalar=false) where T <: Number
+    G = DLRImaginaryTimeGF(T, grid, norb, ξ, scalar)
+
+    t0 = TimeGridPoint(1, -1, BranchPoint(im * 0., 0., kd.imaginary_branch))
+
+    g_τ = Array{ComplexF64, 3}(undef, norb, norb, length(grid.dlr.τ))
+    for (idx, t) in enumerate(G.grid.points)
+        if scalar
+            g_τ[1, 1, idx] = f(t, t0)
+        else
+            g_τ[:, :, idx] = f(t, t0)
+        end
+    end
+
+    G.mat.data[:] = le.tau2dlr(G.grid.dlr, g_τ, axis=3)
+
+    return G
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Make a [`DLRImaginaryTimeGF`](@ref) from a Keldysh.jl AbstractDOS object.
+
+"""
+function DLRImaginaryTimeGF(dos::AbstractDOS, grid::DLRImaginaryTimeGrid)
+    β = grid.contour.β
+    DLRImaginaryTimeGF(grid, 1, fermionic, true) do t1, t2
+        kd.dos2gf(dos, β, t1.bpoint, t2.bpoint)
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Given a scalar-valued imaginary time Green's function ``g(\\tau)``, return its particle-hole
+conjugate ``g(\\beta-\\tau)``.
+"""
+function ph_conj(g::DLRImaginaryTimeGF{T, true})::DLRImaginaryTimeGF{T, true} where {T}
+    g_rev = deepcopy(g)
+
+    dlr = g.grid.dlr
+    g_τ = le.dlr2tau(dlr, g.mat.data, dlr.β .- dlr.τ, axis=3)
+    g_rev.mat.data[:] = le.tau2dlr(dlr, g_τ, axis=3)
+
+    return g_rev
+end
+
+end # module keldysh_dlr

test/Project.toml

@@ -4,6 +4,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Keldysh = "50f2bc7e-5fd7-11e9-13c7-85fb88b4f34f"
 KeldyshED = "675b6b9c-7c2f-11e9-3bf3-dfd4d61640f7"
+Lehmann = "95bf888a-8996-4655-9f35-1c0506bdfefe"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"