Clean low-rank mesolve and add Documentation bibliography

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"

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")

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!
 ```

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}
+}

docs/src/bibliography.md

@@ -0,0 +1,8 @@
+```@meta
+CurrentModule = QuantumToolbox
+```
+
+# [Bibliography](@id doc:Bibliography)
+
+```@bibliography
+```

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)