general Hartree-Fock mean fields

docs/make.jl

@@ -26,6 +26,8 @@ makedocs(;
             "Observables" => "tutorial/observables.md"
             ],
         "Advanced" => [
+            "Non-spatial models" => "advanced/nonspatial.md",
+            "Serializers" => "advanced/serializers.md",
             "Self-consistent mean fields" => "advanced/meanfield.md",
             "Wannier90 imports" => "advanced/wannier90.md"
             ],

docs/src/advanced/nonspatial.md

@@ -0,0 +1,71 @@
+# Non-spatial models
+
+As briefly mentioned when discussing parametric models and modifiers, we have a special syntax that allows models to depend on sites directly, instead of on their spatial location. We call these non-spatial models. A simple example, with an onsite energy proportional to the site **index**
+```julia
+julia> model = @onsite((i; p = 1) --> ind(i) * p)
+ParametricModel: model with 1 term
+  ParametricOnsiteTerm{ParametricFunction{1}}
+    Region            : any
+    Sublattices       : any
+    Cells             : any
+    Coefficient       : 1
+    Argument type     : non-spatial
+    Parameters        : [:p]
+```
+or a modifier that changes a hopping between different cells
+```julia
+julia> modifier = @hopping!((t, i, j; dir = 1) --> (cell(i) - cell(j))[dir])
+HoppingModifier{ParametricFunction{3}}:
+  Region            : any
+  Sublattice pairs  : any
+  Cell distances    : any
+  Hopping range     : Inf
+  Reverse hops      : false
+  Argument type     : non-spatial
+  Parameters        : [:dir]
+```
+
+Note that we use the special syntax `-->` instead of `->`. This indicates that the positional arguments of the function, here called `i` and `j`, are no longer site positions as up to now. These `i, j` are non-spatial arguments, as noted by the `Argument type` property shown above. Instead of a position, a non-spatial argument `i` represent an individual site, whose index is `ind(i)`, its position is `pos(i)` and the cell it occupies on the lattice is `cell(i)`.
+
+Technically `i` is of type `CellSitePos`, which is an internal type not meant for the end user to instantiate. One special property of this type, however, is that it can efficiently index `OrbitalSliceArray`s. We can use this to build a Hamiltonian that depends on an observable, such as a `densitymatrix`. A simple example of a four-site chain with onsite energies shifted by a potential proportional to the local density on each site:
+```julia
+julia> h = LP.linear() |> onsite(2) - hopping(1) |> supercell(4) |> supercell;
+
+julia> h(SA[])
+4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+  2.0+0.0im  -1.0+0.0im       ⋅           ⋅
+ -1.0+0.0im   2.0+0.0im  -1.0+0.0im       ⋅
+      ⋅      -1.0+0.0im   2.0+0.0im  -1.0+0.0im
+      ⋅           ⋅      -1.0+0.0im   2.0+0.0im
+
+julia> g = greenfunction(h, GS.Spectrum());
+
+julia> ρ0 = densitymatrix(g[])(0.5, 0) ## density matrix at chemical potential `µ=0.5` and temperature `kBT = 0`  on all sites
+4×4 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}:
+ 0.138197+0.0im  0.223607+0.0im  0.223607+0.0im  0.138197+0.0im
+ 0.223607+0.0im  0.361803+0.0im  0.361803+0.0im  0.223607+0.0im
+ 0.223607+0.0im  0.361803+0.0im  0.361803+0.0im  0.223607+0.0im
+ 0.138197+0.0im  0.223607+0.0im  0.223607+0.0im  0.138197+0.0im
+
+julia> hρ = h |> @onsite!((o, i; U = 0, ρ = ρ0) --> o + U * ρ[i])
+ParametricHamiltonian{Float64,1,0}: Parametric Hamiltonian on a 0D Lattice in 1D space
+  Bloch harmonics  : 1
+  Harmonic size    : 4 × 4
+  Orbitals         : [1]
+  Element type     : scalar (ComplexF64)
+  Onsites          : 4
+  Hoppings         : 6
+  Coordination     : 1.5
+  Parameters       : [:U, :ρ]
+
+julia> hρ(SA[]; U = 1, ρ = ρ0)
+4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+ 2.1382+0.0im    -1.0+0.0im         ⋅             ⋅
+   -1.0+0.0im  2.3618+0.0im    -1.0+0.0im         ⋅
+        ⋅        -1.0+0.0im  2.3618+0.0im    -1.0+0.0im
+        ⋅             ⋅        -1.0+0.0im  2.1382+0.0im
+```
+Note the `ρ[i]` above. This indexes `ρ` at site `i`. For a multiorbital hamiltonian, this will be a matrix (the local density matrix on each site `i`). Here it is just a number, either ` 0.138197` (sites 1 and 4) or `0.361803` (sites 2 and 3).
+
+!!! tip "Sparse vs dense"
+    The method explained above to build a Hamiltonian `hρ` using `-->` supports all the `SiteSelector` and `HopSelector` functionality of conventional models. Therefore, although the density matrix computed above is dense, its application to the Hamiltonian is sparse: it only touches the onsite matrix elements in this case.

