Non-spatial models (e.g. mean-field models)

docs/make.jl

@@ -23,7 +23,8 @@ makedocs(;
             "Hamiltonians" => "tutorial/hamiltonians.md",
             "Bandstructures" => "tutorial/bandstructures.md",
             "GreenFunctions" => "tutorial/greenfunctions.md",
-            "Observables" => "tutorial/observables.md"
+            "Observables" => "tutorial/observables.md",
+            "Advanced topics" => "tutorial/advanced.md"
             ],
         "Examples" => "examples.md",
         "API" => "api.md",

docs/src/tutorial/advanced.md

@@ -0,0 +1,94 @@
+# Non-spatial models and self-consistent mean-field problems
+
+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> ρ = densitymatrix(g[])(0.5, 0) ## density matrix at chemical potential `µ=0.5` and temperature `kBT = 0`  on all sites
+4×4 OrbitalSliceMatrix{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, ρ) --> 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).
+
+The above provides the tools to implement self-consistent mean field. We just need to find a fixed point `ρf(ρ) = ρ` of the function `ρf` that produces the density matrix of the system.
+
+In the following example we use the FixedPoint.jl package. It provides a simple utility `afps(f, x0)` that computes the solution of `f(x) = x` with initial condition `x0`. The package requires `x0` to be a (real) `AbstractArray`. Note that any other fixed-point search routine that work with `AbstractArray`s should also work.
+```julia
+julia> using FixedPoint
+
+julia> ρf(hρ; µ = 0.5, kBT = 0, U = 0.1) = ρ -> densitymatrix(greenfunction(hρ(; U, ρ), GS.Spectrum())[])(µ, kBT)
+ρf (generic function with 1 method)
+
+julia> (ρsol, ρerror, iters) = @time afps(ρf(hρ; U = 0.4), real(ρ0), tol = 1e-8, ep = 0.95, vel = 0.0); @show iters, ρerror; ρsol
+  0.000627 seconds (1.91 k allocations: 255.664 KiB)
+(iters, ρerror) = (8, 2.0632459629688071e-10)
+4×4 OrbitalSliceMatrix{Matrix{ComplexF64}}:
+ 0.145836+0.0im  0.227266+0.0im  0.227266+0.0im  0.145836+0.0im
+ 0.227266+0.0im  0.354164+0.0im  0.354164+0.0im  0.227266+0.0im
+ 0.227266+0.0im  0.354164+0.0im  0.354164+0.0im  0.227266+0.0im
+ 0.145836+0.0im  0.227266+0.0im  0.227266+0.0im  0.145836+0.0im
+```
+
+!!! note "Bring your own fixed-point solution!"
+    Note that fixed-point calculations can be tricky, and the search algorithm can have a huge impact in convergence (if the problem converges at all!). For this reason, Quantica.jl does not provide built-in fixed-point routines, only the functionality to write functions such as `ρf` above.
+
+!!! tip "Sparse mean fields"
+    The method explained above to build a Hamiltonian 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. Likewise, we could for example use `@hopping!` with a finite `range` to apply a Fock mean field within a finite range. In the future we will support built-in Hartree-Fock model presets with adjustable sparsity.

docs/src/tutorial/glossary.md

@@ -29,6 +29,5 @@ This is a summary of the type of objects you will be studying.
 - **`GreenFunction`**: an `OpenHamiltonian` combined with a `GreenSolver`, which is an algorithm that can in general compute the retarded or advanced Green function at any energy between any subset of sites of the underlying lattice.
   - **`GreenSlice`**: a `GreenFunction` evaluated on a specific set of sites, but at an unspecified energy
   - **`GreenSolution`**: a `GreenFunction` evaluated at a specific energy, but on an unspecified set of sites
-- **`Observable`**: a physical observable that can be expressed in terms of a `GreenFunction`.
-
-  Examples of supported observables include local density of states, current density, transmission probability, conductance and Josephson current
+- **`OrbitalSliceArray`**: an `AbstractArray` that can be indexed with a `SiteSelector`, in addition to the usual scalar indexing. Particular cases are `OrbitalSliceMatrix` and `OrbitalSliceVector`. This is the most common type obtained from `GreenFunction`s and observables obtained from them.
+- **Observables**: Supported observables, obtained from Green functions using various algorithms, include **local density of states**, **density matrices**, **current densities**, **transmission probabilities**, **conductance** and **Josephson currents**

