From 4c7f212ccf1cca80c8108f4b12bdc028586fa436 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Fri, 29 Nov 2024 15:43:56 +0100 Subject: [PATCH 1/2] Change SingleSiteOperator with MultiSiteOperator --- docs/src/resources/api.md | 2 +- docs/src/tutorials/cluster.md | 8 +-- docs/src/tutorials/lowrank.md | 10 ++-- src/spin_lattice.jl | 79 +++++++++++++++++++++++------ test/core-test/low_rank_dynamics.jl | 12 ++--- 5 files changed, 79 insertions(+), 32 deletions(-) diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md index 0b0c37a5..d19c4b86 100644 --- a/docs/src/resources/api.md +++ b/docs/src/resources/api.md @@ -248,7 +248,7 @@ fidelity ```@docs Lattice -SingleSiteOperator +MultiSiteOperator DissipativeIsing ``` diff --git a/docs/src/tutorials/cluster.md b/docs/src/tutorials/cluster.md index e14c86a2..5bd752b9 100644 --- a/docs/src/tutorials/cluster.md +++ b/docs/src/tutorials/cluster.md @@ -11,7 +11,7 @@ We will consider two examples: ### Monte Carlo Quantum Trajectories -Let's consider a 2-dimensional transferse field Ising model with 4x3 spins. The Hamiltonian is given by +Let's consider a 2-dimensional transverse field Ising model with 4x3 spins. The Hamiltonian is given by ```math \hat{H} = \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 \, , @@ -91,9 +91,9 @@ hy = 0.0 hz = 0.0 γ = 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) +Sx = mapreduce(i -> MultiSiteOperator(latt, i=>sigmax()), +, 1:latt.N) +Sy = mapreduce(i -> MultiSiteOperator(latt, i=>sigmay()), +, 1:latt.N) +Sz = mapreduce(i -> MultiSiteOperator(latt, i=>sigmaz()), +, 1:latt.N) H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1) e_ops = [Sx, Sy, Sz] diff --git a/docs/src/tutorials/lowrank.md b/docs/src/tutorials/lowrank.md index 59e292bf..d7084600 100644 --- a/docs/src/tutorials/lowrank.md +++ b/docs/src/tutorials/lowrank.md @@ -47,12 +47,12 @@ Define lr states. Take as initial state all spins up. All other N states are tak i = 1 for j in 1:N_modes global i += 1 - i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), j, latt) * ϕ[1]) + i <= M && (ϕ[i] = MultiSiteOperator(latt, j=>sigmap()) * ϕ[1]) end for k in 1:N_modes-1 for l in k+1:N_modes global i += 1 - i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), k, latt) * SingleSiteOperator(sigmap(), l, latt) * ϕ[1]) + i <= M && (ϕ[i] = MultiSiteOperator(latt, k=>sigmap(), l=>sigmap()) * ϕ[1]) end end for i in i+1:M @@ -85,9 +85,9 @@ hy = 0.0 hz = 0.0 γ = 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) +Sx = mapreduce(i->MultiSiteOperator(latt, i=>sigmax()), +, 1:latt.N) +Sy = mapreduce(i->MultiSiteOperator(latt, i=>sigmay()), +, 1:latt.N) +Sz = mapreduce(i->MultiSiteOperator(latt, i=>sigmaz()), +, 1:latt.N) H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1) e_ops = (Sx, Sy, Sz) diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl index 192422e7..9ce0094a 100644 --- a/src/spin_lattice.jl +++ b/src/spin_lattice.jl @@ -1,4 +1,4 @@ -export Lattice, SingleSiteOperator, DissipativeIsing +export Lattice, MultiSiteOperator, DissipativeIsing @doc raw""" Lattice @@ -15,20 +15,60 @@ end #Definition of many-body operators @doc raw""" - SingleSiteOperator(O::QuantumObject, i::Integer, N::Integer) + MultiSiteOperator(dims::Union{AbstractVector, Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...) -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)}``. +A Julia function for generating a multi-site operator ``\\hat{O} = \\hat{O}_i \\hat{O}_j \\cdots \\hat{O}_k``. The function takes a vector of dimensions `dims` and a list of pairs `pairs` where the first element of the pair is the site index and the second element is the operator acting on that site. + +# Arguments +- `dims::Union{AbstractVector, Tuple}`: A vector of dimensions of the lattice. +- `pairs::Pair{<:Integer,<:QuantumObject}...`: A list of pairs where the first element of the pair is the site index and the second element is the operator acting on that site. + +# Returns +`QuantumObject`: A `QuantumObject` representing the multi-site operator. + +# Example +```jldoctest +julia> op = MultiSiteOperator(Val(8), 5=>sigmax(), 7=>sigmaz()); + +julia> op.dims +8-element SVector{8, Int64} with indices SOneTo(8): + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 +``` """ -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))) +function MultiSiteOperator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...) + sites_unsorted = collect(first.(pairs)) + idxs = sortperm(sites_unsorted) + _sites = sites_unsorted[idxs] + _ops = collect(last.(pairs))[idxs] + _dims = collect(dims) # Use this instead of a Tuple, to avoid type instability when indexing on a slice + + sites, ops = _get_unique_sites_ops(_sites, _ops) + + collect(dims)[sites] == [op.dims[1] for op in ops] || throw(ArgumentError("The dimensions of the operators do not match the dimensions of the lattice.")) + + data = kron(I(prod(_dims[1:sites[1]-1])), ops[1].data) + for i in 2:length(sites) + data = kron(data, I(prod(_dims[sites[i-1]+1:sites[i]-1])), ops[i].data) + end + data = kron(data, I(prod(_dims[sites[end]+1:end]))) + + return QuantumObject(data; type = Operator, dims = dims) +end +function MultiSiteOperator(N::Union{Integer,Val}, pairs::Pair{<:Integer,<:QuantumObject}...) + dims = ntuple(j -> 2, makeVal(N)) + + return MultiSiteOperator(dims, pairs...) +end +function MultiSiteOperator(latt::Lattice, pairs::Pair{<:Integer,<:QuantumObject}...) + return MultiSiteOperator(makeVal(latt.N), pairs...) end -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 periodic_boundary_conditions(i::Integer, N::Integer) = 1 + (i - 1 + N) % N @@ -87,7 +127,7 @@ function DissipativeIsing( boundary_condition::Union{Symbol,Val} = Val(:periodic_bc), order::Integer = 1, ) - S = [SingleSiteOperator(sigmam(), i, latt) for i in 1:latt.N] + S = [MultiSiteOperator(latt, i => sigmam()) for i in 1:latt.N] c_ops = sqrt(γ) .* S op_sum(S, i::CartesianIndex) = @@ -95,19 +135,26 @@ function DissipativeIsing( H = 0 if (Jx != 0 || hx != 0) - S = [SingleSiteOperator(sigmax(), i, latt) for i in 1:latt.N] + S = [MultiSiteOperator(latt, i => sigmax()) 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 || hy != 0) - S = [SingleSiteOperator(sigmay(), i, latt) for i in 1:latt.N] + S = [MultiSiteOperator(latt, i => sigmay()) 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 || hz != 0) - S = [SingleSiteOperator(sigmaz(), i, latt) for i in 1:latt.N] + S = [MultiSiteOperator(latt, i => sigmaz()) 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 + +function _get_unique_sites_ops(sites, ops) + unique_sites = unique(sites) + unique_ops = map(i -> prod(ops[findall(==(i), sites)]), unique_sites) + + return unique_sites, unique_ops +end diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl index 2da7dc68..ef767f30 100644 --- a/test/core-test/low_rank_dynamics.jl +++ b/test/core-test/low_rank_dynamics.jl @@ -15,12 +15,12 @@ i = 1 for j in 1:N_modes i += 1 - i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), j, latt) * ϕ[1]) + i <= M && (ϕ[i] = MultiSiteOperator(latt, j => sigmap()) * ϕ[1]) end for k in 1:N_modes-1 for l in k+1:N_modes i += 1 - i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), k, latt) * SingleSiteOperator(sigmap(), l, latt) * ϕ[1]) + i <= M && (ϕ[i] = MultiSiteOperator(latt, k => sigmap(), l => sigmap()) * ϕ[1]) end end for i in i+1:M @@ -43,11 +43,11 @@ hz = 0.0 γ = 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) + Sx = mapreduce(i -> MultiSiteOperator(latt, i => sigmax()), +, 1:latt.N) + Sy = mapreduce(i -> MultiSiteOperator(latt, i => sigmay()), +, 1:latt.N) + Sz = mapreduce(i -> MultiSiteOperator(latt, i => sigmaz()), +, 1:latt.N) SFxx = mapreduce( - x -> SingleSiteOperator(sigmax(), x[1], latt) * SingleSiteOperator(sigmax(), x[2], latt), + x -> MultiSiteOperator(latt, x[1] => sigmax(), x[2] => sigmax()), +, Iterators.product(1:latt.N, 1:latt.N), ) From e2520aee911d53ba8d9128c663f3ff2ac6833d17 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Fri, 29 Nov 2024 15:47:34 +0100 Subject: [PATCH 2/2] Update changelog --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index daa4b6a1..9c59465b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## Unreleased +- Change `SingleSiteOperator` with the more general `MultiSiteOperator`. ([#324]) + ## [v0.22.0] (2024-11-20) - Change the parameters structure of `sesolve`, `mesolve` and `mcsolve` functions to possibly support automatic differentiation. ([#311]) @@ -33,3 +35,4 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [#311]: https://github.com/qutip/QuantumToolbox.jl/issues/311 [#315]: https://github.com/qutip/QuantumToolbox.jl/issues/315 [#318]: https://github.com/qutip/QuantumToolbox.jl/issues/318 +[#324]: https://github.com/qutip/QuantumToolbox.jl/issues/324