More general Quantica.gaps

docstring

docstring

src/docstrings.jl

@@ -2538,13 +2538,14 @@ true
 deserialize!
 
 """
-    Quantica.gaps(h::Hamiltonian1D{T}, µ = 0; atol = eps(T))
+    Quantica.gaps(h::Hamiltonian1D{T}, µ = 0; atol = eps(T), kw...)
 
 Compute the energy gaps of a 1D Hamiltonian `h` at chemical potential `µ`. The result is a
 `Vector{T}` of the local minima of the `|ϵ(ϕ) - µ|`, where `ϵ(ϕ)` is the energy band closest
 to `µ` and `ϕ ∈ [-π,π]` is the Bloch phase. The `atol` parameter is the absolute tolerance
 used to determine the local minima versus `ϕ`, which are computed using the `Schur` solver
-for 1D Hamiltonians.
+for 1D Hamiltonians. The keywords `kw` are passed to the ArnoldiMethod `partialschur!`
+eigensolver (`kw = (; nev = 1)` by default).
 
 The `LinearMaps` and `ArnoldiMethod` packages must be loaded to enable this functionality.
 
@@ -2562,11 +2563,22 @@ julia> Quantica.gaps(h(U=2))
 ```
 
 # See also:
-    `Quantica.decay_lengths`
+    `Quantica.gap`, `Quantica.decay_lengths`
 
 """
 gaps
 
+"""
+    Quantica.gap(h::Hamiltonian1D{T}, µ = 0; atol = eps(T), kw...)
+
+Compute the minimal gap around `µ`, see `Quantica.gaps`
+
+# See also:
+    `Quantica.gaps`, `Quantica.decay_lengths`
+
+"""
+gap
+
 """
     Quantica.decay_lengths(g::GreenFunctionSchurLead, µ = 0; reverse = false)
     Quantica.decay_lengths(h::AbstractHamiltonian1D, µ = 0; reverse = false)

src/observables.jl

@@ -760,7 +760,7 @@ end
 # gap(::Hamiltonian, µ = 0, ...) -> minimum gap in the bands around µ
 #region
 
-function gaps(h::Hamiltonian{T,<:Any,1}, µ = 0, ϕstore = missing; atol = eps(T)) where {T}
+function gaps(h::Hamiltonian{T,<:Any,1}, µ = 0, ϕstore = missing; atol = eps(T), nev = 1, kw...) where {T}
     g = greenfunction(h, GS.Schur())
     λs = schur_eigvals(g, µ)
     cϕs = λs .= -im .* log.(λs) # saves one allocation
@@ -771,9 +771,11 @@ function gaps(h::Hamiltonian{T,<:Any,1}, µ = 0, ϕstore = missing; atol = eps(T
     rϕs = real.(cϕs)
     iϕs = chopsmall.(abs.(imag.(cϕs)), atol)
     ϕstore === missing || copy!(ϕstore, cϕs)
-    solver = ES.ShiftInvert(ES.ArnoldiMethod(nev = 1), µ)
+    solver = ES.ShiftInvert(ES.ArnoldiMethod(; nev, kw...), µ)
     n = flatsize(h)
-    Δs = [iszero(iϕ) || rank(h(rϕ)-µ*I; tol=atol) < n ? zero(T) : abs(first(energies(spectrum(h, rϕ; solver)))-µ) for (rϕ, iϕ) in zip(rϕs, iϕs)]
+    Δs = [iszero(iϕ) || rank(h(rϕ)-µ*I; tol=atol) < n ?
+        zero(T) :
+        abs(minimum(x->abs(x-µ), energies(spectrum(h, rϕ; solver)))-µ) for (rϕ, iϕ) in zip(rϕs, iϕs)]
     return Δs
 end