Integration path override in densitymatrix and josephson

src/docstrings.jl

@@ -1976,9 +1976,10 @@ said algorithm.
 
 As above, but using a generic algorithm that relies on numerical integration along a contour
 in the complex plane, between points `(ωmin, ωmax)` (or along a polygonal path connecting
-`ωpoints`, a collection of numbers), which should be chosen so as to encompass the full system bandwidth. Keywords
-`quadgk_opts` are passed to the `QuadGK.quadgk` integration routine. See below for
-additiona `opts`.
+`ωpoints`, a collection of numbers), which should be chosen so as to encompass the full
+system bandwidth. When the `ωpoints` are all real, an extra point is added at `ω = µ` to the
+integration path, to better integrate the step in the Fermi function. Keywords `quadgk_opts`
+are passed to the `QuadGK.quadgk` integration routine. See below for additiona `opts`.
 
     densitymatrix(gs::GreenSlice, ωmax::Number; opts...)
 
@@ -1990,9 +1991,16 @@ As above with `ωmin = -ωmax`.
 
 Evaluate the density matrix at chemical potential `μ` and temperature `kBT` (in the same
 units as the Hamiltonian) for the given `g` parameters `params`, if any. The result is given
-as an `OrbitalSliceMatrix`, see its docstring for further details. When the `ωpoints` are
-all real, an extra point is added at `ω = µ` to the integration path, to better integrate
-the step in the Fermi function.
+as an `OrbitalSliceMatrix`, see its docstring for further details.
+
+If the generic integration algorithm is used with complex `ωpoints`, the following form is
+also available:
+
+    ρ(μ = 0, kBT = 0, override = missing; params...)
+
+Here `override` can be a collection of `ωpoints` that will replace the original ones
+provided when defining `ρ`. Alternatively, it can be a function that will be applied to the
+original `ωpoints`. This may be useful when the integration path must depend on `params`.
 
 ## Algorithms and keywords
 
@@ -2035,23 +2043,29 @@ densitymatrix
 For a `gs = g[i::Integer]` slice of the `g::GreenFunction` of a hybrid junction, build a
 `J::Josephson` object representing the equilibrium (static) Josephson current `I_J` flowing
 into `g` through contact `i`, integrated along a polygonal contour connecting `ωpoints` (a
-collection of numbers) in the complex plane. The result of `I_J` is given in units of `qe/h`
-(`q` is the dimensionless carrier charge). `I_J` can be written as ``I_J = Re ∫ dω f(ω)
-j(ω)``, where ``j(ω) = (qe/h) × 2Tr[(ΣʳᵢGʳ - GʳΣʳᵢ)τz]``. Here `f(ω)` is the Fermi function
-with `µ = 0`.
+collection of numbers) in the complex plane. When the `ωpoints` are all real, an extra point
+is added at `ω = 0` to the integration path, to better integrate the step in the Fermi
+function.
+
+The result of `I_J` is given in units of `qe/h` (`q` is the dimensionless carrier charge).
+`I_J` can be written as ``I_J = Re ∫ dω f(ω) j(ω)``, where ``j(ω) = (qe/h) × 2Tr[(ΣʳᵢGʳ -
+GʳΣʳᵢ)τz]``. Here `f(ω)` is the Fermi function with `µ = 0`.
 
     josephson(gs::GreenSlice, ωmax::Real; kw...)
 
 As above, but with `ωpoints = (-ωmax, ωmax)`.
 
 ## Full evaluation
 
-    J(kBT = 0; params...)   # where J::Josephson
+    J(kBT = 0, override = missing; params...)   # where J::Josephson
 
 Evaluate the current `I_J` at chemical potemtial `µ = 0` and temperature `kBT` (in the same
-units as the Hamiltonian) for the given `g` parameters `params`, if any. When the `ωpoints`
-are all real, an extra point is added at `ω = 0` to the integration path, to better
-integrate the step in the Fermi function.
+units as the Hamiltonian) for the given `g` parameters `params`, if any.
+
+It's possible to `override` a complex integration path at evaluation time. In this case
+`override` can be a collection of `ωpoints` that will replace the original ones provided
+when defining `J`. Alternatively, it can be a function that will be applied to the original
+`ωpoints`. This may be useful when the integration path must depend on `params`.
 
 ## Keywords
 

src/observables.jl

@@ -482,13 +482,13 @@ end
 #region ## Constructors ##
 
 # generic fallback (for other solvers)
-(ρ::DensityMatrix)(mu = 0, kBT = 0; params...) = ρ.solver(mu, kBT; params...) |>
-    maybe_OrbitalSliceArray(axes(ρ.gs))
+(ρ::DensityMatrix)(mu = 0, kBT = 0; params...) =
+    ρ.solver(mu, kBT; params...) |> maybe_OrbitalSliceArray(axes(ρ.gs))
 # special case for integrator solver
