Skip to content

Commit

Permalink
Non-spatial models (e.g. mean-field models) (#248)
Browse files Browse the repository at this point in the history
* non-spatial models

* more tests

* test fix

* fix test

* typo spacial => spatial

* change syntax to -->

* cleanups and test fixes

* cleanup, tests and docs

* more docs

* cleanup

* doctest fix
  • Loading branch information
pablosanjose authored Mar 1, 2024
1 parent 2bb98a5 commit 624c354
Show file tree
Hide file tree
Showing 19 changed files with 578 additions and 233 deletions.
3 changes: 2 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
94 changes: 94 additions & 0 deletions docs/src/tutorial/advanced.md
Original file line number Diff line number Diff line change
@@ -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 |> @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> (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((; 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.
5 changes: 2 additions & 3 deletions docs/src/tutorial/glossary.md
Original file line number Diff line number Diff line change
Expand Up @@ -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**
9 changes: 8 additions & 1 deletion docs/src/tutorial/models.md
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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.
Expand Down
7 changes: 3 additions & 4 deletions src/Quantica.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
37 changes: 20 additions & 17 deletions src/apply.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)))

Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)))
Expand Down
Loading

0 comments on commit 624c354

Please sign in to comment.