Forward Mode AD for QobjEvo Hamiltonians with Different Time Dependent Coefficients Fails

[Gavin-Rockwood]

Bug Description

Following the example here, if one adds a QobjEvo term with explicit time dependence, the forward mode AD fails. I'm not sure if this is an issue with QuantumToolbox or something deeper in ForwardDiff or SciMLOperators.

Code to Reproduce the Bug

using QuantumToolbox
using ForwardDiff


function my_f_mesolve_direct(p)
    H = p[1] * a' * a + p[2] * (a + a')
    c_ops = [sqrt(p[3]) * a]
    sol = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
    return real(expect(a' * a, sol.states[end]))
end
println("Example Output: ", ForwardDiff.gradient(my_f_mesolve_direct, params))

function my_f_mesolve_direct_2(p)
    H = p[1] * a' * a + p[2] * (a + a') + QobjEvo(a+a', (p,t) -> 0.0)
    c_ops = [sqrt(p[3]) * a]
    sol = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
    return real(expect(a' * a, sol.states[end]))
end
println("QobjEvo with zero coefficient but no time dependence: ", ForwardDiff.gradient(my_f_mesolve_direct_2, params))

function my_f_mesolve_direct_3(p)
    H = p[1] * a' * a + p[2] * (a + a') + QobjEvo(a+a', (p,t) -> 0*t)
    c_ops = [sqrt(p[3]) * a]
    sol = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
    return real(expect(a' * a, sol.states[end]))
end
println("QobjEvo with infinitesimal time-dependent coefficient: ", ForwardDiff.gradient(my_f_mesolve_direct_3, params))

Code Output

Example Output: [-0.853333333359848, 1.0666666665744058, -0.2133333331248055]
QobjEvo with zero coefficient but no time dependence: [-0.8533333334104487, 1.066666666949228, -0.213333332505229]
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(!Matched::Int16)
   @ Base float.jl:239
  ...


Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ Base ./number.jl:7
  [2] ComplexF64(re::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, im::Int64)
    @ Base ./complex.jl:14
  [3] ComplexF64(x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ Base ./complex.jl:42
  [4] convert(::Type{ComplexF64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ Base ./number.jl:7
  [5] setproperty!(x::ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, f::Symbol, v::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ Base ./Base.jl:52
  [6] update_coefficients!(L::ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, u::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, p::NullParameters, t::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}; kwargs::@Kwargs{})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/7ExCq/src/scalar.jl:201
  [7] update_coefficients!(L::ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, u::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, p::NullParameters, t::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/7ExCq/src/scalar.jl:200
  [8] #_#153
    @ ~/.julia/packages/SciMLOperators/7ExCq/src/basic.jl:414 [inlined]
  [9] ScaledOperator
    @ ~/.julia/packages/SciMLOperators/7ExCq/src/basic.jl:412 [inlined]
 [10] macro expansion
    @ ~/.julia/packages/SciMLOperators/7ExCq/src/basic.jl:650 [inlined]
 [11] (::SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}})(w::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, v::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, u::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, p::NullParameters, t::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}; kwargs::@Kwargs{})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/7ExCq/src/basic.jl:639
 [12] (::ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing})(du::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, u::Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, p::NullParameters, t::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3})
    @ SciMLBase ~/.julia/packages/SciMLBase/YE7xF/src/scimlfunctions.jl:2595
 [13] initialize!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Nothing, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}}, ODESolution{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, 2, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}}, Nothing, Nothing, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}}}, Nothing, ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, true, NullParameters, ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{abstol::Float64, reltol::Float64, save_everystep::Bool, save_end::Bool, saveat::Vector{Float64}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Vector{Vector{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqTsit5.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, OrdinaryDiffEqCore.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Bool, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, DataStructures.FasterForward}, DataStructures.BinaryHeap{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing}, cache::OrdinaryDiffEqTsit5.Tsit5Cache{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False})
    @ OrdinaryDiffEqTsit5 ~/.julia/packages/OrdinaryDiffEqTsit5/o07IY/src/tsit_perform_step.jl:175
 [14] __init(prob::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, true, NullParameters, ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{abstol::Float64, reltol::Float64, save_everystep::Bool, save_end::Bool, saveat::Vector{Float64}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Vector{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_discretes::Bool, save_start::Bool, save_end::Bool, callback::Nothing, dense::Bool, calck::Bool, dt::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, dtmin::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, dtmax::ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Float64, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/qzTLI/src/solve.jl:578
 [15] __init (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/qzTLI/src/solve.jl:11 [inlined]
 [16] __solve(::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, true, NullParameters, ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{abstol::Float64, reltol::Float64, save_everystep::Bool, save_end::Bool, saveat::Vector{Float64}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}; kwargs::@Kwargs{abstol::Float64, reltol::Float64, save_everystep::Bool, save_end::Bool, saveat::Vector{Float64}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/qzTLI/src/solve.jl:6
 [17] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/qzTLI/src/solve.jl:1 [inlined]
 [18] solve_call(_prob::ODEProblem{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}}, Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, true, NullParameters, ODEFunction{true, FullSpecialize, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}, MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{abstol::Float64, reltol::Float64, save_everystep::Bool, save_end::Bool, saveat::Vector{Float64}, alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:127
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:84 [inlined]
 [20] #solve_up#39
    @ ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:563 [inlined]
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:540 [inlined]
 [22] #solve#38
    @ ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:530 [inlined]
 [23] solve
    @ ~/.julia/packages/DiffEqBase/p82Yh/src/solve.jl:520 [inlined]
 [24] #mesolve#160
    @ ~/.julia/packages/QuantumToolbox/nzPVM/src/time_evolution/mesolve.jl:208 [inlined]
 [25] mesolve
    @ ~/.julia/packages/QuantumToolbox/nzPVM/src/time_evolution/mesolve.jl:207 [inlined]
 [26] mesolve(H::QuantumObjectEvolution{Operator, Dimensions{1, Tuple{Space}}, SciMLOperators.AddedOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Tuple{MatrixOperator{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}, SciMLOperators.ScaledOperator{ComplexF64, ScalarOperator{ComplexF64, QuantumToolbox.var"#33#34"{Tuple{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, var"#173#174"}}}, MatrixOperator{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}, SciMLOperators.FilterKwargs{Nothing, Tuple{}}}}}}}, ψ0::QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}, tlist::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, c_ops::Vector{QuantumObject{Operator, Dimensions{1, Tuple{Space}}, SparseArrays.SparseMatrixCSC{Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}, Int64}}}; alg::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, e_ops::Nothing, params::NullParameters, progress_bar::Val{false}, inplace::Val{true}, kwargs::@Kwargs{})
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/nzPVM/src/time_evolution/mesolve.jl:201
 [27] mesolve
    @ ~/.julia/packages/QuantumToolbox/nzPVM/src/time_evolution/mesolve.jl:157 [inlined]
 [28] my_f_mesolve_direct_3(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}})
    @ Main ~/Library/Mobile Documents/iCloud~md~obsidian/Documents/research/02 Projects/ModelBasedTomography/Code/MBT_Autodiff/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X26sZmlsZQ==.jl:24
 [29] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/2nvFM/src/apiutils.jl:24 [inlined]
 [30] vector_mode_gradient(f::typeof(my_f_mesolve_direct_3), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/2nvFM/src/gradient.jl:98
 [31] gradient
    @ ~/.julia/packages/ForwardDiff/2nvFM/src/gradient.jl:20 [inlined]
 [32] gradient(f::typeof(my_f_mesolve_direct_3), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_f_mesolve_direct_3), Float64}, Float64, 3}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/2nvFM/src/gradient.jl:17
 [33] gradient(f::typeof(my_f_mesolve_direct_3), x::Vector{Float64})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/2nvFM/src/gradient.jl:17
 [34] top-level scope
    @ ~/Library/Mobile Documents/iCloud~md~obsidian/Documents/research/02 Projects/ModelBasedTomography/Code/MBT_Autodiff/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X26sZmlsZQ==.jl:27

Expected Behaviour

I would expect that the gradient should still be computed and would be the same as the first 2 examples. I would also note that this is not an issue when I have something like:

H = QobjEvo(p[1] * a' * a + p[2] * (a + a'), (p,t) -> 1+1e-32*t)+QobjEvo(a+a', (p,t) -> 0)

which does work. However,

H = QobjEvo(p[1] * a' * a + p[2] * (a + a'), (p,t) -> 1+1e-32*t)+QobjEvo(a+a', (p,t) -> 0*t)

does not work.

Edit:
Interestingly enough, this also works:

H = QobjEvo(p[2] * (a + a'), (p,t) -> 1+0*t)+QobjEvo(p[1] * a' * a, (p,t) -> t)

It seems that mixing terms that contain the target of the gradent p and those that do not is what causes the issue. For the orignal case, this allows for a working around by just tossing in a random factor of p[1]^0 in order to make ForwardDiff happy.

Your Environment

QuantumToolbox.jl: Quantum Toolbox in Julia
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © QuTiP team 2022 and later.
Current admin team:
    Alberto Mercurio and Yi-Te Huang

Package information:
====================================
Julia              Ver. 1.11.7
QuantumToolbox     Ver. 0.35.0
SciMLOperators     Ver. 1.7.2
LinearSolve        Ver. 3.40.1
OrdinaryDiffEqCore Ver. 1.34.0

System information:
====================================
OS       : macOS (arm64-apple-darwin24.0.0)
CPU      : 14 × Apple M4 Pro
Memory   : 48.0 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-16.0.6 (ORCJIT, apple-m1)
BLAS     : libopenblas64_.dylib (ilp64)
Threads  : 10 (on 10 virtual cores)

Additional Context

No response

[albertomercurio]

This is very odd. Let me investigate.

As a first workaround you could use the inplace keyword argument to the solver. The problem is that it is updating inplace the ScalarOperator, which has a ComplexF64 type and wants to be updated into a Dual object`

function my_f_mesolve_direct_3(p)
    H = p[1] * a' * a + p[2] * (a + a') + QobjEvo(a+a', (p,t) -> 0*t)
    c_ops = [sqrt(p[3]) * a]
    sol = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false), inplace=Val(false))
    return real(expect(a' * a, sol.states[end]))
end
println("QobjEvo with infinitesimal time-dependent coefficient: ", ForwardDiff.gradient(my_f_mesolve_direct_3, params))