all solvers/slicers get the two-argument minimum copy

comments

add test

docs/src/tutorial/bandstructures.md

@@ -18,7 +18,7 @@ States:
 The above destructuring syntax assigns eigenvalues and eigenvectors to `eᵢ` and `ψᵢ`, respectively. The available eigensolvers and their options can be checked in the `EigenSolvers` docstrings.
 
 !!! warning "Arpack solver is not thread-safe"
-    `EigenSolver.Arpack` relies on a Fortran library that is not currently thread-safe. If you launch Julia with multiple threads, they will not be used with this specific solver. Otherwise Arpack would segfault.
+    `EigenSolvers.Arpack` relies on a Fortran library that is not currently thread-safe. If you launch Julia with multiple threads, they will not be used with this specific solver. Otherwise Arpack would segfault.
 
 We define a "bandstructure" of an `h::AbstractHamiltonian` as a linear interpolation of its eigenpairs over a portion of the Brillouin zone, which is discretized with a base mesh of `ϕᵢ` values. At each `ϕᵢ` of the base mesh, the Bloch matrix `h(ϕᵢ)` is diagonalized with `spectrum`. The adjacent eigenpairs `eⱼ(ϕᵢ), ψⱼ(ϕᵢ)` are then connected ("stitched") together into a number of band meshes with vertices `(ϕᵢ..., eⱼ(ϕᵢ))` by maximizing the overlap of adjacent `ψⱼ(ϕᵢ)` (since the bands should be continuuous). Degenerate eigenpairs are collected into a single node of the band mesh.
 

docs/src/tutorial/glossary.md

@@ -22,11 +22,11 @@ This is a summary of the type of objects you will be studying.
   - `ParametricHamiltonian`: like the above, but using a `ParametricModel`, which makes it dependent on a set of free parameters that can be efficiently adjusted.
 
   An `h::AbstractHamiltonian` can be used to produce a Bloch matrix `h(ϕ; params...)` of the same size as the number of orbitals per unit cell, where `ϕ = [ϕᵢ...]` are Bloch phases and `params` are values for the free parameters, if any.
-- **`Spectrum`**: the set of eigenpairs (eigenvalues and corresponding eigenvectors) of a Bloch matrix. It can be computed with a number of `EigenSolver`s.
+- **`Spectrum`**: the set of eigenpairs (eigenvalues and corresponding eigenvectors) of a Bloch matrix. It can be computed with a number of `EigenSolvers`.
 - **`Bandstructure`**: a collection of spectra, evaluated over a discrete mesh (typically a discretization of the Brillouin zone), that is connected to its mesh neighbors into a linearly-interpolated approximation of the `AbstractHamiltonian`'s bandstructure.
 - **`SelfEnergy`**: an operator `Σ(ω)` defined to act on a `LatticeSlice` of an `AbstractHamiltonian` that depends on energy `ω`.
 - **`OpenHamiltonian`**: an `AbstractHamiltonian` combined with a set of `SelfEnergies`
-- **`GreenFunction`**: an `OpenHamiltonian` combined with a `GreenSolver`, which is an algorithm that can in general compute the retarded or advanced Green function at any energy between any subset of sites of the underlying lattice.
+- **`GreenFunction`**: an `OpenHamiltonian` combined with an `AbstractGreenSolver`, which is an algorithm that can in general compute the retarded or advanced Green function at any energy between any subset of sites of the underlying lattice.
   - **`GreenSlice`**: a `GreenFunction` evaluated on a specific set of sites, but at an unspecified energy
   - **`GreenSolution`**: a `GreenFunction` evaluated at a specific energy, but on an unspecified set of sites
 - **`OrbitalSliceArray`**: an `AbstractArray` that can be indexed with a `SiteSelector`, in addition to the usual scalar indexing. Particular cases are `OrbitalSliceMatrix` and `OrbitalSliceVector`. This is the most common type obtained from `GreenFunction`s and observables obtained from them.

docs/src/tutorial/greenfunctions.md

@@ -9,7 +9,7 @@ where `i, j` are orbitals, `H` is the (possibly infinite) Hamiltonian matrix, an
 We split the problem of computing `Gʳᵢⱼ(ω)` of a given `h::AbstractHamiltonian` into four steps:
 
 1. Attach self-energies to `h` using the command `oh = attach(h, args...)`. This produces a new object `oh::OpenHamiltonian` with a number of `Contacts`, numbered `1` to `N`
