Merge pull request #330 from pablosanjose/pointkwargs

Disallow integrand params to be different from integration path params

src/docstrings.jl

@@ -2186,8 +2186,7 @@ josephson
 Return the complex integrand `d::JosephsonIntegrand` whose integral over frequency yields
 the Josephson current, `J(kBT; params...) = Re(∫dx d(x; params...))`, where `ω(x) =
 Quantica.point(x, d)` is a path parametrization over real variable `x` . To evaluate the `d`
-for a given `x` and parameters, use `d(x; params...)`, or `call!(d, x; params...)` for its
-mutating (non-allocating) version.
+for a given `x`, use `d(x)`, or `call!(d, x)` for its mutating (non-allocating) version.
 
     Quantica.integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0, kBT = 0; params...)
 

src/greenfunction.jl

@@ -37,6 +37,9 @@ call!_output(c::Contacts) = selfenergyblocks(c)
 
 ############################################################################################
 # GreenFuntion call! API
+#   `symmetrize` and `post` allow to obtain gij and conj(gji) combinations efficiently.
+#   Useful when computing e.g. the spectral density A = (Gʳ-Gʳ')/2πi on a non-diagonal slice
+#   These kwargs are currently not documented - see specialmatrices.jl
 #region
 
 (g::GreenFunction)(ω; params...) = minimal_callsafe_copy(call!(g, ω; params...))

src/integrator.jl

@@ -157,9 +157,9 @@ options(I::Integrator) = I.quadgk_opts
 
 ## call! ##
 # scalar version
-function call!(I::Integrator{<:Any,Missing}; params...)
+function call!(I::Integrator{<:Any,Missing})
     fx = x -> begin
-        y = call!(I.integrand, x; params...)  # should be a scalar
+        y = call!(I.integrand, x)  # should be a scalar
         I.callback(x, y)
         return y
     end
@@ -169,9 +169,9 @@ function call!(I::Integrator{<:Any,Missing}; params...)
 end
 
 # nonscalar version
-function call!(I::Integrator{<:Any,T}; params...) where {T}
+function call!(I::Integrator{<:Any,T}) where {T}
     fx! = (y, x) -> begin
-        y .= serialize(call!(I.integrand, x; params...))
+        y .= serialize(call!(I.integrand, x))
         I.callback(x, y)
         return nothing
     end
@@ -181,7 +181,7 @@ function call!(I::Integrator{<:Any,T}; params...) where {T}
     return result´
 end
 
-(I::Integrator)(; params...) = copy(call!(I; params...))
+(I::Integrator)() = copy(call!(I))
 
 #endregion
 #endregion

src/observables.jl

@@ -361,11 +361,10 @@ end
 
 # produces integrand_transform!(gs(ω´; omegamap(ω´)..., params...) * f(ω´-mu))
 # with ω´ = path_transform(ω)
-struct DensityMatrixIntegrand{G<:GreenSlice,T,O,P<:AbstractIntegrationPath,PT}
-    gs::G
+struct DensityMatrixIntegrand{T,GF<:Function,P<:AbstractIntegrationPath,PT}
+    gsfunc::GF         # (ω, symmetrize) -> gs(ω; symmetrize, omegamap(ω)..., params...)
     mu::T
     kBT::T
-    omegamap::O        # function: returns changed system parameters for each ω
     path::P            # AbstractIntegrationPath object
     pts::PT            # ω-points that define specific integration path, derived from `path`
 end
@@ -404,7 +403,8 @@ function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegama
     function ifunc(mu, kBT; params...)
         pts = points(path, mu, kBT; params...)
         realpts = realpoints(path, pts)
-        ρd = DensityMatrixIntegrand(gs, T(mu), T(kBT), omegamap, path, pts)
+        gsfunc(ω, symmetrize) = call!(gs, ω; symmetrize, omegamap(ω)..., params...)
+        ρd = DensityMatrixIntegrand(gsfunc, T(mu), T(kBT), path, pts)
         return Integrator(result, ρd, realpts; opts´...)
     end
     return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs)
@@ -431,19 +431,19 @@ post_transform_rho(::AbstractIntegrationPath, _) = identity
     ρ.solver(mu, kBT; params...)
 # special case for integrator solver
 (ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0; params...) =
-    ρ.solver(mu, kBT; params...)(; params...)
+    ρ.solver(mu, kBT; params...)()
 
 (s::DensityMatrixIntegratorSolver)(mu, kBT; params...) =
     s.ifunc(mu, kBT; params...);
 
-(ρi::DensityMatrixIntegrand)(x; params...) = copy(call!(ρi, x; params...))
+(ρi::DensityMatrixIntegrand)(x) = copy(call!(ρi, x))
 
-function call!(ρi::DensityMatrixIntegrand, x; params...)
+function call!(ρi::DensityMatrixIntegrand, x)
     ω = point(x, ρi.path, ρi.pts)
     j = jacobian(x, ρi.path, ρi.pts)
     f = fermi(chopsmall(ω - ρi.mu), inv(ρi.kBT))
     symmetrize = -j*f/(2π*im)