docs/src/advanced/serializers.md

@@ -0,0 +1,63 @@
+
+# Serializers
+
+Serializers are useful to translate between a complex data structure and a stream of plain numbers of a given type. Serialization and deserialization is a common encode/decode operation in programming language.
+
+In Quantica, a `s::Serializer{T}` is an object that takes an `h::AbstractHamiltonian`, a selection of the sites and hoppings to be translated, and an `encoder`/`decoder` pair of functions to translate each element into a portion of the stream. This `s` can then be used to convert the specified elements of `h` into a vector of scalars of type `T` and back, possibly after applying some parameter values. Consider this example from the `serializer` docstring
+```julia
+julia> h1 = LP.linear() |> hopping((r, dr) -> im*dr[1]) - @onsite((r; U = 2) -> U);
+
+julia> as = serializer(Float64, h1; encoder = s -> reim(s), decoder = v -> complex(v[1], v[2]))
+AppliedSerializer : translator between a selection of of matrix elements of an AbstractHamiltonian and a collection of scalars
+  Object            : ParametricHamiltonian
+  Object parameters : [:U]
+  Stream parameter  : :stream
+  Output eltype     : Float64
+  Encoder/Decoder   : Single
+  Length            : 6
+
+julia> v = serialize(as; U = 4)
+6-element Vector{Float64}:
+ -4.0
+  0.0
+ -0.0
+ -1.0
+  0.0
+  1.0
+
+julia> h2 = deserialize!(as, v);
+
+julia> h2 == h1(U = 4)
+true
+
+julia> h3 = hamiltonian(as)
+ParametricHamiltonian{Float64,1,1}: Parametric Hamiltonian on a 1D Lattice in 1D space
+  Bloch harmonics  : 3
+  Harmonic size    : 1 × 1
+  Orbitals         : [1]
+  Element type     : scalar (ComplexF64)
+  Onsites          : 1
+  Hoppings         : 2
+  Coordination     : 2.0
+  Parameters       : [:U, :stream]
+
+julia> h3(stream = v, U = 5) == h1(U = 4)  # stream overwrites the U=5 onsite terms
+true
+```
+The serializer functionality is designed with efficiency in mind. Using the in-place `serialize!`/`deserialize!` pair we can do the encode/decode round trip without allocations
+```
+julia> using BenchmarkTools
+
+julia> v = Vector{Float64}(undef, length(as));
+
+julia> deserialize!(as, serialize!(v, as)) === Quantica.call!(h1, U = 4)
+true
+
+julia> @btime deserialize!($as, serialize!($v, $as));
+  149.737 ns (0 allocations: 0 bytes)
+```
+It also allows powerful compression into relevant degrees of freedom through appropriate use of encoders/decoders, see the `serializer` docstring.
+
+## Serializers of OrbitalSliceArrays
+
+Serialization of `OrbitalSliceArray`s is simpler than for `AbstractHamiltonians`, as there is no need for an intermediate `Serializer` object. To serialize an `m::OrbitalSliceArray` simply do `v = serialize(m)`, or `v = serialize(T::Type, m)` if you want a specific eltype for `v`. To deserialize, just do `m´ = deserialize(m, v)`.

src/apply.jl

@@ -17,7 +17,8 @@ function apply(s::SiteSelector, lat::Lattice{T,E,L}, cells...) where {T,E,L}
     sublats = recursive_push!(Int[], intsublats)
     cells = recursive_push!(SVector{L,Int}[], sanitize_cells(s.cells, Val(L)), cells...)
     unique!(sort!(sublats))
-    unique!(sort!(cells))
+    # we don't sort cells, in case we have received them as an explicit list
+    unique!(cells)
     # isnull: to distinguish in a type-stable way between s.cells === missing and no-selected-cells
     # and the same for sublats
     isnull = (s.cells !== missing && isempty(cells)) || (s.sublats !== missing && isempty(sublats))

src/docstrings.jl

