Merge branch 'main' into multi-traj-sol

Author: ytdHuang

CHANGELOG.md

@@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Improve efficiency of `bloch_redfield_tensor` by avoiding unnecessary conversions. ([#509])
+- Support `SciMLOperators v1.4+`. ([#470])
+- Fix compatibility with `Makie v0.24+`. ([#513])
 - Add `keep_runs_results` option for multi-trajectory solvers to align with `QuTiP`. ([#512])
   - Breaking changes for multi-trajectory solutions:
     - the original fields `expect` and `states` now store the results depend on keyword argument `keep_runs_results` (decide whether to store all trajectories results or not).
@@ -263,6 +266,7 @@ Release date: 2024-11-13
 [#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
 [#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
 [#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
+[#470]: https://github.com/qutip/QuantumToolbox.jl/issues/470
 [#472]: https://github.com/qutip/QuantumToolbox.jl/issues/472
 [#473]: https://github.com/qutip/QuantumToolbox.jl/issues/473
 [#476]: https://github.com/qutip/QuantumToolbox.jl/issues/476
@@ -277,4 +281,6 @@ Release date: 2024-11-13
 [#504]: https://github.com/qutip/QuantumToolbox.jl/issues/504
 [#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506
 [#507]: https://github.com/qutip/QuantumToolbox.jl/issues/507
+[#509]: https://github.com/qutip/QuantumToolbox.jl/issues/509
 [#512]: https://github.com/qutip/QuantumToolbox.jl/issues/512
+[#513]: https://github.com/qutip/QuantumToolbox.jl/issues/513

Project.toml

@@ -56,13 +56,13 @@ KernelAbstractions = "0.9.2"
 LaTeXStrings = "1.2"
 LinearAlgebra = "1"
 LinearSolve = "2, 3"
-Makie = "0.20, 0.21, 0.22, 0.23, 0.24"
+Makie = "0.24"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqTsit5 = "1"
 Pkg = "1"
 Random = "1"
 SciMLBase = "2"
-SciMLOperators = "0.3, 0.4"
+SciMLOperators = "1.4"
 SparseArrays = "1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"

docs/src/users_guide/plotting_the_bloch_sphere.md

@@ -195,8 +195,9 @@ At the end of the last section we saw that the colors and marker shapes of the d
 | `b.sphere_color` | Color of Bloch sphere surface | `0.2` |
 | `b.sphere_alpha` | Transparency of sphere surface | `"#FFDDDD"` |
 | `b.vector_color` | Colors for vectors | `["green", "#CC6600", "blue", "red"]` |
-| `b.vector_width` | Width of vectors | `0.025` |
-| `b.vector_arrowsize` | Scales the size of the arrow head. The first two elements scale the radius (in `x/y` direction) and the last one is the length of the cone. | `[0.07, 0.08, 0.08]` |
+| `b.vector_width` | Width of vectors | `0.02` |
+| `b.vector_tiplength` | Length of vector arrow head | `0.08` |
+| `b.vector_tipradius` | Radius of vector arrow head | `0.05` |
 | `b.view` | Azimuthal and elevation viewing angles in degrees | `[30, 30]` |
 | `b.xlabel` | Labels for x-axis | `[L"x", ""]` (``+x`` and ``-x``) |
 | `b.xlpos` | Positions of x-axis labels | `[1.2, -1.2]` |

docs/src/users_guide/settings.md

@@ -13,7 +13,8 @@ Here, we list out each setting along with the specific functions that will use i
 
 - `tidyup_tol::Float64 = 1e-14` : tolerance for [`tidyup`](@ref) and [`tidyup!`](@ref).
 - `auto_tidyup::Bool = true` : Automatically tidyup during the following situations:
-    * Solving for eigenstates, including [`eigenstates`](@ref), [`eigsolve`](@ref), and [`eigsolve_al`](@ref).
+    * Solving for eigenstates, including [`eigenstates`](@ref), [`eigsolve`](@ref), [`eigsolve_al`](@ref)
+    * Creating [`bloch_redfield_tensor`](@ref) or [`brterm`](@ref), and solving [`brmesolve`](@ref).
 - (to be announced)
 
 ## Change default settings

docs/src/users_guide/time_evolution/brmesolve.md

@@ -162,7 +162,7 @@ H = -Δ/2.0 * sigmax() - ε0/2 * sigmaz()
 
 ohmic_spectrum(ω) = (ω == 0.0) ? γ1 : γ1 / 2 * (ω / (2 * π)) * (ω > 0.0)
 
-R, U = bloch_redfield_tensor(H, [(sigmax(), ohmic_spectrum)])
+R, U = bloch_redfield_tensor(H, ((sigmax(), ohmic_spectrum), ))
 
 R
 ```
@@ -183,6 +183,8 @@ H = - Δ/2.0 * sigmax() - ϵ0/2.0 * sigmaz()
 ohmic_spectrum(ω) = (ω == 0.0) ? γ1 : γ1 / 2 * (ω / (2 * π)) * (ω > 0.0)
 
 a_ops = ((sigmax(), ohmic_spectrum),)
+R = bloch_redfield_tensor(H, a_ops; fock_basis = Val(true))
+
 e_ops = [sigmax(), sigmay(), sigmaz()]
 
 # same initial random ket state in QuTiP doc 
@@ -192,7 +194,7 @@ e_ops = [sigmax(), sigmay(), sigmaz()]
 ])
 
 tlist = LinRange(0, 15.0, 1000)
-sol = brmesolve(H, ψ0, tlist, a_ops, e_ops=e_ops)
+sol = mesolve(R, ψ0, tlist, e_ops=e_ops)
 expt_list = real(sol.expect)
 
 # plot the evolution of state on Bloch sphere

ext/QuantumToolboxMakieExt.jl

@@ -22,7 +22,7 @@ import Makie:
     Sphere,
     lines!,
     scatter!,
-    arrows!,
+    arrows3d!,
     text!,
     Point3f,
     mesh!,
@@ -656,26 +656,17 @@ Scales vectors appropriately and adds `3D` arrow markers.
 function _plot_vectors!(b::Bloch, lscene)
     isempty(b.vectors) && return nothing
 
-    arrow_head_length = b.vector_arrowsize[3]
     for (i, v) in enumerate(b.vectors)
         color = get(b.vector_color, i, RGBAf(0.2, 0.5, 0.8, 0.9))
-        nv = norm(v)
-        (arrow_head_length < nv) || throw(
-            ArgumentError(
-                "The length of vector arrow head (Bloch.vector_arrowsize[3]=$arrow_head_length) should be shorter than vector norm: $nv",
-            ),
-        )
 
-        # multiply by the following factor makes the end point of arrow head represent the actual vector position.
-        vec = (1 - arrow_head_length / nv) * Vec3f(v...)
-        arrows!(
+        arrows3d!(
             lscene,
             [Point3f(0, 0, 0)],
-            [vec],
+            [v],
             color = color,
-            linewidth = b.vector_width,
-            arrowsize = Vec3f(b.vector_arrowsize...),
-            arrowcolor = color,
+            shaftradius = b.vector_width,
+            tiplength = b.vector_tiplength,
+            tipradius = b.vector_tipradius,
             rasterize = 3,
         )
     end

src/QuantumToolbox.jl

@@ -17,6 +17,7 @@ import SciMLBase:
     u_modified!,
     NullParameters,
     ODEFunction,
+    SDEFunction,
     ODEProblem,
     SDEProblem,
     EnsembleProblem,

src/qobj/quantum_object_evo.jl

@@ -512,7 +512,7 @@ function (A::QuantumObjectEvolution)(
         )
     end
 
-    A.data(ψout.data, ψin.data, p, t)
+    A.data(ψout.data, ψin.data, nothing, p, t)
 
     return ψout
 end

src/time_evolution/brmesolve.jl

@@ -37,18 +37,28 @@ function bloch_redfield_tensor(
     U = QuantumObject(rst.vectors, Operator(), H.dimensions)
     sec_cutoff = float(sec_cutoff)
 
-    # in fock basis
-    R0 = liouvillian(H, c_ops)
+    H_new = getVal(fock_basis) ? H : QuantumObject(Diagonal(rst.values), Operator(), H.dimensions)
+    c_ops_new = isnothing(c_ops) ? nothing : map(x -> getVal(fock_basis) ? x : U' * x * U, c_ops)
+    L0 = liouvillian(H_new, c_ops_new)
 
-    # set fock_basis=Val(false) and change basis together at the end
-    R1 = 0
-    isempty(a_ops) || (R1 += mapreduce(x -> _brterm(rst, x[1], x[2], sec_cutoff, Val(false)), +, a_ops))
+    # Check whether we can rotate the terms to the eigenbasis directly in the Hamiltonian space
+    fock_basis_hamiltonian = getVal(fock_basis) && sec_cutoff == -1
 
-    SU = sprepost(U, U') # transformation matrix from eigen basis back to fock basis
+    R = isempty(a_ops) ? 0 : sum(x -> _brterm(rst, x[1], x[2], sec_cutoff, fock_basis_hamiltonian), a_ops)
+
+    # If in fock basis, we need to transform the terms back to the fock basis
+    # Note: we can transform the terms in the Hamiltonian space only if sec_cutoff is -1
+    # otherwise, we need to use the SU superoperator below to transform the entire Liouvillian
+    # at the end, due to the action of M_cut
     if getVal(fock_basis)
-        return R0 + SU * R1 * SU'
+        if fock_basis_hamiltonian
+            return L0 + R # Already rotated in the Hamiltonian space
+        else
+            SU = sprepost(U, U')
+            return L0 + SU * R * SU'
+        end
     else
-        return SU' * R0 * SU + R1, U
+        return L0 + R, U
     end
 end
 
@@ -86,11 +96,21 @@ function brterm(
     fock_basis::Union{Bool,Val} = Val(false),
 )
     rst = eigenstates(H)
-    term = _brterm(rst, a_op, spectra, sec_cutoff, makeVal(fock_basis))
+    U = QuantumObject(rst.vectors, Operator(), H.dimensions)
+
+    # Check whether we can rotate the terms to the eigenbasis directly in the Hamiltonian space
+    fock_basis_hamiltonian = getVal(fock_basis) && sec_cutoff == -1
+
+    term = _brterm(rst, a_op, spectra, sec_cutoff, fock_basis_hamiltonian)
     if getVal(fock_basis)
-        return term
+        if fock_basis_hamiltonian
+            return term # Already rotated in the Hamiltonian space
+        else
+            SU = sprepost(U, U')
+            return SU * term * SU'
+        end
     else
-        return term, Qobj(rst.vectors, Operator(), rst.dimensions)
+        return term, U
     end
 end
 
@@ -99,7 +119,7 @@ function _brterm(
     a_op::T,
     spectra::F,
     sec_cutoff::Real,
-    fock_basis::Union{Val{true},Val{false}},
+    fock_basis_hamiltonian::Union{Bool,Val},
 ) where {T<:QuantumObject{Operator},F<:Function}
     _check_br_spectra(spectra)
 
@@ -110,9 +130,11 @@ function _brterm(
     spectrum = spectra.(skew)
 
     A_mat = U' * a_op.data * U
+    A_mat_spec = A_mat .* spectrum
+    A_mat_spec_t = A_mat .* transpose(spectrum)
 
-    ac_term = (A_mat .* spectrum) * A_mat
-    bd_term = A_mat * (A_mat .* trans(spectrum))
+    ac_term = A_mat_spec * A_mat
+    bd_term = A_mat * A_mat_spec_t
 
     if sec_cutoff != -1
         m_cut = similar(skew)
@@ -124,20 +146,29 @@ function _brterm(
         M_cut = @. abs(vec_skew - vec_skew') < sec_cutoff
     end
 
-    out =
-        0.5 * (
-            + _sprepost(A_mat .* trans(spectrum), A_mat) + _sprepost(A_mat, A_mat .* spectrum) - _spost(ac_term, Id) -
-            _spre(bd_term, Id)
-        )
+    # Rotate the terms to the eigenbasis if possible
+    if getVal(fock_basis_hamiltonian)
+        A_mat = U * A_mat * U'
+        A_mat_spec = U * A_mat_spec * U'
+        A_mat_spec_t = U * A_mat_spec_t * U'
+        ac_term = U * ac_term * U'
+        bd_term = U * bd_term * U'
+    end
+
+    # Remove small values before passing in the Liouville space
+    if settings.auto_tidyup
+        tidyup!(A_mat)
+        tidyup!(A_mat_spec)
+        tidyup!(A_mat_spec_t)
+        tidyup!(ac_term)
+        tidyup!(bd_term)
+    end
+
+    out = (_sprepost(A_mat_spec_t, A_mat) + _sprepost(A_mat, A_mat_spec) - _spost(ac_term, Id) - _spre(bd_term, Id)) / 2
 
     (sec_cutoff != -1) && (out .*= M_cut)
 
-    if getVal(fock_basis)
-        SU = _sprepost(U, U')
-        return QuantumObject(SU * out * SU', SuperOperator(), rst.dimensions)
-    else
-        return QuantumObject(out, SuperOperator(), rst.dimensions)
-    end
+    return QuantumObject(out, SuperOperator(), rst.dimensions)
 end
 
 @doc raw"""

src/time_evolution/time_evolution.jl

@@ -481,40 +481,37 @@ end
 
 A struct to represent the diffusion operator. This is used to perform the diffusion process on N different Wiener processes.
 =#
-struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}} <: AbstractSciMLOperator{T}
+struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}}
     ops::OpType
     function DiffusionOperator(ops::OpType) where {OpType}
         T = mapreduce(eltype, promote_type, ops)
         return new{T,OpType}(ops)
     end
 end
 
-@generated function update_coefficients!(L::DiffusionOperator, u, p, t)
-    ops_types = L.parameters[2].parameters
-    N = length(ops_types)
-    return quote
-        Base.@nexprs $N i -> begin
-            update_coefficients!(L.ops[i], u, p, t)
-        end
-
-        nothing
-    end
-end
-
-@generated function LinearAlgebra.mul!(v::AbstractVecOrMat, L::DiffusionOperator, u::AbstractVecOrMat)
+@generated function (L::DiffusionOperator)(w, v, p, t)
     ops_types = L.parameters[2].parameters
     N = length(ops_types)
     quote
-        M = length(u)
-        S = (size(v, 1), size(v, 2)) # This supports also `v` as a `Vector`
+        M = length(v)
+        S = (size(w, 1), size(w, 2)) # This supports also `w` as a `Vector`
         (S[1] == M && S[2] == $N) || throw(DimensionMismatch("The size of the output vector is incorrect."))
         Base.@nexprs $N i -> begin
-            mul!(@view(v[:, i]), L.ops[i], u)
+            op = L.ops[i]
+            op(@view(w[:, i]), v, v, p, t)
         end
-        return v
+        return w
     end
 end
 
+#TODO: Remove when a new release of SciMLBase.jl >2.104.0 is available
+(f::SDEFunction)(du, u, p, t) =
+    if f.f isa AbstractSciMLOperator
+        f.f(du, u, u, p, t)
+    else
+        f.f(du, u, p, t)
+    end
+
 #######################################
 
 function liouvillian_floquet(