diff --git a/docs/src/tutorial/observables.md b/docs/src/tutorial/observables.md index 9440e739..61b0c601 100644 --- a/docs/src/tutorial/observables.md +++ b/docs/src/tutorial/observables.md @@ -38,16 +38,16 @@ Note that `d[sites...]` produces a vector with the LDOS at sites defined by `sit We can also compute the convolution of the density of states with the Fermi distribution `f(ω)=1/(exp((ω-μ)/kBT) + 1)`, which yields the density matrix in thermal equilibrium, at a given temperature `kBT` and chemical potential `μ`. This is computed with `ρ = densitymatrix(gs, (ωmin, ωmax))`. Here `gs = g[sites...]` is a `GreenSlice`, and `(ωmin, ωmax)` are integration bounds (they should span the full bandwidth of the system). Then, `ρ(µ, kBT = 0; params...)` will yield a matrix over the selected `sites` for a set of model `params`. ```julia julia> ρ = densitymatrix(g[region = RP.circle(1)], (-0.1, 8.1)) -DensityMatrix: density matrix on specified sites using solver of type DensityMatrixIntegratorSolver +DensityMatrix{DensityMatrixIntegratorSolver}: density matrix on specified sites julia> @time ρ(4) - 6.150548 seconds (57.84 k allocations: 5.670 GiB, 1.12% gc time) -5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: - 0.5+0.0im -7.34893e-10-3.94035e-15im 0.204478+1.9366e-14im -7.34889e-10-1.44892e-15im -5.70089e-10+5.48867e-15im - -7.34893e-10+3.94035e-15im 0.5+0.0im 0.200693-2.6646e-14im -5.70089e-10-1.95251e-15im -7.34891e-10-2.13804e-15im - 0.204478-1.9366e-14im 0.200693+2.6646e-14im 0.5+0.0im 0.200693+3.55692e-14im 0.204779-4.27255e-14im - -7.34889e-10+1.44892e-15im -5.70089e-10+1.95251e-15im 0.200693-3.55692e-14im 0.5+0.0im -7.34885e-10-3.49861e-15im - -5.70089e-10-5.48867e-15im -7.34891e-10+2.13804e-15im 0.204779+4.27255e-14im -7.34885e-10+3.49861e-15im 0.5+0.0im + 4.594645 seconds (111.82 k allocations: 4.890 GiB, 2.81% gc time, 0.86% compilation time) +5×5 OrbitalSliceMatrix{ComplexF64,Array}: + 0.5+0.0im -7.35075e-10+3.40256e-15im 0.204478+4.46023e-14im -7.35077e-10-1.13342e-15im -5.70426e-10-2.22213e-15im + -7.35075e-10-3.40256e-15im 0.5+0.0im 0.200693-4.46528e-14im -5.70431e-10+3.53853e-15im -7.35092e-10-4.07992e-16im + 0.204478-4.46023e-14im 0.200693+4.46528e-14im 0.5+0.0im 0.200693+6.7793e-14im 0.204779-6.78156e-14im + -7.35077e-10+1.13342e-15im -5.70431e-10-3.53853e-15im 0.200693-6.7793e-14im 0.5+0.0im -7.351e-10+2.04708e-15im + -5.70426e-10+2.22213e-15im -7.35092e-10+4.07992e-16im 0.204779+6.78156e-14im -7.351e-10-2.04708e-15im 0.5+0.0im ``` Note that the diagonal is `0.5`, indicating half-filling. @@ -55,33 +55,34 @@ Note that the diagonal is `0.5`, indicating half-filling. The default algorithm used here is slow, as it relies on numerical integration in the complex plane. Some GreenSolvers have more efficient implementations. If they exist, they can be accessed by omitting the `(ωmin, ωmax)` argument. For example, using `GS.Spectrum`: ```julia julia> @time g = h |> greenfunction(GS.Spectrum()); - 37.638522 seconds (105 allocations: 2.744 GiB, 0.79% gc time) + 18.249136 seconds (75 allocations: 1.567 GiB, 0.67% gc time) julia> ρ = densitymatrix(g[region = RP.circle(1)]) -DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver +DensityMatrix{DensityMatrixSpectrumSolver}: density matrix on specified sites + +julia> @time ρ(4) # second-run timing + 0.029662 seconds (7 allocations: 688 bytes) +5×5 OrbitalSliceMatrix{ComplexF64,Array}: + 0.5+0.0im -6.6187e-16+0.0im 0.204478+0.0im 2.49658e-15+0.0im -2.6846e-16+0.0im + -6.6187e-16+0.0im 0.5+0.0im 0.200693+0.0im -2.01174e-15+0.0im 1.2853e-15+0.0im + 0.204478+0.0im 0.200693+0.0im 0.5+0.0im 0.200693+0.0im 0.204779+0.0im + 2.49658e-15+0.0im -2.01174e-15+0.0im 0.200693+0.0im 0.5+0.0im 1.58804e-15+0.0im + -2.6846e-16+0.0im 1.2853e-15+0.0im 0.204779+0.0im 1.58804e-15+0.0im 0.5+0.0im -julia> @time ρ(4) - 0.001659 seconds (9 allocations: 430.906 KiB) -5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: - 0.5+0.0im -2.21437e-15+0.0im 0.204478+0.0im 2.67668e-15+0.0im 3.49438e-16+0.0im - -2.21437e-15+0.0im 0.5+0.0im 0.200693+0.0im -1.40057e-15+0.0im -2.92995e-15+0.0im - 0.204478+0.0im 0.200693+0.0im 0.5+0.0im 0.200693+0.0im 0.204779+0.0im - 2.67668e-15+0.0im -1.40057e-15+0.0im 0.200693+0.0im 0.5+0.0im 1.81626e-15+0.0im - 3.49438e-16+0.0im -2.92995e-15+0.0im 0.204779+0.0im 1.81626e-15+0.0im 0.5+0.0im ``` Note, however, that the computation of `g` is much slower in this case, due to the need of a full diagonalization. A better algorithm choice in this case is `GS.KPM`. It requires, however, that we define the region for the density matrix beforehand, as a `nothing` contact. ```julia julia> @time g = h |> attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 10000, bandrange = (0,8))); -Computing moments: 100%|█████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01 - 2.065083 seconds (31.29 k allocations: 11.763 MiB) +Computing moments: 100%|███████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01 + 1.360412 seconds (51.17 k allocations: 11.710 MiB) julia> ρ = densitymatrix(g[1]) -DensityMatrix: density matrix on specified sites with solver of type DensityMatrixKPMSolver +DensityMatrix{DensityMatrixKPMSolver}: density matrix on specified sites julia> @time ρ(4) - 0.006580 seconds (3 allocations: 1.156 KiB) -5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: + 0.004024 seconds (3 allocations: 688 bytes) +5×5 OrbitalSliceMatrix{ComplexF64,Array}: 0.5+0.0im 2.15097e-17+0.0im 0.20456+0.0im 2.15097e-17+0.0im 3.9251e-17+0.0im 2.15097e-17+0.0im 0.5+0.0im 0.200631+0.0im 1.05873e-16+0.0im 1.70531e-18+0.0im 0.20456+0.0im 0.200631+0.0im 0.5+0.0im 0.200631+0.0im 0.20482+0.0im @@ -89,6 +90,9 @@ julia> @time ρ(4) 3.9251e-17+0.0im 1.70531e-18+0.0im 0.20482+0.0im 1.70531e-18+0.0im 0.5+0.0im ``` +!!! note "Alternative integration paths" + The integration algorithm allows many different integration paths that can be adjusted to each problem, see the `Paths` docstring. Another versatile choice is `Paths.radial(ωrate, ϕ)`. This one is called with `ϕ = π/4` when doing `ρ = densitymatrix(gs::GreenSlice, ωrate::Number)`. In the example above this is slightly faster than the `(ωmin, ωmax)` choice, which resorts to `Paths.sawtooth(ωmin, ωmax)`. + ## Current A similar computation can be done to obtain the current density, using `J = current(g(ω), direction = missing)`. This time `J[sᵢ, sⱼ]` yields a sparse matrix of current densities along a given direction for each hopping (or the current norm if `direction = missing`). Passing `J` as a hopping shader yields the equilibrium current in a system. In the above example we can add a magnetic flux to make this current finite @@ -238,14 +242,7 @@ GreenFunction{Float64,2,0}: Green function of a Hamiltonian{Float64,2,0} Coordination : 3.7448 julia> J = josephson(g[1], 4.1) -Integrator: Complex-plane integrator - Integration path : (-4.1 + 1.4901161193847656e-8im, -2.05 + 2.050000014901161im, 0.0 + 1.4901161193847656e-8im) - Integration options : (atol = 1.0e-7,) - Integrand: : - JosephsonIntegrand{Float64} : Equilibrium (dc) Josephson current observable before integration over energy - kBT : 0.0 - Contact : 1 - Number of phase shifts : 0 +Josephson{JosephsonIntegratorSolver}: equilibrium Josephson current at a specific contact julia> qplot(g, children = (; sitecolor = :blue)) ``` @@ -257,37 +254,31 @@ In this case we have chosen to introduce the superconducting leads with a model corresponding to a BCS bulk, but any other self-energy form could be used. We have introduced the phase difference (`phase`) as a model parameter. We can now evaluate the zero-temperature Josephson current simply with ```julia julia> J(phase = 0) --1.974396994480587e-16 +1.992660837638158e-12 julia> J(phase = 0.2) -0.004617597139699372 +0.0046175971391935605 ``` Note that finite temperatures can be taken using the `kBT` keyword argument for `josephson`, see docstring for details. One is often interested in the critical current, which is the maximum of the Josephson current over all phase differences. Quantica.jl can compute the integral over a collection of phase differences simultaneously, which is more efficient that computing them one by one. This is done with ```julia julia> φs = subdiv(0, pi, 11); J = josephson(g[1], 4.1; phases = φs) - Integration path : (-4.1 + 1.4901161193847656e-8im, -2.05 + 2.050000014901161im, 0.0 + 1.4901161193847656e-8im) - Integration options : (atol = 1.0e-7,) - Integrand: : - JosephsonIntegrand{Float64} : Equilibrium (dc) Josephson current observable before integration over energy - kBT : 0.0 - Contact : 1 - Number of phase shifts : 11 +Josephson{JosephsonIntegratorSolver}: equilibrium Josephson current at a specific contact julia> Iφ = J() 11-element Vector{Float64}: - 1.868862401627357e-14 - 0.007231421775452674 - 0.014242855188877 - 0.02081870760779799 - 0.026752065104401878 - 0.031847203848574666 - 0.0359131410974842 - 0.03871895510547465 - 0.039762442694035505 - 0.03680096751905469 - 2.7677727119798235e-14 + -6.361223111882911e-13 + 0.007231421776215144 + 0.01424285518831463 + 0.020818707606469377 + 0.026752065101976884 + 0.031847203846513975 + 0.035913141096514584 + 0.038718955102068034 + 0.03976244268586444 + 0.036800967573567184 + -1.437196514806921e-12 julia> f = Figure(); a = Axis(f[1,1], xlabel = "φ", ylabel = "I [e/h]"); lines!(a, φs, Iφ); scatter!(a, φs, Iφ); f ``` diff --git a/src/docstrings.jl b/src/docstrings.jl index 204e3358..5e02e9c8 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -2050,15 +2050,18 @@ true if `gs` contains user-defined contacts created using `attach(model)` with g ω-dependent models, but note that Quantica will not check for analyticity. Keywords `quadgk_opts` are passed to the `QuadGK.quadgk` integration routine. - densitymatrix(gs::GreenSlice, ωscale::Real; opts..., quadgk_opts...) + densitymatrix(gs::GreenSlice, ωscale::Real; kw...) -As above with `path = Paths.vertical(ωscale)`, which computes the density matrix by -integrating the Green function along a quarter circle from `ω = -Inf` to `ω = µ + im*Inf`, -and then back along a vertical path `ω = µ + im*ωscale*x` from `x = Inf` to `x = 0`. Here -`µ` is the chemical potential and `ωscale` should be some typical energy scale in the system -(the integration time may depend on `ωscale`, but not the result). This vertical integration -path approach is performant but, currently it does not support finite temperatures, unlike -other `AbstractIntegrationPath`s like `Paths.sawtooth`. +As above with `path = Paths.radial(ωscale, π/4)`, which computes the density matrix by +integrating the Green function from `ω = -Inf` to `ω = Inf` along a `ϕ = π/4` radial complex +path, see `Paths` for details. Here `ωscale` should correspond to some typical energy scale +in the system, which dictates the speed at which we integrate the radial paths (the +integration runtime may depend on `ωscale`, but not the result). + + densitymatrix(gs::GreenSlice, ωs::NTuple{N,Real}; kw...) + +As above, but with a `path = Paths.sawtooth(ωs)`, which uses a sawtooth-shaped path touching +points `ωs` on the real axes. Ideally, these values should span the full system bandwidth. ## Full evaluation @@ -2068,9 +2071,6 @@ Evaluate the density matrix at chemical potential `μ` and temperature `kBT` (in 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. -If the generic integration algorithm is used with complex `ωpoints`, the following form is -also available: - ## Algorithms and keywords The generic integration algorithm allows for the following `opts` (see also `josephson`): @@ -2122,13 +2122,16 @@ GʳΣʳᵢ)τz]``. Here `f(ω)` is the Fermi function with `µ = 0`. josephson(gs::GreenSlice, ωscale::Real; kw...) -As above, but with `path = Paths.vertical(ωscale)`, which computes josephson current by -performing the above integral along a quarter circle from `ω = -Inf` to `ω = im*Inf`, and -then back along a vertical path `ω = im*ωscale*x` from `x = Inf` to `x = 0`. Here `ωscale` -should be some typical energy scale in the system, like the superconducting gap (the -integration time may depend on `ωscale`, but not the result). This vertical integration path -approach is performant but, currently it does not support finite temperatures, unlike other -`AbstractIntegrationPath`s like `Paths.sawtooth`. +As above, but with `path = Paths.radial(ωscale, π/4)`, which computes josephson current by +integrating the Green function from `ω = -Inf` to `ω = Inf` along a `ϕ = π/4` radial complex +path, see `Paths` for details. Here `ωscale` should be some typical energy scale in the +system, like the superconducting gap (the integration time may depend on `ωscale`, but not +the result). + + josephson(gs::GreenSlice, ωs::NTuple{N,Real}; kw...) + +As above, but with a `path = Paths.sawtooth(ωs)`, which uses a sawtooth-shaped path touching +points `ωs` on the real axes. Ideally, these values should span the full system bandwidth. ## Full evaluation @@ -2155,21 +2158,21 @@ julia> glead = LP.square() |> hamiltonian(@onsite((; ω = 0) -> 0.0005 * SA[0 1; julia> g0 = LP.square() |> hamiltonian(hopping(SA[1 0; 0 -1]), orbitals = 2) |> supercell(region = r->-2 attach(glead, reverse = true) |> attach(glead) |> greenfunction; -julia> J = josephson(g0[1], 4; omegamap = ω -> (;ω), phases = subdiv(0, pi, 10)) +julia> J = josephson(g0[1], 0.0005; omegamap = ω -> (;ω), phases = subdiv(0, pi, 10)) Josephson: equilibrium Josephson current at a specific contact using solver of type JosephsonIntegratorSolver julia> J(0.0) 10-element Vector{Float64}: - 7.060440509787806e-18 - 0.0008178484258721882 - 0.0016108816082772972 - 0.002355033150366814 - 0.0030277117620820513 - 0.003608482493380227 - 0.004079679643085058 - 0.004426918320990192 - 0.004639358112465513 - 2.2618383948099795e-12 + 1.0130103834038537e-15 + 0.0008178802022977883 + 0.0016109471548779466 + 0.0023551370133118215 + 0.0030278625151304614 + 0.003608696305848759 + 0.00407998998248311 + 0.0044274100715435295 + 0.004640372460465891 + 5.179773035314182e-12 ``` # See also @@ -2178,59 +2181,69 @@ julia> J(0.0) josephson """ - Quantica.integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0) + Quantica.integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0; params...) -Return the complex integrand `d::JosephsonIntegrand` whose integral over frequency yields the -Josephson current, `J(kBT) = Re(∫dω d(ω))`. To evaluate the `d` for a -given `ω` and parameters, use `d(ω; params...)`, or `call!(d, ω; params...)` for its +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. - Quantica.integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0, kBT = 0) + Quantica.integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0, kBT = 0; params...) -Like above for the density matrix `ρ(mu, kBT)`, with `d::DensityMatrixIntegrand`, so that -`ρ(mu, kBT) = -mat_imag(∫dω d(ω))/π`, and `mat_imag(GF::Matrix) = (GF - GF')/2im`. +Like above for the density matrix `ρ`, with `d::DensityMatrixIntegrand`, so that `ρ(mu, kBT; +params...) = -mat_imag(∫dω d(ω; params...))/π`, and `mat_imag(GF::Matrix) = (GF - GF')/2im`. """ integrand """ - Quantica.path(O::Josephson, args...) - Quantica.path(O::DensityMatrix, args...) + Quantica.points(O::Josephson, args...) + Quantica.points(O::DensityMatrix, args...) -Return the vertices of the polygonal integration path used to compute `O(args...)`. +Return the vertices of the integration path used to compute `O(args...)`. """ -path +points """ Paths A Quantica submodule that contains representations of different integration paths in the -complex-ω plane. Available paths are: +complex-ω plane for integrals of the form `∫f(ω)g(ω)dω`, where `f(ω)` is the Fermi +distribution at chemical potential `µ` and temperature `kBT`. Available paths are: - Paths.vertical(ωscale::Real) + Paths.radial(ωscale::Real, ϕ) -A vertical path from a point `µ` on the real axis to `µ + Inf*im`. `ωscale` defines the rate at which the integration variable `x` sweeps this ω path, i.e. `ω = µ+im*ωscale*x`. +A four-segment path from `ω = -Inf` to `ω = Inf`. The path first ascends along an arc of +angle `ϕ` (with `0 <= ϕ < π/2`) and infinite radius. It then converges along a straight line +in the second `ω-µ` quadrant to the chemical potential origin `ω = µ`. Finally it performs +the mirror-symmetric itinerary in the first `ω-µ` quadrant, ending on the real axis at `ω = +Inf`. The radial segments are traversed at a rate dictated by `ωscale`, that should +represent some relevant energy scale in the system for optimal convergence of integrals. At +zero temperature `kBT = 0`, the two last segments don't contribute and are elided. - Paths.sawtooth(ωpoints; slope = 1) + Paths.sawtooth(ωpoints...; slope = 1, imshift = true) A path connecting all points in the `ωpoints` collection (should all be real), but departing between each two into the complex plane and back with the given `slope`, forming a sawtooth path. Note that an extra point `µ` is added to the collection, where `µ` is the chemical -potential. This path allows finite temperature integrals +potential, and the `real(ω)>µ` segments are elided at zero temperatures. If `imshift = true` +all `ωpoints` (of type `T<:Real`) are shifted by `sqrt(eps(T))` to avoid the real axis. - Paths.sawtooth(ωmax::Real; slope = 1) + Paths.sawtooth(ωmax::Real; kw...) As above with `ωpoints = (-ωmax, ωmax)`. - Paths.polygon(ωpoints) + Paths.polygon(ωpoints...) -A ageneral polygonal path connecting `ωpoints`, which can be any collection of real or +A general polygonal path connecting `ωpoints`, which can be any collection of real or complex numbers Paths.polygon(ωfunc::Function) -A ageneral polygonal path connecting `ωpoints = ωfunc(µ, kBT)`. This is useful when the -desired integration path depends on the chemical potential `µ` or temperature `kBT`. +A ageneral polygonal path connecting `ωpoints = ωfunc(µ, kBT; params...)`. This is useful +when the desired integration path depends on the chemical potential `µ`, temperature +`kBT` or other system parameters `params`. # See also `densitymatrix`, `josephson` diff --git a/src/greenfunction.jl b/src/greenfunction.jl index 47f31900..d57a75cc 100644 --- a/src/greenfunction.jl +++ b/src/greenfunction.jl @@ -105,12 +105,16 @@ function Base.getindex(g::GreenSolution, args...; post = identity, symmetrize = end Base.getindex(g::GreenSolution, i::AnyCellOrbitals, j::AnyCellOrbitals) = - slicer(g)[sanitize_cellorbs(i), sanitize_cellorbs(j)] + slicer(g)[sanitize_cellorbs(i, g), sanitize_cellorbs(j, g)] # must ensure that orbindices is not a scalar, to consistently obtain a Matrix +# also, it should not be Colon, as GreenSlicers generally assume an explicit range +sanitize_cellorbs(c::CellOrbitals, g) = sites_to_orbs(c, g) # convert : to range sanitize_cellorbs(c::CellOrbitals) = c -sanitize_cellorbs(c::CellOrbital) = CellOrbitals(cell(c), orbindex(c):orbindex(c)) -sanitize_cellorbs(c::CellOrbitalsGrouped) = CellOrbitals(cell(c), orbindices(c)) +sanitize_cellorbs(c::CellOrbital, _...) = CellOrbitals(cell(c), orbindex(c):orbindex(c)) +sanitize_cellorbs(c::CellOrbitalsGrouped, _...) = CellOrbitals(cell(c), orbindices(c)) +sanitize_cellorbs(::CellOrbitals{<:Any,Colon}) = + internalerror("sanitize_cellorbs: Colon indices leaked!") ## common getindex! shortcut in terms of GreenSlice @@ -125,55 +129,6 @@ orbinds_or_contactinds( c::Union{Colon,Integer,DiagIndices{Colon},DiagIndices{<:Integer}}, _, _) = (r, c) orbinds_or_contactinds(_, _, or, oc) = (or, oc) -## symmetrization functions -# If symmetrize == sym::Number, we do output .= post(sym*gij + conj(sym)*gji'). -# Otherwise just output .= post(gij) - -function maybe_symmetrize!(out, gij, ::Missing; post = identity) - out .= post.(gij) - return out -end - -function maybe_symmetrize!(out, (gij, gji´), sym::Number; post = identity) - out .= post.(sym .* gij .+ conj(sym) .* gji´) - return out -end - -# assumes mat is square -maybe_symmetrize(mat, ::Missing) = mat -maybe_symmetrize(mat, sym::Number) = maybe_symmetrize!(mat, (mat, copy(mat')), sym) - -maybe_symmetrized_view!(output, g, i, j, ::Missing; kw...) = - (maybe_symmetrize!(parent(output), missing, view(g, i, j); kw...); output) - -maybe_symmetrized_view!(output, g, i, j, sym, ; kw...) = - (maybe_symmetrize!(parent(output), (view(g, i, j), view(g, j, i)'), sym; kw...); output) - -maybe_symmetrized_getindex!(output, g, i, j, ::Missing; kw...) = - maybe_symmetrize!(output, g[i, j], missing; kw...) - -function maybe_symmetrized_getindex!(output, g, i, j, sym; kw...) - gij = g[i, j] - gji´ = i == j ? gij' : g[j, i]' - return maybe_symmetrize!(output, (gij, gji´), sym; kw...) -end - -maybe_symmetrized_getindex(g, i, j, ::Missing; post = identity) = post.(g[i, j]) - -function maybe_symmetrized_getindex(g, i, j, sym; kw...) - gij = g[i, j] # careful, this could be non-mutable - if i == j - gij = maybe_symmetrize!(gij, (gij, copy(gij')), sym; kw...) - else - gij = maybe_symmetrize!(gij, (gij, g[j, i]'), sym; kw...) - end - return gij -end - -# in case gij above is non-mutable -maybe_symmetrize!(::StaticArray, (gij, gji´), sym::Number; post = identity) = - post.(sym * gij + conj(sym) * gji´) - ## getindex! - core functions # fastpath for intra and inter-contact, only for g::GreenSolution @@ -184,32 +139,35 @@ getindex!(output, g::GreenSolution, i::Union{Integer,Colon}, j::Union{Integer,Co getindex!(output, g, i::AnyCellOrbitals, j::AnyCellOrbitals; symmetrize = missing, kw...) = maybe_symmetrized_getindex!(output, g, i, j, symmetrize; kw...) -# indexing over several cells -getindex!(output, g, ci::AnyOrbitalSlice, cj::AnyOrbitalSlice; kw...) = - getindex_cells!(output, g, cellsdict(ci), cellsdict(cj); kw...) -getindex!(output, g, ci::AnyOrbitalSlice, cj::AnyCellOrbitals; kw...) = - getindex_cells!(output, g, cellsdict(ci), (cj,); kw...) -getindex!(output, g, ci::AnyCellOrbitals, cj::AnyOrbitalSlice; kw...) = - getindex_cells!(output, g, (ci,), cellsdict(cj); kw...) +# fallback for multicell case +getindex!(output, g, ci::Union{AnyOrbitalSlice,AnyCellOrbitals}, cj::Union{AnyOrbitalSlice,AnyCellOrbitals}; kw...) = + getindex_cells!(output, g, cellinds_iterable_axis(ci), cellinds_iterable_axis(cj); kw...) -function getindex_cells!(output, g, cis, cjs; kw...) +cellinds_iterable_axis(ci::CellIndices) = ((ci,), missing) +cellinds_iterable_axis(ci::AnyOrbitalSlice) = (cellsdict(ci), ci) + +# at this point, either cis or cjs is a multi-cell iterator, so an output view is needed +function getindex_cells!(output, g, (cis, iaxis), (cjs, jaxis); kw...) for ci in cis, cj in cjs - getindex!(output, g, ci, cj; kw...) # will typically call the method below + rows, cols = cell_orbindices(ci, iaxis), cell_orbindices(cj, jaxis) + output´ = view(maybe_orbmat_parent(output), rows, cols) + getindex!(output´, g, ci, cj; kw...) # will call the single-cell method above end return output end -function getindex!(output::OrbitalSliceMatrix, g, i::AnyCellOrbitals, j::AnyCellOrbitals; kw...) - oi, oj = orbaxes(output) - rows, cols = orbrange(oi, cell(i)), orbrange(oj, cell(j)) - v = view(parent(output), rows, cols) - getindex!(v, g, i, j; kw...) - return output -end +# output may be a Matrix (inhomogenous rows/cols) or an OrbitalSliceMatrix (homogeneous) +maybe_orbmat_parent(output::OrbitalSliceMatrix) = parent(output) +maybe_orbmat_parent(output) = output + +# orbital index range in parent(output::OrbitalSliceMatrix) for a given cell +cell_orbindices(ci, axis) = orbrange(axis, cell(ci)) +# if output::Matrix, we take all indices, since output was already cropped along this axis +cell_orbindices(_, ::Missing) = Colon() ## GreenSlicer -> Matrix, dispatches to each solver's implementation -# fallback conversion to CellOrbitals +# simplify CellOrbitalGroups to CellOrbitals. No CellOrbitals{L,Colon} should get here Base.getindex(s::GreenSlicer, i::AnyCellOrbitals, j::AnyCellOrbitals = i; kw...) = getindex(s, sanitize_cellorbs(i), sanitize_cellorbs(j); kw...) @@ -241,8 +199,10 @@ getindex!(output, g, i::OrbitalPairIndices, j::OrbitalPairIndices = i; kw...) = # should also implement a `blockstructure` method so we know how to map sites to indices. #region -getindex_diagonal!(output, g::GreenSolution, i::Union{Colon,Integer}, ker; kw...) = - fill_diagonal!(output, view(g, i, i), contact_kernel_ranges(ker, i, g), ker; post) +function getindex_diagonal!(output, g::GreenSolution, i::Union{Colon,Integer}, ker; symmetrize = missing, kw...) + v = maybe_symmetrized_matrix(view(g, i, i), symmetrize) + return fill_diagonal!(output, v, contact_kernel_ranges(ker, i, g), ker; kw...) +end function getindex_diagonal!(output, g, o::CellOrbitalsGrouped, ker; symmetrize = missing, kw...) rngs = orbital_kernel_ranges(ker, o) @@ -559,19 +519,24 @@ function TMatrixSlicer(g0slicer::GreenSlicer{C}, Σblocks, contactorbs) where {C else nreg = norbitals(contactorbs) g0contacts = zeros(C, nreg, nreg) - off = offsets(contactorbs) - for (j, sj) in enumerate(cellsdict(contactorbs)), (i, si) in enumerate(cellsdict(contactorbs)) - irng = off[i]+1:off[i+1] - jrng = off[j]+1:off[j+1] - g0view = view(g0contacts, irng, jrng) - copy!(g0view, g0slicer[si, sj]) - end + fill_g0contacts!(g0contacts, g0slicer, contactorbs) Σblocks´ = tupleflatten(Σblocks...) tmatrix, gcontacts = t_g_matrices(g0contacts, contactorbs, Σblocks´...) end return TMatrixSlicer(g0slicer, tmatrix, gcontacts, contactorbs) end +function fill_g0contacts!(mat, g0slicer, contactorbs) + off = offsets(contactorbs) + for (j, sj) in enumerate(cellsdict(contactorbs)), (i, si) in enumerate(cellsdict(contactorbs)) + irng = off[i]+1:off[i+1] + jrng = off[j]+1:off[j+1] + g0view = view(mat, irng, jrng) + copy!(g0view, g0slicer[si, sj]) + end + return mat +end + # Takes a precomputed g0contacts (for dummy g0slicer that doesn't implement indexing) function TMatrixSlicer(g0contacts::AbstractMatrix{C}, Σblocks, contactorbs) where {C} tmatrix, gcontacts = t_g_matrices(g0contacts, contactorbs, Σblocks...) diff --git a/src/integrator.jl b/src/integrator.jl index fc79a9ad..581e628e 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -10,54 +10,74 @@ end struct SawtoothPath{T<:Real} <: AbstractIntegrationPath realpts::Vector{T} - pts::Vector{Complex{T}} slope::T + imshift::Bool end # straight from 0 to im*∞, but using a real variable from 0 to inf as parameter # this is because QuadGK doesn't understand complex infinities -struct VerticalPath{T<:AbstractFloat} <: AbstractIntegrationPath +struct RadialPath{T<:AbstractFloat} <: AbstractIntegrationPath rate::T + angle::T + cisinf::Complex{T} end #region ## Constructors ## -SawtoothPath(pts::Vector{T}, slope) where {T<:Real} = - SawtoothPath(pts, complex.(pts), T(slope)) +SawtoothPath(pts::Vector{T}, slope, imshift) where {T<:Real} = + SawtoothPath(pts, T(slope), imshift) + +RadialPath(rate, angle) = RadialPath(promote(float(rate), float(angle))...) + +function RadialPath(rate::T, angle::T) where {T<:AbstractFloat} + 0 <= angle < π/2 || argerror("The radial angle should be in the interval [0, π/2), got $angle") + RadialPath(rate, angle, cis(angle)) +end #endregion #region ## API ## -points(p::PolygonPath, args...) = p.pts -points(p::PolygonPath{<:Function}, mu, kBT) = p.pts(mu, kBT) +points(p::PolygonPath, args...; _...) = p.pts +points(p::PolygonPath{<:Function}, mu, kBT; params...) = p.pts(mu, kBT; params...) + +points(p::RadialPath{T}, mu, kBT; _...) where {T} = iszero(kBT) ? + [mu - T(Inf)*conj(p.cisinf), maybe_imshift(mu)] : + [mu - T(Inf)*conj(p.cisinf), maybe_imshift(mu), mu + T(Inf)*p.cisinf] -points(::VerticalPath{T}, mu, kBT) where {T} = (mu + zero(T)*im, mu + T(Inf)*im) +function points(p::SawtoothPath, mu, kBT; _...) + pts = complex(p.realpts) + if (maximum(real, pts) > real(mu) > minimum(real, pts) && !any(≈(real(mu)), pts)) + # insert µ if it is within extrema and not already included + sort!(push!(pts, mu), by = real) + iszero(kBT) && cut_tail!(x -> real(x) > real(mu), pts) + end + p.imshift && imshift!(pts) + triangular_sawtooth!(pts, p.slope) + return pts +end -realpoints(::VerticalPath{T}, pts) where {T} = (zero(T), T(Inf)) +maybe_imshift(x::T) where {T<:Real} = float(x) + sqrt(eps(float(T)))*im +maybe_imshift(x::Complex) = float(x) + +imshift!(v::AbstractArray{Complex{T}}) where {T} = (v .+= sqrt(eps(T))*im; v) + +# Must be a type-stable tuple, and avoiding 0 in RadialPath +realpoints(p::RadialPath{T}, pts) where {T} = (-T(Inf), T(0), ifelse(length(pts)==2, T(0), T(Inf))) realpoints(::Union{SawtoothPath,PolygonPath}, pts) = 1:length(pts) -point(x, p::VerticalPath, pts) = first(pts) + p.rate*im*x +# note: pts[2] == µ +point(x::Real, p::RadialPath, pts) = pts[2] + p.rate*x*ifelse(x <= 0, conj(p.cisinf), p.cisinf) -function point(x, p::Union{SawtoothPath,PolygonPath}, pts) +function point(x::Real, p::Union{SawtoothPath,PolygonPath}, pts) x0 = floor(Int, x) p0, p1 = pts[x0], pts[ceil(Int, x)] return p0 + (x-x0)*(p1 - p0) end -# the minus sign inverts the VerticalPath, so it is transformed to run from Inf to 0 -jacobian(x, p::VerticalPath, pts) = -p.rate*im -jacobian(x, p::Union{SawtoothPath,PolygonPath}, pts) = pts[ceil(Int, x)] - pts[floor(Int, x)] - -function points(p::SawtoothPath, mu, kBT) - copy!(p.pts, p.realpts) - if !((iszero(kBT) && maximum(real, p.pts) <= mu) || any(≈(mu), p.pts)) - sort!(push!(p.pts, mu), by = real) - iszero(kBT) && cut_tail!(x -> real(x) <= mu, p.pts) - end - triangular_sawtooth!(p.pts, p.slope) - return p.pts -end +# the minus sign inverts the RadialPath, so it is transformed to run from Inf to 0 +jacobian(x::Real, p::RadialPath, pts) = p.rate*ifelse(x <= 0, conj(p.cisinf), p.cisinf) +jacobian(x::Real, p::Union{SawtoothPath,PolygonPath}, pts) = pts[ceil(Int, x)] - pts[floor(Int, x)] function triangular_sawtooth!(pts, slope) for i in 2:length(pts) @@ -78,24 +98,25 @@ end module Paths -using Quantica: Quantica, PolygonPath, SawtoothPath, VerticalPath +using Quantica: Quantica, PolygonPath, SawtoothPath, RadialPath, argerror ## API -sawtooth(ωpoints; kw...) = sawtooth!(float.(ωpoints); kw...) +sawtooth(ω1::Real, ω2::Real, ωs::Real...; kw...) = sawtooth(float.((ω1, ω2, ωs...)); kw...) sawtooth(ωmax::Real; kw...) = sawtooth!([-float(ωmax), float(ωmax)]; kw...) sawtooth(ωpoints::Tuple; kw...) = sawtooth!(collect(float.(ωpoints)); kw...) sawtooth(ωpoints::AbstractVector; kw...) = sawtooth!(float.(ωpoints); kw...) -sawtooth!(pts::AbstractVector{<:AbstractFloat}; slope = 1) = - SawtoothPath(sort!(pts), slope) +sawtooth!(pts::AbstractVector{<:AbstractFloat}; slope = 1, imshift = true) = + SawtoothPath(sort!(pts), slope, imshift) sawtooth!(pts; kw...) = argerror("Path.sawtooth expects a collection of real points, got $pts") -vertical(rate) = VerticalPath(float(rate)) +radial(rate, angle) = RadialPath(rate, angle) polygon(ωmax::Real) = polygon((-ωmax, ωmax)) -polygon(ωpoints) = PolygonPath(ωpoints) +polygon(ω1::Number, ω2::Number, ωs::Number...) = PolygonPath((ω1, ω2, ωs...)) +polygon(ωs) = PolygonPath(ωs) end # Module @@ -148,11 +169,11 @@ function call!(I::Integrator{<:Any,Missing}; params...) end # nonscalar version -function call!(I::Integrator; params...) +function call!(I::Integrator{<:Any,T}; params...) where {T} fx! = (y, x) -> begin y .= serialize(call!(I.integrand, x; params...)) I.callback(x, y) - return y + return nothing end result, err = quadgk!(fx!, serialize(I.result), points(I)...; I.quadgk_opts...) # note: post-processing is not element-wise & can be in-place diff --git a/src/observables.jl b/src/observables.jl index cce695ac..5b5924b4 100644 --- a/src/observables.jl +++ b/src/observables.jl @@ -4,8 +4,8 @@ abstract type IndexableObservable end # any type that can be sliced into (current and ldos) -fermi(ω::C, β = Inf) where {C} = - isinf(β) ? ifelse(real(ω) <= 0, C(1), C(0)) : C(1/(exp(β * ω) + 1)) +fermi(ω::C, β = Inf; atol = sqrt(eps(real(C)))) where {C} = + isinf(β) ? ifelse(abs(ω) < atol, C(0.5), ifelse(real(ω) <= 0, C(1), C(0))) : C(1/(exp(β * ω) + 1)) normal_size(h::AbstractHamiltonian) = normal_size(blockstructure(h)) @@ -391,15 +391,18 @@ densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) = argerror("Dedicated `densitymatrix` algorithm not implemented for $(nameof(typeof(s))). Use generic one instead.") # default integrator solver -densitymatrix(gs::GreenSlice, ωmax::Number; kw...) = - densitymatrix(gs::GreenSlice, Paths.vertical(ωmax); kw...) +densitymatrix(gs::GreenSlice, ωmax::Real; kw...) = + densitymatrix(gs::GreenSlice, Paths.radial(ωmax, π/4); kw...) + +densitymatrix(gs::GreenSlice, ωs::NTuple{<:Any,Real}; kw...) = + densitymatrix(gs::GreenSlice, Paths.sawtooth(ωs); kw...) function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap = Returns((;)), atol = 1e-7, opts...) where {T} result = copy(call!_output(gs)) post = post_transform_rho(path, gs) opts´ = (; post, atol, opts...) - function ifunc(mu, kBT) - pts = points(path, mu, kBT) + 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) return Integrator(result, ρd, realpts; opts´...) @@ -407,9 +410,9 @@ function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegama return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs) end -# we need to add the quarter-circle path arc from -∞ to +im*∞ -function post_transform_rho(::VerticalPath, gs) - arc = gs(0.5*I) +# we need to add the arc path segment from -∞ to ∞ * p.cisinf +function post_transform_rho(p::RadialPath, gs) + arc = gs((p.angle/π)*I) function post(x) x .+= arc return x @@ -428,18 +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...) + ρ.solver(mu, kBT; params...)(; params...) -(s::DensityMatrixIntegratorSolver)(mu, kBT) = - s.ifunc(mu, kBT); +(s::DensityMatrixIntegratorSolver)(mu, kBT; params...) = + s.ifunc(mu, kBT; params...); -(ρi::DensityMatrixIntegrand)(ω; params...) = copy(call!(ρi, ω; params...)) +(ρi::DensityMatrixIntegrand)(x; params...) = copy(call!(ρi, x; params...)) function call!(ρi::DensityMatrixIntegrand, x; params...) ω = point(x, ρi.path, ρi.pts) - symmetrize = -jacobian(x, ρi.path, ρi.pts)/(2π*im) + 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...) - serialize(output) .*= fermi(ω - ρi.mu, inv(ρi.kBT)) return output end @@ -447,9 +451,14 @@ end #region ## API ## -integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = integrand(ρ.solver(mu, kBT)) +integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0; params...) = + integrand(ρ.solver(mu, kBT; params...)) + +points(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0; params...) = + points(integrand(ρ, mu, kBT; params...)) +points(ρ::DensityMatrixIntegrand) = ρ.pts -path(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = path(ρ.solver(mu, kBT)) +point(x, ρi::DensityMatrixIntegrand) = point(x, ρi.path, ρi.pts) temperature(D::DensityMatrixIntegrand) = D.kBT @@ -508,8 +517,11 @@ end #region ## josephson API ## -josephson(gs::GreenSlice{T}, ωmax::Number; kw...) where {T} = - josephson(gs, Paths.vertical(ωmax); kw...) +josephson(gs::GreenSlice, ωmax::Real; kw...) = + josephson(gs, Paths.radial(ωmax, π/4); kw...) + +josephson(gs::GreenSlice, ωs::NTuple{<:Any,Real}; kw...) = + josephson(gs, Paths.sawtooth(ωs); kw...) function josephson(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap = Returns((;)), phases = missing, atol = 1e-7, opts...) where {T} check_nodiag_axes(gs) @@ -522,8 +534,8 @@ function josephson(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap = tauz = tauz_diag.(axes(Σ, 1), normalsize) phases´, traces = sanitize_phases_traces(phases, T) opts´ = (; post = real, atol, opts...) - function ifunc(kBT) - pts = points(path, 0, kBT) + 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, traces, Σfull, Σ, similar(Σ), similar(Σ), similar(Σ), similar(tauz, Complex{T})) @@ -548,20 +560,20 @@ end # generic fallback (for other solvers) (j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...) -# special case for integrator solver +# special case for integrator solver (so we can access integrand etc before integrating) (j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0; params...) = - j.solver(kBT)(; params...) + j.solver(kBT; params...)(; params...) -(s::JosephsonIntegratorSolver)(kBT) = s.ifunc(kBT) +(s::JosephsonIntegratorSolver)(kBT; params...) = s.ifunc(kBT; params...) -(J::JosephsonIntegrand)(ω; params...) = copy(call!(J, ω; params...)) +(J::JosephsonIntegrand)(x; params...) = copy(call!(J, x; params...)) -function call!(J::JosephsonIntegrand, x; params...) - ω = point(x, J.path, J.pts) - gω = call!(J.g, ω; J.omegamap(ω)..., params...) - f = fermi(ω, inv(J.kBT)) - traces = josephson_traces(J, gω, f) - traces .*= jacobian(x, J.path, J.pts) +function call!(Ji::JosephsonIntegrand, x; params...) + ω = point(x, Ji.path, Ji.pts) + gω = call!(Ji.g, ω; Ji.omegamap(ω)..., params...) + f = fermi(ω, inv(Ji.kBT)) + traces = josephson_traces(Ji, gω, f) + traces = mul_scalar_or_array!(traces, jacobian(x, Ji.path, Ji.pts)) return traces end @@ -569,9 +581,12 @@ end #region ## API ## -integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0) = integrand(J.solver(kBT)) +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::JosephsonIntegrand) = J.pts -path(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0) = path(J.solver(kBT)) +point(x, Ji::JosephsonIntegrand) = point(x, Ji.path, Ji.pts) temperature(J::JosephsonIntegrand) = J.kBT diff --git a/src/show.jl b/src/show.jl index 14a1ad5a..5021fd42 100644 --- a/src/show.jl +++ b/src/show.jl @@ -428,26 +428,26 @@ Base.summary(::Transmission) = #endregion -############################################################################################ -# Integrator -#region +# ############################################################################################ +# # Integrator +# #region -function Base.show(io::IO, I::Integrator) - i = get(io, :indent, "") - print(io, i, summary(I), "\n", -"$i Integration path : $(path(I)) -$i Integration options : $(display_namedtuple(options(I))) -$i Integrand : ") - ioindent = IOContext(io, :indent => i * " ") - show(ioindent, integrand(I)) -end +# function Base.show(io::IO, I::Integrator) +# i = get(io, :indent, "") +# print(io, i, summary(I), "\n", +# "$i Integration path : $(path(I)) +# $i Integration options : $(display_namedtuple(options(I))) +# $i Integrand : ") +# ioindent = IOContext(io, :indent => i * " ") +# show(ioindent, integrand(I)) +# end -Base.summary(::Integrator) = "Integrator: Complex-plane integrator" +# Base.summary(::Integrator) = "Integrator: Complex-plane integrator" -display_namedtuple(nt::NamedTuple) = isempty(nt) ? "()" : "$nt" +# display_namedtuple(nt::NamedTuple) = isempty(nt) ? "()" : "$nt" -#endregion +# #endregion ############################################################################################ # Josephson and DensityMatrix integrands diff --git a/src/slices.jl b/src/slices.jl index 9d5f84fe..f40e9c18 100644 --- a/src/slices.jl +++ b/src/slices.jl @@ -359,6 +359,11 @@ sites_to_orbs(c::SparseIndices{<:Any,<:Hamiltonian}, _) = c # unused # sites_to_orbs_nogroups(cs::CellOrbitals, _) = cs +## CellOrbitals{L,Colon} -> CellOrbitals with explicit index range + +sites_to_orbs(cs::CellOrbitals{<:Any,Colon}, g) = + sites_to_orbs_nogroups(sites(cell(cs), :), g) + ## DiagIndices -> DiagIndices function sites_to_orbs(d::DiagIndices, g) @@ -403,14 +408,13 @@ end ## convert CellSites -> CellOrbitalsGrouped -sites_to_orbs(c::CellSites, g, ker...) = - sites_to_orbs(sanitize_cellindices(c, g), blockstructure(g), ker...) +sites_to_orbs(c::CellSites, g) = + sites_to_orbs(sanitize_cellindices(c, g), blockstructure(g)) -function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure, ker...) - os´ = maybe_scalarize(os, ker...) +function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure) sites = siteindices(cs) - groups = _groups(sites, os´) # sites, orbranges - orbinds = _orbinds(sites, groups, os´) + groups = _groups(sites, os) # sites, orbranges + orbinds = _orbinds(sites, groups, os) return CellOrbitalsGrouped(cell(cs), orbinds, Dictionary(groups...)) end diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index 5e87e01a..e0ab9170 100644 --- a/src/solvers/green/bands.jl +++ b/src/solvers/green/bands.jl @@ -888,7 +888,7 @@ minimal_callsafe_copy(s::BandsGreenSlicer, parentham, parentcontacts) = s # it # triggers fast codepath above getindex_diag(gω::GreenSolution{T,<:Any,<:Any,G}, o::CellOrbitalsGrouped, symmetrize) where {T,G<:BandsGreenSlicer{<:Any,Missing}} = - maybe_symmetrize( + maybe_symmetrized_matrix( inf_band_slice(slicer(gω), SparseIndices(o, Missing), SparseIndices(o, Missing)), symmetrize) diff --git a/src/solvers/green/internal.jl b/src/solvers/green/internal.jl index 81d22685..1439f671 100644 --- a/src/solvers/green/internal.jl +++ b/src/solvers/green/internal.jl @@ -9,9 +9,11 @@ struct AppliedModelGreenSolver{F} <: AppliedGreenSolver f::F end -struct ModelGreenSlicer{C,F} <: GreenSlicer{C} +struct ModelGreenSlicer{C<:Complex,L,F} <: GreenSlicer{C} ω::C fω::F + contactorbs::ContactOrbitals{L} + g0contacts::Matrix{C} end ## API ## @@ -24,7 +26,14 @@ end ## GreenFunction API ## -(s::AppliedModelGreenSolver)(ω, args...) = ModelGreenSlicer(ω, s.f(ω)) +function (s::AppliedModelGreenSolver)(ω::C, Σblocks, contactorbs) where {C} + n = norbitals(contactorbs) + fω = s.f(ω) + g0contacts = Matrix{C}(undef, n, n) + slicer = ModelGreenSlicer(ω, fω, contactorbs, g0contacts) + fill_g0contacts!(g0contacts, slicer, contactorbs) + return slicer +end needs_omega_shift(s::AppliedModelGreenSolver) = false @@ -33,4 +42,9 @@ minimal_callsafe_copy(s::Union{ModelGreenSlicer,AppliedModelGreenSolver}, args.. Base.getindex(s::ModelGreenSlicer, is::CellOrbitals, js::CellOrbitals) = [s.fω(i, j) for i in cellorbs(is), j in cellorbs(js)] +Base.view(s::ModelGreenSlicer, i::Integer, j::Integer) = + view(s.g0contacts, contactinds(s.contactorbs, i), contactinds(s.contactorbs, j)) + +Base.view(s::ModelGreenSlicer, ::Colon, ::Colon) = view(s.g0contacts, :, :) + #endregion diff --git a/src/solvers/green/sparselu.jl b/src/solvers/green/sparselu.jl index bdd795f4..f97ad444 100644 --- a/src/solvers/green/sparselu.jl +++ b/src/solvers/green/sparselu.jl @@ -101,6 +101,8 @@ end Base.view(s::SparseLUGreenSlicer, ::Colon, ::Colon) = compute_or_retrieve_green(s, s.unitcindsall, s.unitcindsall, s.source64, s.sourceC) +# Here it is assumed that CellOrbitals has explicit orbindices (i.e. not Colon) +# This is enforced by indexing in greenfunction.jl function Base.view(s::SparseLUGreenSlicer{C}, i::CellOrbitals, j::CellOrbitals) where {C} # we only preallocate if we actually need to call ldiv! below (empty unitg cache) must_call_ldiv! = !isdefined(s, :unitg) diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index aef569d6..fdbca520 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -663,13 +663,14 @@ Broadcast.broadcast_unalias(dest::AbstractOrbitalArray, src::AbstractOrbitalArra # the product. x is a scalar. Inspired by LowRankMatrices.jl #region -struct EigenProduct{T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} <: AbstractMatrix{T} +# U and V may be SubArrays with different index type, so we need MU and MV +struct EigenProduct{T,MU<:AbstractMatrix{Complex{T}},MV<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} <: AbstractMatrix{T} blockstruct::O - U::M - V::M + U::MU + V::MV phi::T x::Complex{T} - function EigenProduct{T,M,O}(blockstruct::O, U::M, V::M, phi::T, x::Complex{T}) where {T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} + function EigenProduct{T,MU,MV,O}(blockstruct::O, U::MU, V::MV, phi::T, x::Complex{T}) where {T,MU<:AbstractMatrix{Complex{T}},MV<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} axes(U, 2) != axes(V, 2) && argerror("U and V must have identical column axis") return new(blockstruct, U, V, phi, x) end @@ -677,8 +678,8 @@ end #region ## Constructors ## -EigenProduct(blockstruct::O, U::M, V::M, phi = zero(T), x = one(Complex{T})) where {T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} = - EigenProduct{T,M,O}(blockstruct, U, V, T(phi), Complex{T}(x)) +EigenProduct(blockstruct::O, U::MU, V::MV, phi = zero(T), x = one(Complex{T})) where {T,MU<:AbstractMatrix{Complex{T}},MV<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} = + EigenProduct{T,MU,MV,O}(blockstruct, U, V, T(phi), Complex{T}(x)) #endregion @@ -726,4 +727,82 @@ Base.:*(P::EigenProduct, x::Number) = x*P LinearAlgebra.tr(P::EigenProduct) = sum(xy -> first(xy) * conj(last(xy)), zip(P.U, P.V)) * P.x #endregion + + +############################################################################################ +## SymmetrizedMatrix +# Given A::AbstractMatrix and sym::Number, build sym*A + (sym*A)' lazily +#region + +struct SymmetrizedMatrix{T,M<:AbstractMatrix{T}} <: AbstractMatrix{T} + mat::M + sym::T + function SymmetrizedMatrix{T,M}(mat::AbstractMatrix{T}, sym::T) where {T,M} + is_square(mat) || argerror("Cannot symmetrized non-square matrix") + return new(mat, sym) + end +end + +SymmetrizedMatrix(mat::M, sym::Number) where {T,M<:AbstractMatrix{T}} = + SymmetrizedMatrix{T,M}(mat, T(sym)) + +# Minimal AbstractArray interface +Base.size(a::SymmetrizedMatrix) = size(a.mat) +Base.IndexStyle(::Type{<:SymmetrizedMatrix}) = IndexCartesian() +Base.@propagate_inbounds Base.getindex(a::SymmetrizedMatrix, i::Int, j::Int) = + a.sym*getindex(a.mat, i, j) + conj(a.sym * getindex(a.mat, j, i)) + +#endregion + +############################################################################################ +## symmetrization functions - see greenfunction.jl +# If symmetrize == sym::Number, we do output .= post(sym*gij + conj(sym)*gji'). +# Otherwise just output .= post(gij) +#region + +function maybe_symmetrize!(out, gij, ::Missing; post = identity) + out .= post.(gij) + return out +end + +function maybe_symmetrize!(out, (gij, gji´), sym::Number; post = identity) + out .= post.(sym .* gij .+ conj(sym) .* gji´) + return out +end + +# see specialmatrices.jl +maybe_symmetrized_matrix(mat::AbstractArray, ::Missing) = mat +maybe_symmetrized_matrix(mat::AbstractArray, sym::Number) = SymmetrizedMatrix(mat, sym) + +maybe_symmetrized_view!(output, g, i, j, ::Missing; kw...) = + (maybe_symmetrize!(parent(output), view(g, i, j), missing; kw...); output) + +maybe_symmetrized_view!(output, g, i, j, sym, ; kw...) = + (maybe_symmetrize!(parent(output), (view(g, i, j), view(g, j, i)'), sym; kw...); output) + +maybe_symmetrized_getindex!(output, g, i, j, ::Missing; kw...) = + maybe_symmetrize!(output, g[i, j], missing; kw...) + +function maybe_symmetrized_getindex!(output, g, i, j, sym; kw...) + gij = g[i, j] + gji´ = i == j ? gij' : g[j, i]' + return maybe_symmetrize!(output, (gij, gji´), sym; kw...) +end + +maybe_symmetrized_getindex(g, i, j, ::Missing; post = identity) = post.(g[i, j]) + +function maybe_symmetrized_getindex(g, i, j, sym; kw...) + gij = g[i, j] # careful, this could be non-mutable + if i == j + gij = maybe_symmetrize!(gij, (gij, copy(gij')), sym; kw...) + else + gij = maybe_symmetrize!(gij, (gij, g[j, i]'), sym; kw...) + end + return gij +end + +# in case gij above is non-mutable +maybe_symmetrize!(::StaticArray, (gij, gji´), sym::Number; post = identity) = + post.(sym * gij + conj(sym) * gji´) + #endregion diff --git a/src/tools.jl b/src/tools.jl index 4aa9e2fd..26a043da 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -71,6 +71,10 @@ chopsmall(xs, atol) = chopsmall.(xs, Ref(atol)) chopsmall(xs) = chopsmall.(xs) chopsmall(xs::UniformScaling, atol...) = I * chopsmall(xs.λ, atol...) +mul_scalar_or_array!(x::Number, factor) = factor * x +mul_scalar_or_array!(x::Tuple, factor) = factor .* x +mul_scalar_or_array!(x::AbstractArray, factor) = (x .*= factor; x) + # Flattens matrix of Matrix{<:Number} into a matrix of Number's function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:Matrix{C}} isempty(ms) && return convert(Matrix{C}, ms) @@ -137,11 +141,11 @@ lengths_to_offsets(v::NTuple{<:Any,Integer}) = (0, cumsum(v)...) lengths_to_offsets(v) = prepend!(cumsum(v), 0) lengths_to_offsets(f::Function, v) = prepend!(accumulate((i,j) -> i + f(j), v; init = 0), 0) -# remove elements of xs after and including the first for which f(x) is false +# remove elements of xs after and including the first for which f(x) is true function cut_tail!(f, xs) i = 1 for outer i in eachindex(xs) - f(xs[i]) || break + f(xs[i]) && break end resize!(xs, i-1) return xs diff --git a/src/types.jl b/src/types.jl index 57aca8e2..46a2e56c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -596,6 +596,8 @@ Base.parent(ls::LatticeSlice) = ls.lat Base.copy(ls::LatticeSlice) = LatticeSlice(ls.lat, copy(ls.cellsdict)) +Base.collect(ci::CellIndices) = CellIndices(ci.cell, collect(ci.inds), ci.type) + #endregion #endregion @@ -2088,8 +2090,8 @@ boundingbox(oh::OpenHamiltonian) = struct ContactOrbitals{L} orbsdict::CellOrbitalsGroupedDict{L} # non-extended orbital indices for merged contact corbsdict::Vector{CellOrbitalsGroupedDict{L}} # same for each contact (alias of Σ's) - contactinds::Vector{Vector{Int}} # orbital indices in orbdict for each contact - offsets::Vector{Int} # orbsdict offset from number of orbs per cell + contactinds::Vector{Vector{Int}} # orbital indices in corbdict for each contact + offsets::Vector{Int} # corbsdict offset from number of orbs per cell end struct Contacts{L,N,S<:NTuple{N,SelfEnergy},O<:OrbitalSliceGrouped} @@ -2260,6 +2262,7 @@ greenindices(g::GreenSlice) = g.rows, g.cols # - OrbitalSliceMatrix for i,j both OrbitalSliceGrouped # - OrbitalSliceMatrix for i,j both SparseIndices of OrbitalSliceGrouped # - Matrix otherwise + similar_slice(::Type{C}, i::OrbitalSliceGrouped, j::OrbitalSliceGrouped) where {C} = OrbitalSliceMatrix{C}(undef, i, j) diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 9814d388..4403458d 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -14,9 +14,8 @@ function testgreen(h, s; kw...) z = zero(SVector{L,Int}) o = Quantica.unitvector(1, SVector{L,Int}) conts = ntuple(identity, ncontacts(h)) - locs = (sites(z, :), sites(z, 2), sites(1:2), sites(o, (1,2)), - CellOrbitals(o, 1), CellOrbitals(z, 1:2), CellOrbitals(z, SA[2,1]), - conts...) + locs = (sites(z, :), sites(z, 2), sites(1:2), sites(o, (1,2)), (; cells = z), + CellOrbitals(o, 1), CellOrbitals(z, :), CellOrbitals(z, SA[1,2]), conts...) for loc in locs, loc´ in locs gs = g[loc, loc´] @test gs isa GreenSlice @@ -333,8 +332,8 @@ end g3 = LP.triangular() |> hamiltonian(hopping(-I) + onsite(1.8I), orbitals = 2) |> supercell(10) |> supercell |> attach(nothing, region = RP.circle(2)) |> attach(nothing, region = RP.circle(2, SA[1,2])) |> greenfunction(GS.KPM(bandrange=(-4,5))) # KPM doesn't support finite temperatures or generic indexing yet, so g3 excluded from this loop - for g in (g1, g2), inds in (diagonal(1), diagonal(:), sitepairs(range = 1)) - ρ0 = densitymatrix(g[inds], 5) + for g in (g1, g2), inds in (diagonal(1), diagonal(:), sitepairs(range = 1)), path in (5, Paths.radial(1,π/6), Paths.sawtooth(5)) + ρ0 = densitymatrix(g[inds], path) ρ = densitymatrix(g[inds]) @test isapprox(ρ0(), ρ(); atol = 1e-7) @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-7) @@ -403,38 +402,38 @@ function testcond(g0; nambu = false) end function testjosephson(g0) - J1 = josephson(g0[1], 4; phases = subdiv(0, pi, 10)) - J2 = josephson(g0[2], 4; phases = subdiv(0, pi, 10)) - J3 = josephson(g0[1], (-4, 0); phases = subdiv(0, pi, 10)) - J4 = josephson(g0[1], range(-4, 0, 3); phases = subdiv(0, pi, 10)) - J5 = josephson(g0[1], (-4,-2+3im,0) .+ sqrt(eps(Float64))*im; phases = subdiv(0, pi, 10)) + J1 = josephson(g0[1], 1; phases = subdiv(0, pi, 10)) + J2 = josephson(g0[2], 1; phases = subdiv(0, pi, 10)) + J3 = josephson(g0[1], (-5, 0); phases = subdiv(0, pi, 10)) + J4 = josephson(g0[1], Paths.sawtooth(range(-5, 0, 3)); phases = subdiv(0, pi, 10)) + p5(µ, T) = ifelse(iszero(T),(-5,-2+3im,0),(-5,-2+3im,0,2+3im,5)) .+ sqrt(eps(Float64))*im + J5 = josephson(g0[1], Paths.polygon(p5); phases = subdiv(0, pi, 10)) j1 = J1() @test all(>=(0), Quantica.chopsmall.(j1)) - @test all(((j1, j2) -> ≈(j1, j2, atol = 0.000001)).(j1, J2())) - @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´ + @test all(((j1, j2) -> ≈(j1, j2, atol = 1e-6)).(j1, J2())) + @test all(((j1, j2) -> ≈(j1, j2, atol = 1e-6)).(j1, J3())) + @test all(((j1, j2) -> ≈(j1, j2, atol = 1e-6)).(j1, J4())) + @test all(((j1, j2) -> ≈(j1, j2, atol = 1e-6)).(j1, J5())) + @test all(((j1, j2) -> ≈(j1, j2, atol = 1e-6)).(J1(0.1), J5(0.1))) + j = Quantica.integrand(J1) - @test Quantica.call!(j, 0.2+0.3im) isa Vector - @test typeof.(Quantica.path.((J1, J2, J3, J4, J5))) == - (Vector{ComplexF64}, Vector{ComplexF64}, NTuple{3, ComplexF64}, Vector{ComplexF64}, NTuple{3, ComplexF64}) + @test Quantica.call!(j, 0.2) isa Vector + # static length of integration points (required for QuadGK.quadgk type-stability) + @test typeof.(Quantica.points.((J1, J2, J3, J4, J5))) == + (Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Tuple{ComplexF64, ComplexF64, ComplexF64}) # integration path logic - J = josephson(g0[1], (-4, -1)) + J = josephson(g0[1], Paths.sawtooth(-4, -1)) @test J() <= eps(Float64) - p1, p2 = Quantica.path(J, 0), Quantica.path(J, 0.2) - @test p1 isa NTuple{3,ComplexF64} # tuple sawtooth - @test p2 isa Vector{ComplexF64} && length(p2) == 5 - J = josephson(g0[1], (-4, 0); phases = subdiv(0, pi, 10)) - p1, p2 = Quantica.path(J, 0), Quantica.path(J, 0.2) - @test p1 isa NTuple{3,ComplexF64} # tuple sawtooth - @test p2 isa NTuple{3,ComplexF64} + p1, p2 = Quantica.points(J, 0), Quantica.points(J, 0.2) + @test p1 isa Vector{ComplexF64} && length(p1) == 3 # tuple sawtooth + @test p2 isa Vector{ComplexF64} && length(p2) == 3 + J = josephson(g0[1], Paths.sawtooth(-4, 0); phases = subdiv(0, pi, 10)) + p1, p2 = Quantica.points(J, 0), Quantica.points(J, 0.2) + @test p1 isa Vector{ComplexF64} && length(p1) == 3 # tuple sawtooth + @test p2 isa Vector{ComplexF64} && length(p2) == 3 @test all(p2 .≈ p1) - J = josephson(g0[1], (-4, 1); phases = subdiv(0, pi, 10)) - p1, p2 = Quantica.path(J, 0), Quantica.path(J, 0.2) + J = josephson(g0[1], Paths.sawtooth(-4, 1); phases = subdiv(0, pi, 10)) + p1, p2 = Quantica.points(J, 0), Quantica.points(J, 0.2) @test p1 isa Vector{ComplexF64} && length(p1) == 3 && maximum(real(p1)) == 0 @test p2 isa Vector{ComplexF64} && length(p2) == 5 && maximum(real(p2)) == 1 end @@ -448,8 +447,9 @@ end @test size(J1(0.2)) == size(J2(0.2)) == (3, 3) @test 2*J1(0.2; B = 0.1) ≈ J2(0.2; B = 0.1) - ρ = densitymatrix(g1[cells = SA[1]], 5) - @test Quantica.path(ρ) isa Vector{ComplexF64} && length(Quantica.path(ρ)) == 3 + ρ = densitymatrix(g1[cells = SA[1]], Paths.sawtooth(5)) + @test length(Quantica.points(ρ)) == 3 + @test Quantica.points(ρ) isa Vector{ComplexF64} @test all(≈(0.5), diag(ρ(0, 0; B=0.3))) # half filling ρ = densitymatrix(g1[cells = SA[1]], 7) @test all(<(0.96), real(diag(ρ(4, 1; B=0.1)))) # thermal depletion @@ -464,12 +464,18 @@ end @test ρ1() ≈ ρ2() atol = 0.00001 @test ρ2() ≈ ρ3() atol = 0.00001 gLU´ = h |> attach(nothing; reg...) |> greenfunction(GS.SparseLU()); - ρ1´ = densitymatrix(gLU´[1], (-3, 3)) + ρ1´ = densitymatrix(gLU´[1], Paths.sawtooth(-3, 3)) @test ρ1() ≈ ρ1´() gSpectrum´ = h |> attach(nothing; reg...) |> greenfunction(GS.Spectrum()); ρ2´ = densitymatrix(gSpectrum´[1]) @test ρ2() ≈ ρ2´() + # parameter-dependent paths + g = LP.linear() |> supercell(3) |> hamiltonian(@onsite((; Δ = 0.1) -> SA[0 Δ; Δ 0]) + hopping(SA[1 0; 0 -1]), orbitals = 2) |> greenfunction + ρ = densitymatrix(g[sites(1:2), sites(1:2)], Paths.polygon((µ,kBT; Δ) -> (-5, -Δ+im, 0, Δ+im, 5))) + @test tr(ρ(; Δ = 0.2)) ≈ 2 + @test all(Quantica.points(ρ; Δ = 0.3) .== (-5, -0.3 + 1.0im, 0, 0.3 + 1.0im, 5)) + ρ = densitymatrix(g2[(; cells = SA[0]), (; cells = SA[1])], 5) ρ0 = ρ(0, 0; B=0.3) @test ρ0 isa OrbitalSliceMatrix @@ -489,9 +495,9 @@ end g0 = LP.square() |> hamiltonian(hopping(1)) |> supercell(region = RP.square(10)) |> attach(glead, reverse = true; region = contact2) |> attach(glead; region = contact1) |> greenfunction; testcond(g0) + glead = LP.square() |> hamiltonian(hopping(I) + onsite(SA[0 1; 1 0]), orbitals = 2) |> supercell((1,0), region = r -> -1 <= r[2] <= 1) |> greenfunction(GS.Schur(boundary = 0)); contact1 = r -> r[1] ≈ 5 && -1 <= r[2] <= 1 contact2 = r -> r[1] ≈ -5 && -1 <= r[2] <= 1 - glead = LP.square() |> hamiltonian(hopping(I) + onsite(SA[0 1; 1 0]), orbitals = 2) |> supercell((1,0), region = r -> -1 <= r[2] <= 1) |> greenfunction(GS.Schur(boundary = 0)); g0 = LP.square() |> hamiltonian(hopping(I), orbitals = 2) |> supercell(region = RP.square(10)) |> attach(glead, reverse = true; region = contact2) |> attach(glead; region = contact1) |> greenfunction; testcond(g0; nambu = true) testjosephson(g0) @@ -501,7 +507,7 @@ end g0 = LP.square() |> hamiltonian(hopping(I), orbitals = 2) |> supercell(region = RP.square(10)) |> attach(glead, reverse = true; region = contact2) |> attach(glead; region = contact1) |> greenfunction; J1 = josephson(g0[1], 4; phases = subdiv(0, pi, 10), omegamap = ω ->(; ω)) J2 = josephson(g0[1], 4; phases = subdiv(0, pi, 10)) - @test Quantica.integrand(J1)(-2 + 0.00001*im) != Quantica.integrand(J2)(-2 + 0.00001*im) + @test Quantica.integrand(J1)(-2) != Quantica.integrand(J2)(-2) # test fermi at zero temperature g = LP.linear() |> hopping(1) |> supercell(3) |> supercell |> greenfunction(GS.Spectrum()) @@ -509,7 +515,7 @@ end # test DensityMatrixSchurSolver g = LP.honeycomb() |> hamiltonian(hopping(I) + @onsite((; w=0) -> SA[w 1; 1 -w], sublats = :A), orbitals = (2,1)) |> supercell((1,-1), region = r->-2 greenfunction(GS.Schur()); - ρ0 = densitymatrix(g[cells = (SA[2],SA[4]), sublats = :A], (-4, 4)) + ρ0 = densitymatrix(g[cells = (SA[2],SA[4]), sublats = :A], Paths.sawtooth(-4, 4)) ρ = densitymatrix(g[cells = (SA[2],SA[4]), sublats = :A]) ρ0sol = ρ(0.2, 0.3) ρsol = ρ(0.2, 0.3) @@ -529,6 +535,23 @@ end ρsol = ρ() @test maximum(abs, ρ0sol - ρsol) < 1e-7 @test typeof(ρ0sol) == typeof(ρsol) + + # off-diagonal rho + is = (; region = r -> 0 <= r[1] <= 4) + js = (; region = r -> 2 <= r[1] <= 6) + ρ0 = densitymatrix(g[is, js]) + ρ1 = densitymatrix(g[is, js], 4) + ρ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) + + is = (; region = r -> 0 <= r[1] <= 4) + js = sites(SA[1], 1:2) + ρ0 = densitymatrix(g[is, js]) + ρ1 = densitymatrix(g[is, js], 4) + ρ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) end @testset "greenfunction aliasing" begin @@ -581,7 +604,7 @@ end ρ12 = m.rho(0.2, 0.3)[sites(1), sites(2)] @test Φ[sites(1), sites(2)] ≈ -1.5 * Q * ρ12 * Q - # nambu + # spinless nambu oh = LP.linear() |> hamiltonian(hopping((r, dr) -> SA[1 sign(dr[1]); -sign(dr[1]) -1]) - onsite(SA[1 0; 0 -1]), orbitals = 2) g = oh |> greenfunction Q = SA[1 0; 0 -1] @@ -593,6 +616,32 @@ end hartree = 0.5*Q*(tr(Q*ρ11)*(1+1/2+1/2)) @test Φ[sites(1), sites(1)] ≈ hartree + fock + # spinful nambu - unrotated + σzτz = SA[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1] + σyτy = SA[0 0 0 -1; 0 0 1 0; 0 1 0 0; -1 0 0 0] + oh = LP.linear() |> hamiltonian(hopping(σzτz) - onsite(0.1 * σyτy), orbitals = 4) + g = oh |> greenfunction + Q = σzτz + m = meanfield(g; selector = (; range = 1), nambu = true, namburotation = false, hartree = r -> 1/(1+norm(r)), fock = 1.5, charge = Q) + Φ = m(0, 0.3) + ρ11 = m.rho(0, 0.3)[sites(1), sites(1)] - SA[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1] + fock = -1.5 * Q * ρ11 * Q + hartree = 0.5*Q*(tr(Q*ρ11)*(1+1/2+1/2)) + @test Φ[sites(1), sites(1)] ≈ hartree + fock + + # spinful nambu - rotated + σzτz = SA[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1] + σ0τx = SA[0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] + oh = LP.linear() |> hamiltonian(hopping(σzτz) - onsite(0.1 * σ0τx), orbitals = 4) + g = oh |> greenfunction + Q = σzτz + m = meanfield(g; selector = (; range = 1), nambu = true, namburotation = true, hartree = r -> 1/(1+norm(r)), fock = 1.5, charge = Q) + Φ = m(0, 0.3) + ρ11 = m.rho(0, 0.3)[sites(1), sites(1)] - SA[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1] + fock = -1.5 * Q * ρ11 * Q + hartree = 0.5*Q*(tr(Q*ρ11)*(1+1/2+1/2)) + @test Φ[sites(1), sites(1)] ≈ hartree + fock + g = LP.linear() |> hopping(1) |> greenfunction pot(r) = 1/norm(r) Φ = meanfield(g; selector = (; range = 10), potential = pot, onsite = 3)()