Skip to content

Commit

Permalink
WIP
Browse files Browse the repository at this point in the history 
fixing tests

fixing tests

fixing tests

test fixes

add and fix more tests

more tests

docs

coverage

fix tests

add tests

revert inhomogeneous slice behavior

fix for CellOrbitals(:) slicing

don't allow CellOrbitals{L,Colon} reach slicers

comment

test fix

generalize API
  • Loading branch information
@pablosanjose
pablosanjose committed Dec 5, 2024
1 parent d586cb0 commit dc72b2e
Show file tree
Hide file tree
Showing 14 changed files with 475 additions and 315 deletions.
93 changes: 42 additions & 51 deletions docs/src/tutorial/observables.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,57 +38,61 @@ Note that `d[sites...]` produces a vector with the LDOS at sites defined by `sit
We can also compute the convolution of the density of states with the Fermi distribution `f(ω)=1/(exp((ω-μ)/kBT) + 1)`, which yields the density matrix in thermal equilibrium, at a given temperature `kBT` and chemical potential `μ`. This is computed with `ρ = densitymatrix(gs, (ωmin, ωmax))`. Here `gs = g[sites...]` is a `GreenSlice`, and `(ωmin, ωmax)` are integration bounds (they should span the full bandwidth of the system). Then, `ρ(µ, kBT = 0; params...)` will yield a matrix over the selected `sites` for a set of model `params`.
```julia
julia> ρ = densitymatrix(g[region = RP.circle(1)], (-0.1, 8.1))
DensityMatrix: density matrix on specified sites using solver of type DensityMatrixIntegratorSolver
DensityMatrix{DensityMatrixIntegratorSolver}: density matrix on specified sites

julia> @time ρ(4)
6.150548 seconds (57.84 k allocations: 5.670 GiB, 1.12% gc time)
5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}:
0.5+0.0im -7.34893e-10-3.94035e-15im 0.204478+1.9366e-14im -7.34889e-10-1.44892e-15im -5.70089e-10+5.48867e-15im
-7.34893e-10+3.94035e-15im 0.5+0.0im 0.200693-2.6646e-14im -5.70089e-10-1.95251e-15im -7.34891e-10-2.13804e-15im
0.204478-1.9366e-14im 0.200693+2.6646e-14im 0.5+0.0im 0.200693+3.55692e-14im 0.204779-4.27255e-14im
-7.34889e-10+1.44892e-15im -5.70089e-10+1.95251e-15im 0.200693-3.55692e-14im 0.5+0.0im -7.34885e-10-3.49861e-15im
-5.70089e-10-5.48867e-15im -7.34891e-10+2.13804e-15im 0.204779+4.27255e-14im -7.34885e-10+3.49861e-15im 0.5+0.0im
4.594645 seconds (111.82 k allocations: 4.890 GiB, 2.81% gc time, 0.86% compilation time)
5×5 OrbitalSliceMatrix{ComplexF64,Array}:
0.5+0.0im -7.35075e-10+3.40256e-15im 0.204478+4.46023e-14im -7.35077e-10-1.13342e-15im -5.70426e-10-2.22213e-15im
-7.35075e-10-3.40256e-15im 0.5+0.0im 0.200693-4.46528e-14im -5.70431e-10+3.53853e-15im -7.35092e-10-4.07992e-16im
0.204478-4.46023e-14im 0.200693+4.46528e-14im 0.5+0.0im 0.200693+6.7793e-14im 0.204779-6.78156e-14im
-7.35077e-10+1.13342e-15im -5.70431e-10-3.53853e-15im 0.200693-6.7793e-14im 0.5+0.0im -7.351e-10+2.04708e-15im
-5.70426e-10+2.22213e-15im -7.35092e-10+4.07992e-16im 0.204779+6.78156e-14im -7.351e-10-2.04708e-15im 0.5+0.0im
```

Note that the diagonal is `0.5`, indicating half-filling.

The default algorithm used here is slow, as it relies on numerical integration in the complex plane. Some GreenSolvers have more efficient implementations. If they exist, they can be accessed by omitting the `(ωmin, ωmax)` argument. For example, using `GS.Spectrum`:
```julia
julia> @time g = h |> greenfunction(GS.Spectrum());
37.638522 seconds (105 allocations: 2.744 GiB, 0.79% gc time)
18.249136 seconds (75 allocations: 1.567 GiB, 0.67% gc time)

