Merge pull request #238 from pablosanjose/typestable_g

Add 32bit support for GreenFunctions

.github/workflows/ci.yml

@@ -5,6 +5,9 @@ on:
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: true
+permissions:
+  actions: write
+  contents: read
 jobs:
   test:
     if: github.event_name == 'push' && contains(toJson(github.event.commits), '[skip test]') == false && contains(toJson(github.event.commits), '[skip tests]') == false
@@ -14,23 +17,24 @@ jobs:
       matrix:
         version:
           - '1'
-          # - 'nightly'
+          - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
         arch:
           - x64
-          - x86
-        exclude:
-          - os: macOS-latest
-            arch: x86
+          # - x86   ## Currently broken on nightly due to OutOfMemoryError()
+        # exclude:
+        #   - os: macOS-latest
+        #     arch: x86
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
       - uses: julia-actions/julia-uploadcodecov@latest

src/bands.jl

@@ -61,8 +61,9 @@ end
 
 bands(rng, rngs...; kw...) = h -> bands(h, rng, rngs...; kw...)
 
-bands(h::Union{Function,AbstractHamiltonian}, rng, rngs...; kw...) =
-    bands(h, mesh(rng, rngs...); kw...)
+bands(h::Function, rng, rngs...; kw...) = bands(h, mesh(rng, rngs...); kw...)
+bands(h::AbstractHamiltonian{T}, rng, rngs::Vararg{Any,L´}; kw...) where {T,L´} =
+    bands(h, convert(Mesh{SVector{L´+1,T}}, mesh(rng, rngs...)); kw...)
 
 bands(h::AbstractHamiltonian{<:Any,<:Any,L}; kw...) where {L} =
     bands(h, default_band_ticks(Val(L))...; kw...)
@@ -660,11 +661,14 @@ end
 function simplex_projector(hkav, verts, vind, εav, mindeg)
     φ = states(verts[vind])
     hproj = φ' * hkav * φ
-    _, P = eigen!(Hermitian(hproj), sortby = ε -> abs(ε - εav))
+    _, P = maybe_eigen!(Hermitian(hproj), sortby = ε -> abs(ε - εav))
     Pthin = view(P, :, 1:mindeg)
     return φ * Pthin
 end
 
+maybe_eigen!(A::AbstractMatrix{<:LinearAlgebra.BlasComplex}; kw...) = eigen!(A; kw...)
+maybe_eigen!(A; kw...) = eigen(A; kw...)    # generic fallback for e.g. Complex16
+
 #endregion
 
 ############################################################################################

src/convert.jl

@@ -15,6 +15,9 @@ Base.convert(::Type{T}, l::CellSites) where T<:CellSites = T(l)
 Base.convert(::Type{T}, l::T) where T<:AbstractHamiltonian = l
 Base.convert(::Type{T}, l::AbstractHamiltonian) where T<:AbstractHamiltonian = T(l)
 
+Base.convert(::Type{T}, l::T) where T<:Mesh = l
+Base.convert(::Type{T}, l::Mesh) where T<:Mesh = T(l)
+
 # Constructors for conversion
 Sublat{T,E}(s::Sublat, name = s.name) where {T<:AbstractFloat,E} =
     Sublat{T,E}([sanitize_SVector(SVector{E,T}, site) for site in sites(s)], name)
@@ -40,6 +43,9 @@ function ParametricHamiltonian{E}(ph::ParametricHamiltonian) where {E}
     return ParametricHamiltonian(hparent, h, ms, ptrs, pars)
 end
 
+Mesh{S}(m::Mesh) where {S} =
+    Mesh(convert(Vector{S}, vertices(m)), neighbors(m), simplices(m))
+
 # We need this to promote different sublats into common dimensionality and type to combine
 # into a lattice
 Base.promote(ss::Sublat{T,E}...) where {T,E} = ss

src/greenfunction.jl

@@ -31,23 +31,30 @@ default_green_solver(::AbstractHamiltonian) = GS.Bands()
 (g::GreenSlice)(; params...) = minimal_callsafe_copy(call!(g; params...))
 (g::GreenSlice)(ω; params...) = copy(call!(g, ω; params...))
 
