skip contact evaluation when building GreenSlice identity

src/greenfunction.jl

@@ -57,7 +57,7 @@ end
 function call!(g::GreenFunction{T}, ω::Complex{T}; params...) where {T}
     gsolver = solver(g)
     contacts´ = contacts(g)
-    Σblocks = call!(contacts´, ω; params...)
+    Σblocks = supports_contacts(gsolver) ? call!(contacts´, ω; params...) : missing
     corbs = contactorbitals(contacts´)
     slicer = build_slicer(gsolver, g, ω, Σblocks, corbs; params...)
     return GreenSolution(g, slicer, Σblocks, corbs)
@@ -79,6 +79,8 @@ retarded_omega(ω::T, s::AppliedGreenSolver) where {T<:Real} =
 # fallback, may be overridden
 needs_omega_shift(s::AppliedGreenSolver) = true
 
+supports_contacts(s::AppliedGreenSolver) = true
+
 #endregion
 
 ############################################################################################

src/observables.jl

@@ -411,13 +411,14 @@ function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegama
 end
 
 # we need to add the arc path segment from -∞ to ∞ * p.cisinf
+# we use the syntax gs(::UniformScaling) to find the identity matrix of our slice, see internal.jl
 function post_transform_rho(p::RadialPath, gs)
     arc = gs((p.angle/π)*I)
-    function post(x)
+    function post!(x)
         x .+= arc
         return x
     end
-    return post
+    return post!
 end
 
 post_transform_rho(::AbstractIntegrationPath, _) = identity

src/solvers/green/internal.jl

@@ -3,6 +3,7 @@
 #   Given a function f(ω) that returns a function fω(::CellOrbital, ::CellOrbital),
 #   implement an AppliedGreenSolver that returns a s::ModelGreenSlicer that returns fω(i,j)
 #   for each single orbital when calling getindex(s, ::CellOrbitals...). view not supported.
+#   Contact self energies are ignored.
 #region
 
 struct AppliedModelGreenSolver{F} <: AppliedGreenSolver
@@ -18,6 +19,7 @@ end
 
 ## API ##
 
+# build an identity matrix over the slice. We use an AppliedModelGreenSolver for this
 function (g::GreenSlice)(x::UniformScaling)
     s = AppliedModelGreenSolver(Returns((i, j) -> ifelse(i == j, x.λ, zero(x.λ))))
     g´ = swap_solver(g, s)
@@ -37,6 +39,8 @@ end
 
 needs_omega_shift(s::AppliedModelGreenSolver) = false
 
+supports_contacts(s::AppliedModelGreenSolver) = false
+
 minimal_callsafe_copy(s::AppliedModelGreenSolver, args...) = s
 
 Base.getindex(s::ModelGreenSlicer, is::CellOrbitals, js::CellOrbitals) =

test/test_greenfunction.jl

@@ -589,6 +589,11 @@ end
     @test real(only(ρ(0, 0; o = -1.5))) < 1e-7
     d = Quantica.integrand(ρ)
     @test_throws MethodError d(2; o = 0.8)
+
+    # test that contacts are ignored when building a GreenSlice identity
+    g = LP.linear() |> supercell(2) |> hopping(1) |> attach(@onsite(ω->error())) |> greenfunction
+    @test_throws ErrorException g(0)[]
+    @test g[sites(1), sites(1:2)](2I) == SA[2 0]  # doesn't throw
 end
 
 @testset "greenfunction aliasing" begin