ZeroMeanField

docs/src/advanced/meanfield.md

@@ -10,14 +10,14 @@ Schematically the process is as follows:
 ```julia
 julia> model_0 = hopping(1); # a possible non-interacting model
 
-julia> model_1 = @onsite((i; Φ = missing) --> Φ[i]);
+julia> model_1 = @onsite((i; Φ = meanfield()) --> Φ[i]);
 
-julia> model_2 = @hopping((i, j; Φ = missing) -> Φ[i, j]);Fock
+julia> model_2 = @hopping((i, j; Φ = meanfield()) -> Φ[i, j]);Fock
 
 julia> h = lat |> hamiltonian(model_0 + model_1 + model_2)
 ```
 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.
+The default value `Φ = meanfield()` is a `ZeroMeanField` object that represents no interactions, so `model_1` and `model_2` vanish by default.
 - 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...)`
 
@@ -74,7 +74,7 @@ With the serializer functionality we can build a version of the fixed-point func
 ```julia
 julia> using SIAMFANLEquations
 
-julia> h = LP.linear() |> supercell(4) |> onsite(r->r[1]) - hopping(1) + @onsite((i; phi = missing) --> phi[i]);
+julia> h = LP.linear() |> supercell(4) |> onsite(r->r[1]) - hopping(1) + @onsite((i; phi = meanfield()) --> phi[i]);
 
 julia> M = meanfield(greenfunction(h); onsite = 1, selector = (; range = 0), fock = nothing)
 MeanField{ComplexF64} : builder of Hartree-Fock mean fields
@@ -125,21 +125,23 @@ Note that the content of `pdata` is passed by `aasol` as a third argument to `f!
 
 ## 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.
+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. This is one way to do it.
+
 ```julia
 julia> using SIAMFANLEquations
 
-julia> h = LP.linear() |> supercell(4) |> supercell |> onsite(1) - hopping(1) + @onsite((i; phi = missing) --> phi[i]);
+julia> h = LP.linear() |> supercell(4) |> supercell |> onsite(1) - hopping(1) + @onsite((i; phi) --> phi[i]);
 
-julia> M´(Φ = missing) = meanfield(greenfunction(h(; phi = Φ), GS.Spectrum()); onsite = 1, selector = (; range = 0), fock = nothing)
+julia> M´(phi = meanfield()) = meanfield(greenfunction(h(; phi), GS.Spectrum()); onsite = 1, selector = (; range = 0), fock = nothing)
+M´ (generic function with 3 methods)
 
 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(Φ))
         return x
-    end;
+           end;
 
 julia> m = 2; x0 = serialize(Φ0); vstore = rand(length(x0), 3m+3);  # order m, initial condition x0, and preallocated space vstore
 

src/docstrings.jl

@@ -2602,6 +2602,12 @@ interaction potentials, and `ρ` is the density matrix evaluated at specific che
 potential and temperature. Also `ν = ifelse(nambu, 0.5, 1.0)`, and `v_F(0) = v_H(0) = U`,
 where `U` is the onsite interaction.
 
+    meanfield()
+
+Build an `M::ZeroMeanField` object that represents a zero-values mean field. It has the
+property that it returns zero no matter how it is indexed (`M[inds...] = 0.0 * I`), so it is
+useful as a default value in a non-spatial model involving mean fields.
+
 ## Keywords
 
 - `potential`: charge-charge potential to use for both Hartree and Fock. Can be a number or a function of position. Default: `1
@@ -2625,26 +2631,34 @@ e.g. `ϕ[sites(2:3), sites(1)]`, see `OrbitalSliceArray` for details.
 # Examples
 
 ```jldoctest
-julia> g = HP.graphene(orbitals = 2, a0 = 1) |> supercell((1,-1)) |> greenfunction;
+julia> model = hopping(I) - @onsite((i; phi = meanfield()) --> phi[i]);
+
+julia> g = LP.honeycomb() |> hamiltonian(model, orbitals = 2) |> supercell((1,-1)) |> greenfunction;
 