-call!(g::GreenFunction, ω::Real; params...) = call!(g, retarded_omega(ω, solver(g)); params...)
+function call!(g::GreenFunction{T}, ω::Real; params...) where {T}
+    ω´ = retarded_omega(real_or_complex_typed(T, ω), solver(g))
+    return call!(g, ω´; params...)
+end
 
-function call!(g::GreenFunction, ω::Complex; params...)
+function call!(g::GreenFunction{T}, ω::Complex; params...) where {T}
+    ω´ = real_or_complex_typed(T, ω)
     h = parent(g)
     contacts´ = contacts(g)
     call!(h; params...)
-    Σblocks = call!(contacts´, ω; params...)
+    Σblocks = call!(contacts´, ω´; params...)
     corbs = contactorbitals(contacts´)
-    slicer = solver(g)(ω, Σblocks, corbs)
+    slicer = solver(g)(ω´, Σblocks, corbs)
     return GreenSolution(g, slicer, Σblocks, corbs)
 end
 
 call!(g::GreenSlice; params...) =
     GreenSlice(call!(greenfunction(g); params...), slicerows(g), slicecols(g))
 
-call!(g::GreenSlice, ω; params...) =
-    call!(greenfunction(g), ω; params...)[slicerows(g), slicecols(g)]
+call!(g::GreenSlice{T}, ω; params...) where {T} =
+    call!(greenfunction(g), real_or_complex_typed(T, ω); params...)[slicerows(g), slicecols(g)]
+
+real_or_complex_typed(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω)
+real_or_complex_typed(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω)
 
 retarded_omega(ω::T, s::AppliedGreenSolver) where {T<:Real} =
     ω + im * sqrt(eps(float(T))) * needs_omega_shift(s)
@@ -434,7 +441,7 @@ Base.view(s::TMatrixSlicer, ::Colon, ::Colon) = view(s.gcontacts, :, :)
 
 function Base.getindex(s::TMatrixSlicer, i::CellOrbitals, j::CellOrbitals)
     g0 = s.g0slicer
-    g0ij = g0[i, j]
+    g0ij = ensure_mutable_matrix(g0[i, j])
     tkk´ = s.tmatrix
     isempty(tkk´) && return g0ij
     k = s.contactorbs
@@ -444,6 +451,9 @@ function Base.getindex(s::TMatrixSlicer, i::CellOrbitals, j::CellOrbitals)
     return gij
 end
 
+ensure_mutable_matrix(m::SMatrix) = Matrix(m)
+ensure_mutable_matrix(m::AbstractMatrix) = m
+
 minimal_callsafe_copy(s::TMatrixSlicer) = TMatrixSlicer(minimal_callsafe_copy(s.g0slicer),
     s.tmatrix, s.gcontacts, s.contactorbs)
 

src/presets/hamiltonians.jl

@@ -6,16 +6,16 @@ module HamiltonianPresets
 
 using Quantica, LinearAlgebra
 
-function graphene(; a0 = 0.246, range = a0/sqrt(3), t0 = 2.7, β = 3, dim = 2, type = Float64, names = (:A, :B), kw...)
+function graphene(; a0 = 0.246, range = neighbors(1), t0 = 2.7, β = 3, dim = 2, type = Float64, names = (:A, :B), kw...)
     lat = LatticePresets.honeycomb(; a0, dim, type, names)
     h = hamiltonian(lat,
-        hopping((r, dr) -> t0 * exp(-β*(sqrt(3) * norm(dr)/a0 - 1)) * I,range = range); kw...)
+        hopping((r, dr) -> t0 * exp(-β*(sqrt(3) * norm(dr)/a0 - 1)) * I, range = range); kw...)
     return h
 end
 
 function twisted_bilayer_graphene(;
     twistindex = 1, twistindices = (twistindex, 1), a0 = 0.246,
-    interlayerdistance = 1.36a0, rangeintralayer = a0/sqrt(3), rangeinterlayer = 4a0/sqrt(3),
+    interlayerdistance = 1.36a0, rangeintralayer = neighbors(1), rangeinterlayer = 4a0/sqrt(3),
     hopintra = 2.70 * I, hopinter = 0.48, modelintra = hopping(hopintra, range = rangeintralayer),
     type = Float64, names = (:Ab, :Bb, :At, :Bt),
     kw...)
@@ -59,4 +59,4 @@ end # module
 
 const HP = HamiltonianPresets
 
