Inverse-free Green's functions

[pablosanjose]

This builds on top of #161

• We replace the QR deflation algorithm with a similar deflation algorithm developed by M. Wimmer, A. Akhmerov et al. that has a similar performance but is more transparent IMO. It also provides the physical retarded/advanced eigenmodes or their invariant subspaces, but does so through an LU inversion instead of a QR factorization. The numerical precision is slightly better, squeezing out one or two more significant digits of the self energy in my tests.
• We introduce a new way to build the Green's function that does not require the computation of the self-energy explicitly, which requires an inverse. This has profound consequences for the case of defective eigenmode matrices, i.e. when there is a bound state or a singularity of any other type. The idea is to replace the self energy by an auxiliary subspace, whose "hamiltonian" and coupling to the unit cell h0 is expressed in terms of the eigenmode invariant subspaces, no inverse needed.
• This "inverse-free" approach is also interesting for future problems involving hybrid systems, such as multiterminal scattering systems and other combination of coupled Hamiltonians, since any reservoir, regardless of its dimensionality, could be coupled to another system through the same mechanism.

The user-facing consequences of this PR are

• No undeflated method exists anymore, so Schur1D() takes no keyword arguments.
• Numerical precision and speed should be improved
• In Deflating Schur algorithm for quasi1D Green functions #161, doing g(ω, 1=>n) used a fastpath that was twice faster than a general g(ω, m=>n). Not anymore. Through careful reuse of Schur factorizations, now everything is faster than the original g(ω, 1=>n), and there is no penalty for the general case.
• A new breaking syntax is introduced: g0 = g(ω) with g::GreensFunction{Schur1DGreensSolver} generates a Schur1DGreensSolution object, instead of defaulting to g(ω, 1=>1). It is the equivalent of Deflating Schur algorithm for quasi1D Green functions #161's g(ω, missing). Then g0[n=>m] computes the actual Greens function matrix. Fastpaths still exist that save a little bit of time when only a specific n,m is needed. The syntax for the fastpath computation is just as before, g(ω, n=>m), with the guarantee that g(ω, n=>m) ≈ g(ω)[n=>m]
• Hamiltonians with intercell couplings beyond nearest neighbors cannot be used to build Green's functions. An error is thrown instructing to enlarge the unit cell with unitcell. This is done to avoid nasty bugs (unitcell is not a block-wise cell expansion, but shuffles sites/sublattices)

Example

julia> using LinearAlgebra, VegaLite, BenchmarkTools, Quantica

julia> h = LatticePresets.square() |>
           hamiltonian(hopping(1)) |>
           unitcell((20, 0), region = r -> -10 <= r[2] <= 10 && norm(r-SA[10,0]) > 5)
Hamiltonian{<:Lattice} : Hamiltonian on a 1D Lattice in 2D space
  Bloch harmonics  : 3 (SparseMatrixCSC, sparse)
  Harmonic size    : 339 × 339
  Orbitals         : ((:a,),)
  Element type     : scalar (ComplexF64)
  Onsites          : 0
  Hoppings         : 1272
  Coordination     : 3.752212389380531

julia> g = greens(h, Schur1D(), boundaries = (0,))
GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method
  Flat matrix size      : 339 × 339
  Flat deflated size    : 21 × 21
  Original element type : scalar (ComplexF64)
  Boundaries            : (0,)

julia> g0 = @btime g(0.2)
  22.541 ms (76 allocations: 9.93 MiB)
Schur1DGreensSolution : Schur1D solution of a GreensFunction at fixed ω
  ω             : 0.2 + 1.4901161193847656e-8im
  Matrix size   : 339 × 339
  Deflated size : 21 × 21
  Element type  : ComplexF64
  Boundary      : 0

julia> ldos = @time [(n,) => -imag.(diag(g0[n=>n])) for n in 1:5];
  0.073697 seconds (66.91 k allocations: 12.890 MiB, 64.73% compilation time)

julia> vlplot(h, ldos, sitesize = DensityShader(), maxthickness = 0.02)

[codecov-commenter]

Codecov Report

Merging #162 (c3f5d63) into master (975a383) will increase coverage by 0.07%.
The diff coverage is 84.01%.

@@            Coverage Diff             @@
##           master     #162      +/-   ##
==========================================
+ Coverage   71.00%   71.08%   +0.07%     
==========================================
  Files          16       17       +1     
  Lines        3501     3493       -8     
==========================================
- Hits         2486     2483       -3     
+ Misses       1015     1010       -5     

Impacted Files	Coverage Δ
src/Quantica.jl	100.00% <ø> (ø)
src/hamiltonian.jl	86.64% <42.85%> (-0.21%)	⬇️
src/greens.jl	80.43% <82.46%> (+3.03%)	⬆️
src/tools.jl	57.33% <94.73%> (-1.60%)	⬇️
src/effective.jl	96.87% <96.87%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 975a383...c3f5d63. Read the comment docs.

[pablosanjose]

Note: this PR does not use the velocity operator to separate advanced from retarded modes. Instead it relies on the traditional approach of adding a small imaginary part to ω. The reason is that the added Green function methods are not only useful for transport observables (which do not depend on a small imag(ω) → 0), but also for spectral observables, for which it is desirable to have a broadening of bound states, given by imag(ω), as they would otherwise be invisible. This approach also works in situations where the velocity operator method is problematic, e.g. when there are flat bands or close to band edges (zero velocity). Note that the user doesn't need to specify the imaginary part, it is automatic when ω is real, unless overridden by using a complex ω. The last commit also provides a useful error when the provided imaginary part is insufficient to resolve advanced and retarded modes.