julia> ρ = densitymatrix(g[region = RP.circle(1)])
DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver
DensityMatrix{DensityMatrixSpectrumSolver}: density matrix on specified sites

julia> @time ρ(4) # second-run timing
0.029662 seconds (7 allocations: 688 bytes)
5×5 OrbitalSliceMatrix{ComplexF64,Array}:
0.5+0.0im -6.6187e-16+0.0im 0.204478+0.0im 2.49658e-15+0.0im -2.6846e-16+0.0im
-6.6187e-16+0.0im 0.5+0.0im 0.200693+0.0im -2.01174e-15+0.0im 1.2853e-15+0.0im
0.204478+0.0im 0.200693+0.0im 0.5+0.0im 0.200693+0.0im 0.204779+0.0im
2.49658e-15+0.0im -2.01174e-15+0.0im 0.200693+0.0im 0.5+0.0im 1.58804e-15+0.0im
-2.6846e-16+0.0im 1.2853e-15+0.0im 0.204779+0.0im 1.58804e-15+0.0im 0.5+0.0im

julia> @time ρ(4)
0.001659 seconds (9 allocations: 430.906 KiB)
5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}:
0.5+0.0im -2.21437e-15+0.0im 0.204478+0.0im 2.67668e-15+0.0im 3.49438e-16+0.0im
-2.21437e-15+0.0im 0.5+0.0im 0.200693+0.0im -1.40057e-15+0.0im -2.92995e-15+0.0im
0.204478+0.0im 0.200693+0.0im 0.5+0.0im 0.200693+0.0im 0.204779+0.0im
2.67668e-15+0.0im -1.40057e-15+0.0im 0.200693+0.0im 0.5+0.0im 1.81626e-15+0.0im
3.49438e-16+0.0im -2.92995e-15+0.0im 0.204779+0.0im 1.81626e-15+0.0im 0.5+0.0im
```

Note, however, that the computation of `g` is much slower in this case, due to the need of a full diagonalization. A better algorithm choice in this case is `GS.KPM`. It requires, however, that we define the region for the density matrix beforehand, as a `nothing` contact.
```julia
julia> @time g = h |> attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 10000, bandrange = (0,8)));
Computing moments: 100%|█████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01
2.065083 seconds (31.29 k allocations: 11.763 MiB)
Computing moments: 100%|███████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01
1.360412 seconds (51.17 k allocations: 11.710 MiB)

julia> ρ = densitymatrix(g[1])
DensityMatrix: density matrix on specified sites with solver of type DensityMatrixKPMSolver
DensityMatrix{DensityMatrixKPMSolver}: density matrix on specified sites

