Skip to content

Commit ecec278

Browse files
ytdHuangytdHuang
committed
remove FFTCorrelation and introduce spectrum_correlation_fft 
1 parent 35878b7 commit ecec278

File tree

7 files changed

+231
-179
lines changed

7 files changed

+231
-179
lines changed

‎benchmarks/correlations_and_spectrum.jl‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,18 @@ function benchmark_correlations_and_spectrum!(SUITE)
99
c_ops = [sqrt(γ * (nth + 1)) * a, sqrt(γ * nth) * a']
1010

1111
ω_l = range(0, 3, length = 1000)
12+
t_l = range(0, 333 * π, length = 1000)
13+
14+
function _calculate_fft_spectrum(H, tlist, c_ops, A, B)
15+
corr = correlation_2op_1t(H, nothing, tlist, c_ops, A, B; progress_bar = Val(false))
16+
ωlist, spec = spectrum_correlation_fft(tlist, corr)
17+
nothing
18+
end
1219

1320
SUITE["Correlations and Spectrum"]["FFT Correlation"] =
14-
@benchmarkable spectrum($H, $ω_l, $(a'), $a, $c_ops, solver = FFTCorrelation(), progress_bar = false)
21+
@benchmarkable _calculate_fft_spectrum($H, $t_l, $c_ops, $(a'), $a)
1522

16-
SUITE["Correlations and Spectrum"]["Exponential Series"] = @benchmarkable spectrum($H, $ω_l, $(a'), $a, $c_ops)
23+
SUITE["Correlations and Spectrum"]["Spectrum"]["Exponential Series"] = @benchmarkable spectrum($H, $ω_l, $c_ops, $(a'), $a)
1724

1825
return nothing
1926
end

‎docs/src/resources/api.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -232,7 +232,9 @@ lr_mesolve
232232
correlation_3op_2t
233233
correlation_2op_2t
234234
correlation_2op_1t
235+
spectrum_correlation_fft
235236
spectrum
237+
ExponentialSeries
236238
```
237239

238240
## [Metrics](@id doc-API:Metrics)

‎src/QuantumToolbox.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@ import DiffEqNoiseProcess: RealWienerProcess!
6060
# other dependencies (in alphabetical order)
6161
import ArrayInterface: allowed_getindex, allowed_setindex!
6262
import Distributed: RemoteChannel
63-
import FFTW: fft, fftshift
63+
import FFTW: fft, ifft, fftfreq, fftshift
6464
import Graphs: connected_components, DiGraph
6565
import IncompleteLU: ilu
6666
import Pkg
@@ -107,6 +107,7 @@ include("time_evolution/time_evolution_dynamical.jl")
107107
# Others
108108
include("permutation.jl")
109109
include("correlations.jl")
110+
include("spectrum.jl")
110111
include("wigner.jl")
111112
include("spin_lattice.jl")
112113
include("arnoldi.jl")

‎src/correlations.jl‎

Lines changed: 62 additions & 171 deletions
Original file line numberDiff line numberDiff line change
@@ -1,242 +1,133 @@
1-
export SpectrumSolver, FFTCorrelation, ExponentialSeries
2-
export correlation_3op_2t, correlation_2op_2t, correlation_2op_1t, spectrum
3-
4-
abstract type SpectrumSolver end
5-
6-
struct FFTCorrelation <: SpectrumSolver end
7-
8-
struct ExponentialSeries{T<:Real,CALC_SS} <: SpectrumSolver
9-
tol::T
10-
ExponentialSeries(tol::T, calc_steadystate::Bool = false) where {T} = new{T,calc_steadystate}(tol)
11-
end
12-
13-
ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(tol, calc_steadystate)
1+
export correlation_3op_2t, correlation_2op_2t, correlation_2op_1t
142

153
@doc raw"""
16-
correlation_3op_2t(H::QuantumObject,
17-
ψ0::QuantumObject,
18-
t_l::AbstractVector,
19-
τ_l::AbstractVector,
4+
correlation_3op_2t(H::AbstractQuantumObject,
5+
ψ0::Union{Nothing,QuantumObject},
6+
tlist::AbstractVector,
7+
τlist::AbstractVector,
8+
c_ops::Union{Nothing,AbstractVector,Tuple},
209
A::QuantumObject,
2110
B::QuantumObject,
22-
C::QuantumObject,
23-
c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
11+
C::QuantumObject;
2412
kwargs...)
2513
26-
Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle``
14+
Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
2715
28-
for a given initial state ``\ket{\psi_0}``.
16+
If the initial state `ψ0` is given as `nothing`, then the [`steadystate`] will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
2917
"""
3018
function correlation_3op_2t(
31-
H::QuantumObject{<:AbstractArray{T1},HOpType},
32-
ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
33-
t_l::AbstractVector,
34-
τ_l::AbstractVector,
35-
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
36-
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
37-
C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
38-
c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
19+
H::AbstractQuantumObject{DataType,HOpType},
20+
ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
21+
tlist::AbstractVector,
22+
τlist::AbstractVector,
23+
c_ops::Union{Nothing,AbstractVector,Tuple},
24+
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
25+
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
26+
C::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject};
3927
kwargs...,
4028
) where {
29+
DataType,
4130
T1,
4231
T2,
4332
T3,
4433
T4,
45-
T5,
4634
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
4735
StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
4836
}
49-
allequal((H.dims, ψ0.dims, A.dims, B.dims, C.dims)) || throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
37+
L = liouvillian(H, c_ops)
38+
if ψ0 isa Nothing
39+
ψ0 = steadystate(L)
40+
end
41+
42+
allequal((L.dims, ψ0.dims, A.dims, B.dims, C.dims)) || throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
5043

51-
kwargs2 = merge((saveat = collect(t_l),), (; kwargs...))
52-
ρt = mesolve(H, ψ0, t_l, c_ops; kwargs2...).states
44+
kwargs2 = merge((saveat = collect(tlist),), (; kwargs...))
45+
ρt = mesolve(L, ψ0, tlist; kwargs2...).states
5346

54-
corr = map((t, ρ) -> mesolve(H, C * ρ * A, τ_l .+ t, c_ops, e_ops = [B]; kwargs...).expect[1, :], t_l, ρt)
47+
corr = map((t, ρ) -> mesolve(L, C * ρ * A, τlist .+ t, e_ops = [B]; kwargs...).expect[1, :], tlist, ρt)
5548

5649
return corr
5750
end
5851

5952
@doc raw"""
60-
correlation_2op_2t(H::QuantumObject,
61-
ψ0::QuantumObject,
62-
t_l::AbstractVector,
63-
τ_l::AbstractVector,
53+
correlation_2op_2t(H::AbstractQuantumObject,
54+
ψ0::Union{Nothing,QuantumObject},
55+
tlist::AbstractVector,
56+
τlist::AbstractVector,
57+
c_ops::Union{Nothing,AbstractVector,Tuple},
6458
A::QuantumObject,
65-
B::QuantumObject,
66-
c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
59+
B::QuantumObject;
6760
reverse::Bool=false,
6861
kwargs...)
6962
70-
Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}``
71-
at different times: ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle``.
63+
Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
64+
65+
If the initial state `ψ0` is given as `nothing`, then the [`steadystate`] will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
7266
7367
When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \right\rangle``.
7468
"""
7569
function correlation_2op_2t(
76-
H::QuantumObject{<:AbstractArray{T1},HOpType},
77-
ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
78-
t_l::AbstractVector,
79-
τ_l::AbstractVector,
80-
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
81-
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
82-
c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
70+
H::AbstractQuantumObject{DataType,HOpType},
71+
ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
72+
tlist::AbstractVector,
73+
τlist::AbstractVector,
74+
c_ops::Union{Nothing,AbstractVector,Tuple},
75+
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
76+
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
8377
reverse::Bool = false,
8478
kwargs...,
8579
) where {
80+
DataType,
8681
T1,
8782
T2,
8883
T3,
89-
T4,
9084
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
9185
StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
9286
}
9387
C = eye(prod(H.dims), dims = H.dims)
9488
if reverse
95-
corr = correlation_3op_2t(H, ψ0, t_l, τ_l, A, B, C, c_ops; kwargs...)
89+
corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, A, B, C; kwargs...)
9690
else
97-
corr = correlation_3op_2t(H, ψ0, t_l, τ_l, C, A, B, c_ops; kwargs...)
91+
corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, C, A, B; kwargs...)
9892
end
9993

10094
return reduce(hcat, corr)
10195
end
10296

10397
@doc raw"""
104-
correlation_2op_1t(H::QuantumObject,
105-
ψ0::QuantumObject,
106-
τ_l::AbstractVector,
98+
correlation_2op_1t(H::AbstractQuantumObject,
99+
ψ0::Union{Nothing,QuantumObject},
100+
τlist::AbstractVector,
101+
c_ops::Union{Nothing,AbstractVector,Tuple},
107102
A::QuantumObject,
108-
B::QuantumObject,
109-
c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
103+
B::QuantumObject;
110104
reverse::Bool=false,
111105
kwargs...)
112106
113-
Returns the one-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` at different times ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle``.
107+
Returns the one-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
108+
109+
If the initial state `ψ0` is given as `nothing`, then the [`steadystate`] will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
114110
115111
When `reverse=true`, the correlation function is calculated as ``\left\langle \hat{A}(0) \hat{B}(\tau) \right\rangle``.
116112
"""
117113
function correlation_2op_1t(
118-
H::QuantumObject{<:AbstractArray{T1},HOpType},
119-
ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
120-
τ_l::AbstractVector,
121-
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
122-
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
123-
c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
114+
H::AbstractQuantumObject{DataType,HOpType},
115+
ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
116+
τlist::AbstractVector,
117+
c_ops::Union{Nothing,AbstractVector,Tuple},
118+
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
119+
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject};
124120
reverse::Bool = false,
125121
kwargs...,
126122
) where {
123+
DataType,
127124
T1,
128125
T2,
129126
T3,
130-
T4,
131127
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
132128
StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
133129
}
134-
corr = correlation_2op_2t(H, ψ0, [0], τ_l, A, B, c_ops; reverse = reverse, kwargs...)
130+
corr = correlation_2op_2t(H, ψ0, [0], τlist, c_ops, A, B; reverse = reverse, kwargs...)
135131

136132
return corr[:, 1]
137133
end
138-
139-
@doc raw"""
140-
spectrum(H::QuantumObject,
141-
ω_list::AbstractVector,
142-
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
143-
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
144-
c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
145-
solver::SpectrumSolver=ExponentialSeries(),
146-
kwargs...)
147-
148-
Returns the emission spectrum
149-
150-
```math
151-
S(\omega) = \int_{-\infty}^\infty \left\langle \hat{A}(\tau) \hat{B}(0) \right\rangle e^{-i \omega \tau} d \tau
152-
```
153-
"""
154-
function spectrum(
155-
H::QuantumObject{MT1,HOpType},
156-
ω_list::AbstractVector,
157-
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
158-
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
159-
c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
160-
solver::SpectrumSolver = ExponentialSeries(),
161-
kwargs...,
162-
) where {
163-
MT1<:AbstractMatrix,
164-
T2,
165-
T3,
166-
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
167-
}
168-
return _spectrum(H, ω_list, A, B, c_ops, solver; kwargs...)
169-
end
170-
171-
function _spectrum(
172-
H::QuantumObject{<:AbstractArray{T1},HOpType},
173-
ω_list::AbstractVector,
174-
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
175-
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
176-
c_ops,
177-
solver::FFTCorrelation;
178-
kwargs...,
179-
) where {T1,T2,T3,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
180-
Nsamples = length(ω_list)
181-
ω_max = abs(maximum(ω_list))
182-
dω = 2 * ω_max / (Nsamples - 1)
183-
ω_l = -ω_max:dω:ω_max
184-
185-
T = 2π / (ω_l[2] - ω_l[1])
186-
τ_l = range(0, T, length = length(ω_l))
187-
188-
ρss = steadystate(H, c_ops)
189-
corr = correlation_2op_1t(H, ρss, τ_l, A, B, c_ops; kwargs...)
190-
191-
S = fftshift(fft(corr)) / length(τ_l)
192-
193-
return ω_l, 2 .* real.(S)
194-
end
195-
196-
function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,true}) where {T}
197-
result = eigen(L)
198-
rates, vecs = result.values, result.vectors
199-
200-
return rates, vecs, steadystate(L).data
201-
end
202-
203-
function _spectrum_get_rates_vecs_ss(L, solver::ExponentialSeries{T,false}) where {T}
204-
result = eigen(L)
205-
rates, vecs = result.values, result.vectors
206-
207-
ss_idx = findmin(abs2, rates)[2]
208-
ρss = vec2mat(@view(vecs[:, ss_idx]))
209-
ρss = (ρss + ρss') / 2
210-
ρss ./= tr(ρss)
211-
212-
return rates, vecs, ρss
213-
end
214-
215-
function _spectrum(
216-
H::QuantumObject{<:AbstractArray{T1},HOpType},
217-
ω_list::AbstractVector,
218-
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
219-
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
220-
c_ops,
221-
solver::ExponentialSeries;
222-
kwargs...,
223-
) where {T1,T2,T3,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
224-
allequal((H.dims, A.dims, B.dims)) || throw(DimensionMismatch("The dimensions of H, A and B must be the same"))
225-
226-
L = liouvillian(H, c_ops)
227-
228-
rates, vecs, ρss = _spectrum_get_rates_vecs_ss(L, solver)
229-
230-
ρ0 = B.data * ρss
231-
v = vecs \ mat2vec(ρ0)
232-
233-
amps = map(i -> v[i] * tr(A.data * vec2mat(@view(vecs[:, i]))), eachindex(rates))
234-
idxs = findall(x -> abs(x) > solver.tol, amps)
235-
amps, rates = amps[idxs], rates[idxs]
236-
237-
# spec = map(ω -> 2 * real(sum(@. amps * (1 / (1im * ω - rates)))), ω_list)
238-
amps_rates = zip(amps, rates)
239-
spec = map(ω -> 2 * real(sum(x -> x[1] / (1im * ω - x[2]), amps_rates)), ω_list)
240-
241-
return spec
242-
end

‎src/deprecated.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,3 +13,7 @@ function deprecated_foo(args...; kwargs...)
1313
error("`deprecated_foo` has been deprecated and will be removed in next major release, please use `new_foo` instead.")
1414
end
1515
=#
16+
17+
export FFTCorrelation
18+
19+
FFTCorrelation() = error("`FFTCorrelation` has been deprecated and will be removed in next major release, please use `spectrum_correlation_fft` to calculate the spectrum with FFT method instead.")

0 commit comments

Comments
 (0)