From 275c5fa5fcc0a219078e5231f2af86fc691b1417 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 5 Dec 2024 15:59:56 +0100 Subject: [PATCH] global quash --- docs/make.jl | 3 + docs/src/advanced/meanfield.md | 3 +- docs/src/tutorial/observables.md | 93 ++++----- src/Quantica.jl | 3 +- src/docstrings.jl | 241 +++++++++++++---------- src/greenfunction.jl | 138 +++++++------ src/integrator.jl | 187 ++++++++++++++++++ src/meanfield.jl | 173 +++++++++++++---- src/observables.jl | 319 ++++++++++--------------------- src/serializer.jl | 17 +- src/show.jl | 52 ++--- src/slices.jl | 16 +- src/solvers/green.jl | 1 + src/solvers/green/bands.jl | 8 +- src/solvers/green/internal.jl | 50 +++++ src/solvers/green/schur.jl | 9 +- src/solvers/green/sparselu.jl | 2 + src/specialmatrices.jl | 259 ++++++++++++++++++------- src/tools.jl | 13 ++ src/types.jl | 70 +++++-- test/test_greenfunction.jl | 151 +++++++++++---- 21 files changed, 1168 insertions(+), 640 deletions(-) create mode 100644 src/integrator.jl create mode 100644 src/solvers/green/internal.jl diff --git a/docs/make.jl b/docs/make.jl index 74ff8337..8d5a52e6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,6 +12,9 @@ makedocs(; prettyurls=get(ENV, "CI", "false") == "true", canonical="https://pablosanjose.github.io/Quantica.jl", assets=["assets/custom.css"], + size_threshold_ignore = [ + "api.md" + ] ), pages=[ "Home" => "index.md", diff --git a/docs/src/advanced/meanfield.md b/docs/src/advanced/meanfield.md index 2402976d..cd3870b9 100644 --- a/docs/src/advanced/meanfield.md +++ b/docs/src/advanced/meanfield.md @@ -77,10 +77,11 @@ julia> using SIAMFANLEquations julia> h = LP.linear() |> supercell(4) |> onsite(r->r[1]) - hopping(1) + @onsite((i; phi = zerofield) --> phi[i]); julia> M = meanfield(greenfunction(h); onsite = 1, selector = (; range = 0), fock = nothing) -MeanField{ComplexF64} : builder of Hartree-Fock mean fields +MeanField{ComplexF64} : builder of Hartree-Fock-Bogoliubov mean fields Charge type : scalar (ComplexF64) Hartree pairs : 4 Mean field pairs : 4 + Nambu : false julia> Φ0 = M(0.0, 0.0); 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 008dc9d3..5e02e9c8 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -716,7 +716,7 @@ julia> h[sites(1), sites(2)] 1.0+0.0im 0.0+0.0im julia> ph = h |> @hopping!((t; p = 3) -> p*t); ph[region = RP.square(1)] -4×4 OrbitalSliceMatrix{ComplexF64,SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}: +4×4 OrbitalSliceMatrix{ComplexF64,SparseMatrixCSC}: 0.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im @@ -1684,7 +1684,7 @@ julia> gω[ss] == gs(0.1; t = 2) true julia> summary(gω[ss]) -"14×14 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}" +"14×14 OrbitalSliceMatrix{ComplexF64,Array}" ``` # See also @@ -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`. @@ -2086,7 +2083,7 @@ The `quadgk_opts` are extra keyword arguments (other than `atol`) to pass on to Currently, the following GreenSolvers implement dedicated densitymatrix algorithms: -- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries are not currently supported. No `opts`. +- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries and non-empty contacts are not currently supported. Assumes Hermitian Hamiltonian. No `opts`. - `GS.Spectrum`: based on summation occupation-weigthed eigenvectors. No `opts`. - `GS.KPM`: based on the Chebyshev expansion of the Fermi function. Currently only works for zero temperature and only supports `nothing` contacts (see `attach`). No `opts`. @@ -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 `ωpoints = (-ωmax, ωmax)`. +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 - 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 """ @@ -2302,7 +2337,7 @@ used e.g. to index another `OrbitalSliceArray` or to inspect the indices of each julia> g = HP.graphene(orbitals = 2) |> supercell((1,-1)) |> greenfunction; julia> d = ldos(g[cells = SA[0]])(2); summary(d) -"8-element OrbitalSliceVector{Float64,Vector{Float64}}" +"8-element OrbitalSliceVector{Float64,Array}" julia> a = only(orbaxes(d)) OrbitalSliceGrouped{Float64,2,1} : collection of subcells of orbitals (grouped by sites) for a 1D lattice in 2D space @@ -2464,6 +2499,10 @@ where `h = Quantica.parent_hamiltonian(as)` is the `AbstractHamiltonian` used to Return an `Array` of the same eltype as `m` that contains all the stored matrix elements of `m`. See `deserialize` for the inverse operation. + serialize(T::Type, m::OrbitalSliceArray) + +Reinterpret `serialize(m)` as a collection with eltype `T` + ## See also `serializer`, `serialize!`, `deserialize`, `deserialize!` @@ -2492,7 +2531,8 @@ serialize(s)` restored (i.e. overwritten). See `serialize` for details. Reconstruct an `OrbitalSliceArray` with the same structure as `m` but with the matrix elements enconded in `v`. This `v` is typically the result of a `serialize` call to a -another similar `m`, but the only requirement is that is has the correct size. +another similar `m`, but the only requirement is that is has the correct size. If `v` has +the wrong eltype, it will be reintepreted to match the eltype of `m`. ## See also `serializer`, `serialize`, `serialize!`, `deserialize!` @@ -2594,10 +2634,10 @@ decay_lengths """ meanfield(g::GreenFunction, args...; kw...) -Build a `M::MeanField` object that can be used to compute the Hartree-Fock mean field `Φ` -between selected sites interacting through a given charge-charge potential. The density -matrix used to build the mean field is obtained with `densitymatrix(g[pair_selection], -args...; kw...)`, see `densitymatrix` for details. +Build a `M::MeanField` object that can be used to compute the Hartree-Fock-Bogoliubov mean +field `Φ` between selected sites interacting through a given charge-charge potential. The +density matrix used to build the mean field is obtained with +`densitymatrix(g[pair_selection], args...; kw...)`, see `densitymatrix` for details. The mean field between site `i` and `j` is defined as `Φᵢⱼ = δᵢⱼ hartreeᵢ + fockᵢⱼ`, where @@ -2611,12 +2651,13 @@ where `U` is the onsite interaction. ## Keywords -- `potential`: charge-charge potential to use for both Hartree and Fock. Can be a number or a function of position. Default: `1 +- `potential`: charge-charge potential to use for both Hartree and Fock. Can be a number or a function of position. Default: `1` - `hartree`: charge-charge potential `v_H` for the Hartree mean field. Can be a number or a function of position. Overrides `potential`. Default: `potential` - `fock`: charge-charge potential `v_F` for the Fock mean field. Can be a number, a function of position or `nothing`. In the latter case all Fock terms (even onsite) will be dropped. Default: `hartree` - `onsite`: charge-charge onsite potential. Overrides both Hartree and Fock potentials for onsite interactions. Default: `hartree(0)` - `charge`: a number (in single-orbital systems) or a matrix (in multi-orbital systems) representing the charge operator on each site. Default: `I` - `nambu::Bool`: specifies whether the model is defined in Nambu space. In such case, `charge` should also be in Nambu space, typically `SA[1 0; 0 -1]` or similar. Default: `false` +- `namburotation::Bool`: if `nambu == true` and spinful systems, specifies whether the spinor basis is `[c↑, c↓, c↓⁺, -c↑⁺]` (`namburotation = true`) or `[c↑, c↓, c↑⁺, c↓⁺]` (`namburotation = false`). Default: `false` - `selector::NamedTuple`: a collection of `hopselector` directives that defines the pairs of sites (`pair_selection` above) that interact through the charge-charge potential. Default: `(; range = 0)` (i.e. onsite) Any additional keywords `kw` are passed to the `densitymatrix` function used to compute the @@ -2626,8 +2667,11 @@ mean field, see above M(µ = 0, kBT = 0; params...) # where M::MeanField -Build an `Φ::OrbitalSliceMatrix` that can be indexed at different sites or pairs of sites, -e.g. `ϕ[sites(2:3), sites(1)]`, see `OrbitalSliceArray` for details. +Build an `Φ::CompressedOrbitalMatrix`, which is a special form of `OrbitalSliceMatrix` that +can be indexed at pairs of individual sites, e.g. `ϕ[sites(2), sites(1)]` to return an +`SMatrix`. This type of matrix is less flexible than `OrbitalSliceMatrix` but is fully +static, and can encode symmetries. Its features are implementation details and are bound to +change. The returned `Φ` is just meant to be used in non-spatial models, see Examples below. # Examples @@ -2637,33 +2681,34 @@ julia> model = hopping(I) - @onsite((i; phi = zerofield) --> phi[i]); # see zer julia> g = LP.honeycomb() |> hamiltonian(model, orbitals = 2) |> supercell((1,-1)) |> greenfunction; julia> M = meanfield(g; selector = (; range = 1), charge = I, potential = 0.05) -MeanField{SMatrix{2, 2, ComplexF64, 4}} : builder of Hartree-Fock mean fields +MeanField{SMatrix{2, 2, ComplexF64, 4}} : builder of Hartree-Fock-Bogoliubov mean fields Charge type : 2 × 2 blocks (ComplexF64) Hartree pairs : 14 Mean field pairs : 28 + Nambu : false julia> phi0 = M(0.2, 0.3); julia> phi0[sites(1), sites(2)] |> Quantica.chopsmall -2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries: - 0.00109527+0.0im ⋅ - ⋅ 0.00109527+0.0im +2×2 SMatrix{2, 2, ComplexF64, 4} with indices SOneTo(2)×SOneTo(2): + 0.00109527+0.0im 0.0+0.0im + 0.0+0.0im 0.00109527+0.0im julia> phi0[sites(1)] |> Quantica.chopsmall -2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries: - 0.296672+0.0im ⋅ - ⋅ 0.296672+0.0im +2×2 SMatrix{2, 2, ComplexF64, 4} with indices SOneTo(2)×SOneTo(2): + 0.296672+0.0im 0.0+0.0im + 0.0+0.0im 0.296672+0.0im julia> phi1 = M(0.2, 0.3; phi = phi0); julia> phi1[sites(1), sites(2)] |> Quantica.chopsmall -2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries: - 0.00307712+0.0im ⋅ - ⋅ 0.00307712+0.0im +2×2 SMatrix{2, 2, ComplexF64, 4} with indices SOneTo(2)×SOneTo(2): + 0.00307712+0.0im 0.0+0.0im + 0.0+0.0im 0.00307712+0.0im ``` # See also - `zerofield`, `densitymatrix`, `OrbitalSliceMatrix` + `zerofield`, `densitymatrix` """ meanfield 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/meanfield.jl b/src/meanfield.jl index 319cb794..bdff0887 100644 --- a/src/meanfield.jl +++ b/src/meanfield.jl @@ -7,14 +7,16 @@ # we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry #region -struct MeanField{B,T,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix} +struct MeanField{B,T,C<:CompressedOrbitalMatrix,S<:SparseMatrixCSC,F<:DensityMatrix} + output::C potHartree::S potFock::S - rhoHartree::H - rhoFock::F + rho::F charge::B + nambu::Bool + is_nambu_rotated::Bool rowcol_ranges::NTuple{2,Vector{UnitRange{Int}}} - onsite_tmp::Vector{Complex{T}} + onsite_tmp::Vector{T} end struct ZeroField end @@ -23,59 +25,120 @@ struct ZeroField end function meanfield(g::GreenFunction{T,E}, args...; potential = Returns(1), hartree = potential, fock = hartree, - onsite = missing, charge = I, nambu::Bool = false, + onsite = missing, charge = I, nambu::Bool = false, namburotation = missing, selector::NamedTuple = (; range = 0), kw...) where {T,E} Vh = sanitize_potential(hartree) Vf = sanitize_potential(fock) - Q = sanitize_charge(charge, blocktype(hamiltonian(g))) U = onsite === missing ? T(Vh(zero(SVector{E,T}))) : T(onsite) - Uf = fock === nothing ? zero(U) : U + Uf = onsite === missing ? T(Vf(zero(SVector{E,T}))) : T(onsite) + Uf = fock === nothing ? zero(Uf) : Uf + B = blocktype(hamiltonian(g)) + is_nambu_rotated´ = sanitize_nambu_rotated(namburotation, nambu, B) + Q = sanitize_charge(charge, B, nambu, is_nambu_rotated´) isempty(boundaries(g)) || argerror("meanfield does not currently support systems with boundaries") isfinite(U) || argerror("Onsite potential must be finite, consider setting `onsite`") - nambu && (!is_square(charge) || iseven(size(charge, 1))) && argerror("Invalid charge matrix for Nambu space") - gsHartree = g[diagonal(; cells = 0, kernel = Q)] - rhoHartree = densitymatrix(gsHartree, args...; kw...) - gsFock = g[sitepairs(; selector..., includeonsite = true)] - rhoFock = densitymatrix(gsFock, args...; kw...) + gs = g[sitepairs(; selector..., includeonsite = true)] + rho = densitymatrix(gs, args...; kw...) lat = lattice(hamiltonian(g)) # The sparse structure of hFock will be inherited by the evaluated mean field. Need onsite. - hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : Vf(dr); selector..., includeonsite = true) + hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : T(Vf(dr)); selector..., includeonsite = true) hHartree = (Uf == U && Vh === Vf) ? hFock : - lat |> hopping((r, dr) -> iszero(dr) ? U : Vh(dr); selector..., includeonsite = true) + lat |> hopping((r, dr) -> iszero(dr) ? U : T(Vh(dr)); selector..., includeonsite = true) + # this drops zeros + potHartree = T.(sum(unflat, harmonics(hHartree))) - potHartree = sum(unflat, harmonics(hHartree)) - nambu && (nonzeros(potHartree) .*= T(1/2)) # to compensate for the factor 2 from nambu trace - - oaxes = orbaxes(call!_output(gsFock)) + oaxes = orbaxes(call!_output(gs)) rowcol_ranges = collect.(orbranges.(oaxes)) - onsite_tmp = similar(diag(parent(call!_output(gsHartree)))) + onsite_tmp = Vector{T}(undef, length(last(rowcol_ranges))) - # build potFock with identical axes as the output of rhoFock + # build potFock with identical axes as the output of rho cells_fock = cells(first(oaxes)) hFock_slice = hFock[(; cells = cells_fock), (; cells = 0)] # this is important for the fast orbrange-based implementation of MeanField evaluation - check_cell_order(hFock_slice, rhoFock) - potFock = parent(hFock_slice) + check_cell_order(hFock_slice, rho) + # this does not drop zeros. Important to keep diagonal zeros + potFock = convert(SparseMatrixCSC{T,Int}, parent(hFock_slice)) + + encoder, decoder = nambu ? NambuEncoderDecoder(is_nambu_rotated´) : (identity, identity) + S = typeof(encoder(zero(Q))) + output = call!_output(g[sitepairs(; selector..., includeonsite = true, kernel = Q)]) + sparse_enc = similar(output, S) + output = CompressedOrbitalMatrix(sparse_enc; encoder, decoder, hermitian = true) - return MeanField(potHartree, potFock, rhoHartree, rhoFock, Q, rowcol_ranges, onsite_tmp) + return MeanField(output, potHartree, potFock, rho, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp) end -sanitize_potential(x::Number) = Returns(x) +sanitize_potential(x::Real) = Returns(x) sanitize_potential(x::Function) = x sanitize_potential(x::Nothing) = Returns(0) -sanitize_potential(_) = argerror("Invalid potential: use a number or a function of position") +sanitize_potential(_) = argerror("Invalid potential: use a real number or a function of position") + +sanitize_nambu_rotated(is_nambu_rotated, nambu, ::Type{<:SMatrix{2,2}}) = false +sanitize_nambu_rotated(is_nambu_rotated, nambu, ::Type{<:SMatrix{4,4}}) = + nambu ? sanitize_nambu_rotated(is_nambu_rotated) : false +sanitize_nambu_rotated(is_nambu_rotated, nambu, B) = + nambu ? nambu_dim_error(B) : false +sanitize_nambu_rotated(::Missing) = + argerror("Must specify `namburotation` (true or false)") +sanitize_nambu_rotated(is_nambu_rotated::Bool) = is_nambu_rotated + +function sanitize_charge(charge, B, nambu, is_rotated) + Q = sanitize_charge(charge, B) + ishermitian(Q) || argerror("Charge $Q should be Hermitian") + nambu && check_nambu(Q, is_rotated) + return Q +end + +sanitize_charge(charge, B) = sanitize_block(B, charge) +sanitize_charge(charge, ::Type{<:SMatrixView}) = meanfield_multiorbital_error() + +check_nambu(Q::S, is_rotated) where {S<:Union{SMatrix{2,2},SMatrix{4,4}}} = + nambu_redundants(Q) ≈ nambu_adjoint_significants(Q, is_rotated) || + nambu_sym_error(Q) +check_nambu(Q::S, is_rotated) where {S<:SMatrix} = nambu_dim_error(S) +check_nambu(Q, is_rotated) = nambu_sym_error(Q) + +nambu_dim_error(::Type{S}) where {N,S<:SMatrix{N,N}} = + argerror("Nambu `meanfield` currently only understand 2×2 and 4×4 Nambu spaces, got $N×$N") +meanfield_multiorbital_error() = + argerror("`meanfield` does not currently support systems with heterogeneous orbitals") +nambu_sym_error(Q) = + argerror("Matrix $Q does not satisfy Nambu symmetry") + +nambu_significants(mat::SMatrix{4,4}) = mat[:, SA[1,2]] +nambu_significants(mat::SMatrix{2,2}) = mat[:, 1] +nambu_redundants(mat::SMatrix{4,4}) = mat[:, SA[3,4]] +nambu_redundants(mat::SMatrix{2,2}) = mat[:, 2] + +nambu_adjoint_significants(mat::SMatrix{N,N}, is_rotated) where {N} = + nambu_adjoint_significants(nambu_significants(mat), is_rotated) + +function nambu_adjoint_significants(lmat::SVector{2}, _) + return -SA[0 1; 1 0] * conj(lmat) +end + +function nambu_adjoint_significants(lmat::SMatrix{4,2}, is_rotated) + if is_rotated + return -SA[0 0 0 -1; 0 0 1 0; 0 1 0 0; -1 0 0 0] * lmat * SA[0 -1; 1 0] + else + return -SA[0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] * lmat + end +end -sanitize_charge(B, t) = sanitize_block(t, B) -sanitize_charge(B, ::Type{<:SMatrixView}) = argerror("meanfield does not currently support systems with heterogeneous orbitals") +function NambuEncoderDecoder(is_nambu_rotated) + encoder = nambu_significants + decoder = smat -> [smat nambu_adjoint_significants(smat, is_nambu_rotated)] + return encoder, decoder +end -function check_cell_order(hFock_slice, rhoFock) +function check_cell_order(hFock_slice, rho) opot = first(orbaxes(hFock_slice)) - orho = first(orbaxes(call!_output(rhoFock.gs))) + orho = first(orbaxes(call!_output(rho.gs))) cells(opot) == cells(orho) || internalerror("meanfield: Cell order mismatch between potential and density matrix") return nothing end @@ -90,35 +153,67 @@ hartree_matrix(m::MeanField) = m.potHartree fock_matrix(m::MeanField) = parent(m.potFock) -function (m::MeanField{B})(args...; chopsmall = true, params...) where {B} +isnambu(m::MeanField) = m.nambu + +is_nambu_rotated(m::MeanField) = m.is_nambu_rotated + +(m::MeanField)(args...; kw...) = copy(call!(m, args...; kw...)) + +function call!(m::MeanField{B}, args...; chopsmall = true, params...) where {B} Q, hartree_pot, fock_pot = m.charge, m.onsite_tmp, m.potFock rowrngs, colrngs = m.rowcol_ranges - trρQ = m.rhoHartree(args...; params...) - mul!(hartree_pot, m.potHartree, diag(parent(trρQ))) - meanfield = m.rhoFock(args...; params...) - mf_parent = parent(meanfield) + check_zero_mu(m, args...) + ρ = m.rho(args...; params...) + trρQ = maybe_nambufy_traces!(diag_real_tr_rho_Q(ρ, Q), m) + mul!(hartree_pot, m.potHartree, trρQ) + meanfield = m.output + meanfield_parent = parent(meanfield) + fill!(nonzeros(meanfield_parent), zero(eltype(meanfield_parent))) if chopsmall - hartree_pot .= Quantica.chopsmall.(hartree_pot) - nonzeros(mf_parent) .= Quantica.chopsmall.(nonzeros(mf_parent)) + hartree_pot .= Quantica.chopsmall.(hartree_pot) # this is a Vector + nzs = nonzeros(parent(ρ)) + nzs .= Quantica.chopsmall.(nzs) end rows, cols, nzs = rowvals(fock_pot), axes(fock_pot, 2), nonzeros(fock_pot) for col in cols + # 1/2 to compensate for the factor 2 from nambu trace viiQ = hartree_pot[col] * Q for ptr in nzrange(fock_pot, col) row = rows[ptr] + row < col && ishermitian(meanfield) && continue # skip upper triangle vij = nzs[ptr] - ρij = view(mf_parent, rowrngs[row], colrngs[col]) + ρij = view(ρ, rowrngs[row], colrngs[col]) vQρijQ = vij * Q * sanitize_block(B, ρij) * Q if row == col - ρij .= viiQ - vQρijQ + meanfield_parent[row, col] = encoder(meanfield)(viiQ - vQρijQ) else - ρij .= -vQρijQ + meanfield_parent[row, col] = encoder(meanfield)(-vQρijQ) end end end return meanfield end +check_zero_mu(_) = nothing +check_zero_mu(m, µ, _...) = !m.nambu || iszero(µ) || + argerror("Tried to evaluate a Nambu mean field at a nonzero µ = $µ") + +# turns tr(ρ*Q) into 0.5*tr((ρ-SA[0 0; 0 I])Q) if nambu +function maybe_nambufy_traces!(traces, m::MeanField{B}) where {B} + if isnambu(m) + shift = tr(m.charge * hole_id(B)) + traces .-= shift + traces .*= 0.5 + end + return traces +end + +diag_real_tr_rho_Q(ρ, Q) = + [real(unsafe_trace_prod(Q, view(parent(ρ), rng, rng))) for rng in siteindexdict(orbaxes(ρ, 2))] + +hole_id(::Type{<:SMatrix{2,2}}) = SA[0 0; 0 1] +hole_id(::Type{<:SMatrix{4,4}}) = SA[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1] +hole_id(S) = nambu_dim_error(S) ## ZeroField diff --git a/src/observables.jl b/src/observables.jl index c5b472a4..5b5924b4 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´ = 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/serializer.jl b/src/serializer.jl index d9df24ad..3e0883fe 100644 --- a/src/serializer.jl +++ b/src/serializer.jl @@ -128,22 +128,27 @@ end ############################################################################################ -# serialize and deserialize for OrbitalSliceArray +# serialize and deserialize for AbstractOrbitalArray # extract underlying data and reconstruct an object using such data #region -serialize(a::OrbitalSliceArray) = serialize_array(parent(a)) +serialize(a::AbstractOrbitalArray) = serialize_array(parent(a)) serialize(a::AbstractArray) = serialize_array(a) +serialize(::Type{T}, a::AbstractArray) where {T} = reinterpret(T, serialize(a)) serialize_array(a::Array) = a serialize_array(a::Diagonal{<:Any,<:Vector}) = a.diag serialize_array(a::SparseMatrixCSC) = nonzeros(a) -deserialize(a, v) = (check_serializer_size(a, v); _deserialize(a, v)) +deserialize(a::AbstractArray{T}, v::AbstractArray) where {T} = + deserialize(a, reinterpret(T, v)) +deserialize(a::AbstractArray{T}, v::AbstractArray{T}) where {T} = + (check_serializer_size(a, v); unsafe_deserialize(a, v)) -_deserialize(a::OrbitalSliceArray, v::AbstractArray) = - OrbitalSliceArray(deserialize_array(parent(a), v), orbaxes(a)) -_deserialize(a::AbstractArray, v::AbstractArray) = deserialize_array(a, v) +unsafe_deserialize(a::AbstractOrbitalArray, v::AbstractArray) = + similar(a, deserialize_array(parent(a), v), orbaxes(a)) +unsafe_deserialize(a::AbstractArray, v::AbstractArray) = + deserialize_array(a, v) deserialize_array(::AbstractArray{<:Any,N}, v::AbstractArray{<:Any,N}) where {N} = v deserialize_array(::Diagonal, v::AbstractVector) = diff --git a/src/show.jl b/src/show.jl index 4deb696b..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 @@ -540,20 +540,19 @@ Base.summary(::CurrentDensitySlice{T}) where {T} = #endregion ############################################################################################ -# OrbitalSliceMatrix +# AbstractOrbitalArray #region # For simplified printing of the array typename -function Base.showarg(io::IO, ::OrbitalSliceMatrix{T,M}, toplevel) where {T,M} +function Base.showarg(io::IO, ::A, toplevel) where {T,M,A<:AbstractOrbitalArray{T,<:Any,M}} toplevel || print(io, "::") - print(io, "OrbitalSliceMatrix{$T,$M}") + print(io, "$(nameof(A)){$T,$(nameof(M))}") end -function Base.showarg(io::IO, ::OrbitalSliceVector{T,M}, toplevel) where {T,M} - toplevel || print(io, "::") - print(io, "OrbitalSliceVector{$T,$M}") -end +Base.nameof(::Type{<:OrbitalSliceMatrix}) = "OrbitalSliceMatrix" +Base.nameof(::Type{<:OrbitalSliceVector}) = "OrbitalSliceVector" +Base.nameof(::Type{<:CompressedOrbitalMatrix}) = "CompressedOrbitalMatrix" #endregion @@ -619,7 +618,7 @@ display_parameters(s::AppliedSerializer{<:Any,<:ParametricHamiltonian}) = string ############################################################################################ -# MeanField and EvaluatedMeanField +# MeanField #region function Base.show(io::IO, s::MeanField{Q}) where {Q} @@ -628,10 +627,13 @@ function Base.show(io::IO, s::MeanField{Q}) where {Q} print(io, i, summary(s), "\n", "$i Charge type : $(displaytype(Q)) $i Hartree pairs : $(nnz(hartree_matrix(s))) -$i Mean field pairs : $(nnz(fock_matrix(s)))") +$i Mean field pairs : $(nnz(fock_matrix(s))) +$i Nambu : $(nambustring(s))") end Base.summary(::MeanField{Q}) where {Q} = - "MeanField{$Q} : builder of Hartree-Fock mean fields" + "MeanField{$Q} : builder of Hartree-Fock-Bogoliubov mean fields" + +nambustring(s) = isnambu(s) ? "true $(is_nambu_rotated(s) ? "(rotated basis)" : "")" : "false" #endregion 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/schur.jl b/src/solvers/green/schur.jl index dfe4a007..de3174c1 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -83,8 +83,8 @@ function nearest_cell_harmonics(h) is_nearest || argerror("Too many or too few harmonics. Perhaps try `supercell` to ensure strictly nearest-cell harmonics.") hm, h0, hp = h[hybrid(-1)], h[hybrid(0)], h[hybrid(1)] - flat(hm) == flat(hp)' || - argerror("The Hamiltonian should have h[1] == h[-1]' to use the Schur solver") + flat(hm) ≈ flat(hp)' || + argerror("The Hamiltonian should have h[1] ≈ h[-1]' to use the Schur solver") return hm, h0, hp end @@ -689,6 +689,7 @@ end # use computed Fermi points to integrate spectral function by segments # returns an AbstractMatrix +# we don't use Integrator, because that is meant for integrals over energy, not momentum function (s::DensityMatrixSchurSolver)(µ, kBT; params...) g = parent(s.gs) λs = propagating_eigvals(g, µ + 0im, 1e-2; params...) @@ -696,7 +697,7 @@ function (s::DensityMatrixSchurSolver)(µ, kBT; params...) xs = sort!(ϕs ./ (2π)) pushfirst!(xs, -0.5) push!(xs, 0.5) - xs = unique!(xs) + xs = [mean(view(xs, rng)) for rng in approxruns(xs)] # elliminate approximate duplicates result = call!_output(s.gs) integrand!(x) = fermi_h!(result, s, 2π * x, µ, inv(kBT); params...) fx! = (y, x) -> (y .= integrand!(x)) @@ -709,7 +710,7 @@ function fermi_h!(result, s, ϕ, µ, β = 0; params...) bs = blockstructure(h) # Similar to spectrum(h, ϕ; params...), but less work (no sort! or sanitization) copy!(s.hmat, call!(h, ϕ; params...)) # sparse to dense - ϵs, psis = eigen!(s.hmat) + ϵs, psis = eigen!(Hermitian(s.hmat)) # special-casing β = Inf with views turns out to be slower fs = (@. ϵs = fermi(ϵs - µ, β)) fpsis = (s.psis .= psis .* transpose(fs)) 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 318467b8..fdbca520 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -492,69 +492,72 @@ checkcontactindices(allcontactinds, hdim) = maximum(allcontactinds) <= hdim || #endregion ############################################################################################ -## OrbitalSliceArray +## AbstractOrbitalArray #region # AbstractArray interface -Base.size(a::OrbitalSliceArray) = size(parent(a)) -Base.iterate(a::OrbitalSliceArray, i...) = iterate(parent(a), i...) -Base.length(a::OrbitalSliceArray) = length(parent(a)) -Base.IndexStyle(::Type{T}) where {M,T<:OrbitalSliceArray{<:Any,<:Any,M}} = IndexStyle(M) -Base.similar(a::OrbitalSliceArray) = OrbitalSliceArray(similar(parent(a)), orbaxes(a)) -Base.similar(a::OrbitalSliceArray, t::Type) = OrbitalSliceArray(similar(parent(a), t), orbaxes(a)) -# doesn't make sense to keep orbaxes in similar with different dimensions. -Base.similar(a::OrbitalSliceArray, dims::Tuple) = similar(parent(a), dims) -Base.copy(a::OrbitalSliceArray) = OrbitalSliceArray(copy(parent(a)), orbaxes(a)) -Base.@propagate_inbounds Base.getindex(a::OrbitalSliceArray, i::Int) = +Base.size(a::AbstractOrbitalArray) = size(parent(a)) +Base.iterate(a::AbstractOrbitalArray, i...) = iterate(parent(a), i...) +Base.length(a::AbstractOrbitalArray) = length(parent(a)) +Base.IndexStyle(::Type{T}) where {M,T<:AbstractOrbitalArray{<:Any,<:Any,M}} = IndexStyle(M) +Base.similar(a::AbstractOrbitalArray) = similar(a, similar(parent(a)), orbaxes(a)) +Base.similar(a::AbstractOrbitalArray, t::Type) = similar(a, similar(parent(a), t), orbaxes(a)) +# doesn't make sense to keep orbaxes in similar with different dimensions, return parent +Base.similar(a::AbstractOrbitalArray, dims::Tuple) = similar(parent(a), dims) +Base.copy(a::AbstractOrbitalArray) = similar(a, copy(parent(a)), orbaxes(a)) +Base.@propagate_inbounds Base.getindex(a::AbstractOrbitalArray, i::Int) = getindex(parent(a), i) -Base.@propagate_inbounds Base.getindex(a::OrbitalSliceArray, I::Vararg{Int, N}) where {N} = +Base.@propagate_inbounds Base.getindex(a::AbstractOrbitalArray, I::Vararg{Int, N}) where {N} = getindex(parent(a), I...) -Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, i::Int) = setindex!(parent(a), v, i) -Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, I::Vararg{Int, N}) where {N} = setindex!(parent(a), v, I...) +Base.@propagate_inbounds Base.setindex!(a::AbstractOrbitalArray, v, i::Int) = setindex!(parent(a), v, i) +Base.@propagate_inbounds Base.setindex!(a::AbstractOrbitalArray, v, I::Vararg{Int, N}) where {N} = setindex!(parent(a), v, I...) -Base.:*(x::Number, a::OrbitalSliceArray) = OrbitalSliceArray(parent(a) * x, orbaxes(a)) -Base.:*(a::OrbitalSliceArray, x::Number) = x * a +Base.:*(x::Number, a::AbstractOrbitalArray) = similar(a, parent(a) * x, orbaxes(a)) +Base.:*(a::AbstractOrbitalArray, x::Number) = x * a -Base.:+(a::OrbitalSliceArray, b::OrbitalSliceArray) = parent(a) + parent(b) +Base.:+(a::AbstractOrbitalArray, b::AbstractOrbitalArray) = parent(a) + parent(b) -Base.:-(a::OrbitalSliceArray, b::OrbitalSliceArray) = parent(a) - parent(b) +Base.:-(a::AbstractOrbitalArray, b::AbstractOrbitalArray) = parent(a) - parent(b) # Additional indexing over sites -Base.getindex(a::OrbitalSliceMatrix; sites...) = getindex(a, siteselector(; sites...)) -Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::NamedTuple = i) = +Base.getindex(a::AbstractOrbitalMatrix; sites...) = getindex(a, siteselector(; sites...)) +Base.getindex(a::AbstractOrbitalMatrix, i::NamedTuple, j::NamedTuple = i) = getindex(a, siteselector(i), siteselector(j)) -Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::SiteSelector) = +Base.getindex(a::AbstractOrbitalMatrix, i::NamedTuple, j::SiteSelector) = getindex(a, siteselector(i), j) -Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::NamedTuple) = +Base.getindex(a::AbstractOrbitalMatrix, i::SiteSelector, j::NamedTuple) = getindex(a, i, siteselector(j)) -# SiteSelector: return a new OrbitalSliceMatrix -function Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::SiteSelector = i) +# SiteSelector: return a new AbstractOrbitalMatrix +function Base.getindex(a::AbstractOrbitalMatrix, i::SiteSelector, j::SiteSelector = i) rowslice, colslice = orbaxes(a) rowslice´, colslice´ = rowslice[i], colslice[j] rows = collect(indexcollection(rowslice, rowslice´)) cols = i === j && rowslice === colslice ? rows : indexcollection(colslice, colslice´) m = parent(a)[rows, cols] - return OrbitalSliceMatrix(m, (rowslice´, colslice´)) + return similar(a, m, (rowslice´, colslice´)) end # CellSites: return an unwrapped Matrix -Base.getindex(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) = +Base.getindex(a::AbstractOrbitalMatrix, i::AnyCellSites, j::AnyCellSites = i) = copy(view(a, i, j)) # CellSitePos: return an unwrapped SMatrix or a scalar, even if out of bounds -Base.getindex(a::OrbitalSliceMatrix, i::C, j::C = i) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} = +Base.getindex(a::AbstractOrbitalMatrix, i::C, j::C = i) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} = checkbounds(Bool, a, i, j) ? sanitize_block(B, view(a, i, j)) : zero(B) -Base.checkbounds(::Type{Bool}, a::OrbitalSliceMatrix, i::CellSitePos, j::CellSitePos) = - checkbounds(Bool, first(orbaxes(a)), i) && checkbounds(Bool, last(orbaxes(a)), j) +Base.checkbounds(::Type{Bool}, a::AbstractOrbitalMatrix, i, j) = + checkbounds_axis_or_orbaxis(a, i, 1) && checkbounds_axis_or_orbaxis(a, j, 2) -function Base.view(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) +Base.view(a::AbstractOrbitalMatrix, i::AnyCellSites, j::AnyCellSites = i) = + view(parent(a), indexcollection(a, i, j)...) + +function indexcollection(a::AbstractOrbitalMatrix, i::AnyCellSites, j::AnyCellSites = i) rowslice, colslice = orbaxes(a) - i´, j´ = apply(i, lattice(rowslice)), apply(j, lattice(colslice)) - rows = indexcollection(rowslice, i´) - cols = j === i && rowslice === colslice ? rows : indexcollection(colslice, j´) - return view(parent(a), rows, cols) + isdiag = j === i && rowslice === colslice + rows = indexcollection(rowslice, apply(i, lattice(rowslice))) + cols = isdiag ? rows : indexcollection(colslice, apply(j, lattice(colslice))) + return rows, cols end # Minimal getindex for OrbitalSliceVector @@ -562,65 +565,96 @@ Base.getindex(a::OrbitalSliceVector, i::AnyCellSites) = copy(view(a, i)) Base.getindex(a::OrbitalSliceVector, i::C) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} = sanitize_block(B, view(a, i)) -function Base.view(a::OrbitalSliceVector, i::AnyCellSites) +Base.view(a::OrbitalSliceVector, i::AnyCellSites) = view(parent(a), indexcollection(a, i)) + +function indexcollection(a::OrbitalSliceVector, i::AnyCellSites) rowslice = only(orbaxes(a)) - i´ = apply(i, lattice(rowslice)) - rows = indexcollection(rowslice, i´) - return view(parent(a), rows) + rows = indexcollection(rowslice, apply(i, lattice(rowslice))) + return rows end LinearAlgebra.diag(a::OrbitalSliceMatrix) = OrbitalSliceVector(diag(parent(a)), (first(orbaxes(a)),)) -LinearAlgebra.norm(a::OrbitalSliceMatrix) = norm(parent(a)) +LinearAlgebra.norm(a::AbstractOrbitalMatrix) = norm(parent(a)) + +## CompressedOrbitalMatrix + +# decoding access +Base.view(a::CompressedOrbitalMatrix{C}, i::AnyCellSites, j::AnyCellSites = i) where {C} = + checkbounds(Bool, a, i, j) ? getindex_decode(a, i, j) : zero(decoder(a)(zero(C))) + +@inline function getindex_decode(a::CompressedOrbitalMatrix, ci, cj) + i, j = only.(indexcollection(a, ci, cj)) + # Only lower-triangle is stored in the hermitian case + if ishermitian(a) + # we find the indices for a[sites(-ni, j), sites(0, i)], + # the conjugate of a[sites(ni, i), sites(0, j)], + ci´, cj´ = sites(-cell(ci), siteindex(cj)), sites(cell(cj), siteindex(ci)) + i´, j´ = only.(indexcollection(a, ci´, cj´)) + val = getindex_decode_hermitian(a, i, j) + val´ = getindex_decode_hermitian(a, i´, j´)' + return 0.5 * (val + val´) + end + return dec(parent(a)[i, j]) +end + +function getindex_decode_hermitian(a::CompressedOrbitalMatrix, i, j) + dec = decoder(a) + if i < j + return dec(parent(a)[j, i])' + elseif i == j + return dec(parent(a)[i, j]) + else + return dec(parent(a)[i, j]) + end +end + +## checkbounds + +# checkbounds axis +checkbounds_axis_or_orbaxis(a::AbstractOrbitalMatrix, i, axis) = + checkbounds(Bool, axes(a, axis), i) + +# checkbounds orbaxis +checkbounds_axis_or_orbaxis(a::AbstractOrbitalMatrix, i::CellIndices, axis) = + checkbounds_orbaxis(orbaxes(a, axis), i) + +function checkbounds_orbaxis(a::LatticeSlice, i::CellIndices) + dict = cellsdict(a) + c = cell(i, a) + if haskey(dict, c) + return siteindex(i) in keys(orbgroups(dict[c])) + else + return false + end +end ## broadcasting # following the manual: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array # disabled for now, since it seems to cause issues with OrbitalSliceMatrix{<:SparseMatrixCSC} -# Broadcast.BroadcastStyle(::Type{O}) where {O<:OrbitalSliceArray} = Broadcast.ArrayStyle{O}() +# Broadcast.BroadcastStyle(::Type{O}) where {O<:AbstractOrbitalArray} = Broadcast.ArrayStyle{O}() # Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{O}}, ::Type{ElType}) where {ElType,O<:OrbitalSliceMatrix{<:Any,<:Array}} = -# OrbitalSliceArray(similar(Array{ElType}, axes(bc)), orbaxes(find_osa(bc))) +# AbstractOrbitalArray(similar(Array{ElType}, axes(bc)), orbaxes(find_osa(bc))) # Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{O}}, ::Type{ElType}) where {ElType,O<:OrbitalSliceMatrix{<:Any,<:SparseMatrixCSC}} = -# OrbitalSliceArray(similar(parent(find_osa(bc)), ElType), orbaxes(find_osa(bc))) +# AbstractOrbitalArray(similar(parent(find_osa(bc)), ElType), orbaxes(find_osa(bc))) # find_osa(bc::Base.Broadcast.Broadcasted) = find_osa(bc.args) # find_osa(args::Tuple) = find_osa(find_osa(args[1]), Base.tail(args)) # find_osa(x) = x # find_osa(::Tuple{}) = nothing -# find_osa(a::OrbitalSliceArray, rest) = a +# find_osa(a::AbstractOrbitalArray, rest) = a # find_osa(::Any, rest) = find_osa(rest) # taken from https://github.com/JuliaArrays/OffsetArrays.jl/blob/756e839563c88faa4ebe4ff971286747863aaff0/src/OffsetArrays.jl#L469 -Base.dataids(A::OrbitalSliceArray) = Base.dataids(parent(A)) -Broadcast.broadcast_unalias(dest::OrbitalSliceArray, src::OrbitalSliceArray) = +Base.dataids(A::AbstractOrbitalArray) = Base.dataids(parent(A)) +Broadcast.broadcast_unalias(dest::AbstractOrbitalArray, src::AbstractOrbitalArray) = parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src) -# ## conversion: If i contains an OrbitalSliceGrouped, wrap it in an OrbitalSliceArray. Otherwise, no-op - -# maybe_OrbitalSliceArray(i) = x -> maybe_OrbitalSliceArray(x, i) - -# maybe_OrbitalSliceArray(x::AbstractMatrix, i) = maybe_OrbitalSliceMatrix(x, i) - -# maybe_OrbitalSliceMatrix(x::Diagonal, i::OrbitalSliceGrouped) = -# maybe_OrbitalSliceMatrix(x, (i, i)) -# maybe_OrbitalSliceMatrix(x, (i, j)::Tuple{OrbitalSliceGrouped,OrbitalSliceGrouped}) = -# OrbitalSliceMatrix(x, (i, j)) -# maybe_OrbitalSliceMatrix(x, i) = x - -## currently unsused -# maybe_OrbitalSliceArray(x::AbstractVector, i) = maybe_OrbitalSliceVector(x, i) -# maybe_OrbitalSliceVector(x, i::DiagIndices{Missing,<:OrbitalSliceGrouped}) = -# maybe_OrbitalSliceVector(x, parent(i)) -# maybe_OrbitalSliceVector(x, i::DiagIndices{<:Any,<:OrbitalSliceGrouped}) = -# maybe_OrbitalSliceVector(x, scalarize(parent(i))) -# maybe_OrbitalSliceVector(x, i::OrbitalSliceGrouped) = OrbitalSliceVector(x, (i,)) -# maybe_OrbitalSliceVector(x, i) = x - - #endregion ############################################################################################ @@ -629,13 +663,14 @@ Broadcast.broadcast_unalias(dest::OrbitalSliceArray, src::OrbitalSliceArray) = # 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 @@ -643,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 @@ -692,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 025b1d80..46a2e56c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -164,6 +164,7 @@ numbertype(::Lattice{T}) where {T} = T zerocell(::Bravais{<:Any,<:Any,L}) where {L} = zero(SVector{L,Int}) zerocell(::Lattice{<:Any,<:Any,L}) where {L} = zero(SVector{L,Int}) +zerocell(::Val{L}) where {L} = zero(SVector{L,Int}) zerocellsites(l::Lattice, i) = sites(zerocell(l), i) @@ -519,6 +520,8 @@ ncells(x::CellIndicesDict) = length(x) ncells(x::CellIndices) = 1 cell(s::CellIndices) = s.cell +cell(s::CellIndices{L}, ::LatticeSlice{<:Any,<:Any,L}) where {L} = s.cell +cell(::CellIndices{0}, ::LatticeSlice{<:Any,<:Any,L}) where {L} = zerocell(Val(L)) cells(l::LatticeSlice) = keys(l.cellsdict) @@ -543,15 +546,6 @@ boundingbox(l::LatticeSlice) = boundingbox(keys(cellsdict(l))) pos(s::CellSitePos) = s.type.r ind(s::CellSitePos) = s.inds -function Base.checkbounds(::Type{Bool}, a::OrbitalSliceGrouped, i::CellSitePos) - dict = cellsdict(a) - if haskey(dict, cell(i)) - return siteindex(i) in keys(orbgroups(dict[cell(i)])) - else - return false - end -end - # iterators sites(l::LatticeSlice) = @@ -564,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 @@ -599,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 @@ -1697,7 +1696,7 @@ LinearAlgebra.ishermitian(h::ParametricHamiltonian) = argerror("`ishermitian(::ParametricHamiltonian)` not supported, as the result can depend on the values of parameters.") copy_lattice(p::ParametricHamiltonian) = ParametricHamiltonian( - copy_lattice(p.hparent), p.h, p.modifiers, p.allptrs, p.allparams) + p.hparent, copy_lattice(p.h), p.modifiers, p.allptrs, p.allparams) copy_harmonics_shallow(p::ParametricHamiltonian) = ParametricHamiltonian( copy_harmonics_shallow(p.hparent), copy_harmonics_shallow(p.h), p.modifiers, p.allptrs, p.allparams) @@ -2091,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} @@ -2247,8 +2246,8 @@ GreenSlice(parent::GreenFunction{T}, rows::GreenIndices, cols::GreenIndices, ::T GreenSlice(parent, rows, cols, similar_slice(C, orbindices(rows), orbindices(cols))) function GreenSlice(parent, rows, cols, type...) - if rows isa DiagIndices || cols isa DiagIndices - rows === cols || argerror("Diagonal indices should be identical for rows and columns") + if rows isa SparseIndices || cols isa SparseIndices + rows === cols || argerror("Sparse indices should be identical for rows and columns") end rows´ = greenindices(rows, parent) cols´ = cols === rows ? rows´ : greenindices(cols, parent) @@ -2263,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) @@ -2332,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)) @@ -2511,7 +2517,12 @@ SparseArrays.nnz(h::BarebonesOperator) = sum(har -> nnz(matrix(har)), harmonics( # if missing, the dimension does not span orbitals but something else #region -struct OrbitalSliceArray{C,N,M<:AbstractArray{C,N},A<:NTuple{N,Union{OrbitalSliceGrouped,Missing}}} <: AbstractArray{C,N} +abstract type AbstractOrbitalArray{C,N,M} <: AbstractArray{C,N} end + +const AbstractOrbitalVector{C,M} = AbstractOrbitalArray{C,1,M} +const AbstractOrbitalMatrix{C,M} = AbstractOrbitalArray{C,2,M} + +struct OrbitalSliceArray{C,N,M<:AbstractArray{C,N},A<:NTuple{N,Union{OrbitalSliceGrouped,Missing}}} <: AbstractOrbitalArray{C,N,M} parent::M orbaxes::A end @@ -2519,6 +2530,15 @@ end const OrbitalSliceVector{C,V,A} = OrbitalSliceArray{C,1,V,A} const OrbitalSliceMatrix{C,M,A} = OrbitalSliceArray{C,2,M,A} +struct CompressedOrbitalMatrix{C<:Union{Number,SArray},M,O<:OrbitalSliceMatrix{C,M},E,D} <: AbstractOrbitalMatrix{C,M} + omat::O + hermitian::Bool + encoder::E + decoder::D +end + +#region ## Contructors ## + OrbitalSliceVector(v::AbstractVector, axes) = OrbitalSliceArray(v, axes) OrbitalSliceMatrix(m::AbstractMatrix, axes) = OrbitalSliceArray(m, axes) @@ -2537,8 +2557,28 @@ function OrbitalSliceMatrix{C}(::UndefInitializer, oi::SparseIndices{<:Hamiltoni return OrbitalSliceArray(mat, (i, j)) end +CompressedOrbitalMatrix(o::OrbitalSliceMatrix; hermitian = false, encoder = identity, decoder = identity) = + CompressedOrbitalMatrix(o, hermitian, encoder, decoder) + +Base.similar(::OrbitalSliceArray, mat, orbaxes) = OrbitalSliceArray(mat, orbaxes) +Base.similar(a::CompressedOrbitalMatrix, mat, orbaxes) = + CompressedOrbitalMatrix(similar(a.omat, mat, orbaxes), a.hermitian, a.encoder, a.decoder) + +#endregion +#region ## API ## + orbaxes(a::OrbitalSliceArray) = a.orbaxes +orbaxes(a::OrbitalSliceArray, i::Integer) = a.orbaxes[i] +orbaxes(a::CompressedOrbitalMatrix, args...) = orbaxes(a.omat, args...) Base.parent(a::OrbitalSliceArray) = a.parent +Base.parent(a::CompressedOrbitalMatrix) = parent(a.omat) + +encoder(a::CompressedOrbitalMatrix) = a.encoder + +decoder(a::CompressedOrbitalMatrix) = a.decoder +LinearAlgebra.ishermitian(a::CompressedOrbitalMatrix) = a.hermitian + +#endregion #endregion diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 5559d2ea..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,9 +535,26 @@ 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 "aliasing" begin +@testset "greenfunction aliasing" begin # Issue #267 g = LP.linear() |> hamiltonian(@hopping((; q = 1) -> q*I), orbitals = 2) |> greenfunction g´ = Quantica.minimal_callsafe_copy(g) @@ -561,15 +584,63 @@ end @test g.contacts.selfenergies[2].solver.hlead[(0,)] === g.contacts.selfenergies[1].solver.hlead[(0,)] @test g.contacts.selfenergies[2].solver.hlead[(1,)] === g.contacts.selfenergies[1].solver.hlead[(-1,)] @test g.contacts.selfenergies[2].solver.hlead[(-1,)] === g.contacts.selfenergies[1].solver.hlead[(1,)] + + # ensure full dealiasing of lattices in attach + model = hopping(SA[1 0; 0 -1]) + @onsite((; µ = 0) -> SA[-µ 0; 0 µ]) + h = LP.linear() |> hamiltonian(model, orbitals = 2) + glead = h |> greenfunction(GS.Schur(boundary = 0)) + g = h |> attach(glead, cells = 1) |> greenfunction(GS.Schur(boundary = 0)); + @test sites(lattice(h)) == [SA[0.0]] end @testset "meanfield" begin - oh = LP.honeycomb() |> hamiltonian(hopping(SA[1 im; -im 1]) - onsite(SA[1 0; 0 -1], sublats = :A), orbitals = 2) |> supercell((1,-1)) + oh = LP.honeycomb() |> hamiltonian(hopping(SA[1 im; -im -1]) - onsite(SA[1 0; 0 -1], sublats = :A), orbitals = 2) |> supercell((1,-1)) g = oh |> greenfunction + Q = SA[0 1; im 0] # nonhermitian charge not allowed + @test_throws ArgumentError meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q) Q = SA[0 -im; im 0] m = meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q) Φ = m(0.2, 0.3) - @test Φ[sites(1), sites(2)] ≈ -1.5 * Q * m.rhoFock(0.2, 0.3)[sites(1), sites(2)] * Q + ρ12 = m.rho(0.2, 0.3)[sites(1), sites(2)] + @test Φ[sites(1), sites(2)] ≈ -1.5 * Q * ρ12 * Q + + # 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] + m = meanfield(g; selector = (; range = 1), nambu = true, hartree = r -> 1/(1+norm(r)), fock = 1.5, charge = Q) + @test_throws ArgumentError m(0.2, 0.3) # µ cannot be nonzero + Φ = m(0, 0.3) + ρ11 = m.rho(0, 0.3)[sites(1), sites(1)] - SA[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 - 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)