docs/src/tutorial/models.md

@@ -120,7 +120,8 @@ ParametricModel: model with 1 term
     Hopping range     : Neighbors(1)
     Reverse hops      : false
     Coefficient       : 1
-    Parameters        : [:B, :t]
+    Argument type     : spatial
+    Parameters        : [:Bz, :t]
 ```
 
 One can linearly combine parametric and non-parametric models freely, omit parameter default values, and use any of the functional argument forms described for `onsite` and `hopping` (although not the constant argument form):
@@ -134,6 +135,7 @@ ParametricModel: model with 2 terms
     Hopping range     : Neighbors(1)
     Reverse hops      : false
     Coefficient       : -4
+    Argument type     : spatial
     Parameters        : [:t]
   OnsiteTerm{Int64}:
     Region            : any
@@ -142,6 +144,10 @@ ParametricModel: model with 2 terms
     Coefficient       : 2
 ```
 
+!!! tip "Non-spatial parametric models with -->"
+    The `->` in the above parametric models `@onsite` and `@hopping`, but also in the modifiers below, can be changed to `-->`. This indicates that the function arguments are no longer treated as site or link positions `r` and `dr`, but as objects `i, j` representing destination and source sites. This allows to address sites directly instead of through their spatial location. See the Mean Field section for further details.
+
+
 ## Modifiers
 
 There is a third model-related functionality known as `OnsiteModifier`s and `HoppingModifier`s. Given a model that defines a set of onsite and hopping amplitudes on a subset of sites and hops, one can define a parameter-dependent modification of a subset of said amplitudes. This is a useful way to introduce a new parameter dependence on an already defined model. Modifiers are built with
@@ -160,6 +166,7 @@ HoppingModifier{ParametricFunction{3}}:
   Cell distances    : any
   Hopping range     : Inf
   Reverse hops      : false