-julia> M = meanfield(g; selector = (; range = 1), charge = I)
+julia> M = meanfield(g; selector = (; range = 1), charge = I, potential = 0.05)
 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> phi = M(0.2, 0.3);
+julia> phi0 = M(0.2, 0.3);
 
-julia> phi[sites(1), sites(2)] |> Quantica.chopsmall
+julia> phi0[sites(1), sites(2)] |> Quantica.chopsmall
 2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
- 0.00239416+0.0im             ⋅
-            ⋅      0.00239416+0.0im
+ 0.00109527+0.0im             ⋅
+            ⋅      0.00109527+0.0im
 
-julia> phi[sites(1)] |> Quantica.chopsmall
+julia> phi0[sites(1)] |> Quantica.chopsmall
 2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
- 5.53838+0.0im          ⋅
-         ⋅      5.53838+0.0im
+ 0.296672+0.0im           ⋅
+          ⋅      0.296672+0.0im
 
+julia> phi1 = M(0.2, 0.3; phi = phi0);
+
+julia> phi1[sites(1), sites(2)] |> Quantica.chopsmall
+2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
+ 0.00307712+0.0im             ⋅
+            ⋅      0.00307712+0.0im
 ```
 """
 meanfield

src/meanfield.jl

@@ -17,6 +17,8 @@ struct MeanField{B,T,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix}
     onsite_tmp::Vector{Complex{T}}
 end
 
+struct ZeroMeanField end
+
 #region ## Constructors ##
 
 function meanfield(g::GreenFunction{T,E}, args...;
@@ -63,10 +65,6 @@ function meanfield(g::GreenFunction{T,E}, args...;
     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)
@@ -92,13 +90,17 @@ hartree_matrix(m::MeanField) = m.potHartree
 
 fock_matrix(m::MeanField) = parent(m.potFock)
 
-function (m::MeanField{B})(args...; params...) where {B}
+function (m::MeanField{B})(args...; chopsmall = true, 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
+    mf_parent = parent(meanfield)
+    if chopsmall
+        hartree_pot .= Quantica.chopsmall.(hartree_pot)
+        nonzeros(mf_parent) .= Quantica.chopsmall.(nonzeros(mf_parent))
+    end
     rows, cols, nzs = rowvals(fock_pot), axes(fock_pot, 2), nonzeros(fock_pot)
     for col in cols
         viiQ = hartree_pot[col] * Q
@@ -117,6 +119,15 @@ function (m::MeanField{B})(args...; params...) where {B}
     return meanfield
 end
 
+
+## ZeroMeanField
+
+meanfield() = ZeroMeanField()
+
+(m::ZeroMeanField)(args...; kw...) = m
+
+Base.getindex(::ZeroMeanField, _...) = 0.0 * I
+
 #endregion
 
 #endregion

src/specialmatrices.jl

@@ -549,9 +549,6 @@ Base.getindex(a::OrbitalSliceMatrix, i::C, j::C = i) where {B,C<:CellSitePos{<:A
 Base.checkbounds(::Type{Bool}, a::OrbitalSliceMatrix, i::CellSitePos, j::CellSitePos) =
     checkbounds(Bool, first(orbaxes(a)), i) && checkbounds(Bool, last(orbaxes(a)), 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))

src/tools.jl

@@ -69,6 +69,7 @@ chopsmall(x::C, atol = sqrt(eps(real(C)))) where {C<:Complex} =
     chopsmall(real(x), atol) + im*chopsmall(imag(x), atol)
 chopsmall(xs, atol) = chopsmall.(xs, Ref(atol))
 chopsmall(xs) = chopsmall.(xs)
+chopsmall(xs::UniformScaling, atol...) = I * chopsmall(xs.λ, atol...)
 
 # Flattens matrix of Matrix{<:Number} into a matrix of Number's
 function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:Matrix{C}}