@@ -2452,6 +2452,16 @@ Construct a `Vector{T}` that encodes a selection of matrix elements of `h(; para
 where `h = Quantica.parent_hamiltonian(as)` is the `AbstractHamiltonian` used to build the
 `AppliedSerializer`, see `serializer`.
 
+    serialize(m::OrbitalSliceArray)
+
+Return an `Array` of the same eltype as `m` that contains all the stored matrix elements of
+`m`. See `deserialize` for the inverse operation.
+
+    serialize(T::Type, m::OrbitalSliceArray)
+
+As the above, but return an `Array` of the specified eltype `T`, with matrix elements
+reinterpreted as `T`.
+
 ## See also
     `serializer`, `serialize!`, `deserialize`, `deserialize!`
 
@@ -2466,7 +2476,6 @@ details.
 
 ## See also
     `serialize`, `serialize!`, `deserialize`, `deserialize!`
-
 """
 serialize!
 
@@ -2477,9 +2486,14 @@ Construct `h(; params...)`, where `h = Quantica.parametric_hamiltonian(as)` is t
 `AbstractHamiltonian` enclosed in `as`, with the matrix elements enconded in `v =
 serialize(s)` restored (i.e. overwritten). See `serialize` for details.
 
+    deserialize(m::OrbitalSliceArray, v)
+
+Reconstruct an `OrbitalSliceArray` with the same structure as `m` but with the matrix
+elements enconded in `v`. This `v` is typically the result of a `serialize` call to a
+another similar `m`, but the only requirement is that is has the correct size.
+
 ## See also
     `serializer`, `serialize`, `serialize!`, `deserialize!`
-
 """
 deserialize
 
@@ -2574,3 +2588,73 @@ julia> Quantica.decay_lengths(h(U=2))
 
 """
 decay_lengths
+
+"""
+    meanfield(g::GreenFunction, args...; kw...)
+
+Build a `M::MeanField` object that can be used to compute the Hartree-Fock mean field `Φ`
+between selected sites interacting through a given charge-charge potential. The density
+matrix used to build the mean field is obtained with `densitymatrix(g[pair_selection],
+args...; kw...)`, see `densitymatrix` for details.
+
+The mean field between site `i` and `j` is defined as `Φᵢⱼ = δᵢⱼ hartreeᵢ + fockᵢⱼ`, where
+
+    hartreeᵢ = ν * Q * Σ_k v_H(r_i-r_k) * tr(ρ[k,k]*Q)
+    fockᵢⱼ  = -v_F(r_i-r_j) * Q * ρ[i,j] * Q
+
+Here `Q` is the charge operator, `v_H` and `v_F` are Hartree and Fock
+interaction potentials, and `ρ` is the density matrix evaluated at specific chemical
+potential and temperature. Also `ν = ifelse(nambu, 0.5, 1.0)`, and `v_F(0) = v_H(0) = U`,
+where `U` is the onsite interaction.
+
+## Keywords
+
+- `potential`: charge-charge potential to use for both Hartree and Fock. Can be a number or a function of position. Default: `1
+- `hartree`: charge-charge potential `v_H` for the Hartree mean field. Can be a number or a function of position. Overrides `potential`. Default: `potential`
+- `fock`: charge-charge potential `v_F` for the Fock mean field. Can be a number, a function of position or `nothing`. In the latter case all Fock terms (even onsite) will be dropped. Default: `hartree`
+- `onsite`: charge-charge onsite potential. Overrides both Hartree and Fock potentials for onsite interactions. Default: `hartree(0)`
+- `charge`: a number (in single-orbital systems) or a matrix (in multi-orbital systems) representing the charge operator on each site. Default: `I`
+- `nambu::Bool`: specifies whether the model is defined in Nambu space. In such case, `charge` should also be in Nambu space, typically `SA[1 0; 0 -1]` or similar. Default: `false`
+- `selector::NamedTuple`: a collection of `hopselector` directives that defines the pairs of sites (`pair_selection` above) that interact through the charge-charge potential. Default: `(; range = 0)` (i.e. onsite)
+
+Any additional keywords `kw` are passed to the `densitymatrix` function used to compute the
+mean field, see above
+
+## Evaluation and Indexing
+
+    M(µ = 0, kBT = 0; params...)    # where M::MeanField
+
+Build an `Φ::OrbitalSliceMatrix` that can be indexed at different sites or pairs of sites,
+e.g. `ϕ[sites(2:3), sites(1)]`, see `OrbitalSliceArray` for details.
+
+# Examples
+
+```jldoctest
+julia> m = meanfield(g; selector = (; range = 10), potential = pot, onsite = 3, fock = nothing)
+MeanField{ComplexF64} : builder of Hartree-Fock mean fields
+  Charge type      : scalar (ComplexF64)
+  Hartree pairs    : 1
+  Mean field pairs : 21
+
+julia> g = HP.graphene(orbitals = 2, a0 = 1) |> supercell((1,-1)) |> greenfunction;
+
+julia> M = meanfield(g; selector = (; range = 1), charge = I)
+MeanField{SMatrix{2, 2, ComplexF64, 4}} : builder of Hartree-Fock mean fields
+  Charge type      : 2 × 2 blocks (ComplexF64)
+  Hartree pairs    : 14
+  Mean field pairs : 28
+
+julia> Φ0 = M(0.2, 0.3);
+
+julia> Φ0[sites(1), sites(2)]
+2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
+   0.00239416+2.1684e-19im         -0.0-0.0im
+ -1.82007e-17+1.02672e-18im  0.00239416+1.13841e-17im
+
+julia> Φ0[sites(1)]
+2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
+      5.53838-2.04976e-18im  -2.06596e-17-1.16987e-18im
+ -2.06596e-17+1.16987e-18im       5.53838-1.53265e-18im
+```
+"""
+meanfield