diff --git a/docs/Project.toml b/docs/Project.toml index 88173fe01..b61136e89 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,4 +2,5 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab" diff --git a/docs/make.jl b/docs/make.jl index d0a631ad6..375ffec53 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,11 +3,14 @@ using QuantumToolbox using Documenter +using DocumenterCitations DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, :(using QuantumToolbox); recursive = true) const DRAFT = false # set `true` to disable cell evaluation +bib = CitationBibliography(joinpath(@__DIR__, "src", "bibliography.bib"), style=:authoryear) + const MathEngine = MathJax3( Dict( :loader => Dict("load" => ["[tex]/physics"]), @@ -58,6 +61,7 @@ const PAGES = [ ], ], "API" => "api.md", + "Bibliography" => "bibliography.md", # "Change Log" => "changelog.md", ] @@ -76,6 +80,7 @@ makedocs(; size_threshold_ignore = ["api.md"], ), draft = DRAFT, + plugins = [bib], ) deploydocs(; repo = "github.com/qutip/QuantumToolbox.jl", devbranch = "main") diff --git a/docs/src/api.md b/docs/src/api.md index dae7f72ae..5ba752d0c 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -182,17 +182,18 @@ mcsolveProblem mcsolveEnsembleProblem ssesolveProblem ssesolveEnsembleProblem -lr_mesolveProblem sesolve mesolve mcsolve ssesolve dfd_mesolve -dsf_mesolve -dsf_mcsolve -lr_mesolve liouvillian liouvillian_generalized +``` + +### [Steady State Solvers](@id doc-API:Steady-State-Solvers) + +```@docs steadystate steadystate_floquet SteadyStateDirectSolver @@ -201,6 +202,21 @@ SteadyStateLinearSolver SteadyStateODESolver ``` +### [Dynamical Shifted Fock method](@id doc-API:Dynamical-Shifted-Fock-method) + +```@docs +dsf_mesolve +dsf_mcsolve +``` + +### [Low-rank time evolution](@id doc-API:Low-rank-time-evolution) + +```@docs +TimeEvolutionLRSol +lr_mesolveProblem +lr_mesolve +``` + ## [Correlations and Spectrum](@id doc-API:Correlations-and-Spectrum) ```@docs @@ -219,6 +235,14 @@ tracedist fidelity ``` +## [Spin Lattice](@id doc-API:Spin-Lattice) + +```@docs +Lattice +SingleSiteOperator +DissipativeIsing +``` + ## [Miscellaneous](@id doc-API:Miscellaneous) ```@docs @@ -243,10 +267,4 @@ PhysicalConstants convert_unit row_major_reshape meshgrid -_calculate_expectation! -_adjM_condition_variational -_adjM_affect! -_adjM_condition_ratio -_pinv! -dBdz! ``` diff --git a/docs/src/bibliography.bib b/docs/src/bibliography.bib new file mode 100644 index 000000000..727766d5f --- /dev/null +++ b/docs/src/bibliography.bib @@ -0,0 +1,14 @@ +@article{gravina2024adaptive, + title = {{Adaptive variational low-rank dynamics for open quantum systems}}, + author = {Gravina, Luca and Savona, Vincenzo}, + journal = {Phys. Rev. Res.}, + volume = {6}, + issue = {2}, + pages = {023072}, + numpages = {18}, + year = {2024}, + month = {Apr}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevResearch.6.023072}, + url = {https://link.aps.org/doi/10.1103/PhysRevResearch.6.023072} +} diff --git a/docs/src/bibliography.md b/docs/src/bibliography.md new file mode 100644 index 000000000..74d56db57 --- /dev/null +++ b/docs/src/bibliography.md @@ -0,0 +1,8 @@ +```@meta +CurrentModule = QuantumToolbox +``` + +# [Bibliography](@id doc:Bibliography) + +```@bibliography +``` diff --git a/docs/src/tutorials/lowrank.md b/docs/src/tutorials/lowrank.md index 030446dd4..0a3c6c86b 100644 --- a/docs/src/tutorials/lowrank.md +++ b/docs/src/tutorials/lowrank.md @@ -1,5 +1,19 @@ # [Low rank master equation](@id doc-tutor:Low-rank-master-equation) +In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explaination of the method, we recommend to read the article [gravina2024adaptive](@cite). + +As a test, we will consider the dissipative Ising model with a transverse field. The Hamiltonian is given by + +```math +\hat{H} = \frac{J_x}{2} \sum_{\langle i,j \rangle} \sigma_i^x \sigma_j^x + \frac{J_y}{2} \sum_{\langle i,j \rangle} \sigma_i^y \sigma_j^y + \frac{J_z}{2} \sum_{\langle i,j \rangle} \sigma_i^z \sigma_j^z - \sum_i h_i \sigma_i^z + h_x \sum_i \sigma_i^x + h_y \sum_i \sigma_i^y + h_z \sum_i \sigma_i^z, +``` + +where the sums are over nearest neighbors, and the collapse operators are given by + +```math +c_i = \sqrt{\gamma} \sigma_i^-. +``` + We start by importing the packages ```@example lowrank @@ -21,42 +35,42 @@ Define lr-space dimensions N_cut = 2 # Number of states of each mode N_modes = latt.N # Number of modes N = N_cut^N_modes # Total number of states -M = Nx * Ny + 1 # Number of states in the LR basis +M = latt.N + 1 # Number of states in the LR basis ``` Define lr states. Take as initial state all spins up. All other N states are taken as those with miniman Hamming distance to the initial state. ```@example lowrank -ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject}}(undef, M) -ϕ[1] = kron(repeat([basis(2, 0)], N_modes)...) +ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,M-1}}(undef, M) +ϕ[1] = kron(fill(basis(2, 1), N_modes)...) -global i = 1 +i = 1 for j in 1:N_modes global i += 1 - i <= M && (ϕ[i] = mb(sp, j, latt) * ϕ[1]) + i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), j, latt) * ϕ[1]) end for k in 1:N_modes-1 for l in k+1:N_modes global i += 1 - i <= M && (ϕ[i] = mb(sp, k, latt) * mb(sp, l, latt) * ϕ[1]) + i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), k, latt) * SingleSiteOperator(sigmap(), l, latt) * ϕ[1]) end end for i in i+1:M ϕ[i] = QuantumObject(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims) normalize!(ϕ[i]) end +nothing # hide ``` Define the initial state ```@example lowrank -z = hcat(broadcast(x -> x.data, ϕ)...) -p0 = 0.0 # Population of the lr states other than the initial state -B = Matrix(Diagonal([1 + 0im; p0 * ones(M - 1)])) +z = hcat(get_data.(ϕ)...) +B = Matrix(Diagonal([1 + 0im; zeros(M - 1)])) S = z' * z # Overlap matrix B = B / tr(S * B) # Normalize B -ρ = QuantumObject(z * B * z', dims = ones(Int, N_modes) * N_cut); # Full density matrix +ρ = QuantumObject(z * B * z', dims = ntuple(i->N_cut, Val(N_modes))); # Full density matrix ``` Define the Hamiltonian and collapse operators @@ -67,26 +81,26 @@ Jx = 0.9 Jy = 1.04 Jz = 1.0 hx = 0.0 +hy = 0.0 +hz = 0.0 γ = 1 -Sx = sum([mb(sx, i, latt) for i in 1:latt.N]) -Sy = sum([mb(sy, i, latt) for i in 1:latt.N]) -Sz = sum([mb(sz, i, latt) for i in 1:latt.N]) -SFxx = sum([mb(sx, i, latt) * mb(sx, j, latt) for i in 1:latt.N for j in 1:latt.N]) +Sx = mapreduce(i->SingleSiteOperator(sigmax(), i, latt), +, 1:latt.N) +Sy = mapreduce(i->SingleSiteOperator(sigmay(), i, latt), +, 1:latt.N) +Sz = mapreduce(i->SingleSiteOperator(sigmaz(), i, latt), +, 1:latt.N) -H, c_ops = TFIM(Jx, Jy, Jz, hx, γ, latt; bc = pbc, order = 1) -e_ops = (Sx, Sy, Sz, SFxx) +H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1) +e_ops = (Sx, Sy, Sz) -tl = LinRange(0, 10, 100); +tl = range(0, 10, 100) +nothing # hide ``` ### Full evolution ```@example lowrank -@time mesol = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]); -A = Matrix(mesol.states[end].data) -λ = eigvals(Hermitian(A)) -Strue = -sum(λ .* log2.(λ)) / latt.N; +sol_me = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]); +Strue = entropy_vn(sol_me.states[end], base=2) / latt.N ``` ### Low Rank evolution @@ -120,26 +134,23 @@ function f_entropy(p, z, B) mul!(C, z, sqrt(B)) mul!(σ, C', C) - λ = eigvals(Hermitian(σ)) - λ = λ[λ.>1e-10] - return -sum(λ .* log2.(λ)) -end; + return entropy_vn(Qobj(Hermitian(σ), type=Operator), base=2) +end ``` Define the options for the low-rank evolution ```@example lowrank -opt = - LRMesolveOptions(err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0); +opt = (err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0); -@time lrsol = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_purity, f_entropy, f_trace), opt = opt); +sol_lr = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_purity, f_entropy, f_trace), opt = opt); ``` Plot the results ```@example lowrank -m_me = real(mesol.expect[3, :]) / Nx / Ny -m_lr = real(lrsol.expvals[3, :]) / Nx / Ny +m_me = real(sol_me.expect[3, :]) / Nx / Ny +m_lr = real(sol_lr.expect[3, :]) / Nx / Ny fig = Figure(size = (500, 350), fontsize = 15) ax = Axis(fig[1, 1], xlabel = L"\gamma t", ylabel = L"M_{z}", xlabelsize = 20, ylabelsize = 20) @@ -148,17 +159,17 @@ lines!(ax, tl, m_me, label = "Fock", linewidth = 2, linestyle = :dash) axislegend(ax, position = :rb) ax2 = Axis(fig[1, 2], xlabel = L"\gamma t", ylabel = "Value", xlabelsize = 20, ylabelsize = 20) -lines!(ax2, tl, 1 .- real(lrsol.funvals[1, :]), label = L"$1-P$", linewidth = 2) +lines!(ax2, tl, 1 .- real(sol_lr.fexpect[1, :]), label = L"$1-P$", linewidth = 2) lines!( ax2, tl, - 1 .- real(lrsol.funvals[3, :]), + 1 .- real(sol_lr.fexpect[3, :]), label = L"$1-\mathrm{Tr}(\rho)$", linewidth = 2, linestyle = :dash, color = :orange, ) -lines!(ax2, tl, real(lrsol.funvals[2, :]) / Nx / Ny, color = :blue, label = L"S", linewidth = 2) +lines!(ax2, tl, real(sol_lr.fexpect[2, :]) / Nx / Ny, color = :blue, label = L"S", linewidth = 2) hlines!(ax2, [Strue], color = :blue, linestyle = :dash, linewidth = 2, label = L"S^{\,\mathrm{true}}_{\mathrm{ss}}") axislegend(ax2, position = :rb) diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl index 26aa48ae0..192422e72 100644 --- a/src/spin_lattice.jl +++ b/src/spin_lattice.jl @@ -1,12 +1,10 @@ -export Lattice, mb, TFIM, nn, sx, sy, sz, sm, sp, pbc, obc +export Lattice, SingleSiteOperator, DissipativeIsing -sx = sigmax() -sy = -sigmay() -sz = -sigmaz() -sm = (sx - 1im * sy) / 2 -sp = (sx + 1im * sy) / 2 +@doc raw""" + Lattice -#Lattice structure +A Julia constructor for a lattice object. The lattice object is used to define the geometry of the lattice. `Nx` and `Ny` are the number of sites in the x and y directions, respectively. `N` is the total number of sites. `lin_idx` is a `LinearIndices` object and `car_idx` is a `CartesianIndices` object, and they are used to efficiently select sites on the lattice. +""" Base.@kwdef struct Lattice{TN<:Integer,TLI<:LinearIndices,TCI<:CartesianIndices} Nx::TN Ny::TN @@ -16,47 +14,100 @@ Base.@kwdef struct Lattice{TN<:Integer,TLI<:LinearIndices,TCI<:CartesianIndices} end #Definition of many-body operators -function mb(s::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, i::Integer, N::Integer) where {T1} - T = s.dims[1] - return QuantumObject(kron(eye(T^(i - 1)), s, eye(T^(N - i))); dims = ntuple(j -> 2, Val(N))) +@doc raw""" + SingleSiteOperator(O::QuantumObject, i::Integer, N::Integer) + +A Julia constructor for a single-site operator. `s` is the operator acting on the site. `i` is the site index, and `N` is the total number of sites. The function returns a `QuantumObject` given by ``\\mathbb{1}^{\\otimes (i - 1)} \\otimes \hat{O} \\otimes \\mathbb{1}^{\\otimes (N - i)}``. +""" +function SingleSiteOperator(O::QuantumObject{DT,OperatorQuantumObject}, i::Integer, N::Integer) where {DT} + T = O.dims[1] + return QuantumObject(kron(eye(T^(i - 1)), O, eye(T^(N - i))); dims = ntuple(j -> 2, Val(N))) end -mb(s::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, i::Integer, latt::Lattice) where {T1} = mb(s, i, latt.N) -mb(s::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, row::Integer, col::Integer, latt::Lattice) where {T1} = - mb(s, latt.idx[row, col], latt.N) -mb(s::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, x::CartesianIndex, latt::Lattice) where {T1} = - mb(s, latt.idx[x], latt.N) +SingleSiteOperator(O::QuantumObject{DT,OperatorQuantumObject}, i::Integer, latt::Lattice) where {DT} = + SingleSiteOperator(O, i, latt.N) +SingleSiteOperator(O::QuantumObject{DT,OperatorQuantumObject}, row::Integer, col::Integer, latt::Lattice) where {DT} = + SingleSiteOperator(O, latt.idx[row, col], latt.N) +SingleSiteOperator(O::QuantumObject{DT,OperatorQuantumObject}, x::CartesianIndex, latt::Lattice) where {DT} = + SingleSiteOperator(O, latt.idx[x], latt.N) #Definition of nearest-neighbour sites on lattice -pbc(i::Integer, N::Integer) = 1 + (i - 1 + N) % N -obc(i::Integer, N::Integer) = (i >= 1 && i <= N) -pbc(i::Vector{Int}, N::Integer) = pbc.(i, N) -obc(i::Vector{Int}, N::Integer) = filter(x -> obc(x, N), i) - -function nn(i::CartesianIndex, latt::Lattice, bc::Function; order::Integer = 1) - row = bc([i[1] + order, i[1] - order], latt.Nx) - col = bc([i[2] + order, i[2] - order], latt.Ny) +periodic_boundary_conditions(i::Integer, N::Integer) = 1 + (i - 1 + N) % N +open_boundary_conditions(i::Integer, N::Integer) = (i >= 1 && i <= N) +periodic_boundary_conditions(i::Vector{Int}, N::Integer) = periodic_boundary_conditions.(i, N) +open_boundary_conditions(i::Vector{Int}, N::Integer) = filter(x -> open_boundary_conditions(x, N), i) + +function nearest_neighbor(i::CartesianIndex, latt::Lattice, ::Val{:periodic_bc}; order::Integer = 1) + row = periodic_boundary_conditions([i[1] + order, i[1] - order], latt.Nx) + col = periodic_boundary_conditions([i[2] + order, i[2] - order], latt.Ny) + return vcat([CartesianIndex(r, i[2]) for r in row], [CartesianIndex(i[1], c) for c in col]) +end + +function nearest_neighbor(i::CartesianIndex, latt::Lattice, ::Val{:open_bc}; order::Integer = 1) + row = periodic_boundary_conditions([i[1] + order, i[1] - order], latt.Nx) + col = periodic_boundary_conditions([i[2] + order, i[2] - order], latt.Ny) return vcat([CartesianIndex(r, i[2]) for r in row], [CartesianIndex(i[1], c) for c in col]) end -function TFIM(Jx::Real, Jy::Real, Jz::Real, hx::Real, γ::Real, latt::Lattice; bc::Function = pbc, order::Integer = 1) - S = [mb(sm, i, latt) for i in 1:latt.N] +@doc """ + DissipativeIsing(Jx::Real, Jy::Real, Jz::Real, hx::Real, hy::Real, hz::Real, γ::Real, latt::Lattice; boundary_condition::Union{Symbol, Val} = Val(:periodic_bc), order::Integer = 1) + +A Julia constructor for a dissipative Ising model. The function returns the Hamiltonian + +```math +\\hat{H} = \\frac{J_x}{2} \\sum_{\\langle i, j \\rangle} \\hat{\\sigma}_i^x \\hat{\\sigma}_j^x + \\frac{J_y}{2} \\sum_{\\langle i, j \\rangle} \\hat{\\sigma}_i^y \\hat{\\sigma}_j^y + \\frac{J_z}{2} \\sum_{\\langle i, j \\rangle} \\hat{\\sigma}_i^z \\hat{\\sigma}_j^z + h_x \\sum_i \\hat{\\sigma}_i^x +``` + +and the collapse operators + +```math +\\hat{c}_i = \\sqrt{\\gamma} \\hat{\\sigma}_i^- +``` + +# Arguments +- `Jx::Real`: The coupling constant in the x-direction. +- `Jy::Real`: The coupling constant in the y-direction. +- `Jz::Real`: The coupling constant in the z-direction. +- `hx::Real`: The magnetic field in the x-direction. +- `hy::Real`: The magnetic field in the y-direction. +- `hz::Real`: The magnetic field in the z-direction. +- `γ::Real`: The local dissipation rate. +- `latt::Lattice`: A [`Lattice`](@ref) object that defines the geometry of the lattice. +- `boundary_condition::Union{Symbol, Val}`: The boundary conditions of the lattice. The possible inputs are `periodic_bc` and `open_bc`, for periodic or open boundary conditions, respectively. The default value is `Val(:periodic_bc)`. +- `order::Integer`: The order of the nearest-neighbour sites. The default value is 1. +""" +function DissipativeIsing( + Jx::Real, + Jy::Real, + Jz::Real, + hx::Real, + hy::Real, + hz::Real, + γ::Real, + latt::Lattice; + boundary_condition::Union{Symbol,Val} = Val(:periodic_bc), + order::Integer = 1, +) + S = [SingleSiteOperator(sigmam(), i, latt) for i in 1:latt.N] c_ops = sqrt(γ) .* S - op_sum(S, i::CartesianIndex) = S[latt.lin_idx[i]] * sum(S[latt.lin_idx[nn(i, latt, bc; order = order)]]) + op_sum(S, i::CartesianIndex) = + S[latt.lin_idx[i]] * sum(S[latt.lin_idx[nearest_neighbor(i, latt, makeVal(boundary_condition); order = order)]]) H = 0 if (Jx != 0 || hx != 0) - S .= [mb(sx, i, latt) for i in 1:latt.N] + S = [SingleSiteOperator(sigmax(), i, latt) for i in 1:latt.N] H += Jx / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx) #/2 because we are double counting H += hx * sum(S) end - if Jy != 0 - S .= [mb(sy, i, latt) for i in 1:latt.N] + if (Jy != 0 || hy != 0) + S = [SingleSiteOperator(sigmay(), i, latt) for i in 1:latt.N] H += Jy / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx) + H += hy * sum(S) end - if Jz != 0 - S .= [mb(sz, i, latt) for i in 1:latt.N] + if (Jz != 0 || hz != 0) + S = [SingleSiteOperator(sigmaz(), i, latt) for i in 1:latt.N] H += Jz / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx) + H += hz * sum(S) end return H, c_ops -end; +end diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl index 4d530bb36..8b71e94d9 100644 --- a/src/time_evolution/lr_mesolve.jl +++ b/src/time_evolution/lr_mesolve.jl @@ -1,56 +1,59 @@ -export lr_mesolve, lr_mesolveProblem, LRTimeEvolutionSol, LRMesolveOptions +export lr_mesolve, lr_mesolveProblem, TimeEvolutionLRSol -#=======================================================# -# STRUCT DEFINITIONS -#=======================================================# - -struct LRTimeEvolutionSol{TT<:Vector{<:Real},TS<:AbstractVector,TE<:Matrix{ComplexF64},TM<:Vector{<:Integer}} +@doc raw""" + struct TimeEvolutionLRSol + +A structure storing the results and some information from solving low-rank master equation time evolution. + +# Fields (Attributes) + +- `times::AbstractVector`: The time list of the evolution. +- `states::Vector{QuantumObject}`: The list of result states. +- `expect::Matrix`: The expectation values corresponding to each time point in `times`. +- `fexpect::Matrix`: The function values at each time point. +- `retcode`: The return code from the solver. +- `alg`: The algorithm which is used during the solving process. +- `abstol::Real`: The absolute tolerance which is used during the solving process. +- `reltol::Real`: The relative tolerance which is used during the solving process. +- `z::Vector{QuantumObject}`: The `z`` matrix of the low-rank algorithm at each time point. +- `B::Vector{QuantumObject}`: The `B` matrix of the low-rank algorithm at each time point. +""" +struct TimeEvolutionLRSol{ + TT<:AbstractVector{<:Real}, + TS<:AbstractVector, + TE<:Matrix{ComplexF64}, + RetT<:Enum, + AlgT<:OrdinaryDiffEqAlgorithm, + AT<:Real, + RT<:Real, + TSZB<:AbstractVector, + TM<:Vector{<:Integer}, +} times::TT - z::TS - B::TS - expvals::TE - funvals::TE + states::TS + expect::TE + fexpect::TE + retcode::RetT + alg::AlgT + abstol::AT + reltol::RT + z::TSZB + B::TSZB M::TM end -struct LRMesolveOptions{AlgType<:OrdinaryDiffEqAlgorithm} - alg::AlgType - progress::Bool - err_max::Real - p0::Real - atol_inv::Real - M_max::Integer - compute_Si::Bool - is_dynamical::Bool - adj_condition::String - Δt::Real -end - -function LRMesolveOptions(; - alg::OrdinaryDiffEqAlgorithm = Tsit5(), - progress::Bool = true, - err_max::Real = 0.0, - p0::Real = 0.0, - atol_inv::Real = 1e-4, - M_max::Integer = typemax(Int), - compute_Si::Bool = true, - _is_dynamical::Bool = err_max > 0, - adj_condition::String = "variational", - Δt::Real = 0.0, +lr_mesolve_options_default = ( + alg = Tsit5(), + progress = true, + err_max = 0.0, + p0 = 0.0, + atol_inv = 1e-4, + M_max = typemax(Int), + compute_Si = true, + is_dynamical = false, + adj_condition = "variational", + Δt = 0.0, ) - return LRMesolveOptions{typeof(alg)}( - alg, - progress, - err_max, - p0, - atol_inv, - M_max, - compute_Si, - _is_dynamical, - adj_condition, - Δt, - ) -end #=======================================================# # ADDITIONAL FUNCTIONS @@ -58,25 +61,22 @@ end select(x::Real, xarr::AbstractArray, retval = false) = retval ? xarr[argmin(abs.(x .- xarr))] : argmin(abs.(x .- xarr)) -@doc raw""" - _pinv!(A, T1, T2; atol::Real=0.0, rtol::Real=(eps(real(float(oneunit(T))))*min(size(A)...))*iszero(atol)) where T - Computes the pseudo-inverse of a matrix A, and stores it in T1. If T2 is provided, it is used as a temporary matrix. - The algorithm is based on the SVD decomposition of A, and is taken from the Julia package LinearAlgebra. - The difference with respect to the original function is that the cutoff is done with a smooth function instead of a step function. - - Parameters - ---------- - A : AbstractMatrix{T} - The matrix to be inverted. - T1 : AbstractMatrix{T} - T2 : AbstractMatrix{T} - Temporary matrices used in the calculation. - atol : Real - Absolute tolerance for the calculation of the pseudo-inverse. - rtol : Real - Relative tolerance for the calculation of the pseudo-inverse. -""" -function _pinv!( +#= + _pinv_smooth!(A, T1, T2; atol::Real=0.0, rtol::Real=(eps(real(float(oneunit(T))))*min(size(A)...))*iszero(atol)) where T + +Computes the pseudo-inverse of a matrix A, and stores it in T1. If T2 is provided, it is used as a temporary matrix. +The algorithm is based on the SVD decomposition of A, and is taken from the Julia package LinearAlgebra. +The difference with respect to the original function is that the cutoff is done with a smooth function instead of a step function. + +# Arguments + + - `A::AbstractMatrix`: The matrix to be inverted. + - `T1::AbstractMatrix`: The matrix where the pseudo-inverse is stored. + - `T2::AbstractMatrix`: A temporary matrix. + - `atol::Real`: The absolute tolerance. + - `rtol::Real`: The relative tolerance. +=# +function _pinv_smooth!( A::AbstractMatrix{T}, T1::AbstractMatrix{T}, T2::AbstractMatrix{T}; @@ -98,22 +98,19 @@ function _pinv!( return mul!(T1, SVD.Vt', T2) end -@doc raw""" - _calculate_expectation!(p,z,B,idx) where T - Calculates the expectation values and function values of the operators and functions in p.e_ops and p.f_ops, respectively, and stores them in p.expvals and p.funvals. - The function is called by the callback _save_affect_lr_mesolve!. - - Parameters - ---------- - p : NamedTuple - The parameters of the problem. - z : AbstractMatrix{T} - The z matrix. - B : AbstractMatrix{T} - The B matrix. - idx : Integer - The index of the current time step. -""" +#= + _calculate_expectation!(p,z,B,idx) + +Calculates the expectation values and function values of the operators and functions in p.e_ops and p.f_ops, respectively, and stores them in p.expvals and p.funvals. +The function is called by the callback _save_affect_lr_mesolve!. + +# Arguments + + - `p::NamedTuple`: The parameters of the problem. + - `z::AbstractMatrix`: The z matrix of the low-rank algorithm. + - `B::AbstractMatrix`: The B matrix of the low-rank algorithm. + - `idx::Integer`: The index of the current time step. +=# function _calculate_expectation!(p, z, B, idx) e_ops = p.e_ops f_ops = p.f_ops @@ -169,20 +166,18 @@ end # CALLBACK FUNCTIONS #=======================================================# -@doc raw""" - _adjM_condition_ratio(u, t, integrator) where T - Condition for the dynamical rank adjustment based on the ratio between the smallest and largest eigenvalues of the density matrix. - The spectrum of the density matrix is calculated efficiently using the properties of the SVD decomposition of the matrix. - - Parameters - ---------- - u : AbstractVector{T} - The current state of the system. - t : Real - The current time. - integrator : ODEIntegrator - The integrator of the problem. -""" +#= + _adjM_condition_ratio(u, t, integrator) + +Condition for the dynamical rank adjustment based on the ratio between the smallest and largest eigenvalues of the density matrix. +The spectrum of the density matrix is calculated efficiently using the properties of the SVD decomposition of the matrix. + +# Arguments + + - `u::AbstractVector`: The current state of the system. + - `t::Real`: The current time. + - `integrator::ODEIntegrator`: The integrator of the problem. +=# function _adjM_condition_ratio(u, t, integrator) ip = integrator.p opt = ip.opt @@ -200,19 +195,17 @@ function _adjM_condition_ratio(u, t, integrator) return (err >= opt.err_max && M < N && M < opt.M_max) end -@doc raw""" - _adjM_condition_variational(u, t, integrator) where T - Condition for the dynamical rank adjustment based on the leakage out of the low-rank manifold. - - Parameters - ---------- - u : AbstractVector{T} - The current state of the system. - t : Real - The current time. - integrator : ODEIntegrator - The integrator of the problem. -""" +#= + _adjM_condition_variational(u, t, integrator) + +Condition for the dynamical rank adjustment based on the leakage out of the low-rank manifold. + +# Arguments + + - `u::AbstractVector`: The current state of the system. + - `t::Real`: The current time. + - `integrator::ODEIntegrator`: The integrator of the problem. +=# function _adjM_condition_variational(u, t, integrator) ip = integrator.p opt = ip.opt @@ -222,16 +215,16 @@ function _adjM_condition_variational(u, t, integrator) return (err >= opt.err_max && M < N && M < opt.M_max) end -@doc raw""" +#= _adjM_affect!(integrator) - Affect function for the dynamical rank adjustment. It increases the rank of the low-rank manifold by one, and updates the matrices accordingly. - If Δt>0, it rewinds the integrator to the previous time step. - Parameters - ---------- - integrator : ODEIntegrator - The integrator of the problem. -""" +Affect function for the dynamical rank adjustment. It increases the rank of the low-rank manifold by one, and updates the matrices accordingly. +If Δt>0, it rewinds the integrator to the previous time step. + +# Arguments + + - `integrator::ODEIntegrator`: The integrator of the problem. +=# function _adjM_affect!(integrator) ip = integrator.p opt = ip.opt @@ -265,8 +258,10 @@ function _adjM_affect!(integrator) ), ) mul!(integrator.p.S, z', z) - !(opt.compute_Si) && - (integrator.p.Si .= _pinv!(Hermitian(integrator.p.S), integrator.temp_MM, integrator.L, atol = opt.atol_inv)) + !(opt.compute_Si) && ( + integrator.p.Si .= + _pinv_smooth!(Hermitian(integrator.p.S), integrator.temp_MM, integrator.L, atol = opt.atol_inv) + ) if Δt > 0 integrator.p = merge(integrator.p, (u_save = copy(integrator.u),)) @@ -286,21 +281,18 @@ end # DYNAMICAL EVOLUTION EQUATIONS #=======================================================# -@doc raw""" - dBdz!(du, u, p, t) where T - Dynamical evolution equations for the low-rank manifold. The function is called by the ODEProblem. - - Parameters - ---------- - du : AbstractVector{T} - The derivative of the state of the system. - u : AbstractVector{T} - The current state of the system. - p : NamedTuple - The parameters of the problem. - t : Real - The current time. -""" +#= + dBdz!(du, u, p, t) + +Dynamical evolution equations for the low-rank manifold. The function is called by the ODEProblem. + +# Arguments + + - `du::AbstractVector`: The derivative of the state vector. + - `u::AbstractVector`: The state vector. + - `p::NamedTuple`: The parameters of the problem. + - `t::Real`: The current time. +=# function dBdz!(du, u, p, t) #NxN H, Γ = p.H, p.Γ @@ -326,8 +318,8 @@ function dBdz!(du, u, p, t) mul!(A0, z, B) # Calculate inverse - opt.compute_Si && (Si .= _pinv!(Hermitian(S), temp_MM, L, atol = opt.atol_inv)) - Bi .= _pinv!(Hermitian(B), temp_MM, L, atol = opt.atol_inv) + opt.compute_Si && (Si .= _pinv_smooth!(Hermitian(S), temp_MM, L, atol = opt.atol_inv)) + Bi .= _pinv_smooth!(Hermitian(B), temp_MM, L, atol = opt.atol_inv) # Calculate the effective Hamiltonian part of L_tilde mul!(dz, H, A0) @@ -360,49 +352,59 @@ function dBdz!(du, u, p, t) return dB .-= temp_MM end +#=======================================================# +# OUTPUT GENNERATION +#=======================================================# + +get_z(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, 1:M*N), N, M) + +get_B(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, (M*N+1):length(u)), M, M) + #=======================================================# # PROBLEM FORMULATION #=======================================================# @doc raw""" - lr_mesolveProblem(H, z, B, tlist, c_ops; e_ops=(), f_ops=(), opt=LRMesolveOptions(), kwargs...) where T - Formulates the ODEproblem for the low-rank time evolution of the system. The function is called by lr_mesolve. - - Parameters - ---------- - H : QuantumObject - The Hamiltonian of the system. - z : AbstractMatrix{T} - The initial z matrix. - B : AbstractMatrix{T} - The initial B matrix. - tlist : AbstractVector{T} - The time steps at which the expectation values and function values are calculated. - c_ops : AbstractVector{QuantumObject} - The jump operators of the system. - e_ops : Tuple{QuantumObject} - The operators whose expectation values are calculated. - f_ops : Tuple{Function} - The functions whose values are calculated. - opt : LRMesolveOptions - The options of the problem. - kwargs : NamedTuple - Additional keyword arguments for the ODEProblem. + lr_mesolveProblem( + H::QuantumObject{DT1,OperatorQuantumObject}, + z::AbstractArray{T2,2}, + B::AbstractArray{T2,2}, + tlist::AbstractVector, + c_ops::Union{AbstractVector,Tuple}=(); + e_ops::Union{AbstractVector,Tuple}=(), + f_ops::Union{AbstractVector,Tuple}=(), + opt::NamedTuple = lr_mesolve_options_default, + kwargs..., + ) + +Formulates the ODEproblem for the low-rank time evolution of the system. The function is called by [`lr_mesolve`](@ref). For more information about the low-rank master equation, see [gravina2024adaptive](@cite). + +# Arguments +- `H::QuantumObject`: The Hamiltonian of the system. +- `z::AbstractArray`: The initial z matrix of the low-rank algorithm. +- `B::AbstractArray`: The initial B matrix of the low-rank algorithm. +- `tlist::AbstractVector`: The time steps at which the expectation values and function values are calculated. +- `c_ops::Union{AbstractVector,Tuple}`: The list of the collapse operators. +- `e_ops::Union{AbstractVector,Tuple}`: The list of the operators for which the expectation values are calculated. +- `f_ops::Union{AbstractVector,Tuple}`: The list of the functions for which the function values are calculated. +- `opt::NamedTuple`: The options of the low-rank master equation. +- `kwargs`: Additional keyword arguments. """ function lr_mesolveProblem( - H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, + H::QuantumObject{DT1,OperatorQuantumObject}, z::AbstractArray{T2,2}, B::AbstractArray{T2,2}, tlist::AbstractVector, - c_ops::AbstractVector = []; - e_ops::Tuple = (), - f_ops::Tuple = (), - opt::LRMesolveOptions{AlgType} = LRMesolveOptions(), + c_ops::Union{AbstractVector,Tuple} = (); + e_ops::Union{AbstractVector,Tuple} = (), + f_ops::Union{AbstractVector,Tuple} = (), + opt::NamedTuple = lr_mesolve_options_default, kwargs..., -) where {T1,T2,AlgType<:OrdinaryDiffEqAlgorithm} +) where {DT1,T2} + Hdims = H.dims # Formulation of problem - H -= 0.5im * sum([Γ' * Γ for Γ in c_ops]) + H -= 0.5im * mapreduce(op -> op' * op, +, c_ops) H = get_data(H) c_ops = get_data.(c_ops) e_ops = get_data.(e_ops) @@ -414,6 +416,11 @@ function lr_mesolveProblem( funvals = Array{ComplexF64}(undef, length(f_ops), length(t_l)) Ml = Array{Int64}(undef, length(t_l)) + opt = merge(lr_mesolve_options_default, opt) + if opt.err_max > 0 + opt = merge(opt, (is_dynamical = true,)) + end + # Initialization of parameters. Scalars represents in order: Tr(S^{-1}L), t0 p = ( N = size(z, 1), @@ -436,6 +443,7 @@ function lr_mesolveProblem( Si = similar(B), u_save = vcat(vec(z), vec(B)), scalars = [0.0, t_l[1]], + Hdims = Hdims, ) mul!(p.S, z', z) @@ -484,43 +492,50 @@ function lr_mesolveProblem( # Initialization of ODEProblem tspan = (t_l[1], t_l[end]) - return ODEProblem(dBdz!, p.u_save, tspan, p; kwargs2...) + return ODEProblem{true}(dBdz!, p.u_save, tspan, p; kwargs2...) end +@doc raw""" + lr_mesolve( + H::QuantumObject{DT1,OperatorQuantumObject}, + z::AbstractArray{T2,2}, + B::AbstractArray{T2,2}, + tlist::AbstractVector, + c_ops::Union{AbstractVector,Tuple}=(); + e_ops::Union{AbstractVector,Tuple}=(), + f_ops::Union{AbstractVector,Tuple}=(), + opt::NamedTuple = lr_mesolve_options_default, + kwargs..., + ) + +Time evolution of an open quantum system using the low-rank master equation. For more information about the low-rank master equation, see [gravina2024adaptive](@cite). + +# Arguments +- `H::QuantumObject`: The Hamiltonian of the system. +- `z::AbstractArray`: The initial z matrix of the low-rank algorithm. +- `B::AbstractArray`: The initial B matrix of the low-rank algorithm. +- `tlist::AbstractVector`: The time steps at which the expectation values and function values are calculated. +- `c_ops::Union{AbstractVector,Tuple}`: The list of the collapse operators. +- `e_ops::Union{AbstractVector,Tuple}`: The list of the operators for which the expectation values are calculated. +- `f_ops::Union{AbstractVector,Tuple}`: The list of the functions for which the function values are calculated. +- `opt::NamedTuple`: The options of the low-rank master equation. +- `kwargs`: Additional keyword arguments. +""" function lr_mesolve( - H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, + H::QuantumObject{DT1,OperatorQuantumObject}, z::AbstractArray{T2,2}, B::AbstractArray{T2,2}, tlist::AbstractVector, - c_ops::AbstractVector = []; - e_ops::Tuple = (), - f_ops::Tuple = (), - opt::LRMesolveOptions{AlgType} = LRMesolveOptions(), + c_ops::Union{AbstractVector,Tuple} = (); + e_ops::Union{AbstractVector,Tuple} = (), + f_ops::Union{AbstractVector,Tuple} = (), + opt::NamedTuple = lr_mesolve_options_default, kwargs..., -) where {T1,T2,AlgType<:OrdinaryDiffEqAlgorithm} +) where {DT1,T2} prob = lr_mesolveProblem(H, z, B, tlist, c_ops; e_ops = e_ops, f_ops = f_ops, opt = opt, kwargs...) return lr_mesolve(prob; kwargs...) end -#=======================================================# -# OUTPUT GENNERATION -#=======================================================# - -get_z(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, 1:M*N), N, M) - -get_B(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, (M*N+1):length(u)), M, M) - -@doc raw""" - lr_mesolve(prob::ODEProblem; kwargs...) - Solves the ODEProblem formulated by lr_mesolveProblem. The function is called by lr_mesolve. - - Parameters - ---------- - prob : ODEProblem - The ODEProblem formulated by lr_mesolveProblem. - kwargs : NamedTuple - Additional keyword arguments for the ODEProblem. -""" function lr_mesolve(prob::ODEProblem; kwargs...) sol = solve(prob, prob.p.opt.alg, tstops = prob.p.times) prob.p.opt.progress && print("\n") @@ -529,13 +544,21 @@ function lr_mesolve(prob::ODEProblem; kwargs...) Ll = length.(sol.u) Ml = @. Int((sqrt(N^2 + 4 * Ll) - N) / 2) - if !haskey(kwargs, :saveat) - Bt = map(x -> get_B(x[1], N, x[2]), zip(sol.u, Ml)) - zt = map(x -> get_z(x[1], N, x[2]), zip(sol.u, Ml)) - else - Bt = get_B(sol.u, N, Ml) - zt = get_z(sol.u, N, Ml) - end - - return LRTimeEvolutionSol(sol.prob.p.times, zt, Bt, prob.p.expvals, prob.p.funvals, prob.p.Ml) + Bt = map(x -> get_B(x[1], N, x[2]), zip(sol.u, Ml)) + zt = map(x -> get_z(x[1], N, x[2]), zip(sol.u, Ml)) + ρt = map(x -> Qobj(x[1] * x[2] * x[1]', type = Operator, dims = prob.p.Hdims), zip(zt, Bt)) + + return TimeEvolutionLRSol( + sol.t, + ρt, + prob.p.expvals, + prob.p.funvals, + sol.retcode, + prob.p.opt.alg, + sol.prob.kwargs[:abstol], + sol.prob.kwargs[:reltol], + zt, + Bt, + Ml, + ) end diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl index 99dd33050..2da7dc686 100644 --- a/test/core-test/low_rank_dynamics.jl +++ b/test/core-test/low_rank_dynamics.jl @@ -2,79 +2,82 @@ # Define lattice Nx, Ny = 2, 3 latt = Lattice(Nx = Nx, Ny = Ny) - N_cut = 2 - N_modes = latt.N - N = N_cut^N_modes - M = Nx * Ny + 1 + ## + N_cut = 2 # Number of states of each mode + N_modes = latt.N # Number of modes + N = N_cut^N_modes # Total number of states + M = latt.N + 1 # Number of states in the LR basis # Define initial state ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,M - 1}}(undef, M) - ϕ[1] = tensor(repeat([basis(2, 0)], N_modes)...) + ϕ[1] = kron(fill(basis(2, 1), N_modes)...) + i = 1 for j in 1:N_modes i += 1 - i <= M && (ϕ[i] = mb(sp, j, latt) * ϕ[1]) + i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), j, latt) * ϕ[1]) end for k in 1:N_modes-1 for l in k+1:N_modes i += 1 - i <= M && (ϕ[i] = mb(sp, k, latt) * mb(sp, l, latt) * ϕ[1]) + i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), k, latt) * SingleSiteOperator(sigmap(), l, latt) * ϕ[1]) end end for i in i+1:M - ϕ[i] = Qobj(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims) + ϕ[i] = QuantumObject(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims) normalize!(ϕ[i]) end - z = hcat(broadcast(x -> x.data, ϕ)...) + + z = hcat(get_data.(ϕ)...) B = Matrix(Diagonal([1 + 0im; zeros(M - 1)])) - S = z' * z - B = B / tr(S * B) - ρ = Qobj(z * B * z', dims = ntuple(i -> 1, Val(N_modes)) .* N_cut) + S = z' * z # Overlap matrix + B = B / tr(S * B) # Normalize B + ρ = Qobj(z * B * z', dims = ntuple(i -> N_cut, Val(N_modes))) # Full density matrix # Define Hamiltonian and collapse operators Jx = 0.9 - Jy = 1.02 + Jy = 1.04 Jz = 1.0 hx = 0.0 + hy = 0.0 + hz = 0.0 γ = 1 - Sz = sum([mb(sz, i, latt) for i in 1:latt.N]) - tl = LinRange(0, 10, 100) - H, c_ops = TFIM(Jx, Jy, Jz, hx, γ, latt; bc = pbc, order = 1) + Sx = mapreduce(i -> SingleSiteOperator(sigmax(), i, latt), +, 1:latt.N) + Sy = mapreduce(i -> SingleSiteOperator(sigmay(), i, latt), +, 1:latt.N) + Sz = mapreduce(i -> SingleSiteOperator(sigmaz(), i, latt), +, 1:latt.N) + SFxx = mapreduce( + x -> SingleSiteOperator(sigmax(), x[1], latt) * SingleSiteOperator(sigmax(), x[2], latt), + +, + Iterators.product(1:latt.N, 1:latt.N), + ) + + H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1) e_ops = (Sz,) + tl = range(0, 10, 100) + # Full solution - mesol = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...], progress_bar = Val(false)) - A = Matrix(mesol.states[end].data) - λ = eigvals(Hermitian(A)) - Strue = -sum(λ .* log2.(λ)) + sol_me = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]) + Strue = entropy_vn(sol_me.states[end], base = 2) / latt.N # Low rank solution function f_entropy(p, z, B) C = p.A0 σ = p.Bi + mul!(C, z, sqrt(B)) mul!(σ, C', C) - λ = eigvals(Hermitian(σ)) - λ = λ[λ.>1e-10] - return -sum(λ .* log2.(λ)) + return entropy_vn(Qobj(Hermitian(σ), type = Operator), base = 2) end - opt = LRMesolveOptions( - err_max = 1e-3, - p0 = 0.0, - atol_inv = 1e-6, - adj_condition = "variational", - Δt = 0.2, - progress = false, - ) - lrsol = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_entropy,), opt = opt) + opt = (err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0, progress = false) + + sol_lr = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_entropy,), opt = opt) # Test - m_me = real(mesol.expect[1, :]) - m_lr = real(lrsol.expvals[1, :]) - @test all(abs.((m_me .- m_lr) ./ m_me) .< 0.1) + S_lr = real(sol_lr.fexpect[1, end]) / latt.N - S_lr = real(lrsol.funvals[1, end]) - @test abs((S_lr - Strue) / Strue) < 0.5 + @test real(sol_me.expect[1, :]) ≈ real(sol_lr.expect[1, :]) atol = 1e-1 + @test S_lr ≈ Strue atol = 1e-1 end