-#endregion
+#endregion

src/sanitizers.jl

@@ -82,17 +82,16 @@ end
 
 #endregion
 
-# ############################################################################################
-# # Block sanitizers
-# #region
+############################################################################################
+# CellIndices sanitizers
+#region
 
-# sanitize_block(S::Type{<:Number}, s, _) = convert(S, first(s))
-# sanitize_block(S::Type{<:SMatrix}, s::SMatrix, size) = sanitize_SMatrix(S, s, size)
-# sanitize_block(::Type{S}, s::Number, size) where {S<:SMatrix} = sanitize_SMatrix(S, S(s*I), size)
-# sanitize_block(::Type{S}, s::UniformScaling, size) where {S<:SMatrix} =
-#     sanitize_SMatrix(S, S(s), size)
+# an inds::Tuple fails some tests because convert(Tuple, ::UnitRange) doesnt exist, but
+# convert(SVector, ::UnitRange) does. Used e.g. in compute_or_retrieve_green @ sparselu.jl
+sanitize_cellindices(inds::Tuple) = SVector(inds)
+sanitize_cellindices(inds) = inds
 
-# #endregion
+#endregion
 
 ############################################################################################
 # Supercell sanitizers
@@ -115,11 +114,11 @@ sanitize_supercell(::Val{L}, v) where {L} =
 # Eigen sanitizers
 #region
 
-sanitize_eigen(ε, Ψ) = Eigen(sorteigs!(ε, Ψ)...)
 sanitize_eigen(ε::AbstractVector, Ψs::AbstractVector{<:AbstractVector}) =
     sanitize_eigen(ε, hcat(Ψs...))
-sanitize_eigen(ε::AbstractVector{<:Real}, Ψ::AbstractMatrix) =
-    sanitize_eigen(complex.(ε), Ψ)
+sanitize_eigen(ε, Ψ) = Eigen(sorteigs!(sanitize_eigen(ε), sanitize_eigen(Ψ))...)
+sanitize_eigen(x::AbstractArray{<:Real}) = complex.(x)
+sanitize_eigen(x::AbstractArray{<:Complex}) = x
 
 function sorteigs!(ϵ::AbstractVector, ψ::AbstractMatrix)
     p = Vector{Int}(undef, length(ϵ))

src/slices.jl

@@ -157,7 +157,7 @@ combine_subcells(c::C, cs::C...) where {C<:CellOrbitals} =
 function combine_subcells(c::C, cs::C...) where {C<:CellOrbitalsGrouped}
     groups´ = merge(orbgroups(c), orbgroups.(cs)...)
     indices´ = union(orbindices(c), orbindices.(cs)...)
-    return CellIndices(cell(c), indices´, OrbitalLikeGrouped(groups´))
+    return CellOrbitalsGrouped(cell(c), indices´, groups´)
 end
 
 #endregion
@@ -299,7 +299,7 @@ sites_to_orbs(c::AnyCellOrbitalsDict, _) = c
 sites_to_orbs(c::AnyCellOrbitals, _) = c
 
 ## convert SiteSlice -> OrbitalSliceGrouped/OrbitalSlice
-Contacts
+
 sites_to_orbs(s::SiteSelector, g) = sites_to_orbs(lattice(g)[s], g)
 sites_to_orbs(kw::NamedTuple, g) = sites_to_orbs(getindex(lattice(g); kw...), g)
 sites_to_orbs(i::Integer, g) = orbslice(selfenergies(contacts(g), i))
@@ -335,13 +335,13 @@ function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure)
     sites = siteindices(cs)
     groups = _groups(sites, os) # sites, orbranges
     orbinds = _orbinds(sites, groups, os)
-    return CellIndices(cell(cs), orbinds, OrbitalLikeGrouped(Dictionary(groups...)))
+    return CellOrbitalsGrouped(cell(cs), orbinds, Dictionary(groups...))
 end
 
 function sites_to_orbs_flat(cs::CellSites, os::OrbitalBlockStructure)
     sites = siteindices(cs)
     orbinds = _orbinds(sites, os)
-    return CellIndices(cell(cs), orbinds, OrbitalLike())
+    return CellOrbitals(cell(cs), orbinds)
 end
 
 _groups(i::Integer, os) = [i], [flatrange(os, i)]