From b9d3f3ffeb40f0236b3d8887d1f9e05628a937f5 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Sat, 1 Feb 2025 22:08:48 +0100 Subject: [PATCH 01/10] refactor build_slicer to access params --- src/Quantica.jl | 6 +- src/greenfunction.jl | 5 +- src/solvers/{eigen.jl => eigensolvers.jl} | 0 src/solvers/green/bands.jl | 2 +- src/solvers/green/internal.jl | 2 +- src/solvers/green/kpm.jl | 2 +- src/solvers/green/schur.jl | 98 +++++++++++++++++-- src/solvers/green/sparselu.jl | 4 +- src/solvers/green/spectrum.jl | 2 +- src/solvers/{green.jl => greensolvers.jl} | 5 +- src/solvers/selfenergy/schur.jl | 1 + .../{selfenergy.jl => selfenergysolvers.jl} | 0 12 files changed, 108 insertions(+), 19 deletions(-) rename src/solvers/{eigen.jl => eigensolvers.jl} (100%) rename src/solvers/{green.jl => greensolvers.jl} (92%) rename src/solvers/{selfenergy.jl => selfenergysolvers.jl} (100%) diff --git a/src/Quantica.jl b/src/Quantica.jl index 657f7562..d8c00843 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -74,9 +74,9 @@ include("sanitizers.jl") # Solvers -include("solvers/eigen.jl") -include("solvers/green.jl") -include("solvers/selfenergy.jl") +include("solvers/eigensolvers.jl") +include("solvers/greensolvers.jl") +include("solvers/selfenergysolvers.jl") # Presets include("presets/regions.jl") diff --git a/src/greenfunction.jl b/src/greenfunction.jl index 2e343a12..a738386c 100644 --- a/src/greenfunction.jl +++ b/src/greenfunction.jl @@ -54,12 +54,11 @@ function call!(g::GreenFunction{T}, ω::T; params...) where {T} end function call!(g::GreenFunction{T}, ω::Complex{T}; params...) where {T} - h = parent(g) # not hamiltonian(h). We want the ParametricHamiltonian if it exists. + gsolver = solver(g) contacts´ = contacts(g) - call!(h; params...) Σblocks = call!(contacts´, ω; params...) corbs = contactorbitals(contacts´) - slicer = solver(g)(ω, Σblocks, corbs) + slicer = build_slicer(gsolver, g, ω, Σblocks, corbs; params...) return GreenSolution(g, slicer, Σblocks, corbs) end diff --git a/src/solvers/eigen.jl b/src/solvers/eigensolvers.jl similarity index 100% rename from src/solvers/eigen.jl rename to src/solvers/eigensolvers.jl diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index e0ab9170..94021f8e 100644 --- a/src/solvers/green/bands.jl +++ b/src/solvers/green/bands.jl @@ -708,7 +708,7 @@ end #region ## call ## -function (s::AppliedBandsGreenSolver)(ω, Σblocks, corbitals) +function build_slicer(s::AppliedBandsGreenSolver, g, ω, Σblocks, corbitals; params...) g0slicer = BandsGreenSlicer(complex(ω), s) gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) return gslicer diff --git a/src/solvers/green/internal.jl b/src/solvers/green/internal.jl index 1439f671..8bc1f12b 100644 --- a/src/solvers/green/internal.jl +++ b/src/solvers/green/internal.jl @@ -26,7 +26,7 @@ end ## GreenFunction API ## -function (s::AppliedModelGreenSolver)(ω::C, Σblocks, contactorbs) where {C} +function build_slicer(s::AppliedModelGreenSolver, g, ω::C, Σblocks, contactorbs; params...) where {C} n = norbitals(contactorbs) fω = s.f(ω) g0contacts = Matrix{C}(undef, n, n) diff --git a/src/solvers/green/kpm.jl b/src/solvers/green/kpm.jl index 4b825c49..89aa7c85 100644 --- a/src/solvers/green/kpm.jl +++ b/src/solvers/green/kpm.jl @@ -195,7 +195,7 @@ end #region ## call ## -function (s::AppliedKPMGreenSolver{T})(ω, Σblocks, corbitals) where {T} +function build_slicer(s::AppliedKPMGreenSolver{T}, g, ω, Σblocks, corbitals; params...) where {T} g0contacts = KPMgreen(s.momenta, ω, s.bandCH) # We rely on TMatrixSlicer to incorporate contact self-energies, and to slice contacts gslicer = TMatrixSlicer(g0contacts, Σblocks, corbitals) diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 4ff3acbc..863a0945 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -1,3 +1,23 @@ +# NOTE: this solver is probably too convoluted and could benefit from a refactor to make it +# more maintainable. The priority here was performance (minimal number of computations) +# The formalism is based on currently unpublished notes (scattering.pdf) + +# Design overview: The G₋₁₋₁, G₁₁ and G∞₀₀ slicers inside SchurGreenSlicer correspond to +# boundary cells in seminf and origin cell in an infinite 1D Hamiltonian, based on a LU +# solver. These are not computed upon creating SchurGreenSlicer, since not all are necessary +# for any given slice. Same for their parent gL, dR and g∞ 0D unitcell GreenFunctions inside +# AppliedSchurGreenSolver. Their respective OpenHamiltonians contain selfenergies with solver +# type SelfEnergySchurSolver which themselves all alias the same SchurFactorsSolver devoted +# to computing objects using the deflated Schur algorithm on which the selfenergies depend. +# When building SchurGreenSlicer with a given ω and params we know for sure that we will +# need the SchurFactorsSolver objects for that ω, but not which of G₋₁₋₁, G₁₁ or G∞₀₀. +# We also don't want to compute the SchurFactorsSolver objects more than once for a given ω +# Hence, we do a call!(fsolver::SchurFactorsSolver, ω) upon constructing SchurGreenSlicer, +# but whenever building the uninitialized slicers G₁₁ etc we do a skipsolve_internal = true +# to avoid recomputing SchurFactorsSolver. Also, since the unitcell Hamiltonians harmonics +# h0, hm, hp in SchurFactorsSolver alias those of hamiltonian(g::GreenFunction), we need to +# do call!(parent(g); params...) upon constructing SchurFactorsSolver too. + ############################################################################################ # SchurFactorsSolver - see scattering.pdf notes for derivations # Auxiliary functions for AppliedSchurGreenSolverSolver @@ -291,7 +311,7 @@ end #endregion ############################################################################################ -# AppliedSchurGreenSolver +# AppliedSchurGreenSolver in 1D #region # Mutable: we delay initialization of some fields until they are first needed (which may be never) @@ -352,7 +372,7 @@ end function apply(s::GS.Schur, h::AbstractHamiltonian1D{T}, contacts::Contacts) where {T} h´ = hamiltonian(h) - fsolver = SchurFactorsSolver(h´, s.shift) + fsolver = SchurFactorsSolver(h´, s.shift) # aliasing of h´ boundary = T(round(only(s.boundary))) ohL, ohR, oh∞, G, G∞ = schur_openhams_types(fsolver, h, boundary) solver = AppliedSchurGreenSolver{G,G∞}(fsolver, boundary, ohL, ohR, oh∞) @@ -401,9 +421,10 @@ function minimal_callsafe_copy(s::AppliedSchurGreenSolver, parentham, _) return s´ end -function (s::AppliedSchurGreenSolver)(ω, Σblocks, corbitals) - # call! fsolver once for all the g's - call!(s.fsolver, ω) +function build_slicer(s::AppliedSchurGreenSolver, g, ω, Σblocks, corbitals; params...) + # overwrites hamiltonian(g) with params whenever parent(g) isa ParametricHamiltonian + # Necessary because its harmonics are aliased in the SchurFactorSolver inside gsolver + call!(parent(g); params...) g0slicer = SchurGreenSlicer(ω, s) gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) return gslicer @@ -430,6 +451,7 @@ end #region # We delay initialization of most fields until they are first needed (which may be never) +# should not have any contacts (we defer to TMatrixSlicer for that) mutable struct SchurGreenSlicer{C,A<:AppliedSchurGreenSolver} <: GreenSlicer{C} ω::C solver::A @@ -437,7 +459,7 @@ mutable struct SchurGreenSlicer{C,A<:AppliedSchurGreenSolver} <: GreenSlicer{C} L::Matrix{C} R::Matrix{C} G₋₁₋₁::SparseLUGreenSlicer{C} # These are independent of parent hamiltonian - G₁₁::SparseLUGreenSlicer{C} # as they only rely on call!_output(fsolver) + G₁₁::SparseLUGreenSlicer{C} # as they only rely on call!_output(solver.fsolver) G∞₀₀::SparseLUGreenSlicer{C} # which is updated after call!(solver.fsolver, ω) L´G∞₀₀::Matrix{C} R´G∞₀₀::Matrix{C} @@ -446,6 +468,9 @@ mutable struct SchurGreenSlicer{C,A<:AppliedSchurGreenSolver} <: GreenSlicer{C} R´G₁₁L::Matrix{C} L´G₋₁₋₁R::Matrix{C} function SchurGreenSlicer{C,A}(ω, solver) where {C,A} + # call! the expensive fsolver only once to compute Schur factors, necessary to + # evaluate unitcell selfenergies in OpenHamiltonians of solver + call!(solver.fsolver, ω) s = new() s.ω = ω s.solver = solver @@ -471,6 +496,8 @@ function Base.getproperty(s::SchurGreenSlicer, f::Symbol) # gs[sites...] will be call!-ing e.g. solver.g∞ with the wrong h0 (the one from # params´...). However, if `gs = g(ω; params...)` a copy was made, so it is safe. if f == :G₋₁₋₁ + # we ran SchurFactorSolver when constructing s, so skipsolve_internal = true + # to avoid running it again s.G₋₁₋₁ = slicer(call!(solver.gL, s.ω; skipsolve_internal = true)) elseif f == :G₁₁ s.G₁₁ = slicer(call!(solver.gR, s.ω; skipsolve_internal = true)) @@ -721,4 +748,63 @@ end #endregion + +############################################################################################ +# AppliedSchurGreenSolver in 2D +#region + +# struct AppliedSchurGreenSolverND{L,S<:AppliedSchurGreenSolver,F} <: AppliedGreenSolver +# solver1D::S +# wrapped_axes::NTuple{L,Int} +# phase_func::F # phase_func(ϕ) = (; param_name = ϕ) +# end + +# const GreenFunctionSchurND{T,L} = GreenFunction{T,<:Any,L,<:AppliedSchurGreenSolverND} + +# # should not have any contacts (we defer to TMatrixSlicer for that) +# struct SchurGreenSlicerND{C,S<:GreenSlicer{C},S<:AppliedSchurGreenSolverND} <: GreenSlicer{C} +# ω::C +# g0slicer1D::S +# solver::S +# end + +# function apply(s::GS.Schur, h::AbstractHamiltonian, contacts::Contacts) +# L = latdim(h) +# wrapped_axes = inds_complement(Val(L), (s.axis,)) +# h´ = @torus(h, wrapped_axes, ϕ_internal) +# phase_func(ϕ_internal) = (; ϕ_internal) +# solver1D = apply(s, h´, contacts) +# return AppliedSchurGreenSolverND(solver1D, wrapped_axes, phase_func) +# end + +# #region ## API ## + +# minimal_callsafe_copy(s::AppliedSchurGreenSolver, args...) = +# AppliedSchurGreenSolverND(minimal_callsafe_copy(s.solver1D, args...), s.wrapped_axes, s.phase_func) + +# needs_omega_shift(s::AppliedSchurGreenSolverND) = needs_omega_shift(s.solver1D) + +# (s::AppliedSchurGreenSolverND)(args...) = SchurGreenSlicerND(s.solver1D(args...), s) + +# # hook in one step above gsolver(ω, Σblocks, corbs), to use params and g - see greenfunction.jl +# function build_slicer(g, gsolver::AppliedSchurGreenSolverND, ω, Σblocks, corbs; params...) +# # overwrites hamiltonian(g) with params whenever parent(g) isa ParametricHamiltonian +# # Necessary because its harmonics are aliased in the SchurFactorSolver inside gsolver +# call!(parent(g); params...) +# g0slicer1D = SchurGreenSlicer(ω, gsolver.solver1D) +# g0slicer = SchurGreenSlicerND(ω, g0slicer1D, s) +# gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) +# return gslicer +# end + +# function Base.getindex(s::SchurGreenSlicerND, i::CellOrbitals, j::CellOrbitals) +# dn = cell(i) - cell(j) +# s1D = s.g0slicer1D +# integrand(ϕ...) = s1D[i,j] +# end + +#endregion + +#endregion + #endregion top diff --git a/src/solvers/green/sparselu.jl b/src/solvers/green/sparselu.jl index f97ad444..2e57159c 100644 --- a/src/solvers/green/sparselu.jl +++ b/src/solvers/green/sparselu.jl @@ -62,7 +62,9 @@ apply(::GS.SparseLU, h::AbstractHamiltonian, cs::Contacts) = #region ## call ## # Σblocks and contactorbitals are not used here, because they are already inside invgreen -function (s::AppliedSparseLUGreenSolver{C})(ω, Σblocks, contactorbitals) where {C} +function build_slicer(s::AppliedSparseLUGreenSolver{C}, g, ω, Σblocks, contactorbitals; params...) where {C} + # We must apply params to hamiltonian(g) because its base harmonic is aliased into invgreen as a MatrixBlock + call!(parent(g); params...) invgreen = s.invgreen nonextrng = orbrange(invgreen) unitcinds = invgreen.unitcinds diff --git a/src/solvers/green/spectrum.jl b/src/solvers/green/spectrum.jl index caad8336..20add2ff 100644 --- a/src/solvers/green/spectrum.jl +++ b/src/solvers/green/spectrum.jl @@ -44,7 +44,7 @@ apply(::GS.Spectrum, h::AbstractHamiltonian, ::Contacts) = #region ## call ## -function (s::AppliedSpectrumGreenSolver)(ω, Σblocks, corbitals) +function build_slicer(s::AppliedSpectrumGreenSolver, g, ω, Σblocks, corbitals; params...) g0slicer = SpectrumGreenSlicer(complex(ω), s) gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) return gslicer diff --git a/src/solvers/green.jl b/src/solvers/greensolvers.jl similarity index 92% rename from src/solvers/green.jl rename to src/solvers/greensolvers.jl index 2a226a4f..6c8df84d 100644 --- a/src/solvers/green.jl +++ b/src/solvers/greensolvers.jl @@ -3,7 +3,7 @@ # All new solver::AbstractGreenSolver must live in the GreenSolvers module, and must implement # - 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 +# - build_slicer(s, g::GreenFunction, ω, Σblocks, ::ContactOrbitals; params...) -> GreenSlicer # - 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 @@ -35,9 +35,10 @@ struct SparseLU <:AbstractGreenSolver end struct Schur{T<:AbstractFloat} <: AbstractGreenSolver shift::T # Tunable parameter in algorithm, see Ω in scattering.pdf boundary::T # Cell index for boundary (float to allow boundary at Inf) + axis::Int # free axis to use for n-dimensional AbstractHamiltonians end -Schur(; shift = 1.0, boundary = Inf) = Schur(shift, float(boundary)) +Schur(; shift = 1.0, boundary = Inf, axis = 1) = Schur(shift, float(boundary), axis) struct KPM{B<:Union{Missing,NTuple{2}},A} <: AbstractGreenSolver order::Int diff --git a/src/solvers/selfenergy/schur.jl b/src/solvers/selfenergy/schur.jl index 13916295..a6ea6147 100644 --- a/src/solvers/selfenergy/schur.jl +++ b/src/solvers/selfenergy/schur.jl @@ -1,6 +1,7 @@ ############################################################################################ # SelfEnergy(h, glead::GreenFunctionSchurEmptyLead; kw...) # Extended self energy solver for deflated ΣL or ΣR Schur factors of lead unitcell +# See notes on solvers/green/schur.jl #region struct SelfEnergySchurSolver{T,B,V<:Union{Missing,Vector{Int}},H<:AbstractHamiltonian} <: ExtendedSelfEnergySolver diff --git a/src/solvers/selfenergy.jl b/src/solvers/selfenergysolvers.jl similarity index 100% rename from src/solvers/selfenergy.jl rename to src/solvers/selfenergysolvers.jl From 165f6b593c3070a7cddfb90b6c549d517d1594fd Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Mon, 3 Feb 2025 16:50:46 +0100 Subject: [PATCH 02/10] SchurGreenSlicer2D with transverse quadgk docstring fix --- src/docstrings.jl | 4 +- src/solvers/green/schur.jl | 144 ++++++++++++++++++-------------- src/solvers/greensolvers.jl | 6 +- src/solvers/selfenergy/schur.jl | 16 ++-- 4 files changed, 95 insertions(+), 75 deletions(-) diff --git a/src/docstrings.jl b/src/docstrings.jl index adbacb11..a6487bbc 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -2396,7 +2396,7 @@ OrbitalSliceGrouped{Float64,2,1} : collection of subcells of orbitals (grouped b Total sites : 4 julia> siteindexdict(a) -4-element Dictionaries.Dictionary{Quantica.CellIndices{1, Int64, Quantica.SiteLike}, UnitRange{Int64}} +4-element Dictionaries.Dictionary{Quantica.CellIndices{1, Int64, Quantica.SiteLike}, UnitRange{Int64}}: CellSites{1,Int64} : 1 site in cell zero Sites : 1 │ 1:2 CellSites{1,Int64} : 1 site in cell zero @@ -2670,7 +2670,7 @@ Compute the minimal gap around `µ`, see `Quantica.gaps` gap """ - Quantica.decay_lengths(g::GreenFunctionSchurLead, µ = 0; reverse = false) + Quantica.decay_lengths(g::GreenFunctionSchurLead1D, µ = 0; reverse = false) Quantica.decay_lengths(h::AbstractHamiltonian1D, µ = 0; reverse = false) Compute the decay lengths of evanescent modes of a 1D `AbstractHamiltonian` `h` or a 1D diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 863a0945..bb9a8fcf 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -101,7 +101,7 @@ function nearest_cell_harmonics(h) dn == SA[0] || dn == SA[1] || dn == SA[-1] end is_nearest || - argerror("Too many or too few harmonics. Perhaps try `supercell` to ensure strictly nearest-cell harmonics.") + argerror("Too many or too few harmonics (need 3, got $(length(harmonics(h)))). 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") @@ -335,8 +335,8 @@ mutable struct AppliedSchurGreenSolver{T,B,O,O∞,G,G∞} <: AppliedGreenSolver end end -const GreenFunctionSchurEmptyLead{T,E} = GreenFunction{T,E,1,<:AppliedSchurGreenSolver,<:Any,<:EmptyContacts} -const GreenFunctionSchurLead{T,E} = GreenFunction{T,E,1,<:AppliedSchurGreenSolver,<:Any,<:Any} +const GreenFunctionSchurEmptyLead1D{T,E} = GreenFunction{T,E,1,<:AppliedSchurGreenSolver,<:Any,<:EmptyContacts} +const GreenFunctionSchurLead1D{T,E} = GreenFunction{T,E,1,<:AppliedSchurGreenSolver,<:Any,<:Any} AppliedSchurGreenSolver{G,G∞}(fsolver::SchurFactorsSolver{T,B}, boundary, ohL::O, ohR::O, oh∞::O∞) where {T,B,O,O∞,G,G∞} = AppliedSchurGreenSolver{T,B,O,O∞,G,G∞}(fsolver, boundary, ohL, ohR, oh∞) @@ -370,7 +370,7 @@ end #region ## apply ## -function apply(s::GS.Schur, h::AbstractHamiltonian1D{T}, contacts::Contacts) where {T} +function apply(s::GS.Schur, h::AbstractHamiltonian1D{T}, contacts) where {T} h´ = hamiltonian(h) fsolver = SchurFactorsSolver(h´, s.shift) # aliasing of h´ boundary = T(round(only(s.boundary))) @@ -397,9 +397,6 @@ function schur_openhams_types(fsolver, h, boundary) return ohL, ohR, oh∞, G, G∞ end -apply(::GS.Schur, h::AbstractHamiltonian, cs::Contacts) = - argerror("Can only use GreenSolvers.Schur with 1D AbstractHamiltonians") - const GFUnit{T,E,H,N,S} = GreenFunction{T,E,0,AppliedSparseLUGreenSolver{Complex{T}},H,Contacts{0,N,S,OrbitalSliceGrouped{T,E,0}}} @@ -423,7 +420,7 @@ end function build_slicer(s::AppliedSchurGreenSolver, g, ω, Σblocks, corbitals; params...) # overwrites hamiltonian(g) with params whenever parent(g) isa ParametricHamiltonian - # Necessary because its harmonics are aliased in the SchurFactorSolver inside gsolver + # Necessary because its harmonics are aliased in the SchurFactorSolver inside s call!(parent(g); params...) g0slicer = SchurGreenSlicer(ω, s) gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) @@ -470,6 +467,7 @@ mutable struct SchurGreenSlicer{C,A<:AppliedSchurGreenSolver} <: GreenSlicer{C} function SchurGreenSlicer{C,A}(ω, solver) where {C,A} # call! the expensive fsolver only once to compute Schur factors, necessary to # evaluate unitcell selfenergies in OpenHamiltonians of solver + # (see skipsolve_internal = true further below) call!(solver.fsolver, ω) s = new() s.ω = ω @@ -492,9 +490,9 @@ function Base.getproperty(s::SchurGreenSlicer, f::Symbol) d = size(s.L, 2) # Issue #268: the result of the following call!'s depends on the current value of h0 # which aliases the parent h. This is only a problem if `s` was obtained through - # `gs = call!(g, ω; params...)`. In that case, doing call!(g, ω; params´...) before - # gs[sites...] will be call!-ing e.g. solver.g∞ with the wrong h0 (the one from - # params´...). However, if `gs = g(ω; params...)` a copy was made, so it is safe. + # `gω = call!(g, ω; params...)`. In that case, doing call!(g, ω; params´...) before + # gω[sites...] will be call!-ing e.g. solver.g∞ with the wrong h0 (the one from + # params´...). However, if `gω = g(ω; params...)` a copy was made, so it is safe. if f == :G₋₁₋₁ # we ran SchurFactorSolver when constructing s, so skipsolve_internal = true # to avoid running it again @@ -646,10 +644,10 @@ end # computes schur_eigenvalues of all lead modes #region -schur_eigvals(g::GreenFunctionSchurLead, ω::Real; params...) = +schur_eigvals(g::GreenFunctionSchurLead1D, ω::Real; params...) = schur_eigvals(g, retarded_omega(ω, solver(g)); params...) -function schur_eigvals(g::GreenFunctionSchurLead, ω::Complex; params...) +function schur_eigvals(g::GreenFunctionSchurLead1D, ω::Complex; params...) h = parent(g) # get the (Parametric)Hamiltonian from g call!(h; params...) # update the (Parametric)Hamiltonian with the params sf = g.solver.fsolver # obtain the SchurFactorSolver that computes the AB pencil @@ -671,7 +669,7 @@ propagating_eigvals(g, ω, margin = 0; params...) = filter!(λ -> abs(λ) ≈ 1, schur_eigvals(g, ω; params...)) : filter!(λ -> 1 - margin < abs(λ) < 1 + margin, schur_eigvals(g, ω; params...)) -function decay_lengths(g::GreenFunctionSchurLead, args...; reverse = false, params...) +function decay_lengths(g::GreenFunctionSchurLead1D, args...; reverse = false, params...) λs = reverse ? advanced_eigvals(g, args...; params...) : retarded_eigvals(g, args...; params...) ls = @. -1/log(abs(λs)) # compute the decay lengths in units of a0 return ls @@ -705,6 +703,7 @@ function densitymatrix(s::AppliedSchurGreenSolver, gs::GreenSlice; atol = 1e-7, hmat = similar_Array(hamiltonian(g)) psis = similar(hmat) orbaxes = orbrows(gs), orbcols(gs) + # possibly override the default quadgk_opts in solver opts = (; atol, quadgk_opts...) solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, opts) return DensityMatrix(solver, gs) @@ -753,55 +752,74 @@ end # AppliedSchurGreenSolver in 2D #region -# struct AppliedSchurGreenSolverND{L,S<:AppliedSchurGreenSolver,F} <: AppliedGreenSolver -# solver1D::S -# wrapped_axes::NTuple{L,Int} -# phase_func::F # phase_func(ϕ) = (; param_name = ϕ) -# end - -# const GreenFunctionSchurND{T,L} = GreenFunction{T,<:Any,L,<:AppliedSchurGreenSolverND} - -# # should not have any contacts (we defer to TMatrixSlicer for that) -# struct SchurGreenSlicerND{C,S<:GreenSlicer{C},S<:AppliedSchurGreenSolverND} <: GreenSlicer{C} -# ω::C -# g0slicer1D::S -# solver::S -# end - -# function apply(s::GS.Schur, h::AbstractHamiltonian, contacts::Contacts) -# L = latdim(h) -# wrapped_axes = inds_complement(Val(L), (s.axis,)) -# h´ = @torus(h, wrapped_axes, ϕ_internal) -# phase_func(ϕ_internal) = (; ϕ_internal) -# solver1D = apply(s, h´, contacts) -# return AppliedSchurGreenSolverND(solver1D, wrapped_axes, phase_func) -# end - -# #region ## API ## - -# minimal_callsafe_copy(s::AppliedSchurGreenSolver, args...) = -# AppliedSchurGreenSolverND(minimal_callsafe_copy(s.solver1D, args...), s.wrapped_axes, s.phase_func) - -# needs_omega_shift(s::AppliedSchurGreenSolverND) = needs_omega_shift(s.solver1D) - -# (s::AppliedSchurGreenSolverND)(args...) = SchurGreenSlicerND(s.solver1D(args...), s) - -# # hook in one step above gsolver(ω, Σblocks, corbs), to use params and g - see greenfunction.jl -# function build_slicer(g, gsolver::AppliedSchurGreenSolverND, ω, Σblocks, corbs; params...) -# # overwrites hamiltonian(g) with params whenever parent(g) isa ParametricHamiltonian -# # Necessary because its harmonics are aliased in the SchurFactorSolver inside gsolver -# call!(parent(g); params...) -# g0slicer1D = SchurGreenSlicer(ω, gsolver.solver1D) -# g0slicer = SchurGreenSlicerND(ω, g0slicer1D, s) -# gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) -# return gslicer -# end - -# function Base.getindex(s::SchurGreenSlicerND, i::CellOrbitals, j::CellOrbitals) -# dn = cell(i) - cell(j) -# s1D = s.g0slicer1D -# integrand(ϕ...) = s1D[i,j] -# end +struct AppliedSchurGreenSolver2D{L,H<:AbstractHamiltonian1D,S<:AppliedSchurGreenSolver,F,P<:NamedTuple} <: AppliedGreenSolver + h1D::H + solver1D::S + axis1D::Int + wrapped_axes::SVector{L,Int} + phase_func::F # phase_func(ϕ) = (; param_name = ϕ) + integrate_opts::P +end + +const GreenFunctionSchur2D{T,L} = GreenFunction{T,<:Any,L,<:AppliedSchurGreenSolver2D} + +# should not have any contacts (we defer to TMatrixSlicer for that) +struct SchurGreenSlicer2D{C,S<:AppliedSchurGreenSolver2D,F<:Function} <: GreenSlicer{C} + ω::C + solver::S + slicer_generator::F +end + +function apply(s::GS.Schur, h::AbstractHamiltonian, _) + L = latdim(h) + L == 2 || + argerror("GreenSolvers.Schur currently only implemented for 1D and 2D AbstractHamiltonians") + axis1D = s.axis + wrapped_axes = sanitize_SVector(inds_complement(Val(L), (axis1D,))) + h1D = @torus(h, wrapped_axes, ϕ_internal) + phase_func(ϕ_internal) = (; ϕ_internal) + # we don't pass contacts to solver1D. They will be applied with T-matrix slicer + solver1D = apply(s, h1D, missing) + return AppliedSchurGreenSolver2D(h1D, solver1D, axis1D, wrapped_axes, phase_func, s.integrate_opts) +end + +#region ## API ## + +function minimal_callsafe_copy(s::AppliedSchurGreenSolver2D, args...) + h1D´ = minimal_callsafe_copy(s.h1D) + solver1D´ = minimal_callsafe_copy(s.solver1D, h1D´, missing) + return AppliedSchurGreenSolver2D(h1D´, solver1D´, s.axis1D, s.wrapped_axes, s.phase_func, s.integrate_opts) +end + +minimal_callsafe_copy(s::SchurGreenSlicer2D, args...) = + SchurGreenSlicer2D(s.ω, minimal_callsafe_copy(s.solver, args...), s.slicer_generator) + +function build_slicer(s::AppliedSchurGreenSolver2D, g, ω, Σblocks, corbitals; params...) + function slicer_generator(s, ϕ_internal) + # updates h1D that is aliased into solver1D.fsolver with the apropriate phases, + # included in params + call!(s.h1D; params..., s.phase_func(ϕ_internal)...) + return SchurGreenSlicer(ω, s.solver1D) + end + g0slicer = SchurGreenSlicer2D(ω, s, slicer_generator) + gslicer = maybe_TMatrixSlicer(g0slicer, Σblocks, corbitals) + return gslicer +end + +function Base.getindex(s::SchurGreenSlicer2D, i::CellOrbitals, j::CellOrbitals) + uas, was = SVector(s.solver.axis1D), s.solver.wrapped_axes + ni, nj = cell(i), cell(j) + i´, j´ = CellOrbitals(ni[uas], orbindices(i)), CellOrbitals(nj[uas], orbindices(j)) + dn = (ni - nj)[was] + function integrand_wrapped(ϕ) + ϕs = sanitize_SVector(ϕ) + s1D = s.slicer_generator(s.solver, ϕs) + gij = s1D[i´, j´] .* cis(-dot(dn, ϕs)) + return gij + end + integral, err = quadgk(integrand_wrapped, -π, π; s.solver.integrate_opts...) + return integral ./ (2π) +end #endregion diff --git a/src/solvers/greensolvers.jl b/src/solvers/greensolvers.jl index 6c8df84d..c753c0ec 100644 --- a/src/solvers/greensolvers.jl +++ b/src/solvers/greensolvers.jl @@ -32,13 +32,15 @@ using Quantica: Quantica, AbstractGreenSolver, I struct SparseLU <:AbstractGreenSolver end -struct Schur{T<:AbstractFloat} <: AbstractGreenSolver +struct Schur{T<:AbstractFloat,O<:NamedTuple} <: AbstractGreenSolver shift::T # Tunable parameter in algorithm, see Ω in scattering.pdf boundary::T # Cell index for boundary (float to allow boundary at Inf) axis::Int # free axis to use for n-dimensional AbstractHamiltonians + integrate_opts::O # quadgk opts for integrals (for hams in more than 1D) end -Schur(; shift = 1.0, boundary = Inf, axis = 1) = Schur(shift, float(boundary), axis) +Schur(; shift = 1.0, boundary = Inf, axis = 1, atol = 1e-7, integrate_opts...) = + Schur(shift, float(boundary), axis, (; atol, integrate_opts...)) struct KPM{B<:Union{Missing,NTuple{2}},A} <: AbstractGreenSolver order::Int diff --git a/src/solvers/selfenergy/schur.jl b/src/solvers/selfenergy/schur.jl index a6ea6147..815c60a0 100644 --- a/src/solvers/selfenergy/schur.jl +++ b/src/solvers/selfenergy/schur.jl @@ -1,5 +1,5 @@ ############################################################################################ -# SelfEnergy(h, glead::GreenFunctionSchurEmptyLead; kw...) +# SelfEnergy(h, glead::GreenFunctionSchurEmptyLead1D; kw...) # Extended self energy solver for deflated ΣL or ΣR Schur factors of lead unitcell # See notes on solvers/green/schur.jl #region @@ -37,7 +37,7 @@ end # semi-infinite lead (possibly by first transforming the lead lattice with `transform`) # and if so, builds the extended Self Energy directly, using the same intercell coupling of # the lead, but using the correct site order of hparent -function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurEmptyLead; reverse = false, transform = missing, kw...) +function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurEmptyLead1D; reverse = false, transform = missing, kw...) lsparent = contact_orbslice(hparent; kw...) schursolver = solver(glead) fsolver = schurfactorsolver(schursolver) @@ -176,7 +176,7 @@ end # through the model coupling. The lead is transformed with `transform` to align it to # hparent. Then we apply the model to the 0D lattice of hparent's selected surface plus the # lead unit cell, and then build an extended self energy -function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurLead, model::AbstractModel; +function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurLead1D, model::AbstractModel; reverse = false, transform = missing, kw...) check_lead_contact_reversal(glead, reverse) schursolver = solver(glead) @@ -262,23 +262,23 @@ function selfenergy_plottables(s::SelfEnergyCouplingSchurSolver, ls::LatticeSlic return (p1, p2, p3) end -check_lead_contact_reversal(::GreenFunctionSchurLead, reverse) = +check_lead_contact_reversal(::GreenFunctionSchurLead1D, reverse) = reverse && @warn "Using `reverse` on a lead with additional contacts. Note that contacts will not be reversed, only the Bravais vectors of the lead" -check_lead_contact_reversal(::GreenFunctionSchurEmptyLead, reverse) = nothing +check_lead_contact_reversal(::GreenFunctionSchurEmptyLead1D, reverse) = nothing #endregion #endregion ############################################################################################ -# SelfEnergy(h, glead::GreenFunctionSchurLead; kw...) +# SelfEnergy(h, glead::GreenFunctionSchurLead1D; kw...) # Regular (Generic) self energy, since Extended is not possible for lead with contacts -# Otherwise equivalent to SelfEnergy(h, glead::GreenFunctionSchurEmptyLead; kw...) +# Otherwise equivalent to SelfEnergy(h, glead::GreenFunctionSchurEmptyLead1D; kw...) # reverse is not allowed since lead contacts cannot be reversed #region -function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurLead; reverse = false, transform = missing, sites...) +function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurLead1D; reverse = false, transform = missing, sites...) blocksizes(blockstructure(glead)) == blocksizes(blockstructure(hparent)) || argerror("The orbital/sublattice structure of parent and lead Hamiltonians do not match") check_lead_contact_reversal(glead, reverse) From beb6edd6bed26d88de69105cb52431091581f6d1 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Mon, 3 Feb 2025 20:29:37 +0100 Subject: [PATCH 03/10] DensityMatrixSchurSolver in 2D test fix --- src/observables.jl | 2 + src/solvers/green/schur.jl | 196 ++++++++++++++++++++++--------------- src/specialmatrices.jl | 17 ++-- src/tools.jl | 1 + test/test_greenfunction.jl | 1 - 5 files changed, 129 insertions(+), 88 deletions(-) diff --git a/src/observables.jl b/src/observables.jl index b34eb24c..5659b3d2 100644 --- a/src/observables.jl +++ b/src/observables.jl @@ -468,6 +468,8 @@ Base.parent(ρ::DensityMatrix) = ρ.gs call!_output(ρ::DensityMatrix) = call!_output(ρ.gs) +solver(ρ::DensityMatrix) = ρ.solver + #endregion #endregion diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index bb9a8fcf..4c05cf07 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -647,10 +647,12 @@ end schur_eigvals(g::GreenFunctionSchurLead1D, ω::Real; params...) = schur_eigvals(g, retarded_omega(ω, solver(g)); params...) -function schur_eigvals(g::GreenFunctionSchurLead1D, ω::Complex; params...) - h = parent(g) # get the (Parametric)Hamiltonian from g +schur_eigvals(g::GreenFunctionSchurLead1D, ω::Complex; params...) = + schur_eigvals((parent(g), g.solver), ω; params...) + +function schur_eigvals((h, solver)::Tuple{AbstractHamiltonian1D,AppliedSchurGreenSolver}, ω::Complex; params...) call!(h; params...) # update the (Parametric)Hamiltonian with the params - sf = g.solver.fsolver # obtain the SchurFactorSolver that computes the AB pencil + sf = solver.fsolver # obtain the SchurFactorSolver that computes the AB pencil update_LR!(sf) # Ensure L and R matrices are updated after updating h update_iG!(sf, ω) # shift inverse G with ω A, B = pencilAB!(sf) # build the pecil @@ -680,74 +682,6 @@ decay_lengths(g::AbstractHamiltonian1D, µ = 0, args...; params...) = #endregion -############################################################################################ -# densitymatrix - dedicated Schur method (integration with Fermi points as segments) -# Does not support non-Nothing contacts -#region - -struct DensityMatrixSchurSolver{T,A,G<:GreenSlice{T},O<:NamedTuple} - gs::G # parent of GreenSlice - orbaxes::A # axes of GreenSlice (orbrows, orbcols) - hmat::Matrix{Complex{T}} # it spans just the unit cell, dense Bloch matrix - psis::Matrix{Complex{T}} # it spans just the unit cell, Eigenstates - opts::O # kwargs for quadgk -end - -## Constructor - -function densitymatrix(s::AppliedSchurGreenSolver, gs::GreenSlice; atol = 1e-7, quadgk_opts...) - isempty(boundaries(s)) || - argerror("Boundaries not implemented for DensityMatrixSchurSolver. Consider using the generic integration solver.") - has_selfenergy(gs) && argerror("The Schur densitymatrix solver currently support only `nothing` contacts") - g = parent(gs) - hmat = similar_Array(hamiltonian(g)) - psis = similar(hmat) - orbaxes = orbrows(gs), orbcols(gs) - # possibly override the default quadgk_opts in solver - opts = (; atol, quadgk_opts...) - solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, opts) - return DensityMatrix(solver, gs) -end - -# API - -## call - -# 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...) - ϕs = @. real(-im*log(λs)) - xs = sort!(ϕs ./ (2π)) - pushfirst!(xs, -0.5) - push!(xs, 0.5) - xs = [mean(view(xs, rng)) for rng in approxruns(xs)] # eliminate approximate duplicates - result = call!_output(s.gs) - integrand!(x) = fermi_h!(result, s, 2π * x, µ, inv(kBT); params...) - fx! = (y, x) -> (y .= integrand!(x)) - quadgk!(fx!, result, xs...; s.opts...) - return result -end - -function fermi_h!(result, s, ϕ, µ, β = 0; params...) - h = hamiltonian(s.gs) - 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!(Hermitian(s.hmat)) - # special-casing β = Inf with views turns out to be slower - fs = (@. ϵs = fermi(ϵs - µ, β)) - fpsis = (s.psis .= psis .* transpose(fs)) - ρcell = EigenProduct(bs, psis, fpsis, ϕ) - getindex!(result, ρcell, s.orbaxes...) - return result -end - -#endregion - - ############################################################################################ # AppliedSchurGreenSolver in 2D #region @@ -755,7 +689,7 @@ end struct AppliedSchurGreenSolver2D{L,H<:AbstractHamiltonian1D,S<:AppliedSchurGreenSolver,F,P<:NamedTuple} <: AppliedGreenSolver h1D::H solver1D::S - axis1D::Int + axis1D::SVector{1,Int} wrapped_axes::SVector{L,Int} phase_func::F # phase_func(ϕ) = (; param_name = ϕ) integrate_opts::P @@ -774,8 +708,8 @@ function apply(s::GS.Schur, h::AbstractHamiltonian, _) L = latdim(h) L == 2 || argerror("GreenSolvers.Schur currently only implemented for 1D and 2D AbstractHamiltonians") - axis1D = s.axis - wrapped_axes = sanitize_SVector(inds_complement(Val(L), (axis1D,))) + axis1D = SVector(s.axis) + wrapped_axes = inds_complement(Val(L), axis1D) h1D = @torus(h, wrapped_axes, ϕ_internal) phase_func(ϕ_internal) = (; ϕ_internal) # we don't pass contacts to solver1D. They will be applied with T-matrix slicer @@ -807,22 +741,126 @@ function build_slicer(s::AppliedSchurGreenSolver2D, g, ω, Σblocks, corbitals; end function Base.getindex(s::SchurGreenSlicer2D, i::CellOrbitals, j::CellOrbitals) - uas, was = SVector(s.solver.axis1D), s.solver.wrapped_axes + uas, was = s.solver.axis1D, s.solver.wrapped_axes ni, nj = cell(i), cell(j) i´, j´ = CellOrbitals(ni[uas], orbindices(i)), CellOrbitals(nj[uas], orbindices(j)) dn = (ni - nj)[was] - function integrand_wrapped(ϕ) - ϕs = sanitize_SVector(ϕ) + function integrand(x) + ϕs = sanitize_SVector(2π * x) s1D = s.slicer_generator(s.solver, ϕs) gij = s1D[i´, j´] .* cis(-dot(dn, ϕs)) return gij end - integral, err = quadgk(integrand_wrapped, -π, π; s.solver.integrate_opts...) - return integral ./ (2π) + integral, err = quadgk(integrand, -0.5, 0.5; s.solver.integrate_opts...) + return integral end +boundaries(s::AppliedSchurGreenSolver2D) = boundaries(s.solver1D) + #endregion #endregion + +############################################################################################ +# densitymatrix - dedicated Schur method (integration with Fermi points as segments) +# 1D and 2D routines. Does not support non-Nothing contacts +#region + +struct DensityMatrixSchurSolver{T,L,A,G<:GreenSlice{T,<:Any,L},O<:NamedTuple} + gs::G # GreenSlice + orbaxes::A # axes of GreenSlice (orbrows, orbcols) + hmat::Matrix{Complex{T}} # it spans just the unit cell, dense Bloch matrix + psis::Matrix{Complex{T}} # it spans just the unit cell, Eigenstates + opts::O # kwargs for quadgk +end + +## Constructor + +function densitymatrix(s::Union{AppliedSchurGreenSolver,AppliedSchurGreenSolver2D}, gs::GreenSlice; atol = 1e-7, quadgk_opts...) + check_no_boundaries_schur(s) + check_no_contacts_schur(gs) + return densitymatrix_schur(gs; atol, quadgk_opts...) +end + +check_no_boundaries_schur(s) = isempty(boundaries(s)) || + argerror("Boundaries not implemented for DensityMatrixSchurSolver. Consider using the generic integration solver.") + +check_no_contacts_schur(gs) = has_selfenergy(gs) && + argerror("The Schur densitymatrix solver currently support only `nothing` contacts") + +function densitymatrix_schur(gs; opts...) + g = parent(gs) + hmat = similar_Array(hamiltonian(g)) + psis = similar(hmat) + orbaxes = orbrows(gs), orbcols(gs) + solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, NamedTuple(opts)) + return DensityMatrix(solver, gs) +end + +# API + +## call + +# 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 +(s::DensityMatrixSchurSolver)(µ, kBT; params...) = integrate_rho_schur(s, µ, kBT; params...) + +function integrate_rho_schur(s::DensityMatrixSchurSolver{<:Any,1}, µ, kBT; params...) + result = call!_output(s.gs) + g = parent(s.gs) + xs = fermi_points_integration_path(g, µ; params...) + integrand!(x) = fermi_h!(result, s, SA[2π * x], µ, inv(kBT); params...) + fx! = (y, x) -> (y .= integrand!(x)) + quadgk!(fx!, result, xs...; s.opts...) + return result +end + +function integrate_rho_schur(s::DensityMatrixSchurSolver{<:Any,2}, µ, kBT; params...) + s2D = solver(s.gs) + result = call!_output(s.gs) + h1D, s1D = (s2D.h1D, s2D.solver1D) + + axisorder = SA[only(s2D.axis1D), only(s2D.wrapped_axes)] + + integrand_inner!(result, (x1, x2)) = + fermi_h!(result, s, 2π * SA[x1, x2][axisorder], µ, inv(kBT); params...) + + function integrand_outer!(result, x2) + xs = fermi_points_integration_path((h1D, s1D), µ; params..., s2D.phase_func(SA[2π*x2])...) + inner! = (y, x1) -> (y .= integrand_inner!(result, (x1, x2))) + quadgk!(inner!, result, xs...; s.opts...) + return result + end + + quadgk!(integrand_outer!, result, -0.5, 0.5; s.opts...) + return result +end + +# this g can be a GreenFunction or a (h1D, s1D) pair +function fermi_points_integration_path(g, µ; params...) + λs = propagating_eigvals(g, µ + 0im, 1e-2; params...) + ϕs = @. real(-im*log(λs)) + xs = sort!(ϕs ./ (2π)) + pushfirst!(xs, -0.5) + push!(xs, 0.5) + xs = [mean(view(xs, rng)) for rng in approxruns(xs)] # eliminate approximate duplicates + return xs +end + +function fermi_h!(result, s, ϕ, µ, β = 0; params...) + h = hamiltonian(s.gs) + 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!(Hermitian(s.hmat)) + # special-casing β = Inf with views turns out to be slower + fs = (@. ϵs = fermi(ϵs - µ, β)) + fpsis = (s.psis .= psis .* transpose(fs)) + ρcell = EigenProduct(bs, psis, fpsis, ϕ) + getindex!(result, ρcell, s.orbaxes...) + return result +end + #endregion top diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index fdbca520..a65366f8 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -659,18 +659,19 @@ Broadcast.broadcast_unalias(dest::AbstractOrbitalArray, src::AbstractOrbitalArra ############################################################################################ ## EigenProduct -# A matrix P = X * U * V' that is stored in terms of matrices U and V without computing +# A matrix P = x * U * V' that is stored in terms of matrices U and V without computing # the product. x is a scalar. Inspired by LowRankMatrices.jl #region # 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} +# phi is a 1D Bloch phase to compute x = cis(-only(dn)*phi) +struct EigenProduct{T,L,MU<:AbstractMatrix{Complex{T}},MV<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} <: AbstractMatrix{T} blockstruct::O U::MU V::MV - phi::T + phi::SVector{L,T} x::Complex{T} - 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} + function EigenProduct{T,L,MU,MV,O}(blockstruct::O, U::MU, V::MV, phi::SVector{L,T}, x::Complex{T}) where {T,L,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 @@ -678,8 +679,8 @@ end #region ## Constructors ## -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)) +EigenProduct(blockstruct::O, U::MU, V::MV, phi::SVector{L} = SA[zero(T)], x = one(Complex{T})) where {T,L,MU<:AbstractMatrix{Complex{T}},MV<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} = + EigenProduct{T,L,MU,MV,O}(blockstruct, U, V, sanitize_SVector(SVector{L,T}, phi), Complex{T}(x)) #endregion @@ -703,7 +704,7 @@ Base.@propagate_inbounds function Base.getindex(P::EigenProduct, i::Integer) end function Base.getindex(P::EigenProduct{T}, ci::AnyCellOrbitals, cj::AnyCellOrbitals) where {T} - phase = isempty(cell(ci)) ? one(T) : cis(only(cell(ci) - cell(cj)) * P.phi) + phase = isempty(cell(ci)) ? one(T) : cis(-dot(cell(ci) - cell(cj), P.phi)) return view(phase * P, orbindices(ci), orbindices(cj)) end @@ -715,7 +716,7 @@ Base.axes(P::EigenProduct) = axes(P.U, 1), axes(P.V, 1) Base.size(P::EigenProduct) = size(P.U, 1), size(P.V, 1) # Base.iterate(a::EigenProduct, i...) = iterate(parent(a), i...) Base.length(P::EigenProduct) = prod(size(P)) -Base.IndexStyle(::Type{T}) where {M,T<:EigenProduct{<:Any,M}} = IndexStyle(M) +Base.IndexStyle(::Type{T}) where {M,T<:EigenProduct{<:Any,<:Any,M}} = IndexStyle(M) Base.similar(P::EigenProduct) = EigenProduct(P.blockstruct, similar(P.U), similar(P.V), P.phi, P.x) Base.similar(P::EigenProduct, t::Type) = EigenProduct(P.blockstruct, similar(P.U, t), similar(P.V, t), P.phi, P.x) Base.copy(P::EigenProduct) = EigenProduct(P.blockstruct, copy(P.U), copy(P.V), P.phi, P.x) diff --git a/src/tools.jl b/src/tools.jl index 729b3c89..92b9c8a9 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -27,6 +27,7 @@ padtuple(t, x, N) = ntuple(i -> i <= length(t) ? t[i] : x, N) end # find indices in range 1:L that are not in inds, assuming inds is sorted and unique +inds_complement(val, inds::SVector) = SVector(inds_complement(val, Tuple(inds))) inds_complement(::Val{L}, inds::NTuple{L´,Integer}) where {L,L´} = _inds_complement(ntuple(zero, Val(L-L´)), 1, 1, inds...) diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 7956718b..69e89ad4 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -139,7 +139,6 @@ end @testset "GreenSolvers applicability" begin h = HP.graphene() - @test_throws ArgumentError greenfunction(h, GS.Schur()) @test_throws ArgumentError greenfunction(h, GS.KPM()) @test_throws ArgumentError greenfunction(h, GS.SparseLU()) @test_throws ArgumentError greenfunction(h, GS.Spectrum()) From 2c00b5aca5714833021465a7c6da16dec62e8a39 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Mon, 3 Feb 2025 21:42:42 +0100 Subject: [PATCH 04/10] test fixes --- src/specialmatrices.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index a65366f8..dc7fb445 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -704,7 +704,8 @@ Base.@propagate_inbounds function Base.getindex(P::EigenProduct, i::Integer) end function Base.getindex(P::EigenProduct{T}, ci::AnyCellOrbitals, cj::AnyCellOrbitals) where {T} - phase = isempty(cell(ci)) ? one(T) : cis(-dot(cell(ci) - cell(cj), P.phi)) + # the sign here is correct + phase = isempty(cell(ci)) ? one(T) : cis(dot(cell(ci) - cell(cj), P.phi)) return view(phase * P, orbindices(ci), orbindices(cj)) end From fe1697d826fa249f8060dda5b3126389d25bcc0b Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Tue, 4 Feb 2025 16:20:59 +0100 Subject: [PATCH 05/10] fermi_points refactor --- src/solvers/green/schur.jl | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 4c05cf07..7c7ee89b 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -839,16 +839,23 @@ function integrate_rho_schur(s::DensityMatrixSchurSolver{<:Any,2}, µ, kBT; para end # this g can be a GreenFunction or a (h1D, s1D) pair -function fermi_points_integration_path(g, µ; params...) +fermi_points_integration_path(g, µ; params...) = + sanitize_integration_path(fermi_points(g, µ; params...)) + +function fermi_points(g, µ; params...) λs = propagating_eigvals(g, µ + 0im, 1e-2; params...) - ϕs = @. real(-im*log(λs)) - xs = sort!(ϕs ./ (2π)) - pushfirst!(xs, -0.5) - push!(xs, 0.5) - xs = [mean(view(xs, rng)) for rng in approxruns(xs)] # eliminate approximate duplicates + xs = @. real(-im*log(λs)) / (2π) return xs end +function sanitize_integration_path(xs, (xmin, xmax) = (-0.5, 0.5)) + sort!(xs) + pushfirst!(xs, xmin) + push!(xs, xmax) + xs´ = [mean(view(xs, rng)) for rng in approxruns(xs)] # eliminate approximate duplicates + return xs´ +end + function fermi_h!(result, s, ϕ, µ, β = 0; params...) h = hamiltonian(s.gs) bs = blockstructure(h) From 4ad4d4e51e0578b9b0dc8822dc7f0c427832c68b Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Tue, 4 Feb 2025 16:29:49 +0100 Subject: [PATCH 06/10] sparse slice support --- src/solvers/green/schur.jl | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 7c7ee89b..0903e384 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -809,32 +809,34 @@ end function integrate_rho_schur(s::DensityMatrixSchurSolver{<:Any,1}, µ, kBT; params...) result = call!_output(s.gs) + data = serialize(result) g = parent(s.gs) xs = fermi_points_integration_path(g, µ; params...) - integrand!(x) = fermi_h!(result, s, SA[2π * x], µ, inv(kBT); params...) + integrand!(x) = fermi_h!(s, SA[2π * x], µ, inv(kBT); params...) fx! = (y, x) -> (y .= integrand!(x)) - quadgk!(fx!, result, xs...; s.opts...) + quadgk!(fx!, data, xs...; s.opts...) return result end function integrate_rho_schur(s::DensityMatrixSchurSolver{<:Any,2}, µ, kBT; params...) s2D = solver(s.gs) result = call!_output(s.gs) + data = serialize(result) h1D, s1D = (s2D.h1D, s2D.solver1D) axisorder = SA[only(s2D.axis1D), only(s2D.wrapped_axes)] - integrand_inner!(result, (x1, x2)) = - fermi_h!(result, s, 2π * SA[x1, x2][axisorder], µ, inv(kBT); params...) + integrand_inner!(x1, x2) = + fermi_h!(s, 2π * SA[x1, x2][axisorder], µ, inv(kBT); params...) - function integrand_outer!(result, x2) + function integrand_outer!(data, x2) xs = fermi_points_integration_path((h1D, s1D), µ; params..., s2D.phase_func(SA[2π*x2])...) - inner! = (y, x1) -> (y .= integrand_inner!(result, (x1, x2))) - quadgk!(inner!, result, xs...; s.opts...) - return result + inner! = (y, x1) -> (y .= integrand_inner!(x1, x2)) + quadgk!(inner!, data, xs...; s.opts...) + return data end - quadgk!(integrand_outer!, result, -0.5, 0.5; s.opts...) + quadgk!(integrand_outer!, data, -0.5, 0.5; s.opts...) return result end @@ -856,7 +858,7 @@ function sanitize_integration_path(xs, (xmin, xmax) = (-0.5, 0.5)) return xs´ end -function fermi_h!(result, s, ϕ, µ, β = 0; params...) +function fermi_h!(s, ϕ, µ, β = 0; params...) h = hamiltonian(s.gs) bs = blockstructure(h) # Similar to spectrum(h, ϕ; params...), but less work (no sort! or sanitization) @@ -866,8 +868,10 @@ function fermi_h!(result, s, ϕ, µ, β = 0; params...) fs = (@. ϵs = fermi(ϵs - µ, β)) fpsis = (s.psis .= psis .* transpose(fs)) ρcell = EigenProduct(bs, psis, fpsis, ϕ) + result = call!_output(s.gs) getindex!(result, ρcell, s.orbaxes...) - return result + data = serialize(result) + return data end #endregion top From a8ec3e7782c7c7d08c5ed6beacfc49002449a56f Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Tue, 4 Feb 2025 17:08:21 +0100 Subject: [PATCH 07/10] tests --- src/docstrings.jl | 4 +++- test/test_greenfunction.jl | 15 ++++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/src/docstrings.jl b/src/docstrings.jl index a6487bbc..97577da6 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -1821,7 +1821,9 @@ Create a selection of site pairs `s::SparseIndices` used to sparsely index into `g::GreenFunction` or `g::GreenSolution`, as `g[s]`. Of the resulting `OrbitalSliceMatrix` only the selected pairs of matrix elements will be computed, leaving the rest as zero (sparse matrix). The sparse matrix spans the minimum number of complete unit cells to -include all site pairs +include all site pairs. + +Tip: if onsite terms are required use `includeonsite = true` as a keyword in `s`. If `kernel = Q` (a matrix instead of `missing`), each of these site blocks `gᵢⱼ` will be replaced by `Tr(kernel * gᵢⱼ)`. diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 69e89ad4..70d4dc6d 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -330,7 +330,7 @@ end g2 = LP.square() |> hopping(1) - onsite(1) |> supercell(3) |> supercell |> attach(nothing, region = RP.circle(2)) |> greenfunction(GS.Spectrum()) 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 + # g3 excluded since KPM doesn't support finite temperatures or generic indexing yet 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]) @@ -338,12 +338,25 @@ end @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-7) @test isapprox(ρ0(0.2, 0.3), ρ(0.2, 0.3)) end + # g2 excluded since it is single-orbital for g in (g1, g3) ρ0 = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])], 5) ρ = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])]) @test isapprox(ρ0(), ρ(); atol = 1e-8) @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-8) end + # 2D Schur solver: path integration is too slow for this solver, not tested here + h = LP.square() |> supercell(1,3) |> hamiltonian(hopping(I), orbitals = 2) |> attach(nothing, region = iszero) + g = h |> greenfunction(GS.Schur()) + for inds in (diagonal(1), diagonal(:), sitepairs(range = 1, includeonsite = true)) + ρ = densitymatrix(g[inds]) + @test Diagonal(ρ(0, 0)) ≈ 0.5*I # half filling kBT == 0 + @test Diagonal(ρ(0, 1)) ≈ 0.5*I # half filling kBT > 0 + end + g´ = h |> greenfunction(GS.Bands(range(0, 2π, 100), range(0, 2π, 100))) + # precision of Bands solver is low, but the results of the two solvers are within 2% + @test maximum(abs, g[](0.2) - g´[](0.2)) < 0.005 + g = HP.graphene(orbitals = 2) |> supercell((1,-1)) |> greenfunction ρ0 = densitymatrix(g[sitepairs(range = 1)], 6) mat = ρ0(0.2, 0.3) From e080924a757e8239892396f95e822689b76e9749 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Tue, 4 Feb 2025 18:02:07 +0100 Subject: [PATCH 08/10] docstrings and test fixes --- docs/src/tutorial/greenfunctions.md | 19 +++++++++--------- src/docstrings.jl | 30 ++++++++++++++++++++--------- src/solvers/green/schur.jl | 4 ++-- test/test_greenfunction.jl | 2 +- 4 files changed, 33 insertions(+), 22 deletions(-) diff --git a/docs/src/tutorial/greenfunctions.md b/docs/src/tutorial/greenfunctions.md index 7de22945..9c6be42c 100644 --- a/docs/src/tutorial/greenfunctions.md +++ b/docs/src/tutorial/greenfunctions.md @@ -63,38 +63,37 @@ The currently implemented `GreenSolvers` (abbreviated as `GS`) are the following - `GS.SparseLU()` - For bounded (`L=0`) AbstractHamiltonians. Default for `L=0`. + For bounded (`L == 0`) AbstractHamiltonians. Default for `L == 0`. Uses a sparse `LU` factorization to compute the inverse of `⟨i|ω - H - Σ(ω)|j⟩`, where `Σ(ω)` is the self-energy from the contacts. - `GS.Spectrum(; spectrum_kw...)` - For bounded (`L=0`) Hamiltonians. This solver does not accept ParametricHamiltonians. Convert to Hamiltonian with `h(; params...)` first. + For bounded (`L == 0`) Hamiltonians. This solver does not accept ParametricHamiltonians. Convert to Hamiltonian with `h(; params...)` first. Uses a diagonalization of `H`, obtained with `spectrum(H; spectrum_kw...)`, to compute the `G⁰ᵢⱼ` using the Lehmann representation `∑ₖ⟨i|ϕₖ⟩⟨ϕₖ|j⟩/(ω - ϵₖ)`. Any eigensolver supported by `spectrum` can be used here. If there are contacts, it dresses `G⁰` using a T-matrix approach, `G = G⁰ + G⁰TG⁰`. - `GS.KPM(order = 100, bandrange = missing, kernel = I)` - For bounded (`L=0`) Hamiltonians, and restricted to sites belonging to contacts (see the section on Contacts). + For bounded (`L == 0`) Hamiltonians, and restricted to sites belonging to contacts (see the section on Contacts). It precomputes the Chebyshev momenta, and incorporates the contact self energy with a T-matrix approach. -- `GS.Schur(boundary = Inf)` +- `GS.Schur(boundary = In, axis = 1, integrate_options...)` - For 1D (`L=1`) AbstractHamiltonians with only nearest-cell coupling. Default for `L=1`. - - Uses a deflating Generalized Schur (QZ) factorization of the generalized eigenvalue problem to compute the unit-cell self energies. - The Dyson equation then yields the Green function between arbitrary unit cells, which is further dressed using a T-matrix approach if the lead has any attached self-energy. + For 1D (`L == 1`) and 2D (`L == 2`) AbstractHamiltonians with only nearest-cell coupling along `axis`. Default for `L=1`. + Uses a deflating Generalized Schur (QZ) factorization of the generalized eigenvalue 1D problem along `axis` to compute the unit-cell self energies. + The Dyson equation then yields the Green function between arbitrary unit cells, which is further dressed using a T-matrix approach if the lead has any attached self-energy. Wavevector along transverse axes in 2D systes are numerically integrated with the QuadGK package using `integrate_options`. - `GS.Bands(bandsargs...; boundary = missing, bandskw...)` - For unbounded (`L>0`) Hamiltonians. + For unbounded (`L > 0`) Hamiltonians. - It precomputes a bandstructure `b = bands(h, bandsargs...; kw..., split = false)` and then uses analytic expressions for the contribution of each subband simplex to the `GreenSolution`. If `boundary = dir => cell_pos`, it takes into account the reflections on an infinite boundary perpendicular to Bravais vector number `dir`, so that all sites with cell index `c[dir] <= cell_pos` are removed. Contacts are incorporated using a T-matrix approach. + It precomputes a bandstructure `b = bands(h, bandsargs...; kw..., split = false)` and then uses analytic expressions for the contribution of each subband simplex to the `GreenSolution`. If `boundary = dir => cell_pos`, it takes into account the reflections on an infinite boundary perpendicular to Bravais vector number `dir`, so that all sites with cell index `c[dir] <= cell_pos` are removed. Contacts are incorporated using a T-matrix approach. Note that this method, while quite general, is approximate, as it uses linear interpolation of bands within each simplex, so it may suffer from precision issues. To retrieve the bands from a `g::GreenFunction` that used the `GS.Bands` solver, we may use `bands(g)`. diff --git a/src/docstrings.jl b/src/docstrings.jl index 97577da6..cb85ae2a 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -1188,13 +1188,14 @@ Currying syntax equivalent to `torus(h, x)`. ## Warning on modifier collisions -For `h::ParametricHamiltonian`s which have modifiers on hoppings, wrapping can be -problematic whenever a modified and an unmodified hopping, or two modified hoppings, get -wrapped into a single hopping in the final `ParametricHamiltonian`. The reason is that the -modifier is applied to the hopping sum. This can be a problem if the modifier is e.g. -position dependent or nonlinear, since then the modified sum may not be equal to the sum of -modified hoppings. Quantica will emit a warning whenever it detects this situation. One -solution in such cases is to use a larger supercell before wrapping. +If two or more hoppings of a `h::ParametricHamiltonian` get stitched into a single one by +`torus`, and any of them depends on parameters, a warning is thrown. The reason is that +these hoppings will be summed, and since the sum is the target of modifiers (because at +least one of the summed hoppings are parameter-dependent), the modifiers will be applied to +the sum, not to the original hoppings before being summed. This is typically not what the +used intends, so the warning should not be ignored. A solution is to use a larger supercell +before calling `torus`. + # Examples @@ -1231,6 +1232,15 @@ Equivalent to `torus(h, phases_or_axes)`, but returning an `n`-dimensional along stitched directions (in addition to the ones specified in `phases_or_axes`, if any). `param_name` can also take any `AbstractArray`, or a `Number` if `n=1`. +## Warning on modifier collisions + +If two or more hoppings get stitched into a single one by `@torus`, and any of them depends +on parameters, a warning is thrown. The reason is that these hoppings will be summed, and +since the sum is the target of modifiers (because at least one of the summed hoppings are +parameter-dependent), the modifiers will be applied to the sum, not to the original hoppings +before being summed. This is typically not what the used intends, so the warning should not +be ignored. A solution is to use a larger supercell before calling `@torus`. + # Examples ```jldoctest @@ -1753,8 +1763,10 @@ possible keyword arguments are - `GS.Spectrum(; spectrum_kw...)` : Diagonalization solver for 0D Hamiltonians using `spectrum(h; spectrum_kw...)` - `spectrum_kw...` : keyword arguments passed on to `spectrum` - This solver does not accept ParametricHamiltonians. Convert to Hamiltonian with `h(; params...)` first. Contact self-energies that depend on parameters are supported. -- `GS.Schur(; boundary = Inf)` : Solver for 1D Hamiltonians based on a deflated, generalized Schur factorization +- `GS.Schur(; boundary = Inf, axis = 1, integrate_opts...)` : Solver for 1D and 2D Hamiltonians based on a deflated, generalized Schur factorization - `boundary` : 1D cell index of a boundary cell, or `Inf` for no boundaries. Equivalent to removing that specific cell from the lattice when computing the Green function. + - If the system is 2D, the wavevector along transverse axis (the one different from the 1D `axis` given in the options) is numerically integrated using QuadGK with options given by `integrate_opts`. + - In 2D systems a warning may be thrown associated to conflicts in torus wrapping which should not be ignored. See `@torus` for details. - `GS.KPM(; order = 100, bandrange = missing, kernel = I)` : Kernel polynomial method solver for 0D Hamiltonians - `order` : order of the expansion in Chebyshev polynomials `Tₙ(h)` of the Hamiltonian `h` (lowest possible order is `n = 0`). - `bandrange` : a `(min_energy, max_energy)::Tuple` interval that encompasses the full band of the Hamiltonian. If `missing`, it is computed automatically, but `using ArnoldiMethod` is required first. @@ -2136,7 +2148,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 and non-empty contacts are not currently supported. Assumes Hermitian Hamiltonian. No `opts`. +- `GS.Schur`: based on numerical integration over Bloch phase (1D and 2D). Boundaries and non-empty contacts are not currently supported. Assumes a Hermitian AbstractHamiltonian. 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`. diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 0903e384..05489765 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -789,12 +789,12 @@ check_no_boundaries_schur(s) = isempty(boundaries(s)) || check_no_contacts_schur(gs) = has_selfenergy(gs) && argerror("The Schur densitymatrix solver currently support only `nothing` contacts") -function densitymatrix_schur(gs; opts...) +function densitymatrix_schur(gs; quadgk_opts...) g = parent(gs) hmat = similar_Array(hamiltonian(g)) psis = similar(hmat) orbaxes = orbrows(gs), orbcols(gs) - solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, NamedTuple(opts)) + solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, NamedTuple(quadgk_opts)) return DensityMatrix(solver, gs) end diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 70d4dc6d..9fb9f563 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -1,5 +1,5 @@ using Quantica: GreenFunction, GreenSlice, GreenSolution, zerocell, CellOrbitals, ncontacts, - solver + solver, Diagonal using ArnoldiMethod # for KPM bandrange From cf62d4f8ec1d0cda07508c35ea9894c505b75383 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Wed, 5 Feb 2025 16:39:51 +0100 Subject: [PATCH 09/10] less minimal_callsafe_copy --- src/greenfunction.jl | 9 ++------- src/solvers/green/bands.jl | 2 -- src/solvers/green/internal.jl | 2 +- src/solvers/green/schur.jl | 17 ----------------- src/solvers/green/sparselu.jl | 8 -------- src/solvers/green/spectrum.jl | 3 --- src/solvers/greensolvers.jl | 3 +-- src/types.jl | 20 ++++++++++++-------- 8 files changed, 16 insertions(+), 48 deletions(-) diff --git a/src/greenfunction.jl b/src/greenfunction.jl index a738386c..9a280c4e 100644 --- a/src/greenfunction.jl +++ b/src/greenfunction.jl @@ -42,7 +42,8 @@ call!_output(c::Contacts) = selfenergyblocks(c) # These kwargs are currently not documented - see specialmatrices.jl #region -(g::GreenFunction)(ω; params...) = minimal_callsafe_copy(call!(g, ω; params...)) +# no need to minimal_callsafe_copy here, because a GreenSolution cannot be call!ed again +(g::GreenFunction)(ω; params...) = call!(g, ω; params...) (g::GreenSlice)(ω; params...) = copy(call!(g, ω; params...)) call!(g::G, ω; params...) where {T,G<:Union{GreenFunction{T},GreenSlice{T}}} = @@ -609,17 +610,11 @@ end ensure_mutable_matrix(m::SMatrix) = Matrix(m) ensure_mutable_matrix(m::AbstractMatrix) = m -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, parentcontacts) = s - #endregion #endregion diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index 94021f8e..ca639cea 100644 --- a/src/solvers/green/bands.jl +++ b/src/solvers/green/bands.jl @@ -877,8 +877,6 @@ 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, parentcontacts) = s # it is read-only - #endregion ############################################################################################ diff --git a/src/solvers/green/internal.jl b/src/solvers/green/internal.jl index 8bc1f12b..e01edca2 100644 --- a/src/solvers/green/internal.jl +++ b/src/solvers/green/internal.jl @@ -37,7 +37,7 @@ end needs_omega_shift(s::AppliedModelGreenSolver) = false -minimal_callsafe_copy(s::Union{ModelGreenSlicer,AppliedModelGreenSolver}, args...) = s +minimal_callsafe_copy(s::AppliedModelGreenSolver, args...) = s Base.getindex(s::ModelGreenSlicer, is::CellOrbitals, js::CellOrbitals) = [s.fω(i, j) for i in cellorbs(is), j in cellorbs(js)] diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 05489765..91e16f61 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -621,20 +621,6 @@ end maybe_SMatrix(G::Matrix, rows::SVector{L}, cols::SVector{L´}) where {L,L´} = SMatrix{L,L´}(G) maybe_SMatrix(G, rows, cols) = G -# TODO: Perhaps too conservative -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)) - isdefined(s, :G₋₁₋₁R) && (s´.G₋₁₋₁R = copy(s.G₋₁₋₁R)) - isdefined(s, :R´G₁₁L) && (s´.R´G₁₁L = copy(s.R´G₁₁L)) - isdefined(s, :L´G₋₁₋₁R) && (s´.L´G₋₁₋₁R = copy(s.L´G₋₁₋₁R)) - return s´ -end #endregion #endregion @@ -725,9 +711,6 @@ function minimal_callsafe_copy(s::AppliedSchurGreenSolver2D, args...) return AppliedSchurGreenSolver2D(h1D´, solver1D´, s.axis1D, s.wrapped_axes, s.phase_func, s.integrate_opts) end -minimal_callsafe_copy(s::SchurGreenSlicer2D, args...) = - SchurGreenSlicer2D(s.ω, minimal_callsafe_copy(s.solver, args...), s.slicer_generator) - function build_slicer(s::AppliedSchurGreenSolver2D, g, ω, Σblocks, corbitals; params...) function slicer_generator(s, ϕ_internal) # updates h1D that is aliased into solver1D.fsolver with the apropriate phases, diff --git a/src/solvers/green/sparselu.jl b/src/solvers/green/sparselu.jl index 2e57159c..9840042a 100644 --- a/src/solvers/green/sparselu.jl +++ b/src/solvers/green/sparselu.jl @@ -146,14 +146,6 @@ similar_source64(s::SparseLUGreenSlicer, j::CellOrbitals) = # getindex must return a Matrix Base.getindex(s::SparseLUGreenSlicer, i::CellOrbitals, j::CellOrbitals) = copy(view(s, i, j)) -# 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, 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 diff --git a/src/solvers/green/spectrum.jl b/src/solvers/green/spectrum.jl index 20add2ff..283d2abd 100644 --- a/src/solvers/green/spectrum.jl +++ b/src/solvers/green/spectrum.jl @@ -24,9 +24,6 @@ spectrum(s::AppliedSpectrumGreenSolver) = s.spectrum minimal_callsafe_copy(s::AppliedSpectrumGreenSolver, parentham, parentcontacts) = AppliedSpectrumGreenSolver(s.spectrum) -minimal_callsafe_copy(s::SpectrumGreenSlicer, parentham, parentcontacts) = - SpectrumGreenSlicer(s.ω, s.solver) - #endregion #region ## apply ## diff --git a/src/solvers/greensolvers.jl b/src/solvers/greensolvers.jl index c753c0ec..ca57dea5 100644 --- a/src/solvers/greensolvers.jl +++ b/src/solvers/greensolvers.jl @@ -11,9 +11,8 @@ # - view(gs, ::Int, ::Int) -> g(ω; kw...) between specific contacts (has error fallback) # - view(gs, ::Colon, ::Colon) -> g(ω; kw...) between all contacts (has error fallback) # - Both of the above are of type `SubArray` -# It must also implement generic slicing, and minimal copying +# It must also implement generic slicing # - gs[i::CellOrbitals, j::CellOrbitals] -> must return a `Matrix` for type stability -# - minimal_callsafe_copy(gs, parentham, parentcontacts) # The user-facing indexing API accepts: # - i::Integer -> Sites of Contact number i # - sites(cell::Tuple, sind::Int)::Subcell -> Single site in a cell diff --git a/src/types.jl b/src/types.jl index 043e9ad0..44f08609 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1171,6 +1171,8 @@ end Base.size(s::SparseMatrixView, i...) = size(s.mat, i...) +# Definition of minimal_callsafe_copy +# if x´ = minimal_callsafe_copy(x), then a call!(x, ...; ...) will not affect x´ in any way minimal_callsafe_copy(s::SparseMatrixView) = SparseMatrixView(view(copy(parent(s.matview)), s.matview.indices...), copy(s.mat), s.ptrs) @@ -1681,6 +1683,8 @@ blocktype(h::ParametricHamiltonian) = blocktype(parent(h)) lattice(h::ParametricHamiltonian) = lattice(hamiltonian(h)) +# Definition of minimal_callsafe_copy +# if x´ = minimal_callsafe_copy(x), then a call!(x, ...; ...) will not affect x´ in any way function minimal_callsafe_copy(p::ParametricHamiltonian) h´ = minimal_callsafe_copy(p.h) modifiers´ = maybe_relink_serializer.(p.modifiers, Ref(h´)) @@ -2420,14 +2424,14 @@ function minimal_callsafe_copy(g::GreenFunction) return GreenFunction(parent´, solver´, contacts´) end -function minimal_callsafe_copy(g::GreenSolution) - parentg´ = minimal_callsafe_copy(g.parent) - 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 +# function minimal_callsafe_copy(g::GreenSolution) +# parentg´ = minimal_callsafe_copy(g.parent) +# 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 minimal_callsafe_copy(g::GreenSlice) = GreenSlice(minimal_callsafe_copy(g.parent), g.rows, g.cols, copy(g.output)) From 61c8c34a5a218809f55dcda7b9c6778573982a4d Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Wed, 5 Feb 2025 16:40:48 +0100 Subject: [PATCH 10/10] remove callsafe copy of GreenSolution --- src/types.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/types.jl b/src/types.jl index 44f08609..63142490 100644 --- a/src/types.jl +++ b/src/types.jl @@ -2424,15 +2424,6 @@ function minimal_callsafe_copy(g::GreenFunction) return GreenFunction(parent´, solver´, contacts´) end -# function minimal_callsafe_copy(g::GreenSolution) -# parentg´ = minimal_callsafe_copy(g.parent) -# 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 - minimal_callsafe_copy(g::GreenSlice) = GreenSlice(minimal_callsafe_copy(g.parent), g.rows, g.cols, copy(g.output))