Merge pull request #332 from pablosanjose/schur

`GreenSolvers.Schur` for 2-dimensional systems

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)`.
 

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")

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,11 @@ 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, callback = Returns(nothing), 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.
+    - `callback` : a function `f(ϕ, z)` for 1D systems or `f(ϕ1, ϕ2, z)` for 2D systems that gets called at each Brillouin zone integration point `ϕ` or `(ϕ1, ϕ2)`, and where `z` is the integrand (an array).
+    - 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`, which is `(; atol = 1e-7, order = 5)` by default.
+    - 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.
@@ -1821,7 +1834,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ᵢⱼ)`.
@@ -2134,7 +2149,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`.
 
@@ -2396,7 +2411,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 +2685,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

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
 

src/observables.jl

@@ -468,6 +468,8 @@ Base.parent(ρ::DensityMatrix) = ρ.gs
 
 call!_output(ρ::DensityMatrix) = call!_output(ρ.gs)
 
+solver(ρ::DensityMatrix) = ρ.solver
+
 #endregion
 
 #endregion

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

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)

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)