-2. Use `g = greenfunction(oh, solver)` to build a `g::GreenFunction` representing `Gʳ` (at arbitrary `ω` and `i,j`), where `oh::OpenHamiltonian` and `solver::GreenSolver` (see `GreenSolvers` below for available solvers)
+2. Use `g = greenfunction(oh, solver)` to build a `g::GreenFunction` representing `Gʳ` (at arbitrary `ω` and `i,j`), where `oh::OpenHamiltonian` and `solver::AbstractGreenSolver` (see `GreenSolvers` below for available solvers)
 3. Evaluate `gω = g(ω; params...)` at fixed energy `ω` and model parameters, which produces a `gω::GreenSolution`
 4. Slice `gω[sᵢ, sⱼ]` or `gω[sᵢ] == gω[sᵢ, sᵢ]` to obtain `Gʳᵢⱼ(ω)` as a flat matrix, where `sᵢ, sⱼ` are either site selectors over sites spanning orbitals `i,j`, integers denoting contacts, `1` to `N`, or `:` denoting all contacts merged together.
 
@@ -48,7 +48,7 @@ julia> gω[cells = 1:2]  # we now ask for the Green function between orbitals in
   -0.2+0.846606im   -0.48-0.113394im    0.44+0.282715im  0.104-0.869285im
 ```
 
-Note that the result is a 4 x 4 matrix, because there are 2 orbitals (one per site) in each of the two unit cells. Note also that the Schur GreenSolver used here allows us to compute the Green function between distant cells with little overhead
+Note that the result is a 4 x 4 matrix, because there are 2 orbitals (one per site) in each of the two unit cells. Note also that the `GreenSolvers.Schur` used here allows us to compute the Green function between distant cells with little overhead
 ```julia
 julia> @time gω[cells = 1:2];
   0.000067 seconds (70 allocations: 6.844 KiB)
@@ -59,7 +59,7 @@ julia> @time gω[cells = (SA[10], SA[100000])];
 
 ## GreenSolvers
 
-The currently implemented `GreenSolver`s (abbreviated as `GS`) are the following
+The currently implemented `GreenSolvers` (abbreviated as `GS`) are the following
 
 - `GS.SparseLU()`
 

src/docstrings.jl

