From 0c4836533ba75ae1c567fdfb3a3c9981601b538d Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 25 Oct 2024 20:18:37 +0200 Subject: [PATCH] ZeroField zerofield fix --- docs/src/advanced/meanfield.md | 18 ++++++++++-------- src/Quantica.jl | 2 +- src/docstrings.jl | 32 +++++++++++++++++++++++--------- src/meanfield.jl | 23 +++++++++++++++++------ src/specialmatrices.jl | 3 --- src/tools.jl | 1 + test/test_hamiltonian.jl | 2 +- 7 files changed, 53 insertions(+), 28 deletions(-) diff --git a/docs/src/advanced/meanfield.md b/docs/src/advanced/meanfield.md index 5d550a16..2402976d 100644 --- a/docs/src/advanced/meanfield.md +++ b/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; Φ = zerofield) --> Φ[i]); -julia> model_2 = @hopping((i, j; Φ = missing) -> Φ[i, j]);Fock +julia> model_2 = @hopping((i, j; Φ = zerofield) -> Φ[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 `Φ = zerofield` is an singleton 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 = zerofield) --> phi[i]); julia> M = meanfield(greenfunction(h); onsite = 1, selector = (; range = 0), fock = nothing) MeanField{ComplexF64} : builder of Hartree-Fock mean fields @@ -125,13 +125,15 @@ 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 = zerofield) = 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); @@ -139,7 +141,7 @@ 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 diff --git a/src/Quantica.jl b/src/Quantica.jl index ca67908b..70794573 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -33,7 +33,7 @@ export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian, conductance, josephson, ldos, current, transmission, densitymatrix, OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes, siteindexdict, serializer, serialize, serialize!, deserialize!, deserialize, - meanfield + meanfield, zerofield export LatticePresets, LP, RegionPresets, RP, HamiltonianPresets, HP, ExternalPresets, EP export EigenSolvers, ES, GreenSolvers, GS diff --git a/src/docstrings.jl b/src/docstrings.jl index 700c2817..f7749742 100644 --- a/src/docstrings.jl +++ b/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. + zerofield + +An sigleton of type `ZeroField` that represents a zero-values field. It has the property +that it returns zero no matter how it is indexed (`zerofield[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 = zerofield) --> 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 diff --git a/src/meanfield.jl b/src/meanfield.jl index 03cd4513..319cb794 100644 --- a/src/meanfield.jl +++ b/src/meanfield.jl @@ -17,6 +17,8 @@ struct MeanField{B,T,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix} onsite_tmp::Vector{Complex{T}} end +struct ZeroField 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 + +## ZeroField + +const zerofield = ZeroField() + +(m::ZeroField)(args...; kw...) = m + +Base.getindex(::ZeroField, _...) = 0.0 * I + #endregion #endregion diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index b13874f9..318467b8 100644 --- a/src/specialmatrices.jl +++ b/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)) diff --git a/src/tools.jl b/src/tools.jl index 748e37de..4adb70d4 100644 --- a/src/tools.jl +++ b/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}} diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 511f78c8..e7dd6049 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -158,7 +158,7 @@ end @test_throws ArgumentError h() # out of bounds indexing - h = LP.linear() |> @hopping((i,j; phi = missing) --> phi[i, j]); + h = LP.linear() |> @hopping((i,j; phi = zerofield) --> phi[i, j]); ϕ = h[cells = 0] @test h(phi = ϕ) isa Hamiltonian end