v1.2.0 (#18)

* Added the Langreth rules

Project.toml

@@ -1,10 +1,11 @@
 name = "KadanoffBaym"
 uuid = "82532805-809c-4ef0-842b-4b00c5e9be5f"
 authors = ["Francisco Meirinhos, Tim Bode"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 
@@ -15,6 +16,7 @@ QuadGK = "2.4.2"
 julia = "1.7"
 
 [extras]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]

src/KadanoffBaym.jl

@@ -1,18 +1,25 @@
 module KadanoffBaym
 
-import OrdinaryDiffEq
+using LinearAlgebra
 using QuadGK
 using AbstractFFTs
-
-export GreenFunction, Symmetrical, SkewHermitian
-export kbsolve!
-export wigner_transform
+import OrdinaryDiffEq
 
 include("utils.jl")
+
 include("gf.jl")
+export GreenFunction, Symmetrical, SkewHermitian
+
 include("vie.jl")
 include("vcabm.jl")
 include("kb.jl")
+export kbsolve!
+
 include("wigner.jl")
+export wigner_transform
+
+include("langreth.jl")
+export TimeOrderedGreenFunction, TimeOrderedConvolution, conv
+export greater, lesser, advanced, retarded
 
 end # module

src/langreth.jl

@@ -0,0 +1,37 @@
+abstract type AbstractTimeOrderedGreenFunction end
+
+"""
+"""
+struct TimeOrderedGreenFunction <: AbstractTimeOrderedGreenFunction
+  L   # Lesser
+  G   # Greater
+
+  TimeOrderedGreenFunction(L, G) = new(L, G)
+end
+
+Base.:*(a::Number, g::TimeOrderedGreenFunction) = TimeOrderedGreenFunction(a * lesser(g), a * greater(g))
+Base.:+(g1::TimeOrderedGreenFunction, g2::TimeOrderedGreenFunction) = TimeOrderedGreenFunction(lesser(g1) + lesser(g2), greater(g1) + greater(g2))
+
+"""
+"""
+struct TimeOrderedConvolution{TA<:AbstractTimeOrderedGreenFunction, TB<:AbstractTimeOrderedGreenFunction} <: AbstractTimeOrderedGreenFunction
+  A::TA
+  B::TB
+  dts::UpperTriangular
+end
+
+"""
+"""
+function conv(A::AbstractTimeOrderedGreenFunction, B::AbstractTimeOrderedGreenFunction, dts)
+  c = TimeOrderedConvolution(A, B, dts)
+  return TimeOrderedGreenFunction(lesser(c), greater(c))
+end
+
+# Langreth's rules
+greater(g::TimeOrderedGreenFunction) = g.G
+lesser(g::TimeOrderedGreenFunction) = g.L
+advanced(g::TimeOrderedGreenFunction) = UpperTriangular(lesser(g) - greater(g))
+retarded(g::TimeOrderedGreenFunction) = adjoint(advanced(g))
+
+greater(g::TimeOrderedConvolution) = (retarded(g.A) .* adjoint(g.dts)) * greater(g.B) + greater(g.A) * (g.dts .* advanced(g.B)) |> skew_hermitify!
+lesser(g::TimeOrderedConvolution) = (retarded(g.A) .* adjoint(g.dts)) * lesser(g.B) + lesser(g.A) * (g.dts .* advanced(g.B)) |> skew_hermitify!

src/utils.jl

@@ -29,3 +29,14 @@ _last2() = throw(ArgumentError("Cannot call last2 on an empty tuple."))
 _last2(v) = throw(ArgumentError("Cannot call last2 on 1-element tuple."))
 @inline _last2(v1, v2) = (v1, v2)
 @inline _last2(v, t...) = _last2(t...)
+
+function skew_hermitify!(x)
+  for i in 1:size(x, 1)
+    for j in 1:(i - 1)
+      x[j,i] = -conj(x[i,j])
+    end
+
+    x[i,i] = eltype(x) <: Real ? 0.0 : im * imag(x[i,i])
+  end
+  return x
+end

test/langreth.jl

@@ -0,0 +1,67 @@
+function integrate(x::AbstractVector, y::AbstractVector)
+  if isone(length(y))
+    return zero(first(y))
+  end
+  @inbounds retval = (x[2] - x[1]) * (y[1] + y[2])
+  @inbounds @fastmath @simd for i in 2:(length(y)-1)
+    retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
+  end
+  return 1 // 2 * retval
+end
+
+function compute(A, B, ts)
+  l = zero(lesser(A))
+  for i in 1:size(l, 1)
+    for j in 1:i
+      l[i,j] += integrate(ts, greater(A)[i, 1:i] .* lesser(B)[1:i, j])
+      l[i,j] -= integrate(ts, lesser(A)[i, 1:j] .* greater(B)[1:j, j])
+      l[i,j] -= integrate(ts[j:i], lesser(A)[i, j:i] .* lesser(B)[j:i, j])
+
+      l[j,i] = -conj(l[i,j])
+    end
+    l[i,i] = eltype(l) <: Real ? 0.0 : im * imag(l[i,i])
+  end
+
+  g = zero(greater(A))
+  for i in 1:size(g, 1)
+    for j in 1:i
+      g[i,j] -= integrate(ts, lesser(A)[i, 1:i] .* greater(B)[1:i, j])
+      g[i,j] += integrate(ts, greater(A)[i, 1:j] .* lesser(B)[1:j, j])
+      g[i,j] += integrate(ts[j:i], greater(A)[i, j:i] .* greater(B)[j:i, j])
+
+      g[j,i] = -conj(g[i,j])
+    end
+    g[i,i] = eltype(g) <: Real ? 0.0 : im * imag(g[i,i])
+  end
+
+  TimeOrderedGreenFunction(l, g)
+end
+
+N = 100
+
+# Suppose a non-equidistant time-grid
+ts = sort(N*rand(N));
+
+# And 2 time-ordered GFs defined on that grid
+a = let
+  x = rand(ComplexF64,N,N)
+  y = rand(ComplexF64,N,N)
+  TimeOrderedGreenFunction(x - x', y - y') # Skew-symmetric L and G components
+end
+b = let
+  x = rand(ComplexF64,N,N)
+  y = rand(ComplexF64,N,N)
+  TimeOrderedGreenFunction(x - x', y - y') # Skew-symmetric L and G components
+end
+
+dts = reduce(hcat, ([KadanoffBaym.calculate_weights(ts[1:i], ones(Int64, i-1), 1e-8, 1e-3); zeros(length(ts)-i)] for i in eachindex(ts))) |> UpperTriangular
+
+⋆(a, b) = conv(a, b, dts)
+
+@testset "Langreth's rules" begin
+  @test greater(a ⋆ b) ≈ greater(compute(a, b, ts))
+  @test lesser(a ⋆ b) ≈ lesser(compute(a, b, ts))
+  @test retarded(a ⋆ b) ≈ retarded(compute(a, b, ts))
+  @test greater(a ⋆ b ⋆ a) ≈ greater(compute(compute(a, b, ts), a, ts))
+  @test retarded(a ⋆ b ⋆ a) ≈ retarded(compute(compute(a, b, ts), a, ts))
+end

test/runtests.jl

@@ -1,4 +1,5 @@
 using KadanoffBaym
+using LinearAlgebra
 using Test
 
 @testset verbose=true "KadanoffBaym.jl" begin
@@ -9,4 +10,8 @@ using Test
   @testset "Solver" begin
     include("kb.jl")
   end
+
+  @testset "Langreth" begin
+    include("langreth.jl")
+  end
 end