Skip to content

Commit 6579bc1

Browse files
albertomercurioalbertomercurio
authored
Improve c_ops handling (qutip#216)
1 parent e76d77d commit 6579bc1

File tree

6 files changed

+63
-86
lines changed

6 files changed

+63
-86
lines changed

‎src/correlations.jl‎

Lines changed: 8 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(t
2020
A::QuantumObject,
2121
B::QuantumObject,
2222
C::QuantumObject,
23-
c_ops::AbstractVector=[];
23+
c_ops::Union{Nothing,AbstractVector}=nothing;
2424
kwargs...)
2525
2626
Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\expval{\hat{A}(t) \hat{B}(t + \tau) \hat{C}(t)}``
@@ -35,7 +35,7 @@ function correlation_3op_2t(
3535
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
3636
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
3737
C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
38-
c_ops::AbstractVector = [];
38+
c_ops::Union{Nothing,AbstractVector} = nothing;
3939
kwargs...,
4040
) where {
4141
T1,
@@ -65,7 +65,7 @@ end
6565
τ_l::AbstractVector,
6666
A::QuantumObject,
6767
B::QuantumObject,
68-
c_ops::AbstractVector=[];
68+
c_ops::Union{Nothing,AbstractVector}=nothing;
6969
reverse::Bool=false,
7070
kwargs...)
7171
@@ -81,7 +81,7 @@ function correlation_2op_2t(
8181
τ_l::AbstractVector,
8282
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
8383
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
84-
c_ops::AbstractVector = [];
84+
c_ops::Union{Nothing,AbstractVector} = nothing;
8585
reverse::Bool = false,
8686
kwargs...,
8787
) where {
@@ -108,7 +108,7 @@ end
108108
τ_l::AbstractVector,
109109
A::QuantumObject,
110110
B::QuantumObject,
111-
c_ops::AbstractVector=[];
111+
c_ops::Union{Nothing,AbstractVector}=nothing;
112112
reverse::Bool=false,
113113
kwargs...)
114114
@@ -122,7 +122,7 @@ function correlation_2op_1t(
122122
τ_l::AbstractVector,
123123
A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
124124
B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
125-
c_ops::AbstractVector = [];
125+
c_ops::Union{Nothing,AbstractVector} = nothing;
126126
reverse::Bool = false,
127127
kwargs...,
128128
) where {
@@ -143,7 +143,7 @@ end
143143
ω_list::AbstractVector,
144144
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
145145
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
146-
c_ops::AbstractVector=[];
146+
c_ops::Union{Nothing,AbstractVector}=nothing;
147147
solver::MySolver=ExponentialSeries(),
148148
kwargs...)
149149
@@ -158,16 +158,14 @@ function spectrum(
158158
ω_list::AbstractVector,
159159
A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
160160
B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
161-
c_ops::Vector{QuantumObject{MT2,COpType}} = Vector{QuantumObject{MT1,HOpType}}([]);
161+
c_ops::Union{Nothing,AbstractVector} = nothing;
162162
solver::MySolver = ExponentialSeries(),
163163
kwargs...,
164164
) where {
165165
MT1<:AbstractMatrix,
166-
MT2<:AbstractMatrix,
167166
T2,
168167
T3,
169168
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
170-
COpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
171169
MySolver<:SpectrumSolver,
172170
}
173171
return _spectrum(H, ω_list, A, B, c_ops, solver; kwargs...)

‎src/qobj/eigsolve.jl‎

Lines changed: 4 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -438,7 +438,7 @@ end
438438

439439
@doc raw"""
440440
eigsolve_al(H::QuantumObject,
441-
T::Real, c_ops::AbstractVector=[];
441+
T::Real, c_ops::Union{Nothing,AbstractVector}=nothing;
442442
alg::OrdinaryDiffEqAlgorithm=Tsit5(),
443443
H_t::Union{Nothing,Function}=nothing,
444444
params::NamedTuple=NamedTuple(),
@@ -454,7 +454,7 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
454454
# Arguments
455455
- `H`: The Hamiltonian (or directly the Liouvillian) of the system.
456456
- `T`: The time at which to evaluate the time evolution
457-
- `c_ops`: A vector of collapse operators
457+
- `c_ops`: A vector of collapse operators. Default is `nothing` meaning the system is closed.
458458
- `alg`: The differential equation solver algorithm
459459
- `H_t`: A function `H_t(t)` that returns the additional term at time `t`
460460
- `params`: A dictionary of additional parameters
@@ -476,7 +476,7 @@ and Floquet open quantum systems. Quantum, 6, 649.
476476
function eigsolve_al(
477477
H::QuantumObject{MT1,HOpType},
478478
T::Real,
479-
c_ops::Vector{QuantumObject{MT2,COpType}} = Vector{QuantumObject{MT1,HOpType}}([]);
479+
c_ops::Union{Nothing,AbstractVector} = nothing;
480480
alg::OrdinaryDiffEqAlgorithm = Tsit5(),
481481
H_t::Union{Nothing,Function} = nothing,
482482
params::NamedTuple = NamedTuple(),
@@ -486,12 +486,7 @@ function eigsolve_al(
486486
maxiter::Int = 200,
487487
eigstol::Real = 1e-6,
488488
kwargs...,
489-
) where {
490-
MT1<:AbstractMatrix,
491-
MT2<:AbstractMatrix,
492-
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
493-
COpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
494-
}
489+
) where {MT1<:AbstractMatrix,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
495490
L = liouvillian(H, c_ops)
496491
prob = mesolveProblem(
497492
L,

‎src/steadystate.jl‎

Lines changed: 9 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ end
4343
H::QuantumObject{MT1,HOpType},
4444
ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
4545
tspan::Real = Inf,
46-
c_ops::Vector{QuantumObject{Tc,COpType}} = QuantumObject{MT1,HOpType}[];
46+
c_ops::Union{Nothing,AbstractVector} = nothing;
4747
solver::SteadyStateODESolver = SteadyStateODESolver(),
4848
reltol::Real = 1.0e-8,
4949
abstol::Real = 1.0e-10,
@@ -68,7 +68,7 @@ or
6868
- `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
6969
- `ψ0::QuantumObject`: The initial state of the system.
7070
- `tspan::Real=Inf`: The final time step for the steady state problem.
71-
- `c_ops::AbstractVector=[]`: The list of the collapse operators.
71+
- `c_ops::Union{Nothing,AbstractVector}=nothing`: The list of the collapse operators.
7272
- `solver::SteadyStateODESolver=SteadyStateODESolver()`: see [`SteadyStateODESolver`](@ref) for more details.
7373
- `reltol::Real=1.0e-8`: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
7474
- `abstol::Real=1.0e-10`: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
@@ -78,18 +78,16 @@ function steadystate(
7878
H::QuantumObject{MT1,HOpType},
7979
ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
8080
tspan::Real = Inf,
81-
c_ops::Vector{QuantumObject{Tc,COpType}} = QuantumObject{MT1,HOpType}[];
81+
c_ops::Union{Nothing,AbstractVector} = nothing;
8282
solver::SteadyStateODESolver = SteadyStateODESolver(),
8383
reltol::Real = 1.0e-8,
8484
abstol::Real = 1.0e-10,
8585
kwargs...,
8686
) where {
8787
MT1<:AbstractMatrix,
8888
T2,
89-
Tc<:AbstractMatrix,
9089
HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
9190
StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
92-
COpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
9391
}
9492
(H.dims != ψ0.dims) && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
9593

@@ -125,22 +123,14 @@ function _steadystate_ode_condition(integrator, abstol, reltol, min_t)
125123
end
126124

127125
function steadystate(
128-
L::QuantumObject{<:AbstractArray{T},SuperOperatorQuantumObject};
126+
H::QuantumObject{<:AbstractArray,OpType},
127+
c_ops::Union{Nothing,AbstractVector} = nothing;
129128
solver::SteadyStateSolver = SteadyStateDirectSolver(),
130129
kwargs...,
131-
) where {T}
132-
return _steadystate(L, solver; kwargs...)
133-
end
134-
135-
function steadystate(
136-
H::QuantumObject{<:AbstractArray{T},OpType},
137-
c_ops::AbstractVector;
138-
solver::SteadyStateSolver = SteadyStateDirectSolver(),
139-
kwargs...,
140-
) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
130+
) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
141131
L = liouvillian(H, c_ops)
142132

143-
return steadystate(L; solver = solver, kwargs...)
133+
return _steadystate(L, solver; kwargs...)
144134
end
145135

146136
function _steadystate(
@@ -232,7 +222,7 @@ end
232222
H_p::QuantumObject{<:AbstractArray,OpType2},
233223
H_m::QuantumObject{<:AbstractArray,OpType3},
234224
ωd::Number,
235-
c_ops::AbstractVector = QuantumObject{MT,OpType1}[];
225+
c_ops::Union{Nothing,AbstractVector} = nothing;
236226
n_max::Integer = 2,
237227
tol::R = 1e-8,
238228
solver::FSolver = SSFloquetLinearSystem,
@@ -310,7 +300,7 @@ function steadystate_floquet(
310300
H_p::QuantumObject{<:AbstractArray,OpType2},
311301
H_m::QuantumObject{<:AbstractArray,OpType3},
312302
ωd::Number,
313-
c_ops::AbstractVector = QuantumObject{MT,OpType1}[];
303+
c_ops::Union{Nothing,AbstractVector} = nothing;
314304
n_max::Integer = 2,
315305
tol::R = 1e-8,
316306
solver::FSolver = SSFloquetLinearSystem(),

‎src/time_evolution/mcsolve.jl‎

Lines changed: 26 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ end
101101
mcsolveProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
102102
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
103103
tlist::AbstractVector,
104-
c_ops::Vector{QuantumObject{Tc, OperatorQuantumObject}}=QuantumObject{Matrix, OperatorQuantumObject}[];
104+
c_ops::Union{Nothing,AbstractVector}=nothing;
105105
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm=Tsit5(),
106106
e_ops::Union{Nothing,AbstractVector}=nothing,
107107
H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
@@ -148,7 +148,7 @@ If the environmental measurements register a quantum jump, the wave function und
148148
- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
149149
- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
150150
- `tlist::AbstractVector`: List of times at which to save the state of the system.
151-
- `c_ops::Vector`: List of collapse operators ``\{\hat{C}_n\}_n``.
151+
- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
152152
- `alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
153153
- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
154154
- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
@@ -170,25 +170,28 @@ If the environmental measurements register a quantum jump, the wave function und
170170
"""
171171
function mcsolveProblem(
172172
H::QuantumObject{MT1,OperatorQuantumObject},
173-
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
173+
ψ0::QuantumObject{<:AbstractArray,KetQuantumObject},
174174
tlist::AbstractVector,
175-
c_ops::Vector{QuantumObject{Tc,OperatorQuantumObject}} = QuantumObject{MT1,OperatorQuantumObject}[];
175+
c_ops::Union{Nothing,AbstractVector} = nothing;
176176
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = Tsit5(),
177177
e_ops::Union{Nothing,AbstractVector} = nothing,
178178
H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
179179
params::NamedTuple = NamedTuple(),
180180
seeds::Union{Nothing,Vector{Int}} = nothing,
181181
jump_callback::TJC = ContinuousLindbladJumpCallback(),
182182
kwargs...,
183-
) where {MT1<:AbstractMatrix,T2,Tc<:AbstractMatrix,TJC<:LindbladJumpCallbackType}
183+
) where {MT1<:AbstractMatrix,TJC<:LindbladJumpCallbackType}
184184
H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
185185