-    output = call!(ρi.gs, ω; symmetrize, ρi.omegamap(ω)..., params...)
+    output = ρi.gsfunc(ω, symmetrize)
     return output
 end
 
@@ -487,13 +487,12 @@ call!_output(ρ::DensityMatrix) = call!_output(ρ.gs)
 #   Keywords opts are passed to quadgk for the integral
 #region
 
-struct JosephsonIntegrand{T<:AbstractFloat,P<:Union{Missing,AbstractArray},O,G<:GreenFunction{T},PA,PT}
-    g::G
+struct JosephsonIntegrand{T<:AbstractFloat,P<:Union{Missing,AbstractArray},GF<:Function,PA,PT}
+    gfunc::GF                   # ω -> g(ω; params...)
     kBT::T
     contactind::Int             # contact index
     tauz::Vector{Int}           # precomputed diagonal of tauz
     phaseshifts::P              # missing or collection of phase shifts to apply
-    omegamap::O                 # function that maps ω to parameters
     path::PA                    # AbstractIntegrationPath
     pts::PT                     # points in actual integration path, derived from `path`
     traces::P                   # preallocated workspace
@@ -537,7 +536,8 @@ function josephson(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap =
     function ifunc(kBT; params...)
         pts = points(path, 0, kBT; params...)
         realpts = realpoints(path, pts)
-        jd = JosephsonIntegrand(g, T(kBT), contact, tauz, phases´, omegamap, path, pts,
+        gfunc(ω) = call!(g, ω; omegamap(ω)..., params...)
+        jd = JosephsonIntegrand(gfunc, T(kBT), contact, tauz, phases´, path, pts,
             traces, Σfull, Σ, similar(Σ), similar(Σ), similar(Σ), similar(tauz, Complex{T}))
         return Integrator(traces, jd, realpts; opts´...)
     end
@@ -562,15 +562,15 @@ end
 (j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...)
 # special case for integrator solver (so we can access integrand etc before integrating)
 (j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0; params...) =
-    j.solver(kBT; params...)(; params...)
+    j.solver(kBT; params...)()
 
 (s::JosephsonIntegratorSolver)(kBT; params...) = s.ifunc(kBT; params...)
 
-(J::JosephsonIntegrand)(x; params...) = copy(call!(J, x; params...))
+(J::JosephsonIntegrand)(x) = copy(call!(J, x))
 
-function call!(Ji::JosephsonIntegrand, x; params...)
+function call!(Ji::JosephsonIntegrand, x)
     ω = point(x, Ji.path, Ji.pts)
-    gω = call!(Ji.g, ω; Ji.omegamap(ω)..., params...)
+    gω = Ji.gfunc(ω)
     f = fermi(ω, inv(Ji.kBT))
     traces = josephson_traces(Ji, gω, f)
     traces = mul_scalar_or_array!(traces, jacobian(x, Ji.path, Ji.pts))
@@ -581,9 +581,11 @@ end
 
 #region ## API ##
 
-integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0; params...) = integrand(J.solver(kBT; params...))
+integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0; params...) =
+    integrand(J.solver(kBT; params...))
 
-points(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0; params...) = points(integrand(J, kBT; params...))
+points(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0; params...) =
+    points(integrand(J, kBT; params...))
 points(J::JosephsonIntegrand) = J.pts
 
 point(x, Ji::JosephsonIntegrand) = point(x, Ji.path, Ji.pts)
@@ -605,7 +607,8 @@ function josephson_traces(J, gω, f)
     return josephson_traces!(J, gr, Σi, f)
 end
 
-josephson_traces!(J::JosephsonIntegrand{<:Any,Missing}, gr, Σi, f) = josephson_one_trace!(J, gr, Σi, f)
+josephson_traces!(J::JosephsonIntegrand{<:Any,Missing}, gr, Σi, f) =
+    josephson_one_trace!(J, gr, Σi, f)
 
 function josephson_traces!(J, gr, Σi, f)
     for (i, phaseshift) in enumerate(J.phaseshifts)

test/test_greenfunction.jl

@@ -552,6 +552,16 @@ end
     ρ2 = densitymatrix(g[is, js], Paths.sawtooth(4))
     @test ρ0() ≈ ρ1() ≈ ρ2()
     @test ρ0(0.1, 0.2) ≈ ρ1(0.1, 0.2) ≈ ρ2(0.1, 0.2)
+
+    # parametric path and system
+    g = LP.linear() |> supercell |> @onsite((; o = 0.5) -> o) |> greenfunction
+    ρ = densitymatrix(g[], Paths.polygon((µ, T; o = 0.5) -> o > µ ? (-1, µ, o + im, 1) : (-1, o + im, µ, 1)))
+    @test real(only(ρ(0, 0; o = -0.5))) ≈ 1
+    @test real(only(ρ(0, 0; o = 0.5))) < 1e-7
+    @test real(only(ρ(0, 0; o = 1.5))) < 1e-7
+    @test real(only(ρ(0, 0; o = -1.5))) < 1e-7
+    d = Quantica.integrand(ρ)
+    @test_throws MethodError d(2; o = 0.8)
 end
 
 @testset "greenfunction aliasing" begin