+  Argument type     : spatial
   Parameters        : [:B]
 ```
 The difference with `model_perierls` is that `model_perierls!` will never add any new hoppings. It will only modify previously existing hoppings in a model. Modifiers are not models themselves, and cannot be summed to other models. They are instead meant to be applied sequentially after applying a model.

src/Quantica.jl

@@ -21,10 +21,9 @@ using Statistics: mean
 using QuadGK
 using SpecialFunctions
 
-export sublat, bravais_matrix, lattice, sites, supercell,
-       hopping, onsite, @onsite, @hopping, @onsite!, @hopping!, plusadjoint, neighbors,
-       siteselector, hopselector, diagonal,
-       hamiltonian,
+export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian,
+       hopping, onsite, @onsite, @hopping, @onsite!, @hopping!, pos, ind, cell,
+       plusadjoint, neighbors, siteselector, hopselector, diagonal,
        unflat, torus, transform, translate, combine,
        spectrum, energies, states, bands, subdiv,
        greenfunction, selfenergy, attach, contact, cellsites,

src/apply.jl

@@ -130,10 +130,10 @@ end
 apply(m::TightbindingModel, latos) = TightbindingModel(apply.(terms(m), Ref(latos)))
 
 apply(t::ParametricOnsiteTerm, lat::Lattice) =
-    ParametricOnsiteTerm(functor(t), apply(selector(t), lat), coefficient(t))
+    ParametricOnsiteTerm(functor(t), apply(selector(t), lat), coefficient(t), is_spatial(t))
 
 apply(t::ParametricHoppingTerm, lat::Lattice) =
-    ParametricHoppingTerm(functor(t), apply(selector(t), lat), coefficient(t))
+    ParametricHoppingTerm(functor(t), apply(selector(t), lat), coefficient(t), is_spatial(t))
 
 apply(m::ParametricModel, lat) = ParametricModel(apply.(terms(m), Ref(lat)))
 
@@ -152,7 +152,8 @@ function apply(m::OnsiteModifier, h::Hamiltonian, shifts = missing)
     asel = apply(sel, lattice(h))
     ptrs = pointers(h, asel, shifts)
     B = blocktype(h)
-    return AppliedOnsiteModifier(sel, B, f, ptrs)
+    spatial = is_spatial(m)
+    return AppliedOnsiteModifier(sel, B, f, ptrs, spatial)
 end
 
 function apply(m::HoppingModifier, h::Hamiltonian, shifts = missing)
@@ -161,14 +162,16 @@ function apply(m::HoppingModifier, h::Hamiltonian, shifts = missing)
     asel = apply(sel, lattice(h))
     ptrs = pointers(h, asel, shifts)
     B = blocktype(h)
-    return AppliedHoppingModifier(sel, B, f, ptrs)
+    spatial = is_spatial(m)
+    return AppliedHoppingModifier(sel, B, f, ptrs, spatial)
 end
 
-function pointers(h::Hamiltonian{T,E}, s::AppliedSiteSelector{T,E}, shifts) where {T,E}
+function pointers(h::Hamiltonian{T,E,L,B}, s::AppliedSiteSelector{T,E,L}, shifts) where {T,E,L,B}
     isempty(cells(s)) || argerror("Cannot constrain cells in an onsite modifier, cell periodicity is assumed.")
-    ptr_r = Tuple{Int,SVector{E,T},Int}[]
+    ptrs = Tuple{Int,SVector{E,T},CellSitePos{T,E,L,B},Int}[]
     lat = lattice(h)
     har0 = first(harmonics(h))
+    dn0 = zerocell(lat)
     umat = unflat(har0)
     rows = rowvals(umat)
     norbs = norbitals(h)
@@ -179,19 +182,20 @@ function pointers(h::Hamiltonian{T,E}, s::AppliedSiteSelector{T,E}, shifts) wher
         r = apply_shift(shifts, r, col)
         if (scol, r) in s
             n = norbs[scol]
-            push!(ptr_r, (p, r, n))
+            sp = CellSitePos(dn0, col, r, B)
+            push!(ptrs, (p, r, sp, n))
         end
     end
-    return ptr_r
+    return ptrs
 end
 
-function pointers(h::Hamiltonian{T,E}, s::AppliedHopSelector{T,E}, shifts) where {T,E}
+function pointers(h::Hamiltonian{T,E,L,B}, s::AppliedHopSelector{T,E,L}, shifts) where {T,E,L,B}
     hars = harmonics(h)
-    ptr_r_dr = [Tuple{Int,SVector{E,T},SVector{E,T},Tuple{Int,Int}}[] for _ in hars]
+    ptrs = [Tuple{Int,SVector{E,T},SVector{E,T},CellSitePos{T,E,L,B},CellSitePos{T,E,L,B},Tuple{Int,Int}}[] for _ in hars]
     lat = lattice(h)
     dn0 = zerocell(lat)
     norbs = norbitals(h)
-    for (har, ptr_r_dr) in zip(hars, ptr_r_dr)
+    for (har, ptrs) in zip(hars, ptrs)
         mh = unflat(har)
         rows = rowvals(mh)
         for scol in sublats(lat), col in siterange(lat, scol), p in nzrange(mh, col)
@@ -205,11 +209,12 @@ function pointers(h::Hamiltonian{T,E}, s::AppliedHopSelector{T,E}, shifts) where
             if (scol => srow, (r, dr), dn) in s
                 ncol = norbs[scol]
                 nrow = norbs[srow]
-                push!(ptr_r_dr, (p, r, dr, (nrow, ncol)))
+                srow, scol = CellSitePos(dn, row, rrow, B), CellSitePos(dn0, col, rcol, B)
+                push!(ptrs, (p, r, dr, srow, scol, (nrow, ncol)))
             end
         end
     end
-    return ptr_r_dr
+    return ptrs
 end
 
 apply_shift(::Missing, r, _) = r
@@ -288,17 +293,15 @@ function apply_map(mapping, hf::Function, ::Type{S}) where {T,S<:SVector{<:Any,T
     return FunctionWrapper{SparseMatrixCSC{Complex{T},Int},Tuple{S}}(sfunc)
 end
 
-
 #endregion
 
-
 ############################################################################################
 # apply CellSites
 #region
 
-apply(c::CellSites{L,Vector{Int}}, ::Lattice{<:Any,<:Any,L}) where {L} = c
+apply(c::AnyCellSite, ::Lattice{<:Any,<:Any,L}) where {L} = c
 
-apply(c::CellSites{L,Int}, ::Lattice{<:Any,<:Any,L}) where {L} = c
+apply(c::CellSites{L,Vector{Int}}, ::Lattice{<:Any,<:Any,L}) where {L} = c
 
 apply(c::CellSites{L,Colon}, l::Lattice{<:Any,<:Any,L}) where {L} =
     CellSites(cell(c), collect(siterange(l)))