From 409d449011eed17e52c9b82d77c1c9d7220cb69c Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 25 Oct 2024 12:47:11 +0200 Subject: [PATCH] general Hartree-Fock mean fields --- docs/make.jl | 2 + docs/src/advanced/meanfield.md | 273 +++++++++++++------------------ docs/src/advanced/nonspatial.md | 71 ++++++++ docs/src/advanced/serializers.md | 63 +++++++ src/apply.jl | 3 +- src/docstrings.jl | 88 +++++++++- src/meanfield.jl | 206 +++++++++++++---------- src/sanitizers.jl | 1 + src/serializer.jl | 23 ++- src/show.jl | 20 +++ src/specialmatrices.jl | 7 +- src/tools.jl | 3 + src/types.jl | 4 + test/test_greenfunction.jl | 76 +++++++-- test/test_show.jl | 2 + 15 files changed, 570 insertions(+), 272 deletions(-) create mode 100644 docs/src/advanced/nonspatial.md create mode 100644 docs/src/advanced/serializers.md diff --git a/docs/make.jl b/docs/make.jl index 931c0194..74ff8337 100644 --- a/docs/make.jl +++ b/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" ], diff --git a/docs/src/advanced/meanfield.md b/docs/src/advanced/meanfield.md index 4b7b9363..1fdf24a3 100644 --- a/docs/src/advanced/meanfield.md +++ b/docs/src/advanced/meanfield.md @@ -1,88 +1,58 @@ -# 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] -``` +# Self-consistent mean-field problems + +Here we show how to solve interacting-electron problems in Quantica, approximated at the mean field level. A mean field is a collection of onsite and hopping terms that are added to a given `h::AbstractHamiltonian`, that depend on the density matrix `ρ`. Since `ρ` itself depends on `h`, this defines a self-consistent problem. -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)`. +If the mean field solution is dubbed `Φ`, the problem consists in finding a fixed point solution to the function `Φ = M(Φ)` for a certain function `M` that takes `Φ`, computes `h` with the added mean field onsite and hopping terms, computes the density matrix, and from that computes the new mean field `Φ`. To attack this problem we will employ non-spatial models and a new `meanfield` constructor. -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: +Schematically the process is as follows: +- We start from an `AbstractHamiltonian` that includes a non-interacting model `model_0` and non-spatial model `model_1 + model_2` with a mean field parameter, e.g. `Φ`, ```julia -julia> h = LP.linear() |> onsite(2) - hopping(1) |> supercell(4) |> supercell; +julia> model_0 = hopping(1); # a possible non-interacting model -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> model_1 = @onsite((i; Φ = missing) --> Φ[i]); -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 +julia> model_2 = @hopping((i, j; Φ = missing) -> Φ[i, j]);Fock + +julia> h = lat |> hamiltonian(model_0 + model_1 + model_2) ``` -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). +Here `model_1` corresponds to Hartree and onsite-Fock mean field terms, while `model_2` corresponds to inter-site Fock terms. +The default value of `Φ = missing` corresponds to a vanishing `model_1 + model_2`, i.e. no interactions. +- We build the `GreenFunction` of `h` with `g = greenfunction(h, solver; kw...)` using the `GreenSolver` of choice +- We construct a `M::MeanField` object using `M = meanfield(g; potential = pot, other_options...)` -!!! 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. Likewise, we could for example use `@hopping!` with a finite `range` to apply a Fock mean field within a finite range. + Here `pot(r)` is the charge-charge interaction potential between electrons. We can also specify `hopselector` directives to define which sites interacts, adding e.g. `selector = (; range = 2)` to `other_options`, to make sites at distance `2` interacting. See `meanfield` docstring for further details. -# Self-consistent mean-field problems +- We evaluate this `M` with `Φ0 = M(µ, kBT; params...)`. + + This computes the density matrix at specific chemical potential `µ` and temperature `kBT`, and for specific parameters of `h` (possibly including `Φ`). Then it computes the appropriate Hartree and Fock terms, and stores them in the returned `Φ0::OrbitalSliceMatrix`, where `Φ0ᵢⱼ = δᵢⱼ hartreeᵢ + fockᵢⱼ`. In normal systems, these terms read + + $$\text{hartree}_i = Q \sum_k v_H(r_i-r_k) \text{tr}(\rho_{kk}Q)$$ + + $$\text{fock}_{ij} = -v_F(r_i-r_j) Q \rho_{ij} Q$$ -The above provides the tools to implement self-consistent mean field problems in Quantica. The problem consists of finding a fixed point `f(h) = h` of the function + where `v_H` and `v_F` are Hartree and Fock charge-charge interaction potentials (by default equal to `pot`), and the charge operator is `Q` (equal to the identity by default, but can be changed to implement e.g. spin-spin interactions). + + When computing `Φ0` we don't specify `Φ` in `params`, so that `Φ0` is evaluated using the non-interacting model, hence its name. + +- The self-consistent condition can be tackled naively by iteration-until-convergence, ```julia -function f(h) - g = greenfunction(h, GS.Spectrum()) - ρ = densitymatrix(g[])(0.5, 0) - return hρ(; ρ) -end +Φ0 = M(µ, kBT; params...) +Φ1 = M(µ, KBT; Φ = Φ0, params...) +Φ2 = M(µ, KBT; Φ = Φ1, params...) +... ``` -that takes a Hamiltonian `h`, computes `ρ` with it, and returns a new Hamiltonian `h = hρ(ρ)`. Its fixed point corresponds to the self-consistent Hamiltonian, and the corresponding `ρ` is the self-consistent mean-field density matrix. +A converged solution `Φ`, if found, should satisfy the fixed-point condition + + Φ_sol ≈ M(µ, KBT; Φ = Φ_sol, params...) -In the following we show how to approach such as problem. We will employ an external library to solve generic fixed-point problems, and show how to make it work with Quantica AbstractHamiltonians efficiently. Many generic fixed-point solver backends exist. Here we use the SIAMFANLEquations.jl package. It provides a simple utility `aasol(f, x0)` that computes the solution of `f(x) = x` with initial condition `x0` using Anderson acceleration. This is an example of how it works to compute the fixed point of function `f(x) = 1 + atan(x)` +Then, the self-consistent Hamiltonian is given by `h(; Φ = Φ_sol, params...)`. + +The key problem is to actually find the fixed point of the `M` function. The naive iteration above is not optimal, and often does not converge. To do a better job we should use a dedicated fixed-point solver. + +## Using an external fixed-point solver + +We now show how to approach such a fixed-point problem. We will employ an external library to solve generic fixed-point problems, and show how to make it work with Quantica `MeanField` objects efficiently. Many generic fixed-point solver backends exist. Here we use the SIAMFANLEquations.jl package. It provides a simple utility `aasol(f, x0)` that computes the solution of `f(x) = x` with initial condition `x0` using Anderson acceleration. This is an example of how it works to compute the fixed point of function `f(x) = 1 + atan(x)` ```julia julia> using SIAMFANLEquations @@ -96,119 +66,96 @@ julia> aasol(f!, x0, m, vstore).solution 2.132267725272908 2.132267725284556 ``` -The package requires as input an in-place version `f!` of the function `f`, and the preallocation of some storage space vstore (see the `aasol` documentation). The package, as [a few others](https://docs.sciml.ai/NonlinearSolve/stable/solvers/fixed_point_solvers/), also requires the variable `x` and the initial condition `x0` to be (real) `AbstractArray` (or a scalar, but we need the former for our use case), hence the broadcast dots above. In our case we will therefore need to translate back and forth from a Hamiltonian `h` to a real vector `x` to pass it to `aasol`. - -## Serializers - -This translation is achieved with Quantica's `serializer` funcionality. 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. 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 self-explanatory 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. +The package requires as input an in-place version `f!` of the function `f`, and the preallocation of some storage space vstore (see the `aasol` documentation). The package, as [a few others](https://docs.sciml.ai/NonlinearSolve/stable/solvers/fixed_point_solvers/), also requires the variable `x` and the initial condition `x0` to be (real) `AbstractArray` (or a scalar, but we need the former for our use case), hence the broadcast dots above. In our case we will therefore need to translate back and forth from an `Φ::OrbitalSliceMatrix` to a real vector `x` to pass it to `aasol`. This translation is achieved using Quantica's `serialize`/`deserialize` funcionality. ## Using Serializers with fixed-point solvers -With the `serializer` functionality we can build a version of the fixed-point function `f` that operates on real vectors. Let's return to our original Hamiltonian +With the serializer functionality we can build a version of the fixed-point function `f` that operates on real vectors. Let's take a 1D Hamiltonian with a sawtooth potential, and build a Hartree mean field (note the `fock = nothing` keyword) ```julia -julia> h = LP.linear() |> onsite(2) - hopping(1) |> supercell(4) |> supercell; +julia> using SIAMFANLEquations -julia> ρ0 = densitymatrix(greenfunction(h, GS.Spectrum())[])(0.5, 0) +julia> h = LP.linear() |> supercell(4) |> onsite(r->r[1]) - hopping(1) + @onsite((i; phi = missing) --> phi[i]); -julia> hρ = h |> @onsite!((o, i; U = 0, ρ = ρ0) --> o + U * ρ[i]); +julia> M = meanfield(greenfunction(h); onsite = 1, selector = (; range = 0), fock = nothing) +MeanField{ComplexF64} : builder of Hartree-Fock mean fields + Charge type : scalar (ComplexF64) + Hartree pairs : 4 + Mean field pairs : 4 -julia> s = serializer(Float64, hρ, siteselector(); encoder = real, decoder = identity) -Serializer{Float64} : encoder/decoder of matrix elements into a collection of scalars - Object : ParametricHamiltonian - Output eltype : Float64 - Encoder/Decoder : Single - Length : 4 +julia> Φ0 = M(0.0, 0.0); -julia> function f!(x, x0, (s, params)) - h = deserialize!(s, x0) - g = greenfunction(h, GS.Spectrum()) - ρ = densitymatrix(g[])(0.5, 0) - serialize!(x, s; ρ = ρ, params...) +julia> function f!(x, x0, (M, Φ0)) + Φ = M(0.0, 0.0; phi = deserialize(Φ0, x0)) + copy!(x, serialize(Float64, Φ)) return x end; ``` Then we can proceed as in the `f(x) = 1 + atan(x)` example ```julia -julia> m = 2; len = length(s); x0 = rand(len); vstore = rand(len, 3m+3); # order m, initial condition x0, and preallocated space vstore - -julia> x = aasol(f!, x0, m, vstore; pdata = (s, (; U = 0.3))).solution -4-element Vector{Float64}: - 2.04319819150754 - 2.1068018084891236 - 2.1068018084820492 - 2.0431981915212862 - -julia> h´ = deserialize!(s, x) -Hamiltonian{Float64,1,0}: Hamiltonian on a 0D Lattice in 1D space - Bloch harmonics : 1 +julia> m = 2; x0 = serialize(Float64, Φ0); vstore = rand(length(x0), 3m+3); # order m, initial condition x0, and preallocated space vstore + +julia> x = aasol(f!, x0, m, vstore; pdata = (M, Φ0)).solution +Iteration terminates on entry to aasol.jl +8-element Vector{Float64}: + 0.565818503090569 + 0.0 + 0.3062161093223203 + 0.0 + 0.0669636234267607 + 0.0 + 0.061001764160349956 + 0.0 + +julia> h´ = h(; phi = deserialize(Φ0, x)) +Hamiltonian{Float64,1,1}: Hamiltonian on a 1D Lattice in 1D space + Bloch harmonics : 3 Harmonic size : 4 × 4 Orbitals : [1] Element type : scalar (ComplexF64) Onsites : 4 - Hoppings : 6 - Coordination : 1.5 + Hoppings : 8 + Coordination : 2.0 julia> h´[()] 4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries: - 2.0432+0.0im -1.0+0.0im ⋅ ⋅ - -1.0+0.0im 2.1068+0.0im -1.0+0.0im ⋅ - ⋅ -1.0+0.0im 2.1068+0.0im -1.0+0.0im - ⋅ ⋅ -1.0+0.0im 2.0432+0.0im + 0.565819+0.0im -1.0+0.0im ⋅ ⋅ + -1.0+0.0im 1.30622+0.0im -1.0+0.0im ⋅ + ⋅ -1.0+0.0im 2.06696+0.0im -1.0+0.0im + ⋅ ⋅ -1.0+0.0im 3.061+0.0im ``` Note that the content of `pdata` is passed by `aasol` as a third argument to `f!`. We use this to pass the serializer `s` and `U` parameter to use. !!! note "Bring your own fixed-point solver!" 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. Numerous packages exist for fixed-point computations in julia. Check [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) for one prominent metapackage. + +## GreenSolvers without support for ParametricHamiltonians + +Some `GreenSolver`'s, like `GS.Spectrum` or `GS.KPM`, do not support `ParametricHamiltonian`s. In such cases, the approach above will fail, since it will not be possible to build `g` before knowing `phi`. In such cases one would need to rebuild the `meanfield` object at each step of the fixed-point solver. +```julia +julia> using SIAMFANLEquations + +julia> h = LP.linear() |> supercell(4) |> supercell |> onsite(1) - hopping(1) + @onsite((i; phi = missing) --> phi[i]); + +julia> M´(Φ = missing) = meanfield(greenfunction(h(; phi = Φ), GS.Spectrum()); onsite = 1, selector = (; range = 0), fock = nothing) + +julia> Φ0 = M´()(0.0, 0.0); + +julia> function f!(x, x0, (M´, Φ0)) + Φ = M´(deserialize(Φ0, x0))(0.0, 0.0) + copy!(x, serialize(Float64, Φ)) + return x + end; + +julia> m = 2; x0 = serialize(Float64, Φ0); vstore = rand(length(x0), 3m+3); # order m, initial condition x0, and preallocated space vstore + +julia> x = aasol(f!, x0, m, vstore; pdata = (M´, Φ0)).solution +8-element Vector{Float64}: + 0.15596283661183488 + 0.0 + 0.3440371633881653 + 0.0 + 0.3440371633881646 + 0.0 + 0.15596283661183474 + 0.0 +``` diff --git a/docs/src/advanced/nonspatial.md b/docs/src/advanced/nonspatial.md new file mode 100644 index 00000000..46b220d3 --- /dev/null +++ b/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. diff --git a/docs/src/advanced/serializers.md b/docs/src/advanced/serializers.md new file mode 100644 index 00000000..7e26f9d7 --- /dev/null +++ b/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)`. diff --git a/src/apply.jl b/src/apply.jl index d81bcd3b..8394ed00 100644 --- a/src/apply.jl +++ b/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)) diff --git a/src/docstrings.jl b/src/docstrings.jl index 75f8d2e1..bbaebabd 100644 --- a/src/docstrings.jl +++ b/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 diff --git a/src/meanfield.jl b/src/meanfield.jl index 6428c455..03cd4513 100644 --- a/src/meanfield.jl +++ b/src/meanfield.jl @@ -1,84 +1,122 @@ -# ############################################################################################ -# # meanfield -# # designed to construct hfield and ffield, such that -# # hfield[i] = ν * Σ_k v_H(r_i-r_k) * tr(ρ[k,k]*Q) -# # ffield[i,j] = v_F(r_i-r_j) * Q * ρ[i,j] * Q -# # where ν = ifelse(nambu, 1/2, 1), and Q is the charge matrix or [q 0; 0 -q] if nambu. -# # we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry -# #region - -# struct MeanField{T,O<:OrbitalSliceMatrix{T},H,F,Q,B} -# Vhmat::SparseMatrixCSC{T,Int} -# Vfmat::O -# hrho::H -# frho::F -# charge::Q -# blocktype::B -# end - -# struct EvaluatedMeanField{H<:OrbitalSliceVector,F<:OrbitalSliceMatrix} -# hfield::H -# ffield::F -# end - -# #region ## Constructors ## - -# function meanfield(g::GreenFunction{T}, args...; -# potential = Returns(0), hartreepotential = potential, fockpotential = hartreepotential, -# charge = I, nambu::Bool = false, selector = (; range = 0), kw...) where {T} - -# Vh = sanitize_potential(hartreepotential) -# Vf = sanitize_potential(fockpotential) - -# B = blocktype(hamiltonian(g)) -# lat = lattice(hamiltonian(g)) -# hFock = lat |> hopping((r, dr) -> Vf(dr); selector..., includeonsite = true) -# hHartree = Vh === Vf ? hFock : -# lat |> hopping((r, dr) -> Vh(dr); selector..., includeonsite = true) - -# ldst, lsrc = ham_to_latslices(hFock) -# odst, osrc = sites_to_orbs(ldst, g), sites_to_orbs(lsrc, g) -# ρhartree = densitymatrix(g[diagonal(osrc, kernel = charge)], args...; kw...) -# ρfock = densitymatrix(g[odst, osrc], args...; kw...) - -# Vhmat = sum(unflat, harmonics(hHartree)) -# nambu && (nonzeros(Vhmat) .*= T(1/2)) -# Vfmat = hFock[ldst, lsrc] - -# return MeanField(Vhmat, Vfmat, ρhartree, ρfock, charge, B) -# end - -# sanitize_potential(x::Number) = Returns(x) -# sanitize_potential(x::Function) = x - -# sanitize_charge(x, nambu) = nambu ? SA[x 0*x; 0*x -x] : x - -# #endregion - -# #region ## API ## - - -# function (m::MeanField{T})(args...; params...) where {T} -# Q, B = m.charge, m.blocktype -# trρQ = m.hrho(args...; params...) -# QvtrρQ = (m.Vhmat * parent(trρQ)) .* Ref(Q) -# hfield = OrbitalSliceVector(QvtrρQ, orbaxes(trρQ)) -# ffield = m.frho(args...; params...) -# if !isempty(ffield) -# oi, oj = orbaxes(ffield) -# for sj in sites(oj, B), si in sites(oi, B) # si, sj are CellSitePos -# vij = T(only(view(m.Vfmat, CellSite(si), CellSite(sj)))) # a scalar -# view(ffield, si, sj) .= vij * (Q * ffield[si, sj] * Q) -# end -# end -# return EvaluatedMeanField(hfield, ffield) -# end - -# Base.view(m::EvaluatedMeanField, i) = view(m.hfield, i) -# Base.view(m::EvaluatedMeanField, i, j) = view(m.ffield, i, j) - -# Base.getindex(m::EvaluatedMeanField, i) = m.hfield[i] -# Base.getindex(m::EvaluatedMeanField, i, j) = m.ffield[i, j] -# #endregion - -# #endregion +############################################################################################ +# meanfield +# designed to construct hartreefield and fockfield, such that +# hartreefield[i] = ν * Q * Σ_k v_H(r_i-r_k) * tr(ρ[k,k]*Q) +# fockfield[i,j] = -v_F(r_i-r_j) * Q * ρ[i,j] * Q +# where ν = ifelse(nambu, 1/2, 1), and Q is the charge matrix or [q 0; 0 -q] if nambu. +# we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry +#region + +struct MeanField{B,T,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix} + potHartree::S + potFock::S + rhoHartree::H + rhoFock::F + charge::B + rowcol_ranges::NTuple{2,Vector{UnitRange{Int}}} + onsite_tmp::Vector{Complex{T}} +end + +#region ## Constructors ## + +function meanfield(g::GreenFunction{T,E}, args...; + potential = Returns(1), hartree = potential, fock = hartree, + onsite = missing, charge = I, nambu::Bool = false, + selector::NamedTuple = (; range = 0), kw...) where {T,E} + + Vh = sanitize_potential(hartree) + Vf = sanitize_potential(fock) + Q = sanitize_charge(charge, blocktype(hamiltonian(g))) + U = onsite === missing ? T(Vh(zero(SVector{E,T}))) : T(onsite) + Uf = fock === nothing ? zero(U) : U + + isempty(boundaries(g)) || argerror("meanfield does not currently support systems with boundaries") + isfinite(U) || argerror("Onsite potential must be finite, consider setting `onsite`") + nambu && (!is_square(charge) || iseven(size(charge, 1))) && argerror("Invalid charge matrix for Nambu space") + + gsHartree = g[diagonal(; cells = 0, kernel = Q)] + rhoHartree = densitymatrix(gsHartree, args...; kw...) + gsFock = g[sitepairs(; selector..., includeonsite = true)] + rhoFock = densitymatrix(gsFock, args...; kw...) + + lat = lattice(hamiltonian(g)) + # The sparse structure of hFock will be inherited by the evaluated mean field. Need onsite. + hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : Vf(dr); selector..., includeonsite = true) + hHartree = (Uf == U && Vh === Vf) ? hFock : + lat |> hopping((r, dr) -> iszero(dr) ? U : Vh(dr); selector..., includeonsite = true) + + potHartree = sum(unflat, harmonics(hHartree)) + nambu && (nonzeros(potHartree) .*= T(1/2)) # to compensate for the factor 2 from nambu trace + + oaxes = orbaxes(call!_output(gsFock)) + rowcol_ranges = collect.(orbranges.(oaxes)) + onsite_tmp = similar(diag(parent(call!_output(gsHartree)))) + + # build potFock with identical axes as the output of rhoFock + cells_fock = cells(first(oaxes)) + hFock_slice = hFock[(; cells = cells_fock), (; cells = 0)] + + # this is important for the fast orbrange-based implementation of MeanField evaluation + check_cell_order(hFock_slice, rhoFock) + potFock = parent(hFock_slice) + + return MeanField(potHartree, potFock, rhoHartree, rhoFock, Q, rowcol_ranges, onsite_tmp) +end + +# currying version +meanfield(; kw...) = g -> meanfield(g; kw...) +meanfield(x; kw...) = g -> meanfield(g, x; kw...) + +sanitize_potential(x::Number) = Returns(x) +sanitize_potential(x::Function) = x +sanitize_potential(x::Nothing) = Returns(0) +sanitize_potential(_) = argerror("Invalid potential: use a number or a function of position") + +sanitize_charge(B, t) = sanitize_block(t, B) +sanitize_charge(B, ::Type{<:SMatrixView}) = argerror("meanfield does not currently support systems with heterogeneous orbitals") + +function check_cell_order(hFock_slice, rhoFock) + opot = first(orbaxes(hFock_slice)) + orho = first(orbaxes(call!_output(rhoFock.gs))) + cells(opot) == cells(orho) || internalerror("meanfield: Cell order mismatch between potential and density matrix") + return nothing +end + +#endregion + +#region ## API ## + +charge(m::MeanField) = m.charge + +hartree_matrix(m::MeanField) = m.potHartree + +fock_matrix(m::MeanField) = parent(m.potFock) + +function (m::MeanField{B})(args...; params...) where {B} + Q, hartree_pot, fock_pot = m.charge, m.onsite_tmp, m.potFock + rowrngs, colrngs = m.rowcol_ranges + trρQ = m.rhoHartree(args...; params...) + mul!(hartree_pot, m.potHartree, diag(parent(trρQ))) + meanfield = m.rhoFock(args...; params...) + mf_parent = parent(meanfield) # should be a sparse matrix + rows, cols, nzs = rowvals(fock_pot), axes(fock_pot, 2), nonzeros(fock_pot) + for col in cols + viiQ = hartree_pot[col] * Q + for ptr in nzrange(fock_pot, col) + row = rows[ptr] + vij = nzs[ptr] + ρij = view(mf_parent, rowrngs[row], colrngs[col]) + vQρijQ = vij * Q * sanitize_block(B, ρij) * Q + if row == col + ρij .= viiQ - vQρijQ + else + ρij .= -vQρijQ + end + end + end + return meanfield +end + +#endregion + +#endregion diff --git a/src/sanitizers.jl b/src/sanitizers.jl index 71619c3c..28d54d95 100644 --- a/src/sanitizers.jl +++ b/src/sanitizers.jl @@ -132,6 +132,7 @@ sanitize_cellindices(c::CellIndices{L}, ::TELtypes{<:Any,<:Any,L´}) where {L,L #region sanitize_block(::Type{C}, a) where {C<:Number} = complex(C)(only(a)) +sanitize_block(::Type{C}, a::UniformScaling) where {C<:Number} = complex(C)(a.λ) sanitize_block(::Type{C}, a) where {C<:SMatrix} = C(a) # here we assume a is already of the correct size and let if fail later otherwise sanitize_block(::Type{C}, a) where {C<:SMatrixView} = a diff --git a/src/serializer.jl b/src/serializer.jl index b25d6c45..d75ab50c 100644 --- a/src/serializer.jl +++ b/src/serializer.jl @@ -135,17 +135,30 @@ end serialize(a::OrbitalSliceArray) = serialize_array(parent(a)) serialize(a::AbstractArray) = serialize_array(a) +serialize(::Type{T}, a::AbstractArray{T}) where {T} = serialize(a) +serialize(::Type{T}, a::AbstractArray) where {T} = reinterpret(T, serialize(a)) serialize_array(a::Array) = a serialize_array(a::Diagonal{<:Any,<:Vector}) = a.diag serialize_array(a::SparseMatrixCSC) = nonzeros(a) -deserialize(a::OrbitalSliceArray, v) = OrbitalSliceArray(deserialize_array(parent(a), v), orbaxes(a)) -deserialize(a::AbstractArray, v) = deserialize_array(a, v) +deserialize(a, v) = (check_serializer_length(a, v); _deserialize(a, v)) -deserialize_array(::A, v::A´) where {N,T,T´,A<:Array{T,N},A´<:Array{T´,N}} = v -deserialize_array(::Diagonal, v::Vector) = Diagonal(v) -deserialize_array(a::SparseMatrixCSC, v::Vector) = SparseMatrixCSC(a.m, a.n, a.colptr, a.rowval, v) +_deserialize(a::OrbitalSliceArray{T}, v::AbstractArray{T}) where {T} = + OrbitalSliceArray(deserialize_array(parent(a), v), orbaxes(a)) +_deserialize(a::AbstractArray{T}, v::AbstractArray{T}) where {T} = deserialize_array(a, v) +_deserialize(a::AbstractArray{T}, v::AbstractArray) where {T} = deserialize(a, reinterpret(T, v)) +deserialize_array(::AbstractArray{<:Any,N}, v::AbstractArray{<:Any,N}) where {N} = v +deserialize_array(::Diagonal, v::AbstractVector) = + Diagonal(convert(Vector, v)) +deserialize_array(a::SparseMatrixCSC, v::AbstractVector) = + SparseMatrixCSC(a.m, a.n, a.colptr, a.rowval, convert(Vector, v)) + +check_serializer_length(a::AbstractArray{T}, v::AbstractArray{T}) where {T} = + length(serialize(a)) == length(v) || argerror("Wrong length of serialized array") + +check_serializer_length(a::AbstractArray, v::AbstractArray{T}) where {T} = + length(serialize(T, a)) == length(v) || argerror("Wrong length of serialized array") #endregion diff --git a/src/show.jl b/src/show.jl index 0d861e65..4deb696b 100644 --- a/src/show.jl +++ b/src/show.jl @@ -614,4 +614,24 @@ display_encdec(::Tuple) = "(Onsite, Hopping pairs)" display_parameters(s::AppliedSerializer{<:Any,<:Hamiltonian}) = "none" display_parameters(s::AppliedSerializer{<:Any,<:ParametricHamiltonian}) = string(parameters(parent_hamiltonian(s))) + +#endregion + + +############################################################################################ +# MeanField and EvaluatedMeanField +#region + +function Base.show(io::IO, s::MeanField{Q}) where {Q} + i = get(io, :indent, "") + ioindent = IOContext(io, :indent => i * " ") + print(io, i, summary(s), "\n", +"$i Charge type : $(displaytype(Q)) +$i Hartree pairs : $(nnz(hartree_matrix(s))) +$i Mean field pairs : $(nnz(fock_matrix(s)))") +end + +Base.summary(::MeanField{Q}) where {Q} = + "MeanField{$Q} : builder of Hartree-Fock mean fields" + #endregion diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index 67333aaf..5cf6d4ef 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -541,12 +541,13 @@ end # CellSites: return an unwrapped Matrix or a scalar Base.getindex(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) = copy(view(a, i, j)) -Base.getindex(a::OrbitalSliceMatrix, i::CellSite, j::CellSite = i) = - only(view(a, i, j)) Base.getindex(a::OrbitalSliceMatrix, i::C, j::C = i) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} = sanitize_block(B, view(a, i, j)) +# it also returns 0I if the matrix is missing (useful for meanfield models) +Base.getindex(::Missing, ::CellSitePos, ::CellSitePos...) = 0I + function Base.view(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) rowslice, colslice = orbaxes(a) i´, j´ = apply(i, lattice(rowslice)), apply(j, lattice(colslice)) @@ -687,5 +688,7 @@ Base.eltype(::EigenProduct{T}) where {T} = Complex{T} Base.:*(x::Number, P::EigenProduct) = EigenProduct(P.blockstruct, P.U, P.V, P.phi, P.x * x) Base.:*(P::EigenProduct, x::Number) = x*P +LinearAlgebra.tr(P::EigenProduct) = sum(xy -> first(xy) * conj(last(xy)), zip(P.U, P.V)) * P.x + #endregion #endregion diff --git a/src/tools.jl b/src/tools.jl index 2554be67..748e37de 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -211,6 +211,9 @@ Base.@propagate_inbounds function swapcols!(a::AbstractMatrix, i, j) end end +is_square(a::AbstractMatrix) = size(a, 1) == size(a, 2) +is_square(a) = false + #endregion ############################################################################################ diff --git a/src/types.jl b/src/types.jl index ee0e91ef..756b2b8f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -431,6 +431,10 @@ CellSitePos(cell, ind, r, B) = CellIndices(sanitize_SVector(Int, cell), ind, Sit # changing only cell CellIndices(cell, c::CellIndices) = CellIndices(cell, c.inds, c.type) +zerocell(i::CellIndices) = CellIndices(zero(i.cell), i.inds, i.type) +zerocell(i::CellIndices, j::CellIndices) = + CellIndices(i.cell-j.cell, i.inds, i.type), CellIndices(zero(j.cell), j.inds, j.type) + # LatticeSlice from an AbstractVector of CellIndices LatticeSlice(lat::Lattice, cs::AbstractVector{<:CellIndices}) = LatticeSlice(lat, cellinds_to_dict(apply.(cs, Ref(lat)))) diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index 6a74f30d..e61ed22d 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -346,6 +346,10 @@ end @test isapprox(ρ0(), ρ(); atol = 1e-8) @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-8) end + g = HP.graphene(orbitals = 2) |> supercell((1,-1)) |> greenfunction + ρ0 = densitymatrix(g[sitepairs(range = 1)], 6) + mat = ρ0(0.2, 0.3) + @test Quantica.nnz(parent(mat)) < size(mat) # no broken sparsity end @testset "OrbitalSliceArray slicing" begin @@ -527,21 +531,6 @@ end @test typeof(ρ0sol) == typeof(ρsol) end -@testset "mean-field models" begin - h1 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = 1) - h2 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = 2) - h3 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = (1,2)) - for h0 in (h1, h2, h3) - ρ0 = densitymatrix(greenfunction(h0, GS.Spectrum())[cells = SA[]])(); - h = h0 |> @onsite!((o, s; ρ = ρ0, t) --> o + t*ρ[s]) - @test diag(h(t = 2)[()]) ≈ 2*diag(ρ0) atol = 0.0000001 - h = h0 |> @hopping!((t, si, sj; ρ = ρ0, α = 2) --> α*ρ[si, sj]) - @test h() isa Quantica.Hamiltonian - diff = (h()[()] - 2ρ0) .* h()[()] - @test iszero(diff) - end -end - @testset "aliasing" begin # Issue #267 g = LP.linear() |> hamiltonian(@hopping((; q = 1) -> q*I), orbitals = 2) |> greenfunction @@ -573,3 +562,60 @@ end @test g.contacts.selfenergies[2].solver.hlead[(1,)] === g.contacts.selfenergies[1].solver.hlead[(-1,)] @test g.contacts.selfenergies[2].solver.hlead[(-1,)] === g.contacts.selfenergies[1].solver.hlead[(1,)] end + +@testset "meanfield" begin + oh = LP.honeycomb() |> hamiltonian(hopping(SA[1 im; -im 1]) - onsite(SA[1 0; 0 -1], sublats = :A), orbitals = 2) |> supercell((1,-1)) + g = oh |> greenfunction + Q = SA[0 -im; im 0] + m = meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q) + Φ = m(0.2, 0.3) + @test Φ[sites(1), sites(2)] ≈ -1.5 * Q * m.rhoFock(0.2, 0.3)[sites(1), sites(2)] * Q + + g = LP.linear() |> hopping(1) |> greenfunction + pot(r) = 1/norm(r) + Φ = meanfield(g; selector = (; range = 10), potential = pot, onsite = 3)() + # onsite fock cancels hartree in the spinless case + @test only(view(Φ, sites(1))) ≈ 0.5 * (0*3 + 2*sum(pot, 1:10)) + # unless we disable fock + Φ = meanfield(g; selector = (; range = 10), potential = pot, onsite = 3, fock = nothing)() + @test only(view(Φ, sites(1))) ≈ 0.5 * (3 + 2*sum(pot, 1:10)) + + g = oh |> greenfunction + @test_throws ArgumentError meanfield(g, nambu = true, charge = I) + g = oh |> greenfunction(GS.Schur(boundary = 0)) + @test_throws ArgumentError meanfield(g) + g = LP.honeycomb() |> hamiltonian(hopping(I), orbitals = (2,1)) |> greenfunction + @test_throws ArgumentError meanfield(g) +end + +@testset begin "OrbitalSliceArray serialization" + g = LP.linear() |> supercell(4) |> supercell |> hamiltonian(hopping(I), orbitals = 2) |> greenfunction + m = g(0.2)[cells = 0] + v = serialize(m) + @test v === parent(m) + @test deserialize(m, v) === m + v = serialize(Float64, m) + @test v isa AbstractMatrix{Float64} + @test deserialize(m, v) === m + m = g(0.2)[diagonal(cells = 0)] + v = serialize(m) + @test v === parent(m).diag + @test deserialize(m, v) === m + v = serialize(Float64, m) + @test v isa AbstractVector{Float64} + @test deserialize(m, v) === m + m = g(0.2)[sitepairs(range = 1)] + v = serialize(m) + @test v === Quantica.nonzeros(parent(m)) + @test deserialize(m, v) === m + v = serialize(Float64, m) + @test v isa AbstractVector{Float64} + @test deserialize(m, v) === m + + m = g(0.2)[(; cells = 0), (; cells = 0, region = r -> r[1] < 2)] + m´ = g(0.2)[cells = 0] + v = serialize(m) + @test_throws ArgumentError deserialize(m´, v) + v = serialize(Float64, m) + @test_throws ArgumentError deserialize(m´, v) +end diff --git a/test/test_show.jl b/test/test_show.jl index df7da98e..83266f67 100644 --- a/test/test_show.jl +++ b/test/test_show.jl @@ -61,4 +61,6 @@ w = EP.wannier90("wannier_test_tb.dat"); @test nothing === show(stdout, w) @test nothing === show(stdout, position(w)) + M = LP.linear() |> hopping(1) |> greenfunction |> meanfield + @test nothing === show(stdout, M) end