diff --git a/src/LindbladVectorizedTensors.jl b/src/LindbladVectorizedTensors.jl index 5cd3d5e..18d691d 100644 --- a/src/LindbladVectorizedTensors.jl +++ b/src/LindbladVectorizedTensors.jl @@ -1,9 +1,5 @@ module LindbladVectorizedTensors -using Base.Filesystem -using CSV -using Combinatorics -using DataFrames using ITensors using IterTools using JSON @@ -12,39 +8,19 @@ using ProgressMeter import ITensors.state, ITensors.op, ITensors.space -export allequal, - chain, +export chain, chop, consecutivepairs, - construct_step_list, - disablegrifqtech, embed_slice, - filenamett, gellmannbasis, gellmannmatrix, - groupresults, - levels, - load_parameters, - oscdimensions, - parse_init_state_osc, - partialtrace, - vec, - vonneumannentropy, - chain_L1_state, - chain_basis_states, - level_subspace_proj, - parse_init_state, - parse_spin_state, - single_ex_state, + jwstring, sitenumber, - jwstring + vec, + vonneumannentropy include("utils.jl") -export spin_current_op_list - -include("current_operators.jl") - include("site_types/spinhalf.jl") include("site_types/vectorized_spinhalf.jl") include("site_types/fermion.jl") @@ -59,19 +35,7 @@ export dissipator_loss, dissipator_gain, dissipator, mixedlindbladplus, mixedlin include("site_types/vectorized_oscillator.jl") -export twositeoperators, - localop, - interactionop, - evolve, - fermioncurrent, - forwardflux, - backwardflux, - dissipator_symmetric, - dissipator_asymmetric, - lindbladian_xy, - hamiltonian_xy, - gkslcommutator, - gkslcommutator_itensor +export gkslcommutator, gkslcommutator_itensor include("operators.jl") diff --git a/src/closure.jl b/src/closure.jl index dcca4a4..78c8ea6 100644 --- a/src/closure.jl +++ b/src/closure.jl @@ -12,12 +12,12 @@ struct Closure innercoupling outercoupling damping - @doc""" - Closure(ω::Vector{<:Real}, γ::Vector{<:Real}, g::Vector{<:Real}, ζ::Vector{<:Complex}) + @doc """ + Closure(ω::Vector{<:Real}, γ::Vector{<:Real}, g::Vector{<:Real}, ζ::Vector{<:Complex}) - Closure is an aggregate type which stores the parameters needed to describe the set of - pseudomodes that make up a Markovian closure. - """ + Closure is an aggregate type which stores the parameters needed to describe the set of + pseudomodes that make up a Markovian closure. + """ function Closure( ω::Vector{<:Real}, γ::Vector{<:Real}, g::Vector{<:Real}, ζ::Vector{<:Complex} ) @@ -633,7 +633,9 @@ function filled_closure_op_adjoint( for (j, site) in enumerate(sitenumbers) # a ρ a† opstring = [repeat(["F⋅ * ⋅F"], site - 1); "A⋅ * ⋅A†"] - ℓ += (gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + ℓ += ( + gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))... + ) # -0.5 (a a† ρ + ρ a a†) = -0.5 (ρ - a† a ρ + ρ a† a) ℓ += 0.5damp(mc, j), "N⋅", site ℓ += 0.5damp(mc, j), "⋅N", site diff --git a/src/current_operators.jl b/src/current_operators.jl deleted file mode 100644 index c1de6e0..0000000 --- a/src/current_operators.jl +++ /dev/null @@ -1,78 +0,0 @@ -# Spin current operators -# ====================== - -function J⁺tag(::SiteType"S=1/2", left_site::Int, i::Int) - # TODO: still useful? - if i == left_site - str = "σx" - elseif i == left_site + 1 - str = "σy" - else - str = "Id" - end - return str -end -function J⁺tag(::SiteType"vecS=1/2", left_site::Int, i::Int) - if i == left_site - str = "vecσx" - elseif i == left_site + 1 - str = "vecσy" - else - str = "vecId" - end - return str -end -function J⁺tag(::SiteType"HvS=1/2", left_site::Int, i::Int) - return J⁺tag(SiteType("vecS=1/2"), left_site, i) -end - -function J⁻tag(::SiteType"S=1/2", left_site::Int, i::Int) - # Just as `J⁺tag`, but for σʸ⊗σˣ - if i == left_site - str = "σy" - elseif i == left_site + 1 - str = "σx" - else - str = "Id" - end - return str -end -function J⁻tag(::SiteType"vecS=1/2", left_site::Int, i::Int) - if i == left_site - str = "vecσy" - elseif i == left_site + 1 - str = "vecσx" - else - str = "vecId" - end - return str -end -function J⁻tag(::SiteType"HvS=1/2", left_site::Int, i::Int) - return J⁻tag(SiteType("vecS=1/2"), left_site, i) -end - -function spin_current_op_list(sites::Vector{Index{Int64}}) - N = length(sites) - # Check if all sites are spin-½ sites. - if all(x -> SiteType("S=1/2") in x, sitetypes.(sites)) - st = SiteType("S=1/2") - MPtype = MPO - elseif all(x -> SiteType("vecS=1/2") in x, sitetypes.(sites)) - st = SiteType("vecS=1/2") - MPtype = MPS - elseif all(x -> SiteType("HvS=1/2") in x, sitetypes.(sites)) - st = SiteType("HvS=1/2") - MPtype = MPS - else - throw( - ArgumentError( - "spin_current_op_list works with SiteTypes " * - "\"S=1/2\", \"vecS=1/2\" or \"HvS=1/2\".", - ), - ) - end - # - J⁺ = [MPtype(sites, [J⁺tag(st, k, i) for i in 1:N]) for k in 1:(N - 1)] - J⁻ = [MPtype(sites, [J⁻tag(st, k, i) for i in 1:N]) for k in 1:(N - 1)] - return -0.5 .* (J⁺ .- J⁻) -end diff --git a/src/operators.jl b/src/operators.jl index ec77769..06571bb 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -1,431 +1,3 @@ -@doc raw""" - twositeoperators(sites::Vector{Index{Int64}}, localcfs::Vector{<:Real}, interactioncfs::Vector{<:Real}) - -Return a list of ITensor operators representing the two-site-gates -decomposition of the Hamiltonian or Lindbladian of the system. - -The element `list[j]` is the operator ``hⱼ,ⱼ₊₁`` (or ``ℓⱼ,ⱼ₊₁``, depending on -the notation) in the expansion -```math -H=\sum_{j=1}^{N}h_{j,j+1}. -``` - -# Arguments -- `sites::Vector{Index}`: the sites making up the system -- `localcfs::Vector{Index}`: coefficients of one-site operators -- `interactioncfs::Vector{Index}`: coefficients of two-site operators - -Index convention: the `j`-th element refers to the ``hⱼ,ⱼ₊₁``/``ℓⱼ,ⱼ₊₁`` term. - -# Usage -This function is extremely specialised to our use-case: the function -automatically chooses some operators representing local and interaction terms, -but their OpNames are currently hardcoded (see `localop` and `interactionop`). -""" -function twositeoperators( - sites::Vector{Index{Int64}}, localcfs::Vector{<:Real}, interactioncfs::Vector{<:Real} -) - list = ITensor[] - # The sites at the edge of the chain need to be treated separately since their - # local coefficients must not be multiplied by 1/2. Instead of writing these - # particular cases out of the loop, I just multiply them beforehand by 2. - localcfs[begin] *= 2 - localcfs[end] *= 2 - for j in 1:(length(sites) - 1) - s1 = sites[j] - s2 = sites[j + 1] - h = - 0.5localcfs[j] * localop(s1) * op("Id", s2) + - 0.5localcfs[j + 1] * op("Id", s1) * localop(s2) + - interactioncfs[j] * interactionop(s1, s2) - push!(list, h) - end - return list -end - -""" - localop(s::Index) - -Return the local ITensor operator associated to `s`'s SiteType. - -# Usage -Currently this is just a way to simplify code, the user doesn't really have a -say in which operator is used, it is just hardcoded. -""" -function localop(s::Index) - if SiteType("S=1/2") ∈ sitetypes(s) - t = 0.5op("σz", s) - elseif SiteType("Osc") ∈ sitetypes(s) - t = op("N", s) - elseif SiteType("vecS=1/2") ∈ sitetypes(s) - t = 0.5im * (op("σz:Id", s) - op("Id:σz", s)) - elseif SiteType("HvS=1/2") ∈ sitetypes(s) - t = -0.5im * (op("σz⋅", s) - op("⋅σz", s)) - elseif SiteType("vecOsc") ∈ sitetypes(s) - t = im * (op("N:Id", s) - op("Id:N", s)) - elseif SiteType("HvOsc") ∈ sitetypes(s) - t = -im * (op("N⋅", s) - op("⋅N", s)) - else - throw(DomainError(s, "SiteType non riconosciuto.")) - end - return t -end - -""" - interactionop(s1::Index, s2::Index) - -Return the ITensor local operator associated to the SiteTypes of `s1` and `s2`. - -# Usage -Currently this is just a way to simplify code, the user doesn't really have a -say in which operator is used, it is just hardcoded. -""" -function interactionop(s1::Index, s2::Index) - if SiteType("S=1/2") ∈ sitetypes(s1) && SiteType("S=1/2") ∈ sitetypes(s2) - t = -0.5 * (op("σ-", s1) * op("σ+", s2) + op("σ+", s1) * op("σ-", s2)) - elseif SiteType("Osc") ∈ sitetypes(s1) && SiteType("S=1/2") ∈ sitetypes(s2) - t = op("X", s1) * op("σx", s2) - elseif SiteType("S=1/2") ∈ sitetypes(s1) && SiteType("Osc") ∈ sitetypes(s2) - t = op("σx", s1) * op("X", s2) - elseif SiteType("Osc") ∈ sitetypes(s1) && SiteType("Osc") ∈ sitetypes(s2) - t = (op("a+", s1) * op("a-", s2) + op("a-", s1) * op("a+", s2)) - # - elseif SiteType("vecS=1/2") ∈ sitetypes(s1) && SiteType("vecS=1/2") ∈ sitetypes(s2) - t = - -0.5im * ( - op("σ-:Id", s1) * op("σ+:Id", s2) + op("σ+:Id", s1) * op("σ-:Id", s2) - - op("Id:σ-", s1) * op("Id:σ+", s2) - op("Id:σ+", s1) * op("Id:σ-", s2) - ) - elseif SiteType("vecOsc") ∈ sitetypes(s1) && SiteType("vecS=1/2") ∈ sitetypes(s2) - t = im * (op("asum:Id", s1) * op("σx:Id", s2) - op("Id:asum", s1) * op("Id:σx", s2)) - elseif SiteType("vecS=1/2") ∈ sitetypes(s1) && SiteType("vecOsc") ∈ sitetypes(s2) - t = im * (op("σx:Id", s1) * op("asum:Id", s2) - op("Id:σx", s1) * op("Id:asum", s2)) - # - elseif SiteType("HvS=1/2") ∈ sitetypes(s1) && SiteType("HvS=1/2") ∈ sitetypes(s2) - t = - 0.5im * ( - op("σ-⋅", s1) * op("σ+⋅", s2) + op("σ+⋅", s1) * op("σ-⋅", s2) - - op("⋅σ-", s1) * op("⋅σ+", s2) - op("⋅σ+", s1) * op("⋅σ-", s2) - ) - elseif SiteType("HvOsc") ∈ sitetypes(s1) && SiteType("HvS=1/2") ∈ sitetypes(s2) - t = -im * (op("asum⋅", s1) * op("σx⋅", s2) - op("⋅asum", s1) * op("⋅σx", s2)) - elseif SiteType("HvS=1/2") ∈ sitetypes(s1) && SiteType("HvOsc") ∈ sitetypes(s2) - t = -im * (op("σx⋅", s1) * op("asum⋅", s2) - op("⋅σx", s1) * op("⋅asum", s2)) - elseif SiteType("HvOsc") ∈ sitetypes(s1) && SiteType("HvOsc") ∈ sitetypes(s2) - # Oscillator-oscillator interaction H = a₁†a₂ + a₂†a₁ that goes in the - # commutator -i[H,ρ] as - # -i (a₁†a₂ ρ + a₂†a₁ ρ - ρ a₁†a₂ - ρ a₂†a₁) - t = - -im * ( - op("a+⋅", s1) * op("a-⋅", s2) + op("a+⋅", s2) * op("a-⋅", s1) - - op("⋅a+", s1) * op("⋅a-", s2) - op("⋅a+", s2) * op("⋅a-", s1) - ) - else - throw(DomainError((s1, s2), "Unrecognised SiteTypes.")) - end - return t -end - -""" - evolve(initialstate, timesteplist, nskip, STorder, linksodd, linkseven, maxerr, maxrank; fout) - -Evolves the initial state in time using the TEBD algorithm. -Return a list of the values of the functions passed within `fout` at each time -instant, together with a list of such instants. - -# Arguments -- `initialstate` → MPS of the initial state -- `timesteplist` → list of the time steps at which the state is evolved -- `nskip` → number of time steps to be skipped between each evaluation - of the output functions -- `STorder` → order of Suzuki-Trotter expansion -- `linksodd` → two-site gates on odd sites -- `linkseven` → two-site gates on even sites -- `maxerr` → maximum truncation error of the MPS during evolution -- `maxrank` → maximum allowed rank of the MPS during evolution -- `fout` → list of functions which take an MPS as their only argument - -The state is evolved through each of the instants specified in `timesteplist`, -but functions in `fout` are evaluated only after `nskip` steps each time. This -could save some computation time if a fine-grained output is not needed. - -# Example -If `fout = [a1, a2]`, the function returns the tuple `(ts, b1, b2)` where -`b1` and `b2` are the result of applying `a1` and `a2` on the state for each -time instant in `ts`, that is `b1[j]` is equal to `a1` applied to the state at -`ts[j]`. -""" -function evolve( - initialstate, timesteplist, nskip, STorder, linksodd, linkseven, maxerr, maxrank; fout -) - @info "Computing Suzuki-Trotter decomposition of the evolution operator." - τ = timesteplist[2] - timesteplist[1] - if STorder == 1 - # Nothing to do if the order is 1. - list = [linksodd(τ), linkseven(τ)] - seq = repeat([1, 2], nskip) - elseif STorder == 2 - # At 2nd order we can group the operators with t/2 (from two consecutive - # steps) in order to save some time and numerical errors. - # Instead of a series of 3*nskip operators between a measurement and the - # next one, we have only 2*nskip+1 of them. - list = [linksodd(0.5τ), linkseven(τ), linksodd(τ)] - seq = [1; repeat([2, 3], nskip - 1); 2; 1] - elseif STorder == 4 - c = (4 - 4^(1 / 3))^(-1) - list = [ - linksodd(c * τ), - linkseven(2c * τ), - linksodd((0.5 - c) * τ), - linkseven((1 - 4c) * τ), - ] - seq = repeat([1, 2, 3, 4, 3, 2, 1], nskip) - # There could be still some margin for improvement: maybe group the first - # and last `linksodd` when nskip > 1. It doesn't make that much of a - # difference, though. - else - throw( - DomainError( - STorder, - "Suzuki-Trotter expansion of order $STorder " * - "is not supported. Only orders 1, 2 and 4 are " * - "available.", - ), - ) - end - - # Applying `list` to a state means evolving it for a time equal to nskip*τ, - # where τ is the integration time step. - @info "Starting time evolution." - state = initialstate - returnvalues = [[f(state) for f in fout]] - - tout = timesteplist[1:nskip:end] - progress = Progress(length(tout), 1, "Simulation in progress", 30) - for _ in timesteplist[(1 + nskip):nskip:end] - for n in seq - state .= apply(list[n], state; cutoff=maxerr, maxdim=maxrank) - end - push!(returnvalues, [f(state) for f in fout]) - next!(progress) - end - - # The list `returnvalues` contains N ≡ length(timesteplist[1:nskip:end]) - # lists as elements; in each one, the j-th element is the result of fout[j] - # when applied to the state at a given time instant. - # We want to reorganise the output, so that the function returns a list `tout` - # with the instant at which the observables are measures, together with as - # many lists as the functions in `fout`: the j-th list will be the result - # of fout[j] applied to the state for each t in `tout`. - outresults = [[returnvalues[k][j] for k in eachindex(tout)] for j in eachindex(fout)] - return tout, outresults... -end - -# Current operators -# ----------------- - -""" - current(sites, leftsite::Int, rightsite::Int) - -Return the current operator from site `leftsite` to `rightside` in `sites`. - -# Usage -The sites (strictly) between the given sites are assigned the ITensor operator -with OpName "Z"; `sites(leftsite)` and `sites(rightsite)` are assigned "plus" -and "minus", and other sites are assigned "Id" operators. All these operators -thus must be correctly defined in order to use this function. -""" -function fermioncurrent(sites, leftsite::Int, rightsite::Int) - # The function doesn't care about the SiteType on which it is applied. - # This allows us to use it also on oscillators, when needed, for example - # to measure the current between the first oscillator sites just beyond - # the edge of the spin chain (it works!). - # The oscillators need to be exactly at the `leftsite` and `rightsite` - # positions, and not between them, otherwise the function doesn't work - # since there is no "Z" operator defined for oscillators. - n = rightsite - leftsite - tags1 = repeat(["Id"], length(sites)) - tags2 = repeat(["Id"], length(sites)) - tags1[leftsite] = "plus" - tags2[leftsite] = "minus" - for l in (leftsite + 1):(rightsite - 1) - tags1[l] = "Z" - tags2[l] = "Z" - end - tags1[rightsite] = "minus" - tags2[rightsite] = "plus" - - if ( - SiteType("S=1/2") ∈ sitetypes(first(sites)) || - SiteType("Osc") ∈ sitetypes(first(sites)) - ) - op = im * (-1)^n * (MPO(sites, tags1) - MPO(sites, tags2)) - elseif ( - SiteType("vecS=1/2") ∈ sitetypes(first(sites)) || - SiteType("HvS=1/2") ∈ sitetypes(first(sites)) || - SiteType("vecOsc") ∈ sitetypes(first(sites)) || - SiteType("HvOsc") ∈ sitetypes(first(sites)) - ) - op = - im * - (-1)^n * - ( - MPS(ComplexF64, sites, string.("vec", tags1)) - - MPS(ComplexF64, sites, string.("vec", tags2)) - ) - # The MPS(sites, states) constructor assumes Float64 as the data type, - # but here we need it to be complex. - else - throw(DomainError(s, "Unrecognised SiteType.")) - end - return op -end - -function forwardflux(sites, leftsite::Int, rightsite::Int) - # Just one of the two terms in j_{k,k'}... - n = rightsite - leftsite - tags1 = repeat(["Id"], length(sites)) - tags1[leftsite] = "plus" - for l in (leftsite + 1):(rightsite - 1) - tags1[l] = "Z" - end - tags1[rightsite] = "minus" - - if ( - SiteType("S=1/2") ∈ sitetypes(first(sites)) || - SiteType("Osc") ∈ sitetypes(first(sites)) - ) - op = im * (-1)^n * (MPO(sites, tags1) - MPO(sites, tags2)) - elseif ( - SiteType("vecS=1/2") ∈ sitetypes(first(sites)) || - SiteType("HvS=1/2") ∈ sitetypes(first(sites)) || - SiteType("vecOsc") ∈ sitetypes(first(sites)) || - SiteType("HvOsc") ∈ sitetypes(first(sites)) - ) - op = im * (-1)^n * MPS(ComplexF64, sites, string.("vec", tags1)) - else - throw(DomainError(s, "Unrecognised SiteType.")) - end - return op -end -function backwardflux(sites, leftsite::Int, rightsite::Int) - # ...and the other one. - n = rightsite - leftsite - tags2 = repeat(["Id"], length(sites)) - tags2[leftsite] = "minus" - for l in (leftsite + 1):(rightsite - 1) - tags2[l] = "Z" - end - tags2[rightsite] = "plus" - - if ( - SiteType("S=1/2") ∈ sitetypes(first(sites)) || - SiteType("Osc") ∈ sitetypes(first(sites)) - ) - op = im * (-1)^n * (MPO(sites, tags1) - MPO(sites, tags2)) - elseif ( - SiteType("vecS=1/2") ∈ sitetypes(first(sites)) || - SiteType("HvS=1/2") ∈ sitetypes(first(sites)) || - SiteType("vecOsc") ∈ sitetypes(first(sites)) || - SiteType("HvOsc") ∈ sitetypes(first(sites)) - ) - op = im * (-1)^n * MPS(ComplexF64, sites, string.("vec", tags2)) - else - throw(DomainError(s, "Unrecognised SiteType.")) - end - return op -end - -""" - hamiltonian_xy(length::Integer, energy::Real, coupling::Real) - -Return the Hamiltonian of an XY spin chain of given `length` and energies. -""" -function hamiltonian_xy(N, energy, coupling) - op = OpSum() - for j in 1:(N - 1) - op += -0.5coupling, "S+", j, "S-", j + 1 - op += -0.5coupling, "S-", j, "S+", j + 1 - end - for j in 1:N - op += energy, "Sz", j - end - return op -end - -""" - lindbladian_xy(length::Integer, spin_energy::Real, spin_coupling::Real) - -Return the Lindbladian operator ``-i[H,—]``, in vectorised form, associated to the -Hamiltonian ``H`` of an XY spin chain of given `length` and energies. -""" -function lindbladian_xy(chain_length::Integer, spin_energy::Real, spin_coupling::Real) - N = chain_length - ε = spin_energy - λ = spin_coupling - - op = OpSum() - for j in 1:(N - 1) - op += +0.5im * λ, "σ-⋅", j, "σ+⋅", j + 1 - op += +0.5im * λ, "σ+⋅", j, "σ-⋅", j + 1 - op += -0.5im * λ, "⋅σ-", j, "⋅σ+", j + 1 - op += -0.5im * λ, "⋅σ+", j, "⋅σ-", j + 1 - end - for j in 1:N - op += -0.5im * ε, "σz⋅", j - op += +0.5im * ε, "⋅σz", j - end - return op -end - -avgn(ε, T) = T == 0 ? 0 : (ℯ^(ε / T) - 1)^(-1) - -""" - dissipator_symmetric(n, excitation_energy, spin_damping_coefficient, temperature) - -Return the dissipator ``D(ρ) ∝ σˣ ρ (σˣ)† - ½ {(σˣ)† σˣ, ρ}`` on site `n`. -""" -function dissipator_symmetric( - n::Integer, excitation_energy::Real, spin_damping_coefficient::Real, temperature::Real -) - # TODO: move the ξ constant out of the function. - ε = excitation_energy - κ = spin_damping_coefficient - T = temperature - ξ = κ * (1 + 2avgn(ε, T)) - - op = OpSum() - op += +ξ, "σx⋅ * ⋅σx", n - op += -ξ, "Id", n - return op -end - -""" - dissipator_asymmetric(n, frequency, damping_coefficient, temperature) - -Return the dissipator ``D(ρ) = γ⁻(σ⁻ρσ⁺ - ½ {σ⁺σ⁻, ρ}) + γ⁺(σ⁺ρσ⁻ - ½ {σ⁻σ⁺, ρ})`` on -site `n`. -""" -function dissipator_asymmetric( - n::Integer, excitation_energy::Real, spin_damping_coefficient::Real, temperature::Real -) - ω = excitation_energy - γ = spin_damping_coefficient - T = temperature - - γ⁺ = γ * avgn(ω, T) - γ⁻ = γ * (1 + avgn(ω, T)) - - op = OpSum() - op += γ⁻, "σ-⋅ * ⋅σ+", n - op += -0.5γ⁻, "σ+⋅ * σ-⋅", n - op += -0.5γ⁻, "⋅σ- * ⋅σ+", n - op += γ⁺, "σ+⋅ * ⋅σ-", n - op += -0.5γ⁺, "σ-⋅ * σ+⋅", n - op += -0.5γ⁺, "⋅σ+ * ⋅σ-", n - return op -end - # Easy construction of the unitary part of the GKSL equation # ---------------------------------------------------------- function makeopsumpairs(args...) diff --git a/src/site_types/fermion.jl b/src/site_types/fermion.jl index 3a457fa..f3b57f9 100644 --- a/src/site_types/fermion.jl +++ b/src/site_types/fermion.jl @@ -3,14 +3,14 @@ function ITensors.op(::OpName"Id", ::SiteType"Fermion") end function ITensors.op!(Op::ITensor, ::OpName"A", ::SiteType"Fermion", s::Index) - return Op[s' => 1, s => 2] = 1.0 + return Op[s' => 1, s => 2] = 1.0 end function ITensors.op!(Op::ITensor, on::OpName"a", st::SiteType"Fermion", s::Index) return ITensors.op!(Op, OpName("A"), st, s) end function ITensors.op!(Op::ITensor, ::OpName"Adag", ::SiteType"Fermion", s::Index) - return Op[s' => 2, s => 1] = 1.0 + return Op[s' => 2, s => 1] = 1.0 end function ITensors.op!(Op::ITensor, on::OpName"adag", st::SiteType"Fermion", s::Index) return ITensors.op!(Op, OpName("Adag"), st, s) diff --git a/src/utils.jl b/src/utils.jl index 8e8eab0..d2a0b58 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,9 +1,9 @@ """ statenamestring(sn::StateName{T}) where T -Return the name of the state "sn" as a string. +Return the name of the state `sn` as a string. """ -function statenamestring(sn::StateName{T}) where T +function statenamestring(sn::StateName{T}) where {T} return string(T) end @@ -29,32 +29,6 @@ function jwstring(; start, stop, op::AbstractString="F") return collect(Iterators.flatten([(op, start + k) for k in 1:(stop - start - 1)])) end -""" - groupresults(obs::AbstractObserver, name::String) - -If the observer returns an array of values, rearrange the results in a big matrix. -""" -function groupresults(obs::AbstractObserver, name::String) - return mapreduce(permutedims, vcat, obs[!, name]) -end - -""" - disablegrifqtech() - -Disable the graphical output if the script is running on Qtech. -""" -function disablegrifqtech() - # The script crashes if it is executed on Qtech, unless we disable - # the graphical output. - if gethostname() == "qtech.fisica.unimi.it" || gethostname() == "qtech2.fisica.unimi.it" - ENV["GKSwstype"] = "100" - @info "Running on remote server. Disabling graphical output." - else - delete!(ENV, "GKSwstype") - # If the key "GKSwstype" doesn't exist then nothing happens. - end -end - # Vectorisation utilities # ======================= # In order to use tr(x'*y) as a tool to extract coefficient the basis must @@ -69,6 +43,7 @@ Compute the vector of coefficients of the matrix `A` wrt the basis `basis`. function vec(A::Matrix, basis::Vector) return [tr(b' * A) for b in basis] end + """ vec(L::Function, basis::Vector) @@ -188,256 +163,6 @@ function canonicalbasis(dim) return [canonicalmatrix(i, j, dim) for (i, j) in [Base.product(1:dim, 1:dim)...]] end -# Spin-half space utilities -# ========================= - -# Chain eigenstate basis -# ---------------------- -# How to measure how much each eigenspace of the number operator (of the -# whole spin chain) "contributes" to a given state ρ? -# We build a projector operator associated to each eigenspace. -# The projector on the m-th level eigenspace can be made by simply adding the -# orthogonal projections on each element of the canonical basis which has m -# "Up" spins and N-m "Down" spins, N being the length of the chain. -# This may result in a very big MPO. If these are needed for more than one -# simulation, calculating them once and for all before the simulations start -# may help to cut down the computation time. -""" - chain_basis_states(n::Int, level::Int) - -Return a list of strings which can be used to build the MPS of all states in -the ``ℂ²ⁿ`` basis that contain `level` "Up" spins. -""" -function chain_basis_states(n::Int, level::Int) - return unique(permutations([ - repeat(["Up"], level) - repeat(["Dn"], n - level) - ])) -end - -""" - level_subspace_proj(sites::Vector{Index{Int64}}, l::Int) - -Return the projector on the subspace with `l` "Up" spins. -""" -function level_subspace_proj(sites::Vector{Index{Int64}}, l::Int) - N = length(sites) - # Check if all sites are spin-½ sites. - if all(x -> SiteType("S=1/2") in x, sitetypes.(sites)) - projs = [ - projector(MPS(sites, names); normalize=false) for - names in chain_basis_states(N, l) - ] - elseif all(x -> SiteType("vecS=1/2") in x, sitetypes.(sites)) || - all(x -> SiteType("HvS=1/2") in x, sitetypes.(sites)) - projs = [MPS(sites, names) for names in chain_basis_states(N, l)] - else - throw( - ArgumentError( - "level_subspace_proj works with SiteTypes " * - "\"S=1/2\", \"vecS=1/2\" or \"HvS=1/2\".", - ), - ) - end - # Somehow return sum(projs) doesn't work… we have to sum manually. - P = projs[1] - for p in projs[2:end] - P = +(P, p; cutoff=1e-10) - end - return P -end - -# First-level chain eigenstates -# ----------------------------- -# It is useful to have at hand the eigenstates of the free chain Hamiltonian -# within the first-level eigenspace of the number operator ([H,N] = 0). -# States |sₖ⟩ with an up spin on site k and a down spin on the others are -# in fact not eigenstates, but they still form a basis for the eigenspace; -# wrt the {sₖ}ₖ (k ∈ {1,…,N}) basis the free chain Hamiltonian is written as -# ⎛ ε λ 0 0 … 0 ⎞ -# ⎜ λ ε λ 0 … 0 ⎟ -# ⎜ 0 λ ε λ … 0 ⎟ -# ⎜ 0 0 λ ε … 0 ⎟ -# ⎜ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⎟ -# ⎝ 0 0 0 0 … ε ⎠ -# whose eigenstates are |vⱼ⟩= ∑ₖ sin(kjπ /(N+1)) |sₖ⟩, con j ∈ {1,…,N}. -# Note that they are not normalised: ‖|vⱼ⟩‖² = (N+1)/2. -""" - single_ex_state(sites::Vector{Index{Int64}}, k::Int) - -Return the MPS of a state with a single excitation on the `k`-th site. -""" -function single_ex_state(sites::Vector{Index{Int64}}, k::Int) - N = length(sites) - if k ∈ 1:N - states = [i == k ? "Up" : "Dn" for i in 1:N] - else - throw( - DomainError( - k, - "Trying to build a state with an excitation localised " * - "at site $k, which does not belong to the chain: please " * - "insert a value between 1 and $N.", - ), - ) - end - return MPS(sites, states) -end - -""" - chain_L1_state(sites::Vector{Index{Int64}}, j::Int) - -Return the `j`-th eigenstate of the chain Hamiltonian in the single-excitation -subspace. -""" -function chain_L1_state(sites::Vector{Index{Int64}}, j::Int) - # FIXME: this seems just wrong. Why not computing directly vec(|vⱼ⟩ ⊗ ⟨vⱼ|), - # without going through the linear combination? - # - # Careful with the coefficients: this isn't |vⱼ⟩ but vec(|vⱼ⟩ ⊗ ⟨vⱼ|), - # so if we want to build it as a linear combination of the |sₖ⟩'s we - # need to square the coefficients. - # Note that vectorised projectors satisfy - # ⟨vec(a⊗aᵀ), vec(b⊗bᵀ)⟩ = tr((a⊗aᵀ)ᵀ b⊗bᵀ) = (aᵀb)². - # We don't expect the norm of this MPS to be 1. What has to be 1, is the - # sum of the inner products of vec(|sₖ⟩ ⊗ ⟨sₖ|) and this vector, for - # k ∈ {1,…,N}. We get ⟨Pₙ|vⱼ⟩ = 1 only if n=j, otherwise it is 0. - N = length(sites) - if !(j in 1:N) - throw( - DomainError( - j, - "Trying to build a chain eigenstate with invalid index " * - "$j: please insert a value between 1 and $N.", - ), - ) - end - states = [ - 2 / (N + 1) * sin(j * k * π / (N + 1))^2 * single_ex_state(sites, k) for k in 1:N - ] - return sum(states) -end - -# Choice of the spin chain's initial state -# ---------------------------------------- -""" - parse_init_state(sites::Vector{Index{Int64}}, state::String) - -Return an MPS representing a particular state of the spin chain, given -by the string `state`: - - * "empty" -> empty state (aka the ground state of the chain Hamiltonian) - * "1locM" -> state with a single excitation at site `M` (``M ∈ {1,…,N}``) - * "1eigM" -> single-excitation eigenstate of the chain Hamiltonian with `M` nodes (``M ∈ {0,…,N-1}``) - -The string is case-insensitive. The length `N` of the chain is computed -from `sites`. -""" -function parse_init_state(sites::Vector{Index{Int64}}, state::String) - state = lowercase(state) - if state == "empty" - v = MPS(sites, "Dn") - elseif occursin(r"^1loc", state) - j = parse(Int, replace(state, "1loc" => "")) - v = single_ex_state(sites, j) - elseif occursin(r"^1eig", state) - j = parse(Int, replace(state, "1eig" => "")) - # The j-th eigenstate has j-1 nodes - v = chain_L1_state(sites, j + 1) - else - throw( - DomainError( - state, - "Unrecognised state: please choose from \"empty\", " * - "\"1locN\" or \"1eigN\".", - ), - ) - end - return v -end - -""" - parse_spin_state(site::Index{Int64}, state::String) - -Return the MPS of a single spin site representing a particular state, given -by the string `state`: - -- "empty", "dn", "down" → spin-down state -- "up" → spin-up state -- "x+" → ``1/√2 ( |+⟩ + |-⟩ )`` state -""" -function parse_spin_state(site::Index{Int64}, state::String) - state = lowercase(state) - if state == "empty" || state == "dn" || state == "down" - v = ITensors.state(site, "Dn") - elseif state == "up" - v = ITensors.state(site, "Up") - elseif state == "x+" - v = 1 / sqrt(2) * (ITensors.state(site, "Up") + ITensors.state(site, "Dn")) - else - throw(DomainError( - state, - "Unrecognised state: please choose from \"empty\", - \"up\", \"down\" or \"x+\".", - )) - end - return MPS([v]) -end - -# Oscillator space utilities -# ========================== - -""" - oscdimensions(N, basedim, decay) - -Compute a decreasing sequence of length `N` where the first two elements are -equal to `basedim` and the following ones are given by -`floor(2 + basedim * ℯ^(-decay * n))`. - -Useful to determine the dimensions of oscillator sites in a TEDOPA chain. -""" -function oscdimensions(length, basedim, decay) - f(j) = 2 + basedim * ℯ^(-decay * j) - return [basedim; basedim; (Int ∘ floor ∘ f).(3:length)] -end - -""" - parse_init_state_osc(site::Index{Int64}, - statename::String; - ) - -Return an MPS representing a particular state of a harmonic oscillator, given -by the string `statename`: - -- "thermal" → thermal equilibrium state -- "fockN" → `N`-th eigenstate of the number operator (element of Fock basis) -- "empty" → alias for "fock0" - -The string is case-insensitive. Other parameters required to build the state -(e.g. frequency, temperature) may be supplied as keyword arguments. -""" -function parse_init_state_osc(site::Index{Int64}, statename::String; kwargs...) - # TODO: maybe remove "init" from title? It is a generic state, after all. - statename = lowercase(statename) - if statename == "thermal" - s = state(site, "ThermEq"; kwargs...) - elseif occursin(r"^fock", statename) - j = parse(Int, replace(statename, "fock" => "")) - s = state(site, "$j") - elseif statename == "empty" - s = state(site, "0") - else - throw( - DomainError( - statename, - "Unrecognised state name; please choose from " * - "\"empty\", \"fockN\" or \"thermal\".", - ), - ) - end - return MPS([s]) -end - # (Von Neumann) Entropy # ===================== """ @@ -482,140 +207,8 @@ function chop(x::Complex; tolerance=1e-10) return Complex(chop(real(x)), chop(imag(x))) end -# Manipulating ITensor objects -# ============================ - -""" - sitetypes(s::Index) - -Return the ITensor tags of `s` as SiteTypes. - -This function is already defined in the ITensor library, but it is not publicly -accessible. -""" -function sitetypes(s::Index) - ts = tags(s) - return SiteType[SiteType(ts.data[n]) for n in 1:length(ts)] -end - -# Reading the parameters of the simulations -# ========================================= - -""" - isjson(filename::String) - -Check if `filename` ends in ".json". - -By design, filenames consisting of only ".json" return `false`. -""" -function isjson(filename::String) - return length(filename) > 5 && filename[(end - 4):end] == ".json" - # We ignore a file whose name is only ".json". -end - -""" - load_parameters(file_list) - -Load the JSON files contained in `file_list` into dictionaries, returning a -list of dictionaries, one for each file. - -If `file_list` is a filename, then the list will contain just one dictionary; -if `file_list` is a directory, every JSON file within it is loaded and a -dictionary is created for each one of them. -""" -function load_parameters(file_list) - if isempty(file_list) - throw(ErrorException("No parameter file provided.")) - end - first_arg = file_list[1] - prev_dir = pwd() - if isdir(first_arg) - # If the first argument is a directory, we read all the JSON files within - # and load them as parameter files; in the end all output will be saved - # in that same directory. We ignore the remaining elements in ARGS. - cd(first_arg) - files = filter(isjson, readdir()) - @info "$first_arg is a directory. Ignoring other command line arguments." - else - # Otherwise, all command line arguments are treated as parameter files. - # Output files will be saved in the pwd() from which the script was - # launched. - files = file_list - end - # Load parameters into dictionaries, one for each file. - parameter_lists = [] - for f in files - open(f) do input - s = read(input, String) - # Add the filename too to the dictionary. - push!(parameter_lists, merge(Dict("filename" => f), JSON.parse(s))) - end - end - cd(prev_dir) - return parameter_lists -end - -function allequal(a) - return all(x -> x == first(a), a) -end - -# Defining the time interval of the simulation -# ============================================ - -""" - construct_step_list(parameters) - -Return a list of time instants at which the time evolution will be evaluated. - -The values run from zero up to ``parameters["simulation_end_time"]``, with a -step size equal to ``parameters["simulation_time_step"]``. -""" -function construct_step_list(parameters) - τ = parameters["simulation_time_step"] - end_time = parameters["simulation_end_time"] - return collect(range(0, end_time; step=τ)) -end - -# Computing expectation values of observables -# =========================================== - -# Function that calculates the eigenvalues of the number operator, given a set -# of projectors on the eigenspaces, and also return their sum (which should -# be 1). -# TODO: these two functions are probably outdated. Anyway they are not -# specific to the occupation levels, so they should have a more general -# description. -function levels(projs::Vector{MPO}, state::MPS) - lev = [real(inner(state, p * state)) for p in projs] - return [lev; sum(lev)] -end -function levels(projs::Vector{MPS}, state::MPS) - lev = [real(inner(p, state)) for p in projs] - return [lev; sum(lev)] -end - # MPS and MPO utilities # ===================== - -""" - vectrace(vecρ::MPS, s::Vector{Index{Int64}}) - -Return the trace of the density matrix ρ which is vectorized and encoded in the -given MPS `vecρ` on sites `s`. -""" -function vectrace(vecρ::MPS, s::Vector{Index{Int64}}) - return dot(MPS("vecId", s), vecρ) -end - -""" - linkdims(m::Union{MPS, MPO}) - -Return a list of the bond dimensions of `m`. -""" -function linkdims(m::Union{MPS,MPO}) - return [ITensors.dim(linkind(m, j)) for j in 1:(length(m) - 1)] -end - """ chain(left::MPS, right::MPS) @@ -675,13 +268,13 @@ function chain(left::MPO, right::MPO) end # Varargs versions - """ chain(a::MPS, b...) Concatenate the given MPSs into a longer MPS, returning their tensor product. """ chain(a::MPS, b...) = chain(a, chain(b...)) + """ chain(a::MPO, b...) @@ -771,24 +364,6 @@ end # Other utilities # =============== - -""" - filenamett(::Dict) - -Extracts the filename from the parameter list, supplied as a Dict, and -formats it in a saw that's safe to use as text in TikZ nodes, plots, etc., -in typewriter font. -""" -function filenamett(d::Dict) - # Get basename and remove the .json extension. - filename = replace(basename(d["filename"]), ".json" => "") - # Sanitise string, escaping underscores. - filename = replace(filename, "_" => "\\_") - # Add \texttt command - filename = raw"\texttt{" * filename * "}" - return filename -end - """ consecutivepairs(v::AbstractVector)