From 4a828d27ad8682656428ffba74c4977954d9d230 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 5 Dec 2024 17:17:37 +0100 Subject: [PATCH] taking integration paths seriously show test --- docs/src/tutorial/observables.md | 93 ++++----- src/Quantica.jl | 3 +- src/docstrings.jl | 185 ++++++++++-------- src/greenfunction.jl | 138 +++++++------ src/integrator.jl | 187 ++++++++++++++++++ src/observables.jl | 319 ++++++++++--------------------- src/show.jl | 30 +-- src/slices.jl | 16 +- src/solvers/green.jl | 1 + src/solvers/green/bands.jl | 8 +- src/solvers/green/internal.jl | 50 +++++ src/solvers/green/sparselu.jl | 2 + src/specialmatrices.jl | 91 ++++++++- src/tools.jl | 13 ++ src/types.jl | 16 +- test/test_greenfunction.jl | 125 ++++++++---- test/test_show.jl | 2 + 17 files changed, 804 insertions(+), 475 deletions(-) create mode 100644 src/integrator.jl create mode 100644 src/solvers/green/internal.jl 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/Quantica.jl b/src/Quantica.jl index 70794573..70e8ea54 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -28,7 +28,7 @@ export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian, plusadjoint, neighbors, siteselector, hopselector, diagonal, sitepairs, unflat, torus, transform, translate, combine, spectrum, energies, states, bands, subdiv, - greenfunction, selfenergy, attach, + greenfunction, selfenergy, attach, Paths, plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults, conductance, josephson, ldos, current, transmission, densitymatrix, OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes, siteindexdict, @@ -63,6 +63,7 @@ include("transform.jl") include("mesh.jl") include("bands.jl") include("greenfunction.jl") +include("integrator.jl") include("observables.jl") include("meanfield.jl") # Plumbing diff --git a/src/docstrings.jl b/src/docstrings.jl index c19475c2..7a919af7 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -2036,23 +2036,32 @@ transmission Compute a `ρ::DensityMatrix` at thermal equilibrium on sites encoded in `gs`. The actual matrix for given system parameters `params`, and for a given chemical potential `mu` and -temperature `kBT` is obtained by calling `ρ(mu = 0, kBT = 0; params...)`. The algorithm used is -specialized for the GreenSolver used, if available. In this case, `opts` are options for +temperature `kBT` is obtained by calling `ρ(mu = 0, kBT = 0; params...)`. The algorithm used +is specialized for the GreenSolver used, if available. In this case, `opts` are options for said algorithm. - densitymatrix(gs::GreenSlice, (ωmin, ωmax); opts..., quadgk_opts...) - densitymatrix(gs::GreenSlice, ωpoints; opts..., quadgk_opts...) + densitymatrix(gs::GreenSlice, path::AbstractIntegrationPath; opts..., quadgk_opts...) -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. 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`. +As above, but using a generic algorithm that relies on numerical integration of the Green +function `gs(ω)` along a contour in the complex-ω plane, defined by the `path` object. See +`Paths` for available paths. Most integration paths employ the Cauchy theorem, and therefore +assume that `gs(ω)` is fully analytic between the path and the real axis. This may not be +true if `gs` contains user-defined contacts created using `attach(model)` with general +ω-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, ωmax::Number; opts...) + densitymatrix(gs::GreenSlice, ωscale::Real; kw...) -As above with `ωmin = -ωmax`. +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 @@ -2062,23 +2071,11 @@ 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: - - ρ(μ = 0, kBT = 0, override_path = 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 The generic integration algorithm allows for the following `opts` (see also `josephson`): - `omegamap`: a function `ω -> (; params...)` that translates `ω` at each point in the integration contour to a set of system parameters. Useful for `ParametricHamiltonians` which include terms `Σ(ω)` that depend on a parameter `ω` (one would then use `omegamap = ω -> (; ω)`). Default: `ω -> (;)`, i.e. no mapped parameters. -- `imshift`: a small imaginary shift to add to the integration contour if all its vertices `ωpoints` are real numbers. Default: `missing` which is equivalent to `sqrt(eps)` for the relevant number type. -- `slope`: if `ωpoints` are all real numbers (typically encompassing the system's bandwidth), the integration contour is transformed into a triangular sawtooth path these points. Between each pair of points the path increases and then decreases linearly with the given `slope`. Default: `1.0`. -- `post`: a function to apply to the result of the integration. Default: `identity`. - `callback`: a function to be called as `callback(x, y)` at each point in the integration, where `x` is the contour point and `y` is the integrand evaluated at that point. Useful for inspection and debugging, e.g. `callback(x, y) = @show x`. Default: `Returns(nothing)`. - `atol`: absolute integration tolerance. The default `1e-7` is chosen to avoid excessive integration times when the current is actually zero. Default `1e-7`. @@ -2107,61 +2104,53 @@ julia> ρ() # with mu = kBT = 0 by default densitymatrix """ - josephson(gs::GreenSlice, ωpoints; opts..., quadgk_opts...) + josephson(gs::GreenSlice, path::AbstractIntegrationPath; opts..., quadgk_opts...) 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. 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. +into `g` through contact `i`. The algorithm relies on the integration of a function of +`gs(ω)` (see below) along a contour in the complex-ω plane defined by `path` (see `Paths` +for available options). Most integration paths employ the Cauchy theorem, and therefore +assume that `gs(ω)` is fully analytic between the path and the real axis. This may not be +true if `gs` contains user-defined contacts created using `attach(model)` with general +ω-dependent models, but note that Quantica will not check for analyticity. Keywords +`quadgk_opts` are passed to the `QuadGK.quadgk` integration routine. 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...) + josephson(gs::GreenSlice, ωscale::Real; kw...) + +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 `ωpoints = (-ωmax, ωmax)`. +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 - J(kBT = 0, override_path = missing; params...) # where J::Josephson + J(kBT = 0; 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. -It's possible to use `override_path` to override a complex integration path at evaluation -time. In this case `override_path` 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 The generic integration algorithm allows for the following `opts` (see also `densitymatrix`): -- `omegamap`: a function `ω -> (; params...)` that translates `ω` at each point in the integration contour to a set of system parameters. Useful for `ParametricHamiltonians` which include terms `Σ(ω)` that depend on a parameter `ω` (one would then use `omegamap = ω -> (; ω)`). Default: `ω -> (;)`, i.e. no mapped parameters. - `phases` : collection of superconducting phase biases to apply to the contact, so as to efficiently compute the full current-phase relation `[I_J(ϕ) for ϕ in phases]`. Note that each phase bias `ϕ` is applied by a `[cis(-ϕ/2)*I 0*I; 0*I cis(ϕ/2)*I]` rotation to the self energy, which is almost free. If `missing`, a single `I_J` is returned. -- `imshift`: a small imaginary shift to add to the integration contour if all its vertices `ωpoints` are real numbers. Default: `missing` which is equivalent to `sqrt(eps)` for the relevant number type. -- `slope`: if `ωpoints`, are all real numbers (typically encompassing the system's bandwidth), the integration contour is transformed into a triangular sawtooth path these points. Between each pair of points the path increases and then decreases linearly with the given `slope`. Default: `1.0`. -- `post`: function to apply to the result of `∫ dω f(ω) j(ω)` to obtain the result, `post = real` by default. +- `omegamap`: a function `ω -> (; params...)` that translates `ω` at each point in the integration contour to a set of system parameters. Useful for `ParametricHamiltonians` which include terms `Σ(ω)` that depend on a parameter `ω` (one would then use `omegamap = ω -> (; ω)`). Default: `ω -> (;)`, i.e. no mapped parameters. - `callback`: a function to be called as `callback(x, y)` at each point in the integration, where `x` is the contour point and `y` is the integrand at that point. Useful for inspection and debugging, e.g. `callback(x, y) = @show x`. Default: `Returns(nothing)`. - `atol`: absolute integration tolerance. The default `1e-7` is chosen to avoid excessive integration times when the current is actually zero. Default `1e-7`. The `quadgk_opts` are extra keyword arguments (other than `atol`) to pass on to the function `QuadGK.quadgk` that is used for the integration. -## Note on analyticity - -A non-zero `slope` parameter (as is the default) moves the integration path into the -upper-half complex-ω plane for increased performance. For this to work it's necessary that -the Green function and it's attached self-energies all be analytic in the upper half-plane -of complex ω. (Technically things will work also with independent analyticity in the -upper-left and upper-right quarter-planes, since the path passes 0 by default). However, no -check of analyticity is performed, so it is up to the user to ensure that. If this is not -possible, consider using `slope = 0`, or choosing a set of `ωpoints` that avoids -non-analyticities and cuts. - # Examples ``` @@ -2169,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 @@ -2192,28 +2181,74 @@ 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) = post(∫dω d(ω))`, with `post = real`. 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` and `post = -Quantica.gf_to_rho!` that computes `-(GF-GF')/(2π*im)` in place for a matrix `GF`. +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 integration path used to compute `O(args...)`. +""" +points -Return the vertices of the polygonal integration path used to compute `O(args...)`. """ -path + Paths + +A Quantica submodule that contains representations of different integration paths in the +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.radial(ωscale::Real, ϕ) + +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, 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, 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; kw...) + +As above with `ωpoints = (-ωmax, ωmax)`. + + Paths.polygon(ωpoints...) + +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; 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` +""" +Paths """ diff --git a/src/greenfunction.jl b/src/greenfunction.jl index 5e9025be..d57a75cc 100644 --- a/src/greenfunction.jl +++ b/src/greenfunction.jl @@ -64,8 +64,8 @@ end call!(g::GreenSlice{T}, ω::T; kw...) where {T} = call!(g, retarded_omega(ω, solver(parent(g))); kw...) -call!(gs::GreenSlice{T}, ω::Complex{T}; post = identity, params...) where {T} = - getindex!(gs, call!(greenfunction(gs), ω; params...); post) +call!(gs::GreenSlice{T}, ω::Complex{T}; post = identity, symmetrize = missing, params...) where {T} = + getindex!(gs, call!(greenfunction(gs), ω; params...); post, symmetrize) real_or_complex_convert(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω) real_or_complex_convert(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω) @@ -98,52 +98,23 @@ Base.view(g::GreenSolution, i::CT, j::CT = i) where {CT<:Union{Integer,Colon}} = # general getindex(::GreenSolution,...) entry point. Relies on GreenSlice constructor # args could be anything: CellSites, CellOrbs, Colon, Integer, DiagIndices... -function Base.getindex(g::GreenSolution, args...; post = identity, kw...) +function Base.getindex(g::GreenSolution, args...; post = identity, symmetrize = missing, kw...) gs = getindex(parent(g), args...; kw...) # get GreenSlice - output = getindex!(gs, g; post) + output = getindex!(gs, g; post, symmetrize) return output 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)) - -## getindex! - preallocated output for getindex - -# fastpath for intra and inter-contact, only for g::GreenSolution -getindex!(output, g::GreenSolution, i::Union{Integer,Colon}, j::Union{Integer,Colon}; post = identity) = - (parent(output) .= post.(view(g, i, j)); output) - -# indexing over single cells. g can be any type that implements g[::AnyCellOrbitals...] -getindex!(output, g, i::AnyCellOrbitals, j::AnyCellOrbitals; post = identity) = - (output .= post.(g[i, j]); output) - -# 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...) - -function getindex_cells!(output, g, cis, cjs; kw...) - for ci in cis, cj in cjs - getindex!(output, g, ci, cj; kw...) # will typically call the method below - 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 +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 @@ -158,9 +129,45 @@ orbinds_or_contactinds( c::Union{Colon,Integer,DiagIndices{Colon},DiagIndices{<:Integer}}, _, _) = (r, c) orbinds_or_contactinds(_, _, or, oc) = (or, oc) +## getindex! - core functions + +# fastpath for intra and inter-contact, only for g::GreenSolution +getindex!(output, g::GreenSolution, i::Union{Integer,Colon}, j::Union{Integer,Colon}; symmetrize = missing, kw...) = + maybe_symmetrized_view!(output, g, i, j, symmetrize; kw...) + +# indexing over single cells. g can be any type that implements g[::AnyCellOrbitals...] +getindex!(output, g, i::AnyCellOrbitals, j::AnyCellOrbitals; symmetrize = missing, kw...) = + maybe_symmetrized_getindex!(output, g, i, j, symmetrize; 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...) + +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 + 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 + +# 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...) @@ -192,45 +199,45 @@ 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; post = identity) = - 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; post = identity) +function getindex_diagonal!(output, g, o::CellOrbitalsGrouped, ker; symmetrize = missing, kw...) rngs = orbital_kernel_ranges(ker, o) - mat = getindex_diag(g, o) - fill_diagonal!(output, mat, rngs, ker; post) + mat = getindex_diag(g, o, symmetrize) + fill_diagonal!(output, mat, rngs, ker; kw...) return output end # g may be a GreenSolution or any other type with a minimal API (e.g. EigenProduct) -function getindex_diagonal!(output, g, i::OrbitalSliceGrouped, ker; post = identity) +function getindex_diagonal!(output, g, i::OrbitalSliceGrouped, ker; symmetrize = missing, kw...) offset = 0 for o in cellsdict(i) rngs = orbital_kernel_ranges(ker, o) - mat = getindex_diag(g, o) - fill_diagonal!(output, mat, rngs, ker; post, offset) + mat = getindex_diag(g, o, symmetrize) + fill_diagonal!(output, mat, rngs, ker; offset, kw...) offset += length(rngs) end return output end -function getindex_sparse!(output, g, hmask::Hamiltonian, ker; post = identity) +function getindex_sparse!(output, g, hmask::Hamiltonian, ker; symmetrize = missing, kw...) bs = blockstructure(g) # used to find blocks of g for nonzeros in hmask oj = sites_to_orbs(sites(zerocell(hmask), :), bs) # full cell for har in harmonics(hmask) oi = CellIndices(dcell(har), oj) - mat_cell = getindex_sparse(g, oi, oj, har) - fill_sparse!(har, mat_cell, bs, ker; post) + mat_cell = getindex_sparse(g, oi, oj, symmetrize) + fill_sparse!(har, mat_cell, bs, ker; kw...) end fill_sparse!(parent(output), hmask) return output end # non-optimized fallback, can be specialized by solvers -# optimized version could compute only site blocks in diagonal -getindex_diag(g, oi) = g[oi, oi] -# optimized version could compute only nonzero blocks from har::Harmonic -getindex_sparse(g, oi, oj, har) = g[oi, oj] +getindex_diag(g, oi, symmetrize) = maybe_symmetrized_getindex(g, oi, oi, symmetrize) +getindex_sparse(g, oi, oj, symmetrize) = maybe_symmetrized_getindex(g, oi, oj, symmetrize) # input index ranges for each output index to fill contact_kernel_ranges(::Missing, o::Colon, gω) = 1:norbitals(contactorbitals(gω)) @@ -512,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 new file mode 100644 index 00000000..581e628e --- /dev/null +++ b/src/integrator.jl @@ -0,0 +1,187 @@ +############################################################################################ +# Integrator paths +#region + +abstract type AbstractIntegrationPath end + +struct PolygonPath{P} <: AbstractIntegrationPath + pts::P # vertices of integration path +end + +struct SawtoothPath{T<:Real} <: AbstractIntegrationPath + realpts::Vector{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 RadialPath{T<:AbstractFloat} <: AbstractIntegrationPath + rate::T + angle::T + cisinf::Complex{T} +end + +#region ## Constructors ## + +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; 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] + +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 + +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) + +# 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::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 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) + push!(pts, 0.5 * (pts[i-1] + pts[i]) + im * slope * 0.5 * (pts[i] - pts[i-1])) + end + sort!(pts, by = real) + return pts +end + +#endregion + +#endregion +#endregion + +############################################################################################ +# Paths module +#region + +module Paths + +using Quantica: Quantica, PolygonPath, SawtoothPath, RadialPath, argerror + +## API + +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, imshift = true) = + SawtoothPath(sort!(pts), slope, imshift) + +sawtooth!(pts; kw...) = argerror("Path.sawtooth expects a collection of real points, got $pts") + +radial(rate, angle) = RadialPath(rate, angle) + +polygon(ωmax::Real) = polygon((-ωmax, ωmax)) +polygon(ω1::Number, ω2::Number, ωs::Number...) = PolygonPath((ω1, ω2, ωs...)) +polygon(ωs) = PolygonPath(ωs) + +end # Module + +#endregion + +############################################################################################ +# Integrator - integrates a function f along a complex path ωcomplex(ω::Real), connecting ωi +# The path is piecewise linear in the form of a triangular sawtooth with a given ± slope +#region + +struct Integrator{I,T,P,O<:NamedTuple,C,F} + integrand::I # call!(integrand, ω::Complex; params...)::Union{Number,Array{Number}} + pts::P # points in the integration interval + result::T # can be missing (for scalar integrand) or a mutable type (nonscalar) + quadgk_opts::O # kwargs for quadgk + callback::C # callback to call at each integration step (callback(ω, i(ω))) + post::F # function to apply to the result of the integration +end + +#region ## Constructor ## + +function Integrator(result, f, path; post = identity, callback = Returns(nothing), quadgk_opts...) + quadgk_opts´ = NamedTuple(quadgk_opts) + return Integrator(f, path, result, quadgk_opts´, callback, post) +end + +Integrator(f, path; kw...) = Integrator(missing, f, path; kw...) + +#endregion + +#region ## API ## + +integrand(I::Integrator) = I.integrand + +points(I::Integrator) = I.pts + +options(I::Integrator) = I.quadgk_opts + +## call! ## +# scalar version +function call!(I::Integrator{<:Any,Missing}; params...) + fx = x -> begin + y = call!(I.integrand, x; params...) # should be a scalar + I.callback(x, y) + return y + end + result, err = quadgk(fx, points(I)...; I.quadgk_opts...) + result´ = I.post(result) + return result´ +end + +# nonscalar version +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 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 + result´ = I.post(unsafe_deserialize(I.result, result)) + return result´ +end + +(I::Integrator)(; params...) = copy(call!(I; params...)) + +#endregion +#endregion diff --git a/src/observables.jl b/src/observables.jl index 74f68912..3969e5f0 100644 --- a/src/observables.jl +++ b/src/observables.jl @@ -51,7 +51,7 @@ check_nodiag_axes(_) = nothing #endregion ############################################################################################ -# Operators +# Operators (wrapped AbstractHamiltonians representing observables other than energy) #region function current(h::AbstractHamiltonian; charge = -I, direction = 1) @@ -63,115 +63,6 @@ end #endregion -############################################################################################ -# Integrator - integrates a function f along a complex path ωcomplex(ω::Real), connecting ωi -# The path is piecewise linear in the form of a triangular sawtooth with a given ± slope -#region - -struct Integrator{I,T,P,O<:NamedTuple,C,F} - integrand::I # call!(integrand, ω::Complex; params...)::Union{Number,Array{Number}} - points::P # a collection of points that form the triangular sawtooth integration path - result::T # can be missing (for scalar integrand) or a mutable type (nonscalar) - quadgk_opts::O # kwargs for quadgk - callback::C # callback to call at each integration step (callback(ω, i(ω))) - post::F # function to apply to the integrand at the end of integration -end - -#region ## Constructor ## - -function Integrator(result, f, pts; imshift = missing, post = identity, slope = 0, callback = Returns(nothing), quadgk_opts...) - imshift´ = imshift === missing ? - sqrt(eps(promote_type(typeof.(float.(real.(pts)))...))) : float(imshift) - sanitize_integration_points(pts) - pts´ = apply_complex_shifts(pts, imshift´, slope) - quadgk_opts´ = NamedTuple(quadgk_opts) - return Integrator(f, pts´, result, quadgk_opts´, callback, post) -end - -Integrator(f, pts; kw...) = Integrator(missing, f, pts; kw...) - -sanitize_integration_points(pts::Vector{<:Real}) = unique!(sort!(pts)) -sanitize_integration_points(pts::Vector) = unique!(pts) -sanitize_integration_points(pts::AbstractRange) = sort(pts) - -# fallback -function sanitize_integration_points(pts) - if promote_type(typeof.(pts)...) <: Real - issorted(pts) || argerror("Real integrated points should in general be sorted, got $pts") - else - allunique(pts) || argerror("Complex integrated points should in general be all unique, got $pts") - end - return pts -end - -# If all points are real, apply sawtooth with slope -apply_complex_shifts(pts::NTuple{<:Any,Real}, imshift, slope) = - iszero(slope) ? pts .+ (im * imshift) : triangular_sawtooth(imshift, slope, pts) -apply_complex_shifts(pts::AbstractVector{<:Real}, imshift, slope) = - iszero(slope) ? (pts .+ (im * imshift)) : triangular_sawtooth(imshift, slope, pts) - -# Otherwise do not modify -apply_complex_shifts(pts, imshift, slope) = pts - -triangular_sawtooth(is, sl, ωs::Tuple) = _triangular_sawtooth(is, sl, (), ωs...) -_triangular_sawtooth(is, sl, ::Tuple{}, ω1, ωs...) = _triangular_sawtooth(is, sl, (ω1 + im * is,), ωs...) -_triangular_sawtooth(is, sl, ωs´, ωn, ωs...) = _triangular_sawtooth(is, sl, - (ωs´..., 0.5 * (real(last(ωs´)) + ωn) + im * (is + sl * 0.5*(ωn - real(last(ωs´)))), ωn + im * is), ωs...) -_triangular_sawtooth(is, sl, ωs´) = ωs´ - -function triangular_sawtooth(imshift, slope, pts::Vector{T}) where {T<:Real} - pts´ = pts .+ (im*imshift) - for i in 2:length(pts) - mid = 0.5 * (pts[i-1] + pts[i]) + im * (imshift + slope * 0.5 * (pts[i] - pts[i-1])) - push!(pts´, mid) - end - sort!(pts´, by = real) - return pts´ -end - -triangular_sawtooth(imshift, slope, pts) = triangular_sawtooth(imshift, slope, [pts...]) - -#endregion - -#region ## API ## - -integrand(I::Integrator) = I.integrand - -path(I::Integrator) = I.points - -options(I::Integrator) = I.quadgk_opts - -## call! ## -# scalar version -function call!(I::Integrator{<:Any,Missing}; params...) - fx = x -> begin - y = call!(I.integrand, x; params...) # should be a scalar - I.callback(x, y) - return y - end - result, err = quadgk(fx, I.points...; I.quadgk_opts...) - result´ = I.post(result) - return result´ -end - -# nonscalar version -function call!(I::Integrator; params...) - fx! = (y, x) -> begin - y .= serialize(call!(I.integrand, x; params...)) - I.callback(x, y) - return y - end - result, err = quadgk!(fx!, serialize(I.result), I.points...; I.quadgk_opts...) - # note: post-processing is not element-wise & can be in-place - result´ = unsafe_deserialize(I.result, I.post(result)) - return result´ -end - -(I::Integrator)(; params...) = copy(call!(I; params...)) - -#endregion -#endregion - ############################################################################################ # ldos: local spectral density # d = ldos(::GreenSolution; kernel = missing) -> d[sites...]::Vector @@ -468,12 +359,15 @@ end # Keywords opts are passed to QuadGK.quadgk for the integral or the algorithm used #region -# this produces gs(ω; params...) * f(ω-mu). Use post = gf_to_rho! after integration -struct DensityMatrixIntegrand{G<:GreenSlice,T,O} +# produces integrand_transform!(gs(ω´; omegamap(ω´)..., params...) * f(ω´-mu)) +# with ω´ = path_transform(ω) +struct DensityMatrixIntegrand{G<:GreenSlice,T,O,P<:AbstractIntegrationPath,PT} gs::G mu::T kBT::T - omegamap::O + 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 # Default solver (integration in complex plane) @@ -486,107 +380,85 @@ struct DensityMatrix{S,G<:GreenSlice} gs::G end -#region ## Constructors ## - -# generic fallback (for other solvers) -(ρ::DensityMatrix)(mu = 0, kBT = 0; params...) = - ρ.solver(mu, kBT; params...) -# special case for integrator solver -(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0, override_path = missing; params...) = - ρ.solver(mu, kBT, override_path)(; params...) - -(s::DensityMatrixIntegratorSolver)(mu, kBT, override_path = missing) = - s.ifunc(mu, kBT, override_path); +#region ## densitymatrix API ## # redirects to specialized method densitymatrix(gs::GreenSlice; kw...) = densitymatrix(solver(parent(gs)), gs::GreenSlice; kw...) -# generic fallback +# generic fallback if no specialized method exists 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; opts...) = densitymatrix(gs, (-ωmax, ωmax); opts...) +densitymatrix(gs::GreenSlice, ωmax::Real; kw...) = + densitymatrix(gs::GreenSlice, Paths.radial(ωmax, π/4); kw...) -function densitymatrix(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), imshift = missing, atol = 1e-7, opts...) where {T} - # check_nodiag_axes(gs) +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)) - opts´ = (; imshift, slope = 1, post = gf_to_rho!, atol, opts...) - ωpoints_vec = collect(promote_type(T, typeof.(ωpoints)...), ωpoints) - function ifunc(mu, kBT, override_path) - ωpoints´ = override_path!(override_path, ωpoints_vec, ωpoints) - ρd = DensityMatrixIntegrand(gs, T(mu), T(kBT), omegamap) - pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, mu, kBT) - return Integrator(result, ρd, pts; opts´...) + post = post_transform_rho(path, gs) + opts´ = (; post, atol, opts...) + 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´...) end return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs) end -# If all pts are real, maybe_insert_mu! inserts mu and orders. pts can be any container. -maybe_insert_mu!(pts´, pts, mu, kBT) = - maybe_insert_mu!(pts´, pts, promote_type(typeof.(pts)...), mu, kBT) - -maybe_insert_mu!(pts´, pts, _, mu, kBT) = pts -function maybe_insert_mu!(pts´, pts, ::Type{<:Real}, mu, kBT) - # union spitting handles this type instability - if (iszero(kBT) && maximum(pts) <= mu) || any(≈(mu), pts) - return pts - else - return maybe_insert_mu!(copyto!(resize!(pts´, length(pts)), pts), mu, kBT) +# 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 end + return post end -function maybe_insert_mu!(pts::AbstractVector{<:Real}, mu, kBT) - sort!(push!(pts, mu)) - # If kBT = 0 it we filter out pts <= mu - # it's unclear whether this is useful. It allocates less, but I don't see any speedup. - iszero(kBT) && filter!(<=(mu), pts) - return pts -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) +post_transform_rho(::AbstractIntegrationPath, _) = identity -function override_path!(pts´, ptsvec::Vector{<:Complex}, pts) - resize!(ptsvec, length(pts)) - return pts´ -end +#endregion -override_path!(override_path, ptsvec, pts) = - argerror("Override of real ωpoints not supported, use complex ωpoints upon construction") +#region ## call API ## -#endregion +# generic fallback (for other solvers) +(ρ::DensityMatrix)(mu = 0, kBT = 0; params...) = + ρ.solver(mu, kBT; params...) +# special case for integrator solver +(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0; params...) = + ρ.solver(mu, kBT; params...)(; params...) -#region ## API ## +(s::DensityMatrixIntegratorSolver)(mu, kBT; params...) = + s.ifunc(mu, kBT; params...); -(gf::DensityMatrixIntegrand)(ω; params...) = copy(call!(gf, ω; params...)) +(ρi::DensityMatrixIntegrand)(x; params...) = copy(call!(ρi, x; params...)) -function call!(gf::DensityMatrixIntegrand, ω; params...) - output = call!(gf.gs, ω; gf.omegamap(ω)..., params...) - serialize(output) .*= fermi(ω - gf.mu, inv(gf.kBT)) +function call!(ρi::DensityMatrixIntegrand, x; params...) + ω = 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...) return output end -function gf_to_rho!(x::AbstractMatrix) - x .= x .- x' - x .*= -1/(2π*im) - return x -end +#endregion -# For diagonal indexing -function gf_to_rho!(x::AbstractVector) - x .= x .- conj.(x) - x .*= -1/(2π*im) - return x -end +#region ## API ## + +integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0; params...) = + integrand(ρ.solver(mu, kBT; params...)) -integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = integrand(ρ.solver(mu, kBT)) +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 @@ -598,6 +470,10 @@ call!_output(ρ::DensityMatrix) = call!_output(ρ.gs) #endregion +#endregion +#endregion + + ############################################################################################ # josephson # The equilibrium (static) Josephson current, in units of qe/h, *from* lead i is given by @@ -611,13 +487,15 @@ 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}} +struct JosephsonIntegrand{T<:AbstractFloat,P<:Union{Missing,AbstractArray},O,G<:GreenFunction{T},PA,PT} g::G 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 Σ::Matrix{Complex{T}} # preallocated workspace, full self-energy ΣggΣ::Matrix{Complex{T}} # preallocated workspace @@ -637,20 +515,15 @@ struct JosephsonIntegratorSolver{I} ifunc::I end -#region ## Constructors ## - -# generic fallback (for other solvers) -(j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...) -# special case for integrator solver -(j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0, override_path = missing; params...) = - j.solver(kBT, override_path)(; params...) +#region ## josephson API ## -(s::JosephsonIntegratorSolver)(kBT, override_path = missing) = s.ifunc(kBT, override_path) +josephson(gs::GreenSlice, ωmax::Real; kw...) = + josephson(gs, Paths.radial(ωmax, π/4); kw...) -josephson(gs::GreenSlice{T}, ωmax::Number; kw...) where {T} = - josephson(gs, (-ωmax, ωmax); kw...) +josephson(gs::GreenSlice, ωs::NTuple{<:Any,Real}; kw...) = + josephson(gs, Paths.sawtooth(ωs); kw...) -function josephson(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), phases = missing, imshift = missing, atol = 1e-7, opts...) where {T} +function josephson(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap = Returns((;)), phases = missing, atol = 1e-7, opts...) where {T} check_nodiag_axes(gs) check_same_contact_slice(gs) contact = rows(gs) @@ -660,21 +533,20 @@ function josephson(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), phases normalsize = normal_size(hamiltonian(g)) tauz = tauz_diag.(axes(Σ, 1), normalsize) 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, override_path) - ωpoints´ = override_path!(override_path, ωpoints_vec, ωpoints) - jd = JosephsonIntegrand(g, T(kBT), contact, tauz, phases´, omegamap, + opts´ = (; post = real, atol, opts...) + 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})) - pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, zero(T), kBT) - return Integrator(traces, jd, pts; opts´...) + return Integrator(traces, jd, realpts; opts´...) end return Josephson(JosephsonIntegratorSolver(ifunc), gs) end sanitize_phases_traces(::Missing, ::Type{T}) where {T} = missing, missing sanitize_phases_traces(phases::Integer, ::Type{T}) where {T} = - sanitize_phases_traces(range(0, π, length = phases), T) + sanitize_phases_traces(range(0, 2π, length = phases), T) function sanitize_phases_traces(phases, ::Type{T}) where {T} phases´ = Complex{T}.(phases) @@ -684,11 +556,37 @@ end #endregion +#region ## call API ## + +# generic fallback (for other solvers) +(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...) + +(s::JosephsonIntegratorSolver)(kBT; params...) = s.ifunc(kBT; params...) + +(J::JosephsonIntegrand)(x; params...) = copy(call!(J, x; params...)) + +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 + +#endregion + #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...)) -path(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0.0) = path(J.solver(kBT)) +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) temperature(J::JosephsonIntegrand) = J.kBT @@ -701,15 +599,6 @@ numphaseshifts(J::JosephsonIntegrand) = numphaseshifts(J.phaseshifts) numphaseshifts(::Missing) = 0 numphaseshifts(phaseshifts) = length(phaseshifts) -(J::JosephsonIntegrand)(ω; params...) = copy(call!(J, ω; params...)) - -function call!(J::JosephsonIntegrand, ω; params...) - gω = call!(J.g, ω; J.omegamap(ω)..., params...) - f = fermi(ω, inv(J.kBT)) - traces = josephson_traces(J, gω, f) - return traces -end - function josephson_traces(J, gω, f) gr = view(gω, J.contactind, J.contactind) Σi = selfenergy!(J.Σ, gω, J.contactind) 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.jl b/src/solvers/green.jl index 0c090ded..2a226a4f 100644 --- a/src/solvers/green.jl +++ b/src/solvers/green.jl @@ -76,3 +76,4 @@ include("green/spectrum.jl") include("green/schur.jl") include("green/kpm.jl") include("green/bands.jl") +include("green/internal.jl") # solvers for internal use only diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index 4ee555ce..e0ab9170 100644 --- a/src/solvers/green/bands.jl +++ b/src/solvers/green/bands.jl @@ -883,12 +883,14 @@ minimal_callsafe_copy(s::BandsGreenSlicer, parentham, parentcontacts) = s # it ############################################################################################ # getindex_diag: optimized cell indexing for BandsGreenSlicer, used by diagonal indexing +# no need to compute full ψ * ψ' if only diagonal is needed #region # triggers fast codepath above -getindex_diag(gω::GreenSolution{T,<:Any,<:Any,G}, o::CellOrbitalsGrouped) where {T,G<:BandsGreenSlicer{<:Any,Missing}} = - inf_band_slice(slicer(gω), SparseIndices(o, Missing), SparseIndices(o, Missing)) +getindex_diag(gω::GreenSolution{T,<:Any,<:Any,G}, o::CellOrbitalsGrouped, symmetrize) where {T,G<:BandsGreenSlicer{<:Any,Missing}} = + maybe_symmetrized_matrix( + inf_band_slice(slicer(gω), SparseIndices(o, Missing), SparseIndices(o, Missing)), + symmetrize) #endregion - #endregion diff --git a/src/solvers/green/internal.jl b/src/solvers/green/internal.jl new file mode 100644 index 00000000..1439f671 --- /dev/null +++ b/src/solvers/green/internal.jl @@ -0,0 +1,50 @@ +############################################################################################ +# AppliedModelGreenSolver - only for internal use +# 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. +#region + +struct AppliedModelGreenSolver{F} <: AppliedGreenSolver + f::F +end + +struct ModelGreenSlicer{C<:Complex,L,F} <: GreenSlicer{C} + ω::C + fω::F + contactorbs::ContactOrbitals{L} + g0contacts::Matrix{C} +end + +## API ## + +function (g::GreenSlice)(x::UniformScaling) + s = AppliedModelGreenSolver(Returns((i, j) -> ifelse(i == j, x.λ, zero(x.λ)))) + g´ = swap_solver(g, s) + return g´(0) +end + +## GreenFunction API ## + +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 + +minimal_callsafe_copy(s::Union{ModelGreenSlicer,AppliedModelGreenSolver}, args...) = s + +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 4adb70d4..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,6 +141,15 @@ 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 true +function cut_tail!(f, xs) + i = 1 + for outer i in eachindex(xs) + f(xs[i]) && break + end + resize!(xs, i-1) + return xs +end # fast tr(A*B) trace_prod(A, B) = (check_sizes(A,B); unsafe_trace_prod(A, B)) diff --git a/src/types.jl b/src/types.jl index 9eb11ee5..46a2e56c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -558,11 +558,14 @@ sites(l::LatticeSlice, blocktype::Type) = cellsites(cs::Union{SiteSlice,OrbitalSliceGrouped}) = cellsites(cs.cellsdict) cellsites(cdict::CellIndicesDict) = (CellSite(cell(cinds), i) for cinds in cdict for i in siteindices(cinds)) +# unused +# cellsites(cs::CellSites) = (CellSite(cell(cs), i) for i in siteindices(cs)) # single-orbital (CellOrbital) iterator cellorbs(co::Union{OrbitalSlice,OrbitalSliceGrouped}) = cellorbs(co.cellsdict) cellorbs(codict::Union{CellOrbitalsDict,CellOrbitalsGroupedDict}) = (CellOrbital(cell(corbs), i) for corbs in codict for i in orbindices(corbs)) +cellorbs(corbs::CellOrbitals) = (CellOrbital(cell(corbs), i) for i in orbindices(corbs)) # orbitals in unit cell for each site orbgroups(s::CellOrbitalsGrouped) = s.type.groups @@ -593,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 @@ -2085,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} @@ -2257,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) @@ -2326,6 +2332,12 @@ latslice(g::GreenFunction, i) = orbslice(g.contacts, i) zerocell(g::Union{GreenFunction,GreenSolution,GreenSlice}) = zerocell(lattice(g)) solver(g::GreenFunction) = g.solver +solver(g::GreenSlice) = g.parent.solver + +swap_solver(g::GreenFunction, s::AppliedGreenSolver) = + GreenFunction(g.parent, s, g.contacts) +swap_solver(g::GreenSlice, s::AppliedGreenSolver) = + GreenSlice(swap_solver(g.parent, s), g.rows, g.cols, similar(g.output)) contacts(g::GreenFunction) = g.contacts contacts(g::Union{GreenSolution,GreenSlice}) = contacts(parent(g)) 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)() diff --git a/test/test_show.jl b/test/test_show.jl index 83266f67..07dfad99 100644 --- a/test/test_show.jl +++ b/test/test_show.jl @@ -39,6 +39,8 @@ @test nothing === show(stdout, sites(SA[1], 2:4)) @test nothing === show(stdout, sites(:)) @test nothing === show(stdout, sites(SA[0], :)) + m = LP.linear() |> hopping(1) |> greenfunction |> meanfield + @test nothing === show(stdout, m) end h = first(hs) @test nothing === show(stdout, serializer(h))