Skip to content

Commit

Permalink
WIP
Browse files Browse the repository at this point in the history
  • Loading branch information
@pablosanjose
pablosanjose committed Dec 3, 2024
1 parent d586cb0 commit 0cf324a
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 72 deletions.
62 changes: 32 additions & 30 deletions src/docstrings.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2052,13 +2052,11 @@ true if `gs` contains user-defined contacts created using `attach(model)` with g
densitymatrix(gs::GreenSlice, ωscale::Real; opts..., quadgk_opts...)
As above with `path = Paths.vertical(ωscale)`, which computes the density matrix by
integrating the Green function along a quarter circle from `ω = -Inf` to `ω = µ + im*Inf`,
and then back along a vertical path `ω = µ + im*ωscale*x` from `x = Inf` to `x = 0`. Here
`µ` is the chemical potential and `ωscale` should be some typical energy scale in the system
(the integration time may depend on `ωscale`, but not the result). This vertical integration
path approach is performant but, currently it does not support finite temperatures, unlike
other `AbstractIntegrationPath`s like `Paths.sawtooth`.
As above with `path = Paths.radial(ωscale, π/4)`, which computes the density matrix by
integrating the Green function from `ω = -Inf` to `ω = Inf` along a `ϕ = π/4` radial complex
path, see `Paths` for details. Here `ωscale` should correspond to some typical energy scale
in the system, which dictates the speed at which we integrate the radial paths (the
integration runtime may depend on `ωscale`, but not the result).
## Full evaluation
Expand All @@ -2068,9 +2066,6 @@ Evaluate the density matrix at chemical potential `μ` and temperature `kBT` (in
units as the Hamiltonian) for the given `g` parameters `params`, if any. The result is given
as an `OrbitalSliceMatrix`, see its docstring for further details.
If the generic integration algorithm is used with complex `ωpoints`, the following form is
also available:
## Algorithms and keywords
The generic integration algorithm allows for the following `opts` (see also `josephson`):
Expand Down Expand Up @@ -2155,21 +2150,21 @@ julia> glead = LP.square() |> hamiltonian(@onsite((; ω = 0) -> 0.0005 * SA[0 1;
julia> g0 = LP.square() |> hamiltonian(hopping(SA[1 0; 0 -1]), orbitals = 2) |> supercell(region = r->-2<r[2]<2 && r[1]≈0) |> attach(glead, reverse = true) |> attach(glead) |> greenfunction;
julia> J = josephson(g0[1], 4; omegamap = ω -> (;ω), phases = subdiv(0, pi, 10))
julia> J = josephson(g0[1], 0.0005; omegamap = ω -> (;ω), phases = subdiv(0, pi, 10))
Josephson: equilibrium Josephson current at a specific contact using solver of type JosephsonIntegratorSolver
julia> J(0.0)
10-element Vector{Float64}:
7.060440509787806e-18
0.0008178484258721882
0.0016108816082772972
0.002355033150366814
0.0030277117620820513
0.003608482493380227
0.004079679643085058
0.004426918320990192
0.004639358112465513
2.2618383948099795e-12
1.0130103834038537e-15
0.0008178802022977883
0.0016109471548779466
0.0023551370133118215
0.0030278625151304614
0.003608696305848759
0.00407998998248311
0.0044274100715435295
0.004640372460465891
5.179773035314182e-12
```
# See also
Expand All @@ -2194,37 +2189,44 @@ Like above for the density matrix `ρ(mu, kBT)`, with `d::DensityMatrixIntegrand
integrand

"""
Quantica.path(O::Josephson, args...)
Quantica.path(O::DensityMatrix, args...)
Quantica.points(O::Josephson, args...)
Quantica.points(O::DensityMatrix, args...)
Return the vertices of the polygonal integration path used to compute `O(args...)`.
Return the vertices of the integration path used to compute `O(args...)`.
"""
path
points

"""
Paths
A Quantica submodule that contains representations of different integration paths in the
complex-ω plane. Available paths are:
complex-ω plane for integrals of the form `∫f(ω)g(ω)dω`, where `f(ω)` is the Fermi
distribution at chemical potential `µ` and temperature `kBT`. Available paths are:
Paths.vertical(ωscale::Real)
Paths.radial(ωscale::Real, ϕ)
A vertical path from a point `µ` on the real axis to `µ + Inf*im`. `ωscale` defines the rate at which the integration variable `x` sweeps this ω path, i.e. `ω = µ+im*ωscale*x`.
A four-segment path from `ω = -Inf` to `ω = Inf`. The path first ascends along an arc of
angle `ϕ` (with `0 <= ϕ < π/2`) and infinite radius. It then converges along a straight line
in the second `ω-µ` quadrant to the chemical potential origin `ω = µ`. Finally it performs
the mirror-symmetric itinerary in the first `ω-µ` quadrant, ending on the real axis at `ω =
Inf`. The radial segments are traversed at a rate dictated by `ωscale`, that should
represent some relevant energy scale in the system for optimal convergence of integrals. At
zero temperature `kBT = 0`, the two last segments don't contribute and are elided.
Paths.sawtooth(ωpoints; slope = 1)
A path connecting all points in the `ωpoints` collection (should all be real), but departing
between each two into the complex plane and back with the given `slope`, forming a sawtooth
path. Note that an extra point `µ` is added to the collection, where `µ` is the chemical
potential. This path allows finite temperature integrals
potential, and the `real(ω)>µ` segments are elided at zero temperatures.
Paths.sawtooth(ωmax::Real; slope = 1)
As above with `ωpoints = (-ωmax, ωmax)`.
Paths.polygon(ωpoints)
A ageneral polygonal path connecting `ωpoints`, which can be any collection of real or
A general polygonal path connecting `ωpoints`, which can be any collection of real or
complex numbers
Paths.polygon(ωfunc::Function)
Expand Down
61 changes: 41 additions & 20 deletions src/integrator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,39 +16,43 @@ end

# straight from 0 to im*∞, but using a real variable from 0 to inf as parameter
# this is because QuadGK doesn't understand complex infinities
struct VerticalPath{T<:AbstractFloat} <: AbstractIntegrationPath
struct RadialPath{T<:AbstractFloat} <: AbstractIntegrationPath
rate::T
angle::T
cisinf::Complex{T}
pts::Vector{Complex{T}}
end

#region ## Constructors ##

SawtoothPath(pts::Vector{T}, slope) where {T<:Real} =
SawtoothPath(pts, complex.(pts), T(slope))

RadialPath(rate, angle) = RadialPath(promote(float(rate), float(angle))...)

function RadialPath(rate::T, angle::T) where {T<:AbstractFloat}
0 <= angle < π/2 || argerror("The radial angle should be in the interval [0, π/2), got $angle")
RadialPath(rate, angle, cis(angle), Complex{T}[])
end

#endregion

#region ## API ##

points(p::PolygonPath, args...) = p.pts
points(p::PolygonPath{<:Function}, mu, kBT) = p.pts(mu, kBT)

points(::VerticalPath{T}, mu, kBT) where {T} = (mu + zero(T)*im, mu + T(Inf)*im)

realpoints(::VerticalPath{T}, pts) where {T} = (zero(T), T(Inf))
realpoints(::Union{SawtoothPath,PolygonPath}, pts) = 1:length(pts)

point(x, p::VerticalPath, pts) = first(pts) + p.rate*im*x

function point(x, p::Union{SawtoothPath,PolygonPath}, pts)
x0 = floor(Int, x)
p0, p1 = pts[x0], pts[ceil(Int, x)]
return p0 + (x-x0)*(p1 - p0)
function points(p::RadialPath{T}, mu, kBT) where {T}
if iszero(kBT)
resize!(p.pts, 2)
copyto!(p.pts, (mu - T(Inf)*conj(p.cisinf), mu + zero(T)*im))
else
resize!(p.pts, 3)
copyto!(p.pts, (mu - T(Inf)*conj(p.cisinf), mu + zero(T)*im, mu + T(Inf)*p.cisinf))
end
return p.pts
end

# the minus sign inverts the VerticalPath, so it is transformed to run from Inf to 0
jacobian(x, p::VerticalPath, pts) = -p.rate*im
jacobian(x, p::Union{SawtoothPath,PolygonPath}, pts) = pts[ceil(Int, x)] - pts[floor(Int, x)]

function points(p::SawtoothPath, mu, kBT)
copy!(p.pts, p.realpts)
if !((iszero(kBT) && maximum(real, p.pts) <= mu) || any(≈(mu), p.pts))
Expand All @@ -59,6 +63,23 @@ function points(p::SawtoothPath, mu, kBT)
return p.pts
end

# Must be a type-stable tuple, and avoiding 0 in RadialPath
realpoints(p::RadialPath{T}, pts) where {T} = (-T(Inf), T(0), ifelse(length(pts)==2, T(0), T(Inf)))
realpoints(::Union{SawtoothPath,PolygonPath}, pts) = 1:length(pts)

# note: pts[2] == µ
point(x, p::RadialPath, pts) = pts[2] + p.rate*x*ifelse(x <= 0, conj(p.cisinf), p.cisinf)

function point(x, p::Union{SawtoothPath,PolygonPath}, pts)
x0 = floor(Int, x)
p0, p1 = pts[x0], pts[ceil(Int, x)]
return p0 + (x-x0)*(p1 - p0)
end

# the minus sign inverts the RadialPath, so it is transformed to run from Inf to 0
jacobian(x, p::RadialPath, pts) = p.rate*ifelse(x <= 0, conj(p.cisinf), p.cisinf)
jacobian(x, p::Union{SawtoothPath,PolygonPath}, pts) = pts[ceil(Int, x)] - pts[floor(Int, x)]

function triangular_sawtooth!(pts, slope)
for i in 2:length(pts)
push!(pts, 0.5 * (pts[i-1] + pts[i]) + im * slope * 0.5 * (pts[i] - pts[i-1]))
Expand All @@ -78,7 +99,7 @@ end

module Paths

using Quantica: Quantica, PolygonPath, SawtoothPath, VerticalPath
using Quantica: Quantica, PolygonPath, SawtoothPath, RadialPath

## API

Expand All @@ -92,7 +113,7 @@ sawtooth!(pts::AbstractVector{<:AbstractFloat}; slope = 1) =

sawtooth!(pts; kw...) = argerror("Path.sawtooth expects a collection of real points, got $pts")

vertical(rate) = VerticalPath(float(rate))
radial(rate, angle) = RadialPath(rate, angle)

polygon(ωmax::Real) = polygon((-ωmax, ωmax))
polygon(ωpoints) = PolygonPath(ωpoints)
Expand Down Expand Up @@ -148,11 +169,11 @@ function call!(I::Integrator{<:Any,Missing}; params...)
end

# nonscalar version
function call!(I::Integrator; params...)
function call!(I::Integrator{<:Any,T}; params...) where {T}
fx! = (y, x) -> begin
y .= serialize(call!(I.integrand, x; params...))
I.callback(x, y)
return y
return nothing
end
result, err = quadgk!(fx!, serialize(I.result), points(I)...; I.quadgk_opts...)
# note: post-processing is not element-wise & can be in-place
Expand Down
17 changes: 10 additions & 7 deletions src/observables.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,7 +392,7 @@ densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) =

# default integrator solver
densitymatrix(gs::GreenSlice, ωmax::Number; kw...) =
densitymatrix(gs::GreenSlice, Paths.vertical(ωmax); kw...)
densitymatrix(gs::GreenSlice, Paths.radial(ωmax, π/4); kw...)

function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegamap = Returns((;)), atol = 1e-7, opts...) where {T}
result = copy(call!_output(gs))
Expand All @@ -407,9 +407,9 @@ function densitymatrix(gs::GreenSlice{T}, path::AbstractIntegrationPath; omegama
return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs)
end

# we need to add the quarter-circle path arc from -∞ to +im*∞
function post_transform_rho(::VerticalPath, gs)
arc = gs(0.5*I)
# we need to add the arc path segment from -∞ to ∞ * p.cisinf
function post_transform_rho(p::RadialPath, gs)
arc = gs((p.angle/π)*I)
function post(x)
x .+= arc
return x
Expand Down Expand Up @@ -437,9 +437,10 @@ post_transform_rho(::AbstractIntegrationPath, _) = identity

function call!(ρi::DensityMatrixIntegrand, x; params...)
ω = point(x, ρi.path, ρi.pts)
symmetrize = -jacobian(x, ρi.path, ρi.pts)/(2π*im)
j = jacobian(x, ρi.path, ρi.pts)
f = fermi(chopsmall(ω - ρi.mu), inv(ρi.kBT))
symmetrize = -j*f/(2π*im)
output = call!(ρi.gs, ω; symmetrize, ρi.omegamap(ω)..., params...)
serialize(output) .*= fermi(ω - ρi.mu, inv(ρi.kBT))
return output
end

Expand All @@ -449,7 +450,9 @@ end

integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = integrand(ρ.solver(mu, kBT))

path(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = path(ρ.solver(mu, kBT))
points(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = points(ρ.solver(mu, kBT))

points(d::DensityMatrixIntegrand) = d.pts

temperature(D::DensityMatrixIntegrand) = D.kBT

Expand Down
30 changes: 15 additions & 15 deletions src/show.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -428,26 +428,26 @@ Base.summary(::Transmission) =

#endregion

############################################################################################
# Integrator
#region
# ############################################################################################
# # Integrator
# #region

function Base.show(io::IO, I::Integrator)
i = get(io, :indent, "")
print(io, i, summary(I), "\n",
"$i Integration path : $(path(I))
$i Integration options : $(display_namedtuple(options(I)))
$i Integrand : ")
ioindent = IOContext(io, :indent => i * " ")
show(ioindent, integrand(I))
end
# function Base.show(io::IO, I::Integrator)
# i = get(io, :indent, "")
# print(io, i, summary(I), "\n",
# "$i Integration path : $(path(I))
# $i Integration options : $(display_namedtuple(options(I)))
# $i Integrand : ")
# ioindent = IOContext(io, :indent => i * " ")
# show(ioindent, integrand(I))
# end

Base.summary(::Integrator) = "Integrator: Complex-plane integrator"
# Base.summary(::Integrator) = "Integrator: Complex-plane integrator"


display_namedtuple(nt::NamedTuple) = isempty(nt) ? "()" : "$nt"
# display_namedtuple(nt::NamedTuple) = isempty(nt) ? "()" : "$nt"

#endregion
# #endregion

############################################################################################
# Josephson and DensityMatrix integrands
Expand Down

0 comments on commit 0cf324a

Please sign in to comment.