julia> @time ρ(4)
0.006580 seconds (3 allocations: 1.156 KiB)
5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}:
0.004024 seconds (3 allocations: 688 bytes)
5×5 OrbitalSliceMatrix{ComplexF64,Array}:
0.5+0.0im 2.15097e-17+0.0im 0.20456+0.0im 2.15097e-17+0.0im 3.9251e-17+0.0im
2.15097e-17+0.0im 0.5+0.0im 0.200631+0.0im 1.05873e-16+0.0im 1.70531e-18+0.0im
0.20456+0.0im 0.200631+0.0im 0.5+0.0im 0.200631+0.0im 0.20482+0.0im
2.15097e-17+0.0im 1.05873e-16+0.0im 0.200631+0.0im 0.5+0.0im 1.70531e-18+0.0im
3.9251e-17+0.0im 1.70531e-18+0.0im 0.20482+0.0im 1.70531e-18+0.0im 0.5+0.0im
```

!!! note "Alternative integration paths"
The integration algorithm allows many different integration paths that can be adjusted to each problem, see the `Paths` docstring. Another versatile choice is `Paths.radial(ωrate, ϕ)`. This one is called with `ϕ = π/4` when doing `ρ = densitymatrix(gs::GreenSlice, ωrate::Number)`. In the example above this is slightly faster than the `(ωmin, ωmax)` choice, which resorts to `Paths.sawtooth(ωmin, ωmax)`.

## Current

A similar computation can be done to obtain the current density, using `J = current(g(ω), direction = missing)`. This time `J[sᵢ, sⱼ]` yields a sparse matrix of current densities along a given direction for each hopping (or the current norm if `direction = missing`). Passing `J` as a hopping shader yields the equilibrium current in a system. In the above example we can add a magnetic flux to make this current finite
Expand Down Expand Up @@ -238,14 +242,7 @@ GreenFunction{Float64,2,0}: Green function of a Hamiltonian{Float64,2,0}
Coordination : 3.7448

julia> J = josephson(g[1], 4.1)
Integrator: Complex-plane integrator
Integration path : (-4.1 + 1.4901161193847656e-8im, -2.05 + 2.050000014901161im, 0.0 + 1.4901161193847656e-8im)
Integration options : (atol = 1.0e-7,)
Integrand: :
JosephsonIntegrand{Float64} : Equilibrium (dc) Josephson current observable before integration over energy
kBT : 0.0
Contact : 1
Number of phase shifts : 0
Josephson{JosephsonIntegratorSolver}: equilibrium Josephson current at a specific contact

julia> qplot(g, children = (; sitecolor = :blue))
```
Expand All @@ -257,37 +254,31 @@ In this case we have chosen to introduce the superconducting leads with a model
corresponding to a BCS bulk, but any other self-energy form could be used. We have introduced the phase difference (`phase`) as a model parameter. We can now evaluate the zero-temperature Josephson current simply with
```julia
julia> J(phase = 0)
-1.974396994480587e-16
1.992660837638158e-12

julia> J(phase = 0.2)
0.004617597139699372
0.0046175971391935605
```
Note that finite temperatures can be taken using the `kBT` keyword argument for `josephson`, see docstring for details.