@@ -1490,7 +1490,7 @@ Simple selector of sites with given `site_indices` in a given cell at `cell_inde
 cellsites
 
 """
-    greenfunction(h::Union{AbstractHamiltonian,OpenHamiltonian}, solver::GreenSolver)
+    greenfunction(h::Union{AbstractHamiltonian,OpenHamiltonian}, solver::AbstractGreenSolver)
 
 Build a `g::GreenFunction` of Hamiltonian `h` using `solver`. See `GreenSolvers` for
 available solvers. If `solver` is not provided, a default solver is chosen automatically

src/greenfunction.jl

@@ -469,15 +469,16 @@ end
 ensure_mutable_matrix(m::SMatrix) = Matrix(m)
 ensure_mutable_matrix(m::AbstractMatrix) = m
 
-minimal_callsafe_copy(s::TMatrixSlicer, parentham) = TMatrixSlicer(
-    minimal_callsafe_copy(s.g0slicer, parentham), s.tmatrix, s.gcontacts, s.contactorbs)
+minimal_callsafe_copy(s::TMatrixSlicer, parentham, parentcontacts) = TMatrixSlicer(
+    minimal_callsafe_copy(s.g0slicer, parentham, parentcontacts),
+    s.tmatrix, s.gcontacts, s.contactorbs)
 
 Base.view(::NothingSlicer, i::Union{Integer,Colon}...) =
     internalerror("view(::NothingSlicer): unreachable reached")
 
 Base.getindex(::NothingSlicer, i::CellOrbitals...) = argerror("Slicer does not support generic indexing")
 
-minimal_callsafe_copy(s::NothingSlicer, parentham) = s
+minimal_callsafe_copy(s::NothingSlicer, parentham, parentcontacts) = s
 
 #endregion
 #endregion

src/solvers/green.jl

@@ -4,15 +4,15 @@
 #     - apply(solver, h::AbstractHamiltonian, c::Contacts) -> AppliedGreenSolver
 #   All new s::AppliedGreenSolver must implement (with Σblock a [possibly nested] tuple of MatrixBlock's)
 #      - s(ω, Σblocks, ::ContactOrbitals) -> GreenSlicer
-#      - minimal_callsafe_copy(s, parentham)  # injects any aliases from parentham
+#      - minimal_callsafe_copy(s, parentham, parentcontacts)  # injects aliases from parent
 #      - optional: needs_omega_shift(s) (has a `true` default fallback)
 #   A gs::GreenSlicer's allows to compute G[gi, gi´]::AbstractMatrix for indices gi
 #   To do this, it must implement contact slicing (unless it relies on TMatrixSlicer)
 #      - view(gs, ::Int, ::Int) -> g(ω; kw...) between specific contacts (has error fallback)
 #      - view(gs, ::Colon, ::Colon) -> g(ω; kw...) between all contacts (has error fallback)
 #   It must also implement generic slicing, and minimal copying
 #      - gs[i::CellOrbitals, j::CellOrbitals] -> must return a Matrix for type stability
-#      - minimal_callsafe_copy(gs, parentham)
+#      - minimal_callsafe_copy(gs, parentham, parentcontacts)
 #   The user-facing indexing API accepts:
 #      - i::Integer -> Sites of Contact number i
 #      - cellsites(cell::Tuple, sind::Int)::Subcell -> Single site in a cell
@@ -22,6 +22,7 @@
 #      - sel::SiteSelector ~ NamedTuple -> forms a LatticeSlice
 #   Optional: to properly plot boundaries, an ::AbstractGreenSolver may also implement
 #      - boundaries(s::AbstractGreenSolver) -> collection of (dir => cell)::Pair{Int,Int}
+#   Aliasing: Green solvers may only alias fields from the parent Hamiltonian and Contacts
 ############################################################################################
 
 module GreenSolvers

src/solvers/green/bands.jl

@@ -653,7 +653,7 @@ end
 #region ## API ##
 
 # Parent hamiltonian needs to be non-parametric, so no need to alias
-minimal_callsafe_copy(s::AppliedBandsGreenSolver, parentham) = s
+minimal_callsafe_copy(s::AppliedBandsGreenSolver, parentham, parentcontacts) = s
 
 needs_omega_shift(s::AppliedBandsGreenSolver) = false
 
@@ -667,7 +667,7 @@ boundaries(s::AppliedBandsGreenSolver) = (s.boundaryorbs.boundary,)
 #region ## apply ##
 
 function apply(s::GS.Bands,  h::AbstractHamiltonian{T,<:Any,L}, cs::Contacts) where {T,L}
-    L == 0 && argerror("Cannot use GreenSolver.Bands with 0D AbstractHamiltonians")
+    L == 0 && argerror("Cannot use GreenSolvers.Bands with 0D AbstractHamiltonians")
     ticks = s.bandsargs
     kw = s.bandskw
     b = bands(h, ticks...; kw..., projectors = true, split = false)
@@ -889,7 +889,7 @@ view_or_copy(ψ, rows::Union{Colon,AbstractRange}, cols::Union{Colon,AbstractRan
     view(ψ, rows, cols)
 view_or_copy(ψ, rows, cols) = ψ[rows, cols]
 
-minimal_callsafe_copy(s::BandsGreenSlicer, parentham) = s  # it is read-only
+minimal_callsafe_copy(s::BandsGreenSlicer, parentham, parentcontacts) = s  # it is read-only
 
 #endregion
 

src/solvers/green/kpm.jl

@@ -141,7 +141,7 @@ end
 moments(s::AppliedKPMGreenSolver) = s.moments
 
 # Parent ham needs to be non-parametric, so no need to alias
-minimal_callsafe_copy(s::AppliedKPMGreenSolver, parentham) = s
+minimal_callsafe_copy(s::AppliedKPMGreenSolver, parentham, parentcontacts) = s
 
 needs_omega_shift(s::AppliedKPMGreenSolver) = false
 

src/solvers/green/schur.jl

@@ -207,6 +207,8 @@ end
 checkmodes(whichmodes) = sum(whichmodes) == length(whichmodes) ÷ 2 ||
     argerror("Cannot differentiate retarded from advanced modes. Consider increasing imag(ω) or check that your Hamiltonian is Hermitian")
 
+# this does not include a parentcontacts because there are none
+# (SchurFactorsSolver is not an AppliedGreenSolver, so it may have different API)
 function minimal_callsafe_copy(s::SchurFactorsSolver, parentham)
     hm´, h0´, hp´ = nearest_cell_harmonics(parentham)
     s´ = SchurFactorsSolver(s.shift, hm´, h0´, hp´, s.l_leq_r, copy(s.iG),
@@ -373,12 +375,10 @@ function schur_openhams_types(fsolver, h, boundary)
 end
 
 apply(::GS.Schur, h::AbstractHamiltonian, cs::Contacts) =
-    argerror("Can only use GreenSolver.Schur with 1D AbstractHamiltonians")
-
-const Contacts0D{T,E,N,S} = Contacts{0,N,S,OrbitalSliceGrouped{T,E,0}}
+    argerror("Can only use GreenSolvers.Schur with 1D AbstractHamiltonians")
 
 const GFUnit{T,E,H,N,S} =
-    GreenFunction{T,E,0,AppliedSparseLUGreenSolver{Complex{T},Contacts0D{T,E,N,S}},H,Contacts0D{T,E,N,S}}
+    GreenFunction{T,E,0,AppliedSparseLUGreenSolver{Complex{T}},H,Contacts{0,N,S,OrbitalSliceGrouped{T,E,0}}}
 
 green_type(::H,::S) where {T,E,H<:AbstractHamiltonian{T,E},S} =
     GFUnit{T,E,H,1,Tuple{S}}
@@ -389,7 +389,7 @@ green_type(::H,::S1,::S2) where {T,E,H<:AbstractHamiltonian{T,E},S1,S2} =
 
 #region ## call API ##
 
-function minimal_callsafe_copy(s::AppliedSchurGreenSolver, parentham)
+function minimal_callsafe_copy(s::AppliedSchurGreenSolver, parentham, _)
     fsolver´ = minimal_callsafe_copy(s.fsolver, parentham)
     ohL´, ohR´, oh∞´, G, G∞ = schur_openhams_types(fsolver´, parentham, s.boundary)
     s´ = AppliedSchurGreenSolver{G,G∞}(fsolver´, s.boundary, ohL´, ohR´, oh∞´)
@@ -594,11 +594,11 @@ maybe_SMatrix(G::Matrix, rows::SVector{L}, cols::SVector{L´}) where {L,L´} = S
 maybe_SMatrix(G, rows, cols) = G
 
 # TODO: Perhaps too conservative
-function minimal_callsafe_copy(s::SchurGreenSlicer, parentham)
-    s´ = SchurGreenSlicer(s.ω, minimal_callsafe_copy(s.solver, parentham))
-    isdefined(s, :G₋₁₋₁)    && (s´.G₋₁₋₁    = minimal_callsafe_copy(s.G₋₁₋₁))
-    isdefined(s, :G₁₁)      && (s´.G₁₁      = minimal_callsafe_copy(s.G₁₁))
-    isdefined(s, :G∞₀₀)     && (s´.G∞₀₀     = minimal_callsafe_copy(s.G∞₀₀))
+function minimal_callsafe_copy(s::SchurGreenSlicer, parentham, parentcontacts)
+    s´ = SchurGreenSlicer(s.ω, minimal_callsafe_copy(s.solver, parentham, parentcontacts))
+    isdefined(s, :G₋₁₋₁)    && (s´.G₋₁₋₁    = minimal_callsafe_copy(s.G₋₁₋₁, parentham, parentcontacts))
+    isdefined(s, :G₁₁)      && (s´.G₁₁      = minimal_callsafe_copy(s.G₁₁, parentham, parentcontacts))
+    isdefined(s, :G∞₀₀)     && (s´.G∞₀₀     = minimal_callsafe_copy(s.G∞₀₀, parentham, parentcontacts))
     isdefined(s, :L´G∞₀₀)   && (s´.L´G∞₀₀   = copy(s.L´G∞₀₀))
     isdefined(s, :R´G∞₀₀)   && (s´.R´G∞₀₀   = copy(s.R´G∞₀₀))
     isdefined(s, :G₁₁L)     && (s´.G₁₁L     = copy(s.G₁₁L))

src/solvers/green/sparselu.jl

@@ -4,9 +4,9 @@
 #   ization, possibly using inverse-free self-energies (using extended sites).
 #region
 
-struct AppliedSparseLUGreenSolver{C,CT} <: AppliedGreenSolver
-    invgreen::InverseGreenBlockSparse{C}
-    contacts::CT
+# invgreen aliases contacts (they are not implemented through a TMatrixSlicer)
+struct AppliedSparseLUGreenSolver{C} <: AppliedGreenSolver
+    invgreen::InverseGreenBlockSparse{C}    # aliases parent contacts
 end
 
 mutable struct SparseLUGreenSlicer{C} <:GreenSlicer{C}
@@ -40,10 +40,9 @@ unitcellinds_contacts(s::SparseLUGreenSlicer, i::Integer) =
         argerror("Cannot access contact $i, there are $(length(s.unitcinds)) contacts")
 unitcellinds_contacts_merged(s::SparseLUGreenSlicer) = s.unitcindsall
 
-function minimal_callsafe_copy(s::AppliedSparseLUGreenSolver, parentham)
-    contacts´ = minimal_callsafe_copy(s.contacts)
-    invgreen´ = inverse_green(parentham, contacts´)
-    return AppliedSparseLUGreenSolver(invgreen´, contacts´)
+function minimal_callsafe_copy(s::AppliedSparseLUGreenSolver, parentham, parentcontacts)
+    invgreen´ = inverse_green(parentham, parentcontacts)
+    return AppliedSparseLUGreenSolver(invgreen´)
 end
 
 #endregion
@@ -52,11 +51,11 @@ end
 
 function apply(::GS.SparseLU, h::AbstractHamiltonian0D, cs::Contacts)
     invgreen = inverse_green(h, cs)
-    return AppliedSparseLUGreenSolver(invgreen, cs)
+    return AppliedSparseLUGreenSolver(invgreen)
 end
 
 apply(::GS.SparseLU, h::AbstractHamiltonian, cs::Contacts) =
-    argerror("Can only use GreenSolver.SparseLU with 0D AbstractHamiltonians")
+    argerror("Can only use GreenSolvers.SparseLU with 0D AbstractHamiltonians")
 
 #endregion
 
@@ -146,13 +145,12 @@ Base.getindex(s::SparseLUGreenSlicer, i::CellOrbitals, j::CellOrbitals) = copy(v
 # the lazy unitg field only aliases source64 or a copy of it. It is not necessary to
 # maintain the alias, as this is just a prealloc for a full-cell slice. We don't even need
 # to copy it, since once it is computed, it is never modified, only read
-function minimal_callsafe_copy(s::SparseLUGreenSlicer{C}, parentham) where {C}
+function minimal_callsafe_copy(s::SparseLUGreenSlicer{C}, parentham, parentcontacts) where {C}
     s´ = SparseLUGreenSlicer{C}(s.fact, s.nonextrng, s.unitcinds, s.unitcindsall, copy(s.source64))
     isdefined(s, :unitg) && (s´.unitg = s.unitg)
     return s´
 end
 
-
 #endregion
 
 #endregion

src/solvers/green/spectrum.jl

@@ -21,10 +21,11 @@ end
 spectrum(s::AppliedSpectrumGreenSolver) = s.spectrum
 
 # Parent ham needs to be non-parametric, so no need to alias
-minimal_callsafe_copy(s::AppliedSpectrumGreenSolver, parentham) =
+minimal_callsafe_copy(s::AppliedSpectrumGreenSolver, parentham, parentcontacts) =
     AppliedSpectrumGreenSolver(s.spectrum)
 
-minimal_callsafe_copy(s::SpectrumGreenSlicer, parentham) = SpectrumGreenSlicer(s.ω, s.solver)
+minimal_callsafe_copy(s::SpectrumGreenSlicer, parentham, parentcontacts) =
+    SpectrumGreenSlicer(s.ω, s.solver)
 
 #endregion
 

src/solvers/selfenergy.jl

@@ -5,14 +5,15 @@
 #     - call!(s::ExtendedSelfEnergySolver, ω; params...) -> (Vᵣₑ, gₑₑ⁻¹, Vₑᵣ) AbsMats
 #         With the extended case, the equivalent Σreg reads Σreg = VᵣₑgₑₑVₑᵣ
 #     - call!_output(s::AbstractSelfEnergySolver) -> object returned by call!(s, ω; kw...)
-#     - minimal_callsafe_copy(s::AbstractSelfEnergySolver)
+#     - minimal_callsafe_copy(s::AbstractSelfEnergySolver)  # no aliasing
 #   These AbstractMatrices are flat, defined on the LatticeSlice in parent SelfEnergy
 #       Note: `params` are only needed in cases where s adds new parameters that must be
 #       applied (e.g. SelfEnergyModelSolver). Otherwise one must assume that any parent
 #       ParametricHamiltonian to GreenFunction has already been call!-ed before calling s.
 #   Optional: AbstractSelfEnergySolver's can also implement `selfenergy_plottables`
 #     - selfenergy_plottables(s::AbstractSelfEnergySolver, parent_latslice)
 #       -> collection of tuples to be passed to plotlattice!(axis, tup...) for visualization
+#   Aliasing: AbstractSelfEnergySolver's are not allowed to alias anything from outside
 ############################################################################################
 
 ############################################################################################

src/solvers/selfenergy/schur.jl

@@ -231,7 +231,7 @@ function minimal_callsafe_copy(s::SelfEnergyCouplingSchurSolver)
     gunit´ = minimal_callsafe_copy(s.gunit)
     s´ = SelfEnergyCouplingSchurSolver(gunit´, hcoupling´,
         minimal_callsafe_copy(s.V´, hcoupling´),
-        inverse_green(solver(gunit´)),
+        inverse_green(solver(gunit´)),      # alias of the invgreen from the solver
         minimal_callsafe_copy(s.V, hcoupling´))
     return s´
 end

src/types.jl

@@ -1273,10 +1273,6 @@ function update!(s::InverseGreenBlockSparse, ω)
     return update!(bsm)
 end
 
-minimal_callsafe_copy(s::InverseGreenBlockSparse, parentham, contacts) =
-    InverseGreenBlockSparse(minimal_callsafe_copy(s.mat, parentham, contacts), s.nonextrng, s.unitcinds,
-    s.unitcindsall, copy(s.source))
-
 #endregion
 #endregion
 
@@ -2186,15 +2182,17 @@ copy_lattice(g::GreenSlice) = GreenSlice(
 
 function minimal_callsafe_copy(g::GreenFunction)
     parent´ = minimal_callsafe_copy(g.parent)
-    solver´ = minimal_callsafe_copy(g.solver, parent´)
     contacts´ = minimal_callsafe_copy(g.contacts)
+    solver´ = minimal_callsafe_copy(g.solver, parent´, contacts´)
     return GreenFunction(parent´, solver´, contacts´)
 end
 
 function minimal_callsafe_copy(g::GreenSolution)
     parentg´ = minimal_callsafe_copy(g.parent)
-    g´ = GreenSolution(parentg´,
-        minimal_callsafe_copy(g.slicer, hamiltonian(parentg´)), g.contactΣs, g.contactorbs)
+    parentham = hamiltonian(parentg´)
+    parentcontacts = contacts(parentg´)
+    slicer´ = minimal_callsafe_copy(g.slicer, parentham, parentcontacts)
+    g´ = GreenSolution(parentg´, slicer´, g.contactΣs, g.contactorbs)
     return g´
 end
 

test/test_greenfunction.jl

@@ -414,4 +414,11 @@ end
     gs = Quantica.call!(g, 0.4, q = 0.1)
     parent(g)(q = 2)
     @test_broken gs[cells = 0] == m
+    # Ensure that g.solver.invgreen alias of parent contacts is not broken by minimal_callsafe_copy
+    glead = LP.linear() |> @onsite((; o = 1) -> o) - hopping(1) |> greenfunction(GS.Schur(boundary = 0));
+    g = LP.linear() |> -hopping(1) |> supercell |> attach(glead) |> greenfunction;
+    g´ = Quantica.minimal_callsafe_copy(g);
+    @test g´(0, o = 0)[] == g(0, o = 0)[]
+    @test g´(0, o = 1)[] == g(0, o = 1)[]
+
 end