Skip to content

Commit

Permalink
Merge branch 'nambu3' into nambu2
Browse files Browse the repository at this point in the history
  • Loading branch information
pablosanjose committed Nov 21, 2024
2 parents 7ca44cd + 1b01340 commit a170a74
Show file tree
Hide file tree
Showing 4 changed files with 29 additions and 25 deletions.
2 changes: 1 addition & 1 deletion src/docstrings.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2086,7 +2086,7 @@ The `quadgk_opts` are extra keyword arguments (other than `atol`) to pass on to
Currently, the following GreenSolvers implement dedicated densitymatrix algorithms:
- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries are not currently supported. No `opts`.
- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries and non-empty contacts are not currently supported. Assumes Hermitian Hamiltonian. No `opts`.
- `GS.Spectrum`: based on summation occupation-weigthed eigenvectors. No `opts`.
- `GS.KPM`: based on the Chebyshev expansion of the Fermi function. Currently only works for zero temperature and only supports `nothing` contacts (see `attach`). No `opts`.
Expand Down
43 changes: 22 additions & 21 deletions src/meanfield.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,17 +7,16 @@
# we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry
#region

struct MeanField{B,T,C<:CompressedOrbitalMatrix,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix}
struct MeanField{B,T,C<:CompressedOrbitalMatrix,S<:SparseMatrixCSC,F<:DensityMatrix}
output::C
potHartree::S
potFock::S
rhoHartree::H
rhoFock::F
rho::F
charge::B
nambu::Bool
is_nambu_rotated::Bool
rowcol_ranges::NTuple{2,Vector{UnitRange{Int}}}
onsite_tmp::Vector{Complex{T}}
onsite_tmp::Vector{T}
end

struct ZeroField end
Expand All @@ -41,38 +40,36 @@ function meanfield(g::GreenFunction{T,E}, args...;
isempty(boundaries(g)) || argerror("meanfield does not currently support systems with boundaries")
isfinite(U) || argerror("Onsite potential must be finite, consider setting `onsite`")

gsHartree = g[diagonal(; cells = 0, kernel = Q)]
rhoHartree = densitymatrix(gsHartree, args...; kw...)
gsFock = g[sitepairs(; selector..., includeonsite = true)]
rhoFock = densitymatrix(gsFock, args...; kw...)
gs = g[sitepairs(; selector..., includeonsite = true)]
rho = densitymatrix(gs, 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 : T(Vf(dr)); selector..., includeonsite = true)
hHartree = (Uf == U && Vh === Vf) ? hFock :
lat |> hopping((r, dr) -> iszero(dr) ? U : T(Vh(dr)); selector..., includeonsite = true)

potHartree = sum(unflat, harmonics(hHartree))
potHartree = T.(sum(unflat, harmonics(hHartree)))

oaxes = orbaxes(call!_output(gsFock))
oaxes = orbaxes(call!_output(gs))
rowcol_ranges = collect.(orbranges.(oaxes))
onsite_tmp = similar(diag(parent(call!_output(gsHartree))))
onsite_tmp = Vector{T}(undef, length(last(rowcol_ranges)))

# build potFock with identical axes as the output of rhoFock
# build potFock with identical axes as the output of rho
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)
check_cell_order(hFock_slice, rho)
potFock = T.(parent(hFock_slice))

encoder, decoder = nambu ? NambuEncoderDecoder(is_nambu_rotated´) : (identity, identity)
S = typeof(encoder(zero(Q)))
output = call!_output(g[sitepairs(; selector..., includeonsite = true, kernel = Q)])
sparse_enc = similar(output, S)
output = CompressedOrbitalMatrix(sparse_enc; encoder, decoder, hermitian = true)

return MeanField(output, potHartree, potFock, rhoHartree, rhoFock, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp)
return MeanField(output, potHartree, potFock, rho, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp)
end

sanitize_potential(x::Real) = Returns(x)
Expand All @@ -91,6 +88,7 @@ sanitize_nambu_rotated(is_nambu_rotated::Bool) = is_nambu_rotated

function sanitize_charge(charge, B, nambu, is_rotated)
Q = sanitize_charge(charge, B)
ishermitian(Q) || argerror("Charge $Q should be Hermitian")
nambu && check_nambu(Q, is_rotated)
return Q
end
Expand Down Expand Up @@ -137,9 +135,9 @@ function NambuEncoderDecoder(is_nambu_rotated)
return encoder, decoder
end

