From 34490686615de0f6b8fb0dbfc1d0ff97a32543b0 Mon Sep 17 00:00:00 2001 From: Davide Ferracin Date: Tue, 12 Mar 2024 00:43:53 +0200 Subject: [PATCH] Remove old, unused, too-specific functions Remove a bunch of functions that were sitting in PseudomodesTTEDOPA since its creation but were never used, and would much probably never be used in the future, such as functions for creating initial states for spin chains, plot utilities or generic-sounding "local operator" functions that actually returned very specific operators based on the SiteType of the Index. --- src/LindbladVectorizedTensors.jl | 46 +--- src/closure.jl | 14 +- src/current_operators.jl | 78 ------ src/operators.jl | 428 ------------------------------ src/site_types/fermion.jl | 4 +- src/utils.jl | 433 +------------------------------ 6 files changed, 19 insertions(+), 984 deletions(-) delete mode 100644 src/current_operators.jl 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)