186186
haskey(kwargs, :save_idxs) &&
187187
throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
188188

189+
c_ops isa Nothing &&
190+
throw(ArgumentError("The list of collapse operators must be provided. Use sesolveProblem instead."))
191+
189192
t_l = convert(Vector{Float64}, tlist) # Convert it into Float64 to avoid type instabilities for OrdinaryDiffEq.jl
190193

191-
H_eff = H - T2(0.5im) * mapreduce(op -> op' * op, +, c_ops)
194+
H_eff = H - 1im * mapreduce(op -> op' * op, +, c_ops) / 2
192195

193196
if e_ops isa Nothing
194197
expvals = Array{ComplexF64}(undef, 0, length(t_l))
@@ -254,15 +257,15 @@ function mcsolveProblem(
254257
end
255258

256259
function mcsolveProblem(
257-
H_eff::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
258-
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
260+
H_eff::QuantumObject{<:AbstractArray,OperatorQuantumObject},
261+
ψ0::QuantumObject{<:AbstractArray,KetQuantumObject},
259262
t_l::AbstractVector,
260263
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm,
261264
H_t::Union{Nothing,Function,TimeDependentOperatorSum},
262265
params::NamedTuple,
263266
jump_callback::ContinuousLindbladJumpCallback;
264267
kwargs...,
265-
) where {T1,T2}
268+
)
266269
cb1 = ContinuousCallback(
267270
LindbladJumpContinuousCondition,
268271
LindbladJumpAffect!,
@@ -283,9 +286,9 @@ end
283286
mcsolveEnsembleProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
284287
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
285288
tlist::AbstractVector,
286-
c_ops::Vector{QuantumObject{Tc, OperatorQuantumObject}}=QuantumObject{Matrix, OperatorQuantumObject}[];
289+
c_ops::Union{Nothing,AbstractVector}=nothing;
287290
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm=Tsit5(),
288-
e_ops::Vector{QuantumObject{Te, OperatorQuantumObject}}=QuantumObject{Matrix, OperatorQuantumObject}[],
291+
e_ops::Union{Nothing,AbstractVector}=nothing,
289292
H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
290293
params::NamedTuple=NamedTuple(),
291294
jump_callback::TJC=ContinuousLindbladJumpCallback(),
@@ -332,9 +335,9 @@ If the environmental measurements register a quantum jump, the wave function und
332335
- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
333336
- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
334337
- `tlist::AbstractVector`: List of times at which to save the state of the system.
335-
- `c_ops::Vector`: List of collapse operators ``\{\hat{C}_n\}_n``.
338+
- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
336339
- `alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
337-
- `e_ops::Vector`: List of operators for which to calculate expectation values.
340+
- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
338341
- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
339342
- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
340343
- `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -358,17 +361,17 @@ function mcsolveEnsembleProblem(
358361
H::QuantumObject{MT1,OperatorQuantumObject},
359362
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
360363
tlist::AbstractVector,
361-
c_ops::Vector{QuantumObject{Tc,OperatorQuantumObject}} = QuantumObject{MT1,OperatorQuantumObject}[];
364+
c_ops::Union{Nothing,AbstractVector} = nothing;
362365
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = Tsit5(),
363-
e_ops::Vector{QuantumObject{Te,OperatorQuantumObject}} = QuantumObject{MT1,OperatorQuantumObject}[],
366+
e_ops::Union{Nothing,AbstractVector} = nothing,
364367
H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
365368
params::NamedTuple = NamedTuple(),
366369
jump_callback::TJC = ContinuousLindbladJumpCallback(),
367370
seeds::Union{Nothing,Vector{Int}} = nothing,
368371
prob_func::Function = _mcsolve_prob_func,
369372
output_func::Function = _mcsolve_output_func,
370373
kwargs...,
371-
) where {MT1<:AbstractMatrix,T2,Tc<:AbstractMatrix,Te<:AbstractMatrix,TJC<:LindbladJumpCallbackType}
374+
) where {MT1<:AbstractMatrix,T2,TJC<:LindbladJumpCallbackType}
372375
prob_mc = mcsolveProblem(
373376
H,
374377
ψ0,
@@ -392,9 +395,9 @@ end
392395
mcsolve(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
393396
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
394397
tlist::AbstractVector,
395-
c_ops::Vector{QuantumObject{Tc, OperatorQuantumObject}}=QuantumObject{Matrix, OperatorQuantumObject}[];
398+
c_ops::Union{Nothing,AbstractVector}=nothing;
396399
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm=Tsit5(),
397-
e_ops::Vector{QuantumObject{Te, OperatorQuantumObject}}=QuantumObject{Matrix, OperatorQuantumObject}[],
400+
e_ops::Union{Nothing,AbstractVector}=nothing,
398401
H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
399402
params::NamedTuple=NamedTuple(),
400403
n_traj::Int=1,
@@ -441,9 +444,9 @@ If the environmental measurements register a quantum jump, the wave function und
441444
- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
442445
- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
443446
- `tlist::AbstractVector`: List of times at which to save the state of the system.
444-
- `c_ops::Vector`: List of collapse operators ``\{\hat{C}_n\}_n``.
447+
- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
445448
- `alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
446-
- `e_ops::Vector`: List of operators for which to calculate expectation values.
449+
- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
447450
- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
448451
- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
449452
- `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -470,9 +473,9 @@ function mcsolve(
470473
H::QuantumObject{MT1,OperatorQuantumObject},
471474
ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
472475
tlist::AbstractVector,
473-
c_ops::Vector{QuantumObject{Tc,OperatorQuantumObject}} = QuantumObject{MT1,OperatorQuantumObject}[];
476+
c_ops::Union{Nothing,AbstractVector} = nothing;
474477
alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = Tsit5(),
475-
e_ops::Vector{QuantumObject{Te,OperatorQuantumObject}} = QuantumObject{MT1,OperatorQuantumObject}[],
478+
e_ops::Union{Nothing,AbstractVector} = nothing,
476479
H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
477480
params::NamedTuple = NamedTuple(),
478481
seeds::Union{Nothing,Vector{Int}} = nothing,
@@ -482,7 +485,7 @@ function mcsolve(
482485
prob_func::Function = _mcsolve_prob_func,
483486
output_func::Function = _mcsolve_output_func,
484487
kwargs...,
485-
) where {MT1<:AbstractMatrix,T2,Tc<:AbstractMatrix,Te<:AbstractMatrix,TJC<:LindbladJumpCallbackType}
488+
) where {MT1<:AbstractMatrix,T2,TJC<:LindbladJumpCallbackType}
486489
if !isnothing(seeds) && length(seeds) != n_traj
487490
throw(ArgumentError("Length of seeds must match n_traj ($n_traj), but got $(length(seeds))"))
488491
end

0 commit comments

Comments
 (0)