function check_cell_order(hFock_slice, rhoFock)
function check_cell_order(hFock_slice, rho)
opot = first(orbaxes(hFock_slice))
orho = first(orbaxes(call!_output(rhoFock.gs)))
orho = first(orbaxes(call!_output(rho.gs)))
cells(opot) == cells(orho) || internalerror("meanfield: Cell order mismatch between potential and density matrix")
return nothing
end
Expand All @@ -164,15 +162,15 @@ function call!(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
check_zero_mu(m, args...)
trρQ = maybe_nambufy_traces!(diag(parent(m.rhoHartree(args...; params...))), m)
ρ = m.rho(args...; params...)
trρQ = maybe_nambufy_traces!(diag_real_tr_rho_Q(ρ, Q), m)
mul!(hartree_pot, m.potHartree, trρQ)
ρFock = m.rhoFock(args...; params...)
meanfield = m.output
meanfield_parent = parent(meanfield)
fill!(nonzeros(meanfield_parent), zero(eltype(meanfield_parent)))
if chopsmall
hartree_pot .= Quantica.chopsmall.(hartree_pot) # this is a Vector
nzs = nonzeros(parent(ρFock))
nzs = nonzeros(parent(ρ))
nzs .= Quantica.chopsmall.(nzs)
end
rows, cols, nzs = rowvals(fock_pot), axes(fock_pot, 2), nonzeros(fock_pot)
Expand All @@ -183,7 +181,7 @@ function call!(m::MeanField{B}, args...; chopsmall = true, params...) where {B}
row = rows[ptr]
row < col && ishermitian(meanfield) && continue # skip upper triangle
vij = nzs[ptr]
ρij = view(ρFock, rowrngs[row], colrngs[col])
ρij = view(ρ, rowrngs[row], colrngs[col])
vQρijQ = vij * Q * sanitize_block(B, ρij) * Q
if row == col
meanfield_parent[row, col] = encoder(meanfield)(viiQ - vQρijQ)
Expand All @@ -209,6 +207,9 @@ function maybe_nambufy_traces!(traces, m::MeanField{B}) where {B}
return traces
end

diag_real_tr_rho_Q(ρ, Q) =
[real(unsafe_trace_prod(Q, view(ρ, rng, rng))) for rng in siteindexdict(orbaxes(ρ, 2))]

hole_id(::Type{<:SMatrix{2,2}}) = SA[0 0; 0 1]
hole_id(::Type{<:SMatrix{4,4}}) = SA[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]
hole_id(S) = nambu_dim_error(S)
Expand Down
3 changes: 2 additions & 1 deletion src/solvers/green/schur.jl
Original file line number Diff line number Diff line change
Expand Up @@ -689,6 +689,7 @@ end

# use computed Fermi points to integrate spectral function by segments
# returns an AbstractMatrix
# we don't use Integrator, because that is meant for integrals over energy, not momentum
function (s::DensityMatrixSchurSolver)(µ, kBT; params...)
g = parent(s.gs)
λs = propagating_eigvals(g, µ + 0im, 1e-2; params...)
Expand All @@ -709,7 +710,7 @@ function fermi_h!(result, s, ϕ, µ, β = 0; params...)
bs = blockstructure(h)
# Similar to spectrum(h, ϕ; params...), but less work (no sort! or sanitization)
copy!(s.hmat, call!(h, ϕ; params...)) # sparse to dense
ϵs, psis = eigen!(s.hmat)
ϵs, psis = eigen!(Hermitian(s.hmat))
# special-casing β = Inf with views turns out to be slower
fs = (@. ϵs = fermi(ϵs - µ, β))
fpsis = (s.psis .= psis .* transpose(fs))
Expand Down
6 changes: 4 additions & 2 deletions test/test_greenfunction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -566,10 +566,12 @@ 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 1; im 0] # nonhermitian charge not allowed
@test_throws ArgumentError meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q)
Q = SA[0 -im; im 0]
m = meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q)
Φ = m(0.2, 0.3)
ρ12 = m.rhoFock(0.2, 0.3)[sites(1), sites(2)]
ρ12 = m.rho(0.2, 0.3)[sites(1), sites(2)]
@test Φ[sites(1), sites(2)] -1.5 * Q * ρ12 * Q

# nambu
Expand All @@ -579,7 +581,7 @@ end
m = meanfield(g; selector = (; range = 1), nambu = true, hartree = r -> 1/(1+norm(r)), fock = 1.5, charge = Q)
@test_throws ArgumentError m(0.2, 0.3) # µ cannot be nonzero
Φ = m(0, 0.3)
ρ11 = m.rhoFock(0, 0.3)[sites(1), sites(1)] - SA[0 0; 0 1]
ρ11 = m.rho(0, 0.3)[sites(1), sites(1)] - SA[0 0; 0 1]
fock = -1.5 * Q * ρ11 * Q
hartree = 0.5*Q*(tr(Q*ρ11)*(1+1/2+1/2))
@test Φ[sites(1), sites(1)] hartree + fock
Expand Down

0 comments on commit a170a74

Please sign in to comment.