diff --git a/Project.toml b/Project.toml index 6cd6f76b..29a029e1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "QuantumOpticsBase" uuid = "4f57444f-1401-5e15-980d-4471b28d5678" -version = "0.4.17" +version = "0.4.18" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/manybody.jl b/src/manybody.jl index a4c8a9bf..dfb7f25a 100644 --- a/src/manybody.jl +++ b/src/manybody.jl @@ -1,3 +1,32 @@ +struct SortedVector{T, OT} <: AbstractVector{T} + sortedvector::Vector{T} + ord::OT + function SortedVector(occ::AbstractVector{T}, ord::OT=Base.Order.Forward) where {T, OT} + if issorted(occ, order=ord) + new{T, OT}(occ, ord) + else + new{T, OT}(sort(occ, order=ord), ord) + end + end +end +Base.:(==)(sv1::SortedVector, sv2::SortedVector) = sv1.sortedvector == sv2.sortedvector +Base.size(sv::SortedVector) = (length(sv.sortedvector),) +Base.@propagate_inbounds function Base.getindex(sv::SortedVector, i::Int) + @boundscheck !checkbounds(Bool, sv.sortedvector, i) && throw(BoundsError(sv, i)) + return sv.sortedvector[i] +end +Base.union(sv1::SortedVector{T}, svs::SortedVector{T}...) where {T} = + SortedVector(union(sv1.sortedvector, (occ.sortedvector for occ in svs)...)) + +# Special methods for fast operator construction +allocate_buffer(sv::AbstractVector) = similar(first(sv)) +function state_index(sv::SortedVector, state) + ret = searchsortedfirst(sv.sortedvector, state, order = sv.ord) + ret == length(sv) + 1 && return nothing + return sv.sortedvector[ret] == state ? ret : nothing +end +state_index(sv::AbstractVector, state) = findfirst(==(state), sv) + """ ManyBodyBasis(b, occupations) @@ -7,18 +36,18 @@ The basis has to know the associated one-body basis `b` and which occupation sta should be included. The occupations_hash is used to speed up checking if two many-body bases are equal. """ -struct ManyBodyBasis{S,B,H,UT} <: Basis - shape::S +struct ManyBodyBasis{B,O,UT} <: Basis + shape::Int onebodybasis::B - occupations::Vector{S} + occupations::O occupations_hash::UT - - function ManyBodyBasis{S,B,H}(onebodybasis::B, occupations::Vector{S}) where {S,B,H} - @assert isa(H, UInt) - new{S,B,H,typeof(H)}([length(occupations)], onebodybasis, occupations, hash(hash.(occupations))) + function ManyBodyBasis{B,O}(onebodybasis::B, occupations::O) where {B,O<:AbstractVector{Vector{Int}}} + h = hash(hash.(occupations)) + new{B,O,typeof(h)}(length(occupations), onebodybasis, occupations, h) end end -ManyBodyBasis(onebodybasis::B, occupations::Vector{S}) where {B,S} = ManyBodyBasis{S,B,hash(hash.(occupations))}(onebodybasis,occupations) +ManyBodyBasis(onebodybasis::B, occupations::O) where {B,O} = ManyBodyBasis{B,O}(onebodybasis, occupations) +ManyBodyBasis(onebodybasis::B, occupations::Vector{T}) where {B,T} = ManyBodyBasis(onebodybasis, SortedVector(occupations)) """ fermionstates(Nmodes, Nparticles) @@ -28,8 +57,8 @@ Generate all fermionic occupation states for N-particles in M-modes. `Nparticles` can be a vector to define a Hilbert space with variable particle number. """ -fermionstates(Nmodes::Int, Nparticles::Int) = _distribute_fermions(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]) -fermionstates(Nmodes::Int, Nparticles::Vector{Int}) = vcat([fermionstates(Nmodes, N) for N in Nparticles]...) +fermionstates(Nmodes::Int, Nparticles::Int) = SortedVector(_distribute_fermions(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]), Base.Reverse) +fermionstates(Nmodes::Int, Nparticles::Vector{Int}) = union((fermionstates(Nmodes, N) for N in Nparticles)...) fermionstates(onebodybasis::Basis, Nparticles) = fermionstates(length(onebodybasis), Nparticles) """ @@ -40,11 +69,11 @@ Generate all bosonic occupation states for N-particles in M-modes. `Nparticles` can be a vector to define a Hilbert space with variable particle number. """ -bosonstates(Nmodes::Int, Nparticles::Int) = _distribute_bosons(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]) -bosonstates(Nmodes::Int, Nparticles::Vector{Int}) = vcat([bosonstates(Nmodes, N) for N in Nparticles]...) +bosonstates(Nmodes::Int, Nparticles::Int) = SortedVector(_distribute_bosons(Nparticles, Nmodes, 1, zeros(Int, Nmodes), Vector{Int}[]), Base.Reverse) +bosonstates(Nmodes::Int, Nparticles::Vector{Int}) = union((bosonstates(Nmodes, N) for N in Nparticles)...) bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis), Nparticles) -==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash==b2.occupations_hash && b1.onebodybasis==b2.onebodybasis +==(b1::ManyBodyBasis, b2::ManyBodyBasis) = b1.occupations_hash == b2.occupations_hash && b1.onebodybasis == b2.onebodybasis """ basisstate([T=ComplexF64,] b::ManyBodyBasis, occupation::Vector) @@ -52,51 +81,36 @@ bosonstates(onebodybasis::Basis, Nparticles) = bosonstates(length(onebodybasis), Return a ket state where the system is in the state specified by the given occupation numbers. """ -function basisstate(::Type{T}, basis::ManyBodyBasis, occupation::Vector) where T - index = findfirst(isequal(occupation), basis.occupations) +function basisstate(::Type{T}, basis::ManyBodyBasis, occupation::Vector) where {T} + index = state_index(basis.occupations, occupation) if isa(index, Nothing) throw(ArgumentError("Occupation not included in many-body basis.")) end basisstate(T, basis, index) end -function isnonzero(occ1, occ2, index) - for i=1:length(occ1) - if i == index - if occ1[i] != occ2[i] + 1 - return false - end - else - if occ1[i] != occ2[i] - return false - end - end - end - true -end - """ create([T=ComplexF64,] b::ManyBodyBasis, index) Creation operator for the i-th mode of the many-body basis `b`. """ -function create(::Type{T}, b::ManyBodyBasis, index) where T +function create(::Type{T}, b::ManyBodyBasis, index) where {T} Is = Int[] Js = Int[] Vs = T[] # <{m}_i| at |{m}_j> - for i=1:length(b) - occ_i = b.occupations[i] + buffer = allocate_buffer(b.occupations) + for (i, occ_i) in enumerate(b.occupations) if occ_i[index] == 0 continue end - for j=1:length(b) - if isnonzero(occ_i, b.occupations[j], index) - push!(Is, i) - push!(Js, j) - push!(Vs, sqrt(occ_i[index])) - end - end + copyto!(buffer, occ_i) + buffer[index] -= 1 + j = state_index(b.occupations, buffer) + j === nothing && continue + push!(Is, i) + push!(Js, j) + push!(Vs, sqrt(occ_i[index])) end return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b))) end @@ -107,23 +121,23 @@ create(b::ManyBodyBasis, index) = create(ComplexF64, b, index) Annihilation operator for the i-th mode of the many-body basis `b`. """ -function destroy(::Type{T}, b::ManyBodyBasis, index) where T +function destroy(::Type{T}, b::ManyBodyBasis, index) where {T} Is = Int[] Js = Int[] Vs = T[] + buffer = allocate_buffer(b.occupations) # <{m}_j| a |{m}_i> - for i=1:length(b) - occ_i = b.occupations[i] + for (i, occ_i) in enumerate(b.occupations) if occ_i[index] == 0 continue end - for j=1:length(b) - if isnonzero(occ_i, b.occupations[j], index) - push!(Is, j) - push!(Js, i) - push!(Vs, sqrt(occ_i[index])) - end - end + copyto!(buffer, occ_i) + buffer[index] -= 1 + j = state_index(b.occupations, buffer) + j === nothing && continue + push!(Is, j) + push!(Js, i) + push!(Vs, sqrt(occ_i[index])) end return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b))) end @@ -134,8 +148,8 @@ destroy(b::ManyBodyBasis, index) = destroy(ComplexF64, b, index) Particle number operator for the i-th mode of the many-body basis `b`. """ -function number(::Type{T}, b::ManyBodyBasis, index) where T - diagonaloperator(b, T[b.occupations[i][index] for i in 1:length(b)]) +function number(::Type{T}, b::ManyBodyBasis, index) where {T} + diagonaloperator(b, T[occ[index] for occ in b.occupations]) end number(b::ManyBodyBasis, index) = number(ComplexF64, b, index) @@ -144,57 +158,33 @@ number(b::ManyBodyBasis, index) = number(ComplexF64, b, index) Total particle number operator. """ -function number(::Type{T}, b::ManyBodyBasis) where T - diagonaloperator(b, T[sum(b.occupations[i]) for i in 1:length(b)]) +function number(::Type{T}, b::ManyBodyBasis) where {T} + diagonaloperator(b, T[sum(occ) for occ in b.occupations]) end number(b::ManyBodyBasis) = number(ComplexF64, b) -function isnonzero(occ1, occ2, index1, index2) - for i=1:length(occ1) - if i == index1 && i == index2 - if occ1[i] != occ2[i] - return false - end - elseif i == index1 - if occ1[i] != occ2[i] + 1 - return false - end - elseif i == index2 - if occ1[i] != occ2[i] - 1 - return false - end - else - if occ1[i] != occ2[i] - return false - end - end - end - true -end - """ transition([T=ComplexF64,] b::ManyBodyBasis, to, from) Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes. + +Note that `to` and `from` can be collections of indices. The resulting operator in this case +will be equal to ``a^\\dagger_{to_1} a^\\dagger_{to_2} \\ldots a_{from_2} a_{from_1}``. """ -function transition(::Type{T}, b::ManyBodyBasis, to, from) where T +function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T} Is = Int[] Js = Int[] Vs = T[] + buffer = allocate_buffer(b.occupations) # <{m}_j| at_to a_from |{m}_i> - for i=1:length(b) - occ_i = b.occupations[i] - if occ_i[from] == 0 - continue - end - for j=1:length(b) - occ_j = b.occupations[j] - if isnonzero(occ_j, occ_i, to, from) - push!(Is, j) - push!(Js, i) - push!(Vs, sqrt(occ_i[from]) * sqrt(occ_j[to])) - end - end + for (i, occ_i) in enumerate(b.occupations) + C = state_transition!(buffer, occ_i, from, to) + C === nothing && continue + j = state_index(b.occupations, buffer) + j === nothing && continue + push!(Is, j) + push!(Js, i) + push!(Vs, C) end return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b))) end @@ -229,7 +219,7 @@ different modes of the N-particle basis. function manybodyoperator(basis::ManyBodyBasis, op) @assert op.basis_l == op.basis_r if op.basis_l == basis.onebodybasis - result = manybodyoperator_1(basis, op) + result = manybodyoperator_1(basis, op) elseif op.basis_l == basis.onebodybasis ⊗ basis.onebodybasis result = manybodyoperator_2(basis, op) else @@ -239,56 +229,59 @@ function manybodyoperator(basis::ManyBodyBasis, op) end function manybodyoperator_1(basis::ManyBodyBasis, op::Operator) - N = length(basis) S = length(basis.onebodybasis) result = DenseOperator(basis) - @inbounds for n=1:N, m=1:N - for j=1:S, i=1:S - C = coefficient(basis.occupations[m], basis.occupations[n], i, j) - if C != 0. - result.data[m,n] += C*op.data[i,j] - end + buffer = allocate_buffer(basis.occupations) + @inbounds for j = 1:S, i = 1:S + value = op.data[i, j] + iszero(value) && continue + for (m, occ) in enumerate(basis.occupations) + C = state_transition!(buffer, occ, i, j) + C === nothing && continue + n = state_index(basis.occupations, buffer) + n === nothing && continue + result.data[m, n] += C * value end end return result end -manybodyoperator_1(basis::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyoperator_1(basis,dagger(op))) +manybodyoperator_1(basis::ManyBodyBasis, op::AdjointOperator) = dagger(manybodyoperator_1(basis, dagger(op))) function manybodyoperator_1(basis::ManyBodyBasis, op::SparseOpPureType) N = length(basis) - M = op.data Is = Int[] Js = Int[] Vs = ComplexF64[] - @inbounds for colindex = 1:M.n - for i=M.colptr[colindex]:M.colptr[colindex+1]-1 - row = M.rowval[i] - value = M.nzval[i] - for m=1:N, n=1:N - C = coefficient(basis.occupations[m], basis.occupations[n], row, colindex) - if C != 0. - push!(Is, m) - push!(Js, n) - push!(Vs, C * value) - end - end + buffer = allocate_buffer(basis.occupations) + @inbounds for (row, column, value) in zip(findnz(op.data)...) + for (m, occ) in enumerate(basis.occupations) + C = state_transition!(buffer, occ, row, column) + C === nothing && continue + n = state_index(basis.occupations, buffer) + n === nothing && continue + push!(Is, m) + push!(Js, n) + push!(Vs, C * value) end end return SparseOperator(basis, sparse(Is, Js, Vs, N, N)) end function manybodyoperator_2(basis::ManyBodyBasis, op::Operator) - N = length(basis) S = length(basis.onebodybasis) @assert S^2 == length(op.basis_l) - @assert S^2 == length(op.basis_r) result = DenseOperator(basis) op_data = reshape(op.data, S, S, S, S) - occupations = basis.occupations - @inbounds for m=1:N, n=1:N - for l=1:S, k=1:S, j=1:S, i=1:S - C = coefficient(occupations[m], occupations[n], (i, j), (k, l)) - result.data[m,n] += C*op_data[i, j, k, l] + buffer = allocate_buffer(basis.occupations) + @inbounds for l = 1:S, k = 1:S, j = 1:S, i = 1:S + value = op_data[i, j, k, l] + iszero(value) && continue + for (m, occ) in enumerate(basis.occupations) + C = state_transition!(buffer, occ, (i, j), (k, l)) + C === nothing && continue + n = state_index(basis.occupations, buffer) + n === nothing && continue + result.data[m, n] += C * value end end return result @@ -300,21 +293,17 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOpType) Is = Int[] Js = Int[] Vs = ComplexF64[] - occupations = basis.occupations - rows = rowvals(op.data) - values = nonzeros(op.data) - @inbounds for column=1:S^2, j in nzrange(op.data, column) - row = rows[j] - value = values[j] - for m=1:N, n=1:N - # println("row:", row, " column:"column, ind_left) - index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2 + row]) - C = coefficient(occupations[m], occupations[n], index[1:2], index[3:4]) - if C!=0. - push!(Is, m) - push!(Js, n) - push!(Vs, C * value) - end + buffer = allocate_buffer(basis.occupations) + @inbounds for (row, column, value) in zip(findnz(op.data)...) + for (m, occ) in enumerate(basis.occupations) + index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2+row]) + C = state_transition!(buffer, occ, index[1:2], index[3:4]) + C === nothing && continue + n = state_index(basis.occupations, buffer) + n === nothing && continue + push!(Is, m) + push!(Js, n) + push!(Vs, C * value) end end return SparseOperator(basis, sparse(Is, Js, Vs, N, N)) @@ -327,154 +316,101 @@ end Expectation value of the one-body operator `op` in respect to the many-body `state`. """ -function onebodyexpect(op::AbstractOperator, state::Ket) - @assert isa(state.basis, ManyBodyBasis) +function onebodyexpect(op::AbstractOperator, state::Union{Ket,AbstractOperator}) + bas = basis(state) + @assert bas isa ManyBodyBasis @assert op.basis_l == op.basis_r - if state.basis.onebodybasis == op.basis_l - result = onebodyexpect_1(op, state) - # Not yet implemented: - # elseif state.basis.basis ⊗ state.basis.basis == op.basis_l - # result = onebodyexpect_2(op, state) + if bas.onebodybasis == op.basis_l + return onebodyexpect_1(op, state) + elseif bas.onebodybasis ⊗ bas.onebodybasis == op.basis_l + # Not yet implemented + throw(ArgumentError("`onebodyexpect` is not implemented for two-body states yet")) else throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis of the state.")) end - result -end - -function onebodyexpect(op::AbstractOperator, state::AbstractOperator) - @assert op.basis_l == op.basis_r - @assert state.basis_l == state.basis_r - @assert isa(state.basis_l, ManyBodyBasis) - if state.basis_l.onebodybasis == op.basis_l - result = onebodyexpect_1(op, state) - # Not yet implemented - # elseif state.basis.basis ⊗ state.basis.basis == op.basis_l - # result = onebodyexpect_2(op, state) - else - throw(ArgumentError("The basis of the given operator has to either be equal to b or b ⊗ b where b is the 1st quantization basis associated to the nparticle basis of the state.")) - end - result -end -onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state=states] - -function onebodyexpect_1(op::Operator, state::Ket) - N = length(state.basis) - S = length(state.basis.onebodybasis) - result = complex(0.) - occupations = state.basis.occupations - for m=1:N, n=1:N - value = conj(state.data[m])*state.data[n] - for i=1:S, j=1:S - C = coefficient(occupations[m], occupations[n], i, j) - if C != 0. - result += C*op.data[i,j]*value - end - end - end - result end -function onebodyexpect_1(op::Operator, state::Operator) - N = length(state.basis_l) - S = length(state.basis_l.onebodybasis) - result = complex(zero(promote_type(eltype(op),eltype(state)))) - occupations = state.basis_l.occupations - @inbounds for s=1:N, t=1:N - value = state.data[t,s] - for i=1:S, j=1:S - C = coefficient(occupations[s], occupations[t], i, j) - if !iszero(C) - result += C*op.data[i,j]*value - end +onebodyexpect(op::AbstractOperator, states::Vector) = [onebodyexpect(op, state) for state = states] + +get_value(state::Ket, m, n) = conj(state.data[m]) * state.data[n] +get_value(state::Operator, m, n) = state.data[n, m] +function onebodyexpect_1(op::Operator, state) + b = basis(state) + occupations = b.occupations + S = length(b.onebodybasis) + buffer = allocate_buffer(occupations) + result = complex(0.0) + for i = 1:S, j = 1:S + value = op.data[i, j] + iszero(value) && continue + for (m, occ) in enumerate(occupations) + C = state_transition!(buffer, occ, i, j) + C === nothing && continue + n = state_index(occupations, buffer) + n === nothing && continue + result += C * value * get_value(state, m, n) end end result end -function onebodyexpect_1(op::SparseOpPureType, state::Ket) - N = length(state.basis) - S = length(state.basis.onebodybasis) - result = complex(0.) - occupations = state.basis.occupations - M = op.data - @inbounds for colindex = 1:M.n - for i=M.colptr[colindex]:M.colptr[colindex+1]-1 - row = M.rowval[i] - value = M.nzval[i] - for m=1:N, n=1:N - C = coefficient(occupations[m], occupations[n], row, colindex) - if C != 0. - result += C*value*conj(state.data[m])*state.data[n] - end - end +function onebodyexpect_1(op::SparseOpPureType, state) + b = basis(state) + occupations = b.occupations + buffer = allocate_buffer(occupations) + result = complex(0.0) + @inbounds for (row, column, value) in zip(findnz(op.data)...) + for (m, occ) in enumerate(occupations) + C = state_transition!(buffer, occ, row, column) + C === nothing && continue + n = state_index(occupations, buffer) + n === nothing && continue + result += C * value * get_value(state, m, n) end end result end -function onebodyexpect_1(op::SparseOpPureType, state::Operator) - N = length(state.basis_l) - S = length(state.basis_l.onebodybasis) - result = complex(0.) - occupations = state.basis_l.occupations - M = op.data - @inbounds for colindex = 1:M.n - for i=M.colptr[colindex]:M.colptr[colindex+1]-1 - row = M.rowval[i] - value = M.nzval[i] - for s=1:N, t=1:N - C = coefficient(occupations[s], occupations[t], row, colindex) - if C != 0. - result += C*value*state.data[t,s] - end - end - end +Base.@propagate_inbounds function state_transition!(buffer, occ_m, at_indices, a_indices) + any(==(0), (occ_m[m] for m in at_indices)) && return nothing + result = 1 + copyto!(buffer, occ_m) + for i in at_indices + result *= buffer[i] + result == 0 && return nothing + buffer[i] -= 1 end - result -end - - -""" -Calculate the matrix element <{m}|at_1...at_n a_1...a_n|{n}>. -""" -Base.@propagate_inbounds function coefficient(occ_m, occ_n, at_indices, a_indices) - any(==(0), (occ_m[m] for m in at_indices)) && return 0. - any(==(0), (occ_n[n] for n in a_indices)) && return 0. - C = prod(√, (occ_m[m] for m in at_indices)) * prod(√, (occ_n[n] for n in a_indices)) - for i in 1:length(occ_m) - vm = occ_m[i] - vn = occ_n[i] - i in at_indices && (vm -= 1) - i in a_indices && (vn -= 1) - vm != vn && return zero(C) + for i in a_indices + buffer[i] += 1 + result *= buffer[i] end - return C + return √result end function _distribute_bosons(Nparticles, Nmodes, index, occupations, results) - if index==Nmodes + if index == Nmodes occupations[index] = Nparticles push!(results, copy(occupations)) else - for n=Nparticles:-1:0 + for n = Nparticles:-1:0 occupations[index] = n - _distribute_bosons(Nparticles-n, Nmodes, index+1, occupations, results) + _distribute_bosons(Nparticles - n, Nmodes, index + 1, occupations, results) end end return results end function _distribute_fermions(Nparticles, Nmodes, index, occupations, results) - if (Nmodes-index)+1 x^2, sv1) == [1, 4, 9, 16, 25] +@test map(x -> x^2, sv2) == [25, 16, 9, 4, 1] +@test findfirst(iseven, sv1) == 2 +@test findfirst(iseven, sv2) == 2 + # Test state creation Nmodes = 5 b = GenericBasis(Nmodes) @@ -105,6 +118,10 @@ t32 = transition(b_mb, 3, 2) @test 1e-12 > D(t32*basisstate(b_mb, [0, 2, 1, 0]), 2*basisstate(b_mb, [0, 1, 2, 0])) @test 1e-12 > D(number(b_mb, 2), transition(b_mb, 2, 2)) +t22_33 = transition(b_mb, (2, 2), (3, 3)) +d3 = destroy(b_mb, 3) +c2 = create(b_mb, 2) +@test t22_33 ≈ c2 * c2 * d3 * d3 # Single particle operator in second quantization b_single = GenericBasis(Nmodes)