-(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0; params...) = ρ.solver(mu, kBT)(; params...) |>
-    maybe_OrbitalSliceArray(axes(ρ.gs))
+(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0, override = missing; params...) =
+    ρ.solver(mu, kBT, override)(; params...) |> maybe_OrbitalSliceArray(axes(ρ.gs))
 
-(s::DensityMatrixIntegratorSolver)(mu, kBT) = s.ifunc(mu, kBT);
+(s::DensityMatrixIntegratorSolver)(mu, kBT, override = missing) = s.ifunc(mu, kBT, override);
 
 # redirects to specialized method
 densitymatrix(gs::GreenSlice; kw...) =
@@ -506,8 +506,12 @@ function densitymatrix(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), ims
     result = similar_Matrix(gs)
     opts´ = (; imshift, slope = 1, post = gf_to_rho!, atol, opts...)
     ωpoints_vec = collect(promote_type(T, typeof.(ωpoints)...), ωpoints)
-    ifunc(mu, kBT) = Integrator(result, DensityMatrixIntegrand(gs, T(mu), T(kBT), omegamap),
-        maybe_insert_mu!(ωpoints_vec, ωpoints, mu, kBT); opts´...)
+    function ifunc(mu, kBT, override)
+        ωpoints´ = override_path!(override, ωpoints_vec, ωpoints)
+        ρd = DensityMatrixIntegrand(gs, T(mu), T(kBT), omegamap)
+        pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, zero(T), kBT)
+        return Integrator(result, ρd, pts; opts´...)
+    end
     return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs)
 end
 
@@ -535,6 +539,18 @@ end
 
 maybe_insert_mu!(pts, mu, kBT) = pts
 
+override_path!(::Missing, ptsvec, pts) = pts
+override_path!(::Missing, ptsvec::Vector{<:Complex}, pts) = pts
+override_path!(f::Function, ptsvec::Vector{<:Complex}, pts) = f.(pts)
+
+function override_path!(pts´, ptsvec::Vector{<:Complex}, pts)
+    resize!(ptsvec, length(pts))
+    return pts´
+end
+
+override_path!(override, ptsvec, pts) =
+    argerror("Override of real ωpoints not supported, use complex ωpoints upon construction")
+
 #endregion
 
 #region ## API ##
@@ -608,9 +624,10 @@ end
 # generic fallback (for other solvers)
 (j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...)
 # special case for integrator solver
-(j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0; params...) = j.solver(kBT)(; params...)
+(j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0, override = missing; params...) =
+    j.solver(kBT, override)(; params...)
 
-(s::JosephsonIntegratorSolver)(kBT) = s.ifunc(kBT)
+(s::JosephsonIntegratorSolver)(kBT, override = missing) = s.ifunc(kBT, override)
 
 josephson(gs::GreenSlice{T}, ωmax::Number; kw...) where {T} =
     josephson(gs, (-ωmax, ωmax); kw...)
@@ -627,10 +644,11 @@ function josephson(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), phases
     phases´, traces = sanitize_phases_traces(phases, T)
     opts´ = (; imshift, slope = 1, post = real, atol, opts...)
     ωpoints_vec = collect(promote_type(T, typeof.(ωpoints)...), ωpoints)
-    function ifunc(kBT)
+    function ifunc(kBT, override)
+        ωpoints´ = override_path!(override, ωpoints_vec, ωpoints)
         jd = JosephsonIntegrand(g, T(kBT), contact, tauz, phases´, omegamap,
             traces, Σfull, Σ, similar(Σ), similar(Σ), similar(Σ), similar(tauz, Complex{T}))
-        pts = maybe_insert_mu!(ωpoints_vec, ωpoints, zero(T), kBT)
+        pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, zero(T), kBT)
         return Integrator(traces, jd, pts; opts´...)
     end
     return Josephson(JosephsonIntegratorSolver(ifunc), gs)

test/test_greenfunction.jl

@@ -355,6 +355,10 @@ function testjosephson(g0)
     @test all(((j1, j2) -> ≈(j1, j2, atol = 0.000001)).(j1, J3()))
     @test all(((j1, j2) -> ≈(j1, j2, atol = 0.000001)).(j1, J4()))
     @test all(((j1, j2) -> ≈(j1, j2, atol = 0.000001)).(j1, J5()))
+    # path override
+    j5 = J5(0.2)
+    j5´ = J5(0.2, (-4,-2+0.1im,0) .+ sqrt(eps(Float64))*im)
+    @test j5 != j5´ && j5 ≈ j5´
     j = Quantica.integrand(J1)
     @test Quantica.call!(j, 0.2+0.3im) isa Vector
     @test typeof.(Quantica.path.((J1, J2, J3, J4, J5))) ==