One is often interested in the critical current, which is the maximum of the Josephson current over all phase differences. Quantica.jl can compute the integral over a collection of phase differences simultaneously, which is more efficient that computing them one by one. This is done with
```julia
julia> φs = subdiv(0, pi, 11); J = josephson(g[1], 4.1; phases = φs)
Integration path : (-4.1 + 1.4901161193847656e-8im, -2.05 + 2.050000014901161im, 0.0 + 1.4901161193847656e-8im)
Integration options : (atol = 1.0e-7,)
Integrand: :
JosephsonIntegrand{Float64} : Equilibrium (dc) Josephson current observable before integration over energy
kBT : 0.0
Contact : 1
Number of phase shifts : 11
Josephson{JosephsonIntegratorSolver}: equilibrium Josephson current at a specific contact

julia> Iφ = J()
11-element Vector{Float64}:
1.868862401627357e-14
0.007231421775452674
0.014242855188877
0.02081870760779799
0.026752065104401878
0.031847203848574666
0.0359131410974842
0.03871895510547465
0.039762442694035505
0.03680096751905469
2.7677727119798235e-14
-6.361223111882911e-13
0.007231421776215144
0.01424285518831463
0.020818707606469377
0.026752065101976884
0.031847203846513975
0.035913141096514584
0.038718955102068034
0.03976244268586444
0.036800967573567184
-1.437196514806921e-12

julia> f = Figure(); a = Axis(f[1,1], xlabel = "φ", ylabel = "I [e/h]"); lines!(a, φs, Iφ); scatter!(a, φs, Iφ); f
```
Expand Down
113 changes: 63 additions & 50 deletions src/docstrings.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2050,15 +2050,18 @@ true if `gs` contains user-defined contacts created using `attach(model)` with g
ω-dependent models, but note that Quantica will not check for analyticity. Keywords
`quadgk_opts` are passed to the `QuadGK.quadgk` integration routine.
densitymatrix(gs::GreenSlice, ωscale::Real; opts..., quadgk_opts...)
densitymatrix(gs::GreenSlice, ωscale::Real; kw...)
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).
densitymatrix(gs::GreenSlice, ωs::NTuple{N,Real}; kw...)
As above, but with a `path = Paths.sawtooth(ωs)`, which uses a sawtooth-shaped path touching
points `ωs` on the real axes. Ideally, these values should span the full system bandwidth.
## Full evaluation
Expand All @@ -2068,9 +2071,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 @@ -2122,13 +2122,16 @@ GʳΣʳᵢ)τz]``. Here `f(ω)` is the Fermi function with `µ = 0`.
josephson(gs::GreenSlice, ωscale::Real; kw...)
As above, but with `path = Paths.vertical(ωscale)`, which computes josephson current by
performing the above integral 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 `ωscale`
should be some typical energy scale in the system, like the superconducting gap (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, but with `path = Paths.radial(ωscale, π/4)`, which computes josephson current by
integrating the Green function from `ω = -Inf` to `ω = Inf` along a `ϕ = π/4` radial complex
path, see `Paths` for details. Here `ωscale` should be some typical energy scale in the
system, like the superconducting gap (the integration time may depend on `ωscale`, but not
the result).
josephson(gs::GreenSlice, ωs::NTuple{N,Real}; kw...)
As above, but with a `path = Paths.sawtooth(ωs)`, which uses a sawtooth-shaped path touching
points `ωs` on the real axes. Ideally, these values should span the full system bandwidth.
## Full evaluation
Expand All @@ -2155,21 +2158,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 @@ -2178,59 +2181,69 @@ julia> J(0.0)
josephson

"""
Quantica.integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0)
Quantica.integrand(J::Josephson{<:JosephsonIntegratorSolver}, kBT = 0; params...)
Return the complex integrand `d::JosephsonIntegrand` whose integral over frequency yields the
Josephson current, `J(kBT) = Re(∫dω d(ω))`. To evaluate the `d` for a
given `ω` and parameters, use `d(ω; params...)`, or `call!(d, ω; params...)` for its
Return the complex integrand `d::JosephsonIntegrand` whose integral over frequency yields
the Josephson current, `J(kBT; params...) = Re(∫dx d(x; params...))`, where `ω(x) =
Quantica.point(x, d)` is a path parametrization over real variable `x` . To evaluate the `d`
for a given `x` and parameters, use `d(x; params...)`, or `call!(d, x; params...)` for its
mutating (non-allocating) version.
Quantica.integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0, kBT = 0)
Quantica.integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0, kBT = 0; params...)
Like above for the density matrix `ρ(mu, kBT)`, with `d::DensityMatrixIntegrand`, so that
`ρ(mu, kBT) = -mat_imag(∫dω d(ω))/π`, and `mat_imag(GF::Matrix) = (GF - GF')/2im`.
Like above for the density matrix `ρ`, with `d::DensityMatrixIntegrand`, so that `ρ(mu, kBT;
params...) = -mat_imag(∫dω d(ω; params...))/π`, and `mat_imag(GF::Matrix) = (GF - GF')/2im`.
"""
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)
Paths.sawtooth(ωpoints...; slope = 1, imshift = true)
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. If `imshift = true`
all `ωpoints` (of type `T<:Real`) are shifted by `sqrt(eps(T))` to avoid the real axis.
Paths.sawtooth(ωmax::Real; slope = 1)
Paths.sawtooth(ωmax::Real; kw...)
As above with `ωpoints = (-ωmax, ωmax)`.
Paths.polygon(ωpoints)
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)
A ageneral polygonal path connecting `ωpoints = ωfunc(µ, kBT)`. This is useful when the
desired integration path depends on the chemical potential `µ` or temperature `kBT`.
A ageneral polygonal path connecting `ωpoints = ωfunc(µ, kBT; params...)`. This is useful
when the desired integration path depends on the chemical potential `µ`, temperature
`kBT` or other system parameters `params`.
# See also
`densitymatrix`, `josephson`
Expand Down
Loading

0 comments on commit dc72b2e

Please sign in to comment.