From dc4ce88e02a0596b4101ca399124924959921a94 Mon Sep 17 00:00:00 2001 From: Davide Ferracin Date: Mon, 11 Mar 2024 13:57:05 +0200 Subject: [PATCH] Move files over from PseudomodesTTEDOPA --- Project.toml | 11 + docs/Project.toml | 2 + docs/make.jl | 2 + docs/src/gksl_example.md | 83 ++ docs/src/included_site_types.md | 186 ++++ docs/src/index.md | 54 +- src/LindbladVectorizedTensors.jl | 79 +- src/closure.jl | 739 ++++++++++++++++ src/current_operators.jl | 78 ++ src/exchange_interaction.jl | 279 ++++++ src/operators.jl | 535 ++++++++++++ src/site_types/electron.jl | 41 + src/site_types/fermion.jl | 23 + src/site_types/fermionic_dot3.jl | 152 ++++ src/site_types/oscillator.jl | 68 ++ src/site_types/spinhalf.jl | 17 + src/site_types/vectorized_electron.jl | 96 +++ src/site_types/vectorized_fermion.jl | 113 +++ src/site_types/vectorized_fermionic_dot3.jl | 147 ++++ src/site_types/vectorized_oscillator.jl | 294 +++++++ src/site_types/vectorized_spinhalf.jl | 158 ++++ src/spin_chain.jl | 220 +++++ src/utils.jl | 886 ++++++++++++++++++++ 23 files changed, 4254 insertions(+), 9 deletions(-) create mode 100644 docs/src/gksl_example.md create mode 100644 docs/src/included_site_types.md create mode 100644 src/closure.jl create mode 100644 src/current_operators.jl create mode 100644 src/exchange_interaction.jl create mode 100644 src/operators.jl create mode 100644 src/site_types/electron.jl create mode 100644 src/site_types/fermion.jl create mode 100644 src/site_types/fermionic_dot3.jl create mode 100644 src/site_types/oscillator.jl create mode 100644 src/site_types/spinhalf.jl create mode 100644 src/site_types/vectorized_electron.jl create mode 100644 src/site_types/vectorized_fermion.jl create mode 100644 src/site_types/vectorized_fermionic_dot3.jl create mode 100644 src/site_types/vectorized_oscillator.jl create mode 100644 src/site_types/vectorized_spinhalf.jl create mode 100644 src/spin_chain.jl create mode 100644 src/utils.jl diff --git a/Project.toml b/Project.toml index 5b96b7d..b9fc55c 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,19 @@ authors = ["Davide Ferracin and contributors"] version = "0.0.0" [compat] +ITensors = "0.3.19" julia = "1.6" +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" + [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/Project.toml b/docs/Project.toml index 9a12789..a466965 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" LindbladVectorizedTensors = "6455c26f-8985-4381-925b-8a93509e9dbd" +HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" diff --git a/docs/make.jl b/docs/make.jl index 309037f..a76af27 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,10 +1,12 @@ using LindbladVectorizedTensors +using ITensors using Documenter DocMeta.setdocmeta!(LindbladVectorizedTensors, :DocTestSetup, :(using LindbladVectorizedTensors); recursive=true) makedocs(; modules=[LindbladVectorizedTensors], + checkdocs = :none, authors="Davide Ferracin and contributors", sitename="LindbladVectorizedTensors.jl", format=Documenter.HTML(; diff --git a/docs/src/gksl_example.md b/docs/src/gksl_example.md new file mode 100644 index 0000000..6b53660 --- /dev/null +++ b/docs/src/gksl_example.md @@ -0,0 +1,83 @@ +# Application: GKSL equation + +With the states and operators defined by this package, we can define a Lindbladian +operator acting on density matrices in a GKSL equation. +Let's take, as an example, a system of ``N=3`` harmonic oscillators described by the +Hamiltonian + +```math +H= +\sum_{j=1}^{N}\omega_j a^\dagger_j a_j + +\sum_{j=1}^{N-1}\lambda_j (a^\dagger_j a_{j+1} + a^\dagger_{j+1} a_j) +``` + +and where each oscillator is damped according to the dissipator operator + +```math +\mathcal{D}(\rho)= +\sum_{j=1}^{N} \gamma_j \bigl( a_j \rho a^\dagger_j - \tfrac12 \rho a^\dagger_j a_j - +\tfrac12 a^\dagger_j a_j \rho \bigr). +``` + +First, create the site indices. In this example, we will use an arbitrary value for the +number of levels of the oscillators: + +```julia +N = 3 +s = siteinds("vOsc", 3; dim=5) +``` + +Create the OpSum object representing the unitary part of the evolution, ``-i[H,{-}]``: + +```julia +ℓ = OpSum() +for j in 1:N + ℓ += ω[j] * gkslcommutator("N", j) +end +for j in 1:(N - 1) + ℓ += + λ[j] * ( + gkslcommutator("Adag", site1, "A", site2) + + gkslcommutator("A", site1, "Adag", site2) + ) +end +``` + +Finally, add the dissipator: + +```julia +for j in 1:N + ℓ += γ, "A⋅ * ⋅Adag", n + ℓ += -γ / 2, "N⋅", n + ℓ += -γ / 2, "⋅N", n +end +``` + +This OpSum object can now, for example, be cast into an MPO and used to evolve a state with +the TDVP algorithm. Here is how ``L`` can be applied to a state ``ρ``: + +```julia +ρ = MPS([state("ThermEq", s[j]; temperature=1.0, frequency=5.0) for j in 1:N]) +L = MPO(ℓ, s) +apply(L, ρ) +``` + +Computing the expectation value of an operator ``A`` on a state ``ρ`` means +computing ``\operatorname{tr}(Aρ)``: with this package, this is done by contracting the +ITensor state representing the adjoint of the vectorized operator with the ITensor state +representing ``ρ``. +For example, if ``A`` is the number operator on site 3, i.e ``I\otimes I\otimes N``, +then we need the following code. + +```julia +vec_n = MPS(s, [i == 3 ? "vN" : "vId" for i in 1:N]) +result = dot(vec_n, ρ) +``` + +As a particular case, the trace of the state can be obtained by "measuring" the identity +operator: + +```julia +vec_id = MPS(s, "vId") +trace = dot(vec_id, ρ) +``` diff --git a/docs/src/included_site_types.md b/docs/src/included_site_types.md new file mode 100644 index 0000000..dc445b0 --- /dev/null +++ b/docs/src/included_site_types.md @@ -0,0 +1,186 @@ +# SiteTypes included with PseudomodesTTEDOPA + +This package extends the SiteType collection provided by ITensor, by defining types for +(vectorized) density matrices associated to built-in and operators acting on them. + +For each of these site types, each operator from the "original", non-vectorized ITensor site +type is promoted as a pre- or post-multiplication operator. +In other words, if there exists an operator "A" for the "Fermion" type, then you will +find operators "A⋅" and "⋅A" (the dot here is a `\cdot`) available for the "vFermion" +site type. + +Moreover, the `gkslcommutator` function allows you to easily define the unitary term (the +commutator) in the GKSL equation: given a list of operator names and site indices, +it returns an OpSum object representing the ``x\mapsto -i[A,x]`` map, for example +```julia-repl +julia> gkslcommutator("A", 1) +sum( + 0.0 - 1.0im A⋅(1,) + 0.0 + 1.0im ⋅A(1,) +) + +julia> gkslcommutator("A", 1, "B", 3) +sum( + 0.0 - 1.0im A⋅(1,) B⋅(3,) + 0.0 + 1.0im ⋅A(1,) ⋅B(3,) +) +``` + +Operators are also vectorized, with the aim of computing expectation values: if +``\{a_i\}_{i=1}^{d}`` and ``\{r_i\}_{i=1}^{d}`` are, respectively, the coordinates of +the operator ``A`` and the state ``\rho`` with respect to the Gell-Mann basis, then the +expectation value ``\operatorname{tr}(\rho A)`` is given by ``\sum_{i=1}^{d}\overline{a_i}r_i``. +Below you can find some operators transformed this way: they are actually "states" in the +ITensor formalism, even if they represent operators. +Their operator names are usually obtained by prefixing a `v` to the original ITensor name. + +## "vS=1/2" SiteType + +Site indices with the "vS=1/2" site type represent the density matrix of a ``S=1/2`` spin, +namely its coefficients with respect to the basis of 2x2 Gell-Mann matrices. + +Making a single "vS=1/2" site or collection of N "vS=1/2" sites +``` +s = siteind("vS=1/2") +sites = siteinds("vS=1/2", N) +``` + + +#### "vS=1/2" states + +The available state names for "vS=1/2" sites are: +- `"Up"` spin in the ``|{\uparrow}\rangle\langle{\uparrow}|`` state +- `"Dn"` spin in the ``|{\downarrow}\rangle\langle{\downarrow}|`` state + +#### "vS=1/2" vectorized operators + +Vectorized operators or associated with "vS=1/2" sites can be made using the `state` +function, for example +``` +Sx4 = state("vSx", s, 4) +N3 = state("vN", sites[3]) +``` + +Spin operators: +- `"vId"` Identity operator ``I_2`` +- `"vσx"` Pauli x matrix ``\sigma^x`` +- `"vσy"` Pauli y matrix ``\sigma^y`` +- `"vσz"` Pauli z matrix ``\sigma^z`` +- `"vSx"` Spin x operator ``S^x = \frac{1}{2} \sigma_x`` +- `"vSy"` Spin y operator ``S^y = \frac{1}{2} \sigma_y`` +- `"vSz"` Spin z operator ``S^z = \frac{1}{2} \sigma_z`` +- `"vN"` Number operator ``N = \frac{1}{2} (I_2+\sigma_z)`` + + +## "Osc" SiteType + +A SiteType representing a harmonic oscillator. It inherits definitions from the "Qudit" +type. + +Available keyword arguments for customization: +- `dim` (default: 2): dimension of the index (number of oscillator levels) +For example: +``` +sites = siteinds("Osc", N; dm=3) +``` + +#### "vOsc" Operators + +Operators associated with "Osc" sites: + +Single-oscillator operators: +- `"A"` (aliases: `"a"`) annihilation operator +- `"Adag"` (aliases: `"adag"`, `"a†"`) creation operator +- `"Asum"` equal to ``a+a^\dagger`` +- `"X"` equal to ``\frac{1}{\sqrt{2}}(a^\dagger+a)`` +- `"Y"` equal to ``\frac{i}{\sqrt{2}}(a^\dagger-a)`` +- `"N"` (aliases: `"n"`) number operator + + +## "vOsc" SiteType + +Available keyword arguments for customization: +- `dim` (default: 2): dimension of the index (number of oscillator levels) +For example: +``` +sites = siteinds("vOsc", N; dim=4) +``` + +#### "vOsc" states + +States associated with "vOsc" sites. +- "`n`", where `n` is an integer between `0` and `dim-1`, gives the Fock state + ``|n\rangle\langle n|`` +- `"ThermEq"`, with additional parameters `temperature` and `frequency`, gives the thermal + equilibrium (Gibbs) state ``\operatorname{tr}(\exp(-\frac{\omega}{T} +N))\exp(-\frac{\omega}{T} N)`` +- `"X⋅ThermEq"`, with additional parameters `temperature` and `frequency`, gives the thermal + equilibrium state (as in the previous point) multiplied by ``X=\frac{1}{\sqrt{2}}(a^\dagger+a)`` on the left (this is not actually a state, but it is useful when computing the correlation function of the ``X`` operator) + +Example: +```julia +s = siteind("vOsc", N; dim=4) +rho_eq = state("ThermEq", s; temperature=1.0, frequency=5.0) +fock_st = state("3", s) +``` + + +#### "vOsc" Operators + +Vectorized operators associated with "vOsc" sites. + +Single-oscillator operators (see the operator for "Osc" sites for their meaning): +- `"vA"` +- `"vAdag"` +- `"vN"` +- `"vX"` +- `"vY"` +- `"vId"` + + +## "vFermion" SiteType + +Site indices with the "vFermion" SiteType represent spinless fermion sites with the states +``|0\rangle``, ``|1\rangle``, corresponding to zero fermions or one fermion. + +#### "vFermion" states + +The available state names for "vFermion" sites are: +- `"Emp"` unoccupied fermion site +- `"Occ"` occupied fermion site + +#### "vFermion" operators + +Vectorized operators associated with "vFermion" sites: + +- `"vId"` Identity operator +- `"vN"` Density operator +- `"vA"` (aliases: `"va"`) Fermion annihilation operator +- `"vAdag"` (aliases: `"vadag"`, `"vA†"`, `"va†"`) Fermion creation operator + +Note that these creation and annihilation operators do not include Jordan-Wigner strings. + + +## "vElectron" SiteType + +The states of site indices with the "vElectron" SiteType correspond to +``|0\rangle``, ``|{\uparrow}\rangle``, ``|{\downarrow}\rangle``, ``|{\uparrow\downarrow}\rangle``. + +#### "vElectron" states + +The available state names for "vElectron" sites are: +- `"Emp"` unoccupied electron site +- `"Up"` electron site occupied with one up electron +- `"Dn"` electron site occupied with one down electron +- `"UpDn"` electron site occupied with two electrons (one up, one down) + +#### "vElectron" operators + +Vectorized operators associated with "vElectron" sites: + +- `"vId"` Identity operator +- `"vNtot"` Total density operator +- `"vNup"` Up density operator +- `"vNdn"` Down density operator +- `"vNupNdn"` Product of `n↑` and `n↓` + diff --git a/docs/src/index.md b/docs/src/index.md index 59a6c7d..0fd5b13 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,14 +1,52 @@ -```@meta -CurrentModule = LindbladVectorizedTensors -``` +# PseudomodesTTEDOPA.jl + +Documentation for PseudomodesTTEDOPA.jl -# LindbladVectorizedTensors +### Spin chain operators +```@docs +exchange_interaction +exchange_interaction_adjoint +exchange_interaction′ +spin_chain +spin_chain_adjoint +spin_chain′ +``` -Documentation for [LindbladVectorizedTensors](https://github.com/phaerrax/LindbladVectorizedTensors.jl). +### Pseudomode operators +```@docs +dissipator_loss +dissipator_gain +dissipator +mixedlindbladplus +mixedlindbladminus +``` -```@index +### Markovian closure operators +```@docs +Closure +``` +```@docs +closure +``` +```@docs +length(mc::Closure) +freqs(mc::Closure) +innercoups(mc::Closure) +outercoups(mc::Closure) +damps(mc::Closure) +freq(mc::Closure, j::Int) +innercoup(mc::Closure, j::Int) +outercoup(mc::Closure, j::Int) +damp(mc::Closure, j::Int) ``` -```@autodocs -Modules = [LindbladVectorizedTensors] +```@docs +closure_op(mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int) +closure_op_adjoint( + mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int, gradefactor::Int +) +closure_op′ +filled_closure_op +filled_closure_op_adjoint +filled_closure_op′ ``` diff --git a/src/LindbladVectorizedTensors.jl b/src/LindbladVectorizedTensors.jl index 7eba3fa..5cd3d5e 100644 --- a/src/LindbladVectorizedTensors.jl +++ b/src/LindbladVectorizedTensors.jl @@ -1,5 +1,82 @@ module LindbladVectorizedTensors -# Write your package code here. +using Base.Filesystem +using CSV +using Combinatorics +using DataFrames +using ITensors +using IterTools +using JSON +using LinearAlgebra +using ProgressMeter + +import ITensors.state, ITensors.op, ITensors.space + +export allequal, + 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, + sitenumber, + jwstring + +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") +include("site_types/vectorized_fermion.jl") +include("site_types/oscillator.jl") +include("site_types/electron.jl") +include("site_types/vectorized_electron.jl") +include("site_types/fermionic_dot3.jl") +include("site_types/vectorized_fermionic_dot3.jl") + +export dissipator_loss, dissipator_gain, dissipator, mixedlindbladplus, mixedlindbladminus + +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 + +include("operators.jl") + +include("spin_chain.jl") +include("exchange_interaction.jl") +include("closure.jl") end diff --git a/src/closure.jl b/src/closure.jl new file mode 100644 index 0000000..dcca4a4 --- /dev/null +++ b/src/closure.jl @@ -0,0 +1,739 @@ +export Closure, closure +export freq, innercoup, outercoup, damp, freqs, innercoups, outercoups, damps +export closure_op, closure_op_adjoint, closure_op′ +export filled_closure_op, filled_closure_op_adjoint, filled_closure_op′ + +""" +A Closure object stores the parameters that define a Markovian closure, providing different +ways to input them. +""" +struct Closure + frequency + innercoupling + outercoupling + damping + @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. + """ + function Closure( + ω::Vector{<:Real}, γ::Vector{<:Real}, g::Vector{<:Real}, ζ::Vector{<:Complex} + ) + if (length(ω) - 1 != length(g) || length(ω) != length(γ) || length(ω) != length(ζ)) + error("Lengths of input parameters do not match.") + end + return new(ω, g, ζ, γ) + end +end + +""" + closure( + Ω::Real, K::Real, + α::Matrix{<:Real}, β::Matrix{<:Real}, w::Matrix{<:Real} + ) + +Construct a Closure object with asymptotic frequency `Ω` and coupling `K`, with the +universal parameter set given by `α`, `β` and `w`, each of which is assumed to be a +two-column matrix with the real parts of the parameters in the first column and the +imaginary parts in the second one. +""" +function closure(Ω::Real, K::Real, α::Matrix{<:Real}, β::Matrix{<:Real}, w::Matrix{<:Real}) + return closure( + Ω, K, α[:, 1] .+ im .* α[:, 2], β[:, 1] .+ im .* β[:, 2], w[:, 1] .+ im .* w[:, 2] + ) +end + +""" + closure( + Ω::Real, K::Real, + α::Vector{<:Complex}, β::Vector{<:Complex}, w::Vector{<:Complex} + ) + +Construct a Closure object with asymptotic frequency `Ω` and coupling `K`, with the +universal parameter set given by `α`, `β` and `w`, which are assumed to be vectors of +complex numbers. +""" +function closure( + Ω::Real, K::Real, α::Vector{<:Complex}, β::Vector{<:Complex}, w::Vector{<:Complex} +) + frequency = @. Ω - 2K * imag(α) + damping = @. -4K * real(α) + innercoupling = @. -2K * imag(β) + outercoupling = @. K * w + return Closure(frequency, damping, innercoupling, outercoupling) +end + +""" + length(mc::Closure) + +Return the length of a Markovian closure, i.e. the number of pseudomodes. +""" +Base.length(mc::Closure) = Base.length(mc.frequency) + +""" + freqs(mc::Closure) + +Return the list of frequencies of the pseudomodes in the given Markovian closure. +""" +freqs(mc::Closure) = mc.frequency + +""" + innercoups(mc::Closure) + +Return the list of inner couplings of the pseudomodes in the Markovian closure `mc`, +i.e. the coupling constant between the pseudomodes. +""" +innercoups(mc::Closure) = mc.innercoupling + +""" + outercoups(mc::Closure) + +Return the list of outer couplings of the pseudomodes in the Markovian closure `mc`, +i.e. the coupling constants with the last point of the TEDOPA chain. +""" +outercoups(mc::Closure) = mc.outercoupling + +""" + damps(mc::Closure) + +Return the list of damping coefficients of the pseudomodes in the Markovian closure `mc`. +""" +damps(mc::Closure) = mc.damping + +""" + freq(mc::Closure, j::Int) + +Return the frequency of the `j`-th pseudomode of the Markovian closure `mc`. +""" +freq(mc::Closure, j::Int) = mc.frequency[j] + +""" + innercoup(mc::Closure, j::Int) + +Return the coupling coefficient between the `j`-th and the `j+1`-th pseudomodes of the +Markovian closure `mc`. +""" +innercoup(mc::Closure, j::Int) = mc.innercoupling[j] + +""" + outercoup(mc::Closure, j::Int) + +Return the coupling coefficient of the `j`-th pseudomode in the Markovian closure `mc`, +with the last point of the TEDOPA chain. +""" +outercoup(mc::Closure, j::Int) = mc.outercoupling[j] + +""" + damp(mc::Closure, j::Int) + +Return the damping coefficient of the `j`-th pseudomode in the Markovian closure `mc`. +""" +damp(mc::Closure, j::Int) = mc.damping[j] + +############################################################################################ + +closure_op(::SiteType, ::Closure, ::Vector{<:Int}, ::Int) = nothing + +""" + closure_op(mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int) + +Return an OpSum object encoding the Markovian closure operators with parameters given by +`mc`, on sites `sites`, linked to the main TEDOPA/thermofield chain on site +`chain_edge_site`. +This closure replaces a chain starting from an empty state. +""" +function closure_op(mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int) + @assert length(mc) == length(sites) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = closure_op(st, mc, sitenumber.(sites), chain_edge_site) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"closure_op\" function not found for Index tags $(tags(sites[1]))" + ), + ) +end + +function closure_op( + st::SiteType"vFermion", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + ℓ += ( + outercoup(mc, j) * gkslcommutator("A†", chain_edge_site, jws..., "A", site) + + conj(outercoup(mc, j)) * + gkslcommutator("A", chain_edge_site, jws..., "A†", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a ρ a† + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "A⋅ * ⋅A†"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + # -0.5 (a† a ρ + ρ a† a) + ℓ += -0.5damp(mc, j), "N⋅", site + ℓ += -0.5damp(mc, j), "⋅N", site + end + return ℓ +end + +function closure_op( + st::SiteType"vS=1/2", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + ℓ += ( + outercoup(mc, j) * gkslcommutator("σ+", chain_edge_site, "σ-", site) + + conj(outercoup(mc, j)) * gkslcommutator("σ-", chain_edge_site, "σ+", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a ρ a† + ℓ += damp(mc, j), "σ-⋅ * ⋅σ+", site + # -0.5 (a† a ρ + ρ a† a) + ℓ += -0.5damp(mc, j), "N⋅", site + ℓ += -0.5damp(mc, j), "⋅N", site + end + return ℓ +end + +function closure_op( + st::SiteType"vElectron", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + # c↑ᵢ† c↑ᵢ₊ₙ = a↑ᵢ† Fᵢ Fᵢ₊₁ ⋯ Fᵢ₊ₙ₋₁ a↑ᵢ₊ₙ + # c↑ᵢ₊ₙ† c↑ᵢ = -a↑ᵢ Fᵢ Fᵢ₊₁ ⋯ Fᵢ₊ₙ₋₁ a↑ᵢ₊ₙ† + # c↓ᵢ† c↓ᵢ₊ₙ = a↓ᵢ† Fᵢ₊₁ Fᵢ₊₂ ⋯ Fᵢ₊ₙ a↓ᵢ₊ₙ + # c↓ᵢ₊ₙ† c↓ᵢ = -a↓ᵢ Fᵢ₊₁ Fᵢ₊₂ ⋯ Fᵢ₊ₙ a↓ᵢ₊ₙ† + + jws = jwstring(; start=chain_edge_site, stop=site) + # ζⱼ c↑₀† c↑ⱼ (0 = chain edge, j = pseudomode) + ℓ += + outercoup(mc, j) * gkslcommutator("Aup†F", chain_edge_site, jws..., "Aup", site) + # conj(ζⱼ) c↑ⱼ† c↑₀ + ℓ += + -conj(outercoup(mc, j)) * + gkslcommutator("AupF", chain_edge_site, jws..., "Aup†", site) + # ζⱼ c↓₀† c↓ⱼ + ℓ += + outercoup(mc, j) * gkslcommutator("Adn†", chain_edge_site, jws..., "FAdn", site) + # conj(ζⱼ) c↓ⱼ† c↓₀ + ℓ += + -conj(outercoup(mc, j)) * + gkslcommutator("Adn", chain_edge_site, jws..., "FAdn†", site) + end + + # Dissipative part + for (j, site) in enumerate(sitenumbers) + # c↑ₖ ρ c↑ₖ† = F₁ ⋯ Fₖ₋₁ a↑ₖ ρ a↑ₖ† Fₖ₋₁ ⋯ F₁ + # Remember that: + # • Fⱼ = (1 - 2 N↑ₖ) (1 - 2 N↓ₖ); + # • Fⱼ and aₛ,ₖ commute only on different sites; + # • {a↓ₖ, a↓ₖ†} = {a↑ₖ, a↑ₖ†} = 1; + # • Fₖ anticommutes with a↓ₖ, a↓ₖ†, a↑ₖ and a↑ₖ†. + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Aup⋅ * ⋅Aup†"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + # c↓ₖ ρ c↓ₖ† = F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ ρ a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "FAdn⋅ * ⋅Adn†F"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + + # -½ (c↑ₖ† c↑ₖ ρ + ρ c↑ₖ† c↑ₖ) = -½ (a↑ₖ† a↑ₖ ρ + ρ a↑ₖ† a↑ₖ) + ℓ += -0.5damp(mc, j), "Nup⋅", site + ℓ += -0.5damp(mc, j), "⋅Nup", site + # -½ (c↓ₖ† c↓ₖ ρ + ρ c↓ₖ† c↓ₖ) = -½ (a↓ₖ† a↓ₖ ρ + ρ a↓ₖ† a↓ₖ) + ℓ += -0.5damp(mc, j), "Ndn⋅", site + ℓ += -0.5damp(mc, j), "⋅Ndn", site + end + return ℓ +end + +############################################################################################ + +closure_op_adjoint(::SiteType, ::Closure, ::Vector{<:Int}, ::Int, ::Int) = nothing + +""" + closure_op_adjoint( + mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int, gradefactor::Int + ) + +Return an OpSum object encoding the adjoint of the Markovian closure operators with +parameters given by `mc`, on sites `sites`, linked to the main TEDOPA/thermofield chain on +site `chain_edge_site`. The integer `gradefactor` is the parity (1: even, -1: odd) of the +operator subject to the time evolution. +This closure replaces a chain starting from an empty state. +""" +function closure_op_adjoint( + mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int, gradefactor::Int +) + @assert length(mc) == length(sites) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = closure_op_adjoint(st, mc, sitenumber.(sites), chain_edge_site, gradefactor) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"closure_op\" function not found for Index tags $(tags(sites[1]))" + ), + ) +end + +function closure_op_adjoint( + st::SiteType"vFermion", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = spin_chain_adjoint(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + ℓ += -( + outercoup(mc, j) * gkslcommutator("A†", chain_edge_site, jws..., "A", site) + + conj(outercoup(mc, j)) * + gkslcommutator("A", chain_edge_site, jws..., "A†", site) + ) + end + + 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)))... + ) + # -0.5 (a† a ρ + ρ a† a) + ℓ += -0.5damp(mc, j), "N⋅", site + ℓ += -0.5damp(mc, j), "⋅N", site + end + return ℓ +end + +function closure_op_adjoint( + st::SiteType"vS=1/2", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = spin_chain_adjoint(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + ℓ += -( + outercoup(mc, j) * gkslcommutator("σ+", chain_edge_site, "σ-", site) + + conj(outercoup(mc, j)) * gkslcommutator("σ-", chain_edge_site, "σ+", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a† ρ a + ℓ += gradefactor * damp(mc, j), "σ+⋅ * ⋅σ-", site + # -0.5 (a† a ρ + ρ a† a) + ℓ += -0.5damp(mc, j), "N⋅", site + ℓ += -0.5damp(mc, j), "⋅N", site + end + return ℓ +end + +function closure_op_adjoint( + st::SiteType"vElectron", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = spin_chain_adjoint(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + # c↑ᵢ† c↑ᵢ₊ₙ = a↑ᵢ† Fᵢ Fᵢ₊₁ ⋯ Fᵢ₊ₙ₋₁ a↑ᵢ₊ₙ + # c↑ᵢ₊ₙ† c↑ᵢ = -a↑ᵢ Fᵢ Fᵢ₊₁ ⋯ Fᵢ₊ₙ₋₁ a↑ᵢ₊ₙ† + # c↓ᵢ† c↓ᵢ₊ₙ = a↓ᵢ† Fᵢ₊₁ Fᵢ₊₂ ⋯ Fᵢ₊ₙ a↓ᵢ₊ₙ + # c↓ᵢ₊ₙ† c↓ᵢ = -a↓ᵢ Fᵢ₊₁ Fᵢ₊₂ ⋯ Fᵢ₊ₙ a↓ᵢ₊ₙ† + + jws = jwstring(; start=chain_edge_site, stop=site) + # ζⱼ c↑₀† c↑ⱼ (0 = chain edge, j = pseudomode) + ℓ += + -outercoup(mc, j) * + gkslcommutator("Aup†F", chain_edge_site, jws..., "Aup", site) + # conj(ζⱼ) c↑ⱼ† c↑₀ + ℓ += + conj(outercoup(mc, j)) * + gkslcommutator("AupF", chain_edge_site, jws..., "Aup†", site) + # ζⱼ c↓₀† c↓ⱼ + ℓ += + -outercoup(mc, j) * + gkslcommutator("Adn†", chain_edge_site, jws..., "FAdn", site) + # conj(ζⱼ) c↓ⱼ† c↓₀ + ℓ += + conj(outercoup(mc, j)) * + gkslcommutator("Adn", chain_edge_site, jws..., "FAdn†", site) + end + + # Dissipative part + for (j, site) in enumerate(sitenumbers) + # Remember that: + # • Fⱼ = (1 - 2 N↑ₖ) (1 - 2 N↓ₖ); + # • Fⱼ and aₛ,ₖ commute only on different sites; + # • {a↓ₖ, a↓ₖ†} = {a↑ₖ, a↑ₖ†} = 1; + # • Fₖ anticommutes with a↓ₖ, a↓ₖ†, a↑ₖ and a↑ₖ†. + + # c↑ₖ† X c↑ₖ = a↑ₖ† Fₖ₋₁ ⋯ F₁ X F₁ ⋯ Fₖ₋₁ a↑ₖ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Aup†⋅ * ⋅Aup"] + ℓ += ( + gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))... + ) + # c↓ₖ† X c↓ₖ = a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ X F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Adn†F⋅ * ⋅FAdn"] + ℓ += ( + gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))... + ) + + # -½ (c↑ₖ† c↑ₖ ρ + ρ c↑ₖ† c↑ₖ) = -½ (a↑ₖ† a↑ₖ ρ + ρ a↑ₖ† a↑ₖ) + ℓ += -0.5damp(mc, j), "Nup⋅", site + ℓ += -0.5damp(mc, j), "⋅Nup", site + # -½ (c↓ₖ† c↓ₖ ρ + ρ c↓ₖ† c↓ₖ) = -½ (a↓ₖ† a↓ₖ ρ + ρ a↓ₖ† a↓ₖ) + ℓ += -0.5damp(mc, j), "Ndn⋅", site + ℓ += -0.5damp(mc, j), "⋅Ndn", site + end + return ℓ +end + +const closure_op′ = closure_op_adjoint + +############################################################################################ + +filled_closure_op(::SiteType, ::Closure, ::Vector{<:Int}, ::Int) = nothing + +""" + filled_closure_op(mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int) + +Return an OpSum object encoding the Markovian closure operators with parameters given by +`mc`, on sites `sites`, linked to the main TEDOPA/thermofield chain on site +`chain_edge_site`. +This closure replaces a chain starting from a completely filled state. +""" +function filled_closure_op(mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int) + @assert length(mc) == length(sites) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = filled_closure_op(st, mc, sitenumber.(sites), chain_edge_site) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"filled_closure_op\" function not found for " * + "Index tags $(tags(sites[1]))", + ), + ) +end + +function filled_closure_op( + st::SiteType"vFermion", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + ℓ += ( + outercoup(mc, j) * gkslcommutator("A†", chain_edge_site, jws..., "A", site) + + conj(outercoup(mc, j)) * + gkslcommutator("A", chain_edge_site, jws..., "A†", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a ρ a† + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "A†⋅ * ⋅A"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + # -0.5 (a a† ρ + ρ a a†) + ℓ += 0.5damp(mc, j), "N⋅", site + ℓ += 0.5damp(mc, j), "⋅N", site + ℓ += -damp(mc, j), "Id", site + end + return ℓ +end + +function filled_closure_op( + st::SiteType"vS=1/2", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + ℓ += ( + outercoup(mc, j) * gkslcommutator("σ+", chain_edge_site, "σ-", site) + + conj(outercoup(mc, j)) * gkslcommutator("σ-", chain_edge_site, "σ+", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a† ρ a + ℓ += damp(mc, j), "σ+⋅ * ⋅σ-", site + # -0.5 (a a† ρ + ρ a a†) + ℓ += 0.5damp(mc, j), "N⋅", site + ℓ += 0.5damp(mc, j), "⋅N", site + ℓ += -damp(mc, j), "Id", site + end + return ℓ +end + +function filled_closure_op( + st::SiteType"vElectron", mc::Closure, sitenumbers::Vector{<:Int}, chain_edge_site::Int +) + ℓ = spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + # ζⱼ c↑₀† c↑ⱼ (0 = chain edge, j = pseudomode) + ℓ += + outercoup(mc, j) * gkslcommutator("Aup†F", chain_edge_site, jws..., "Aup", site) + # conj(ζⱼ) c↑ⱼ† c↑₀ + ℓ += + -conj(outercoup(mc, j)) * + gkslcommutator("AupF", chain_edge_site, jws..., "Aup†", site) + # ζⱼ c↓₀† c↓ⱼ + ℓ += + outercoup(mc, j) * gkslcommutator("Adn†", chain_edge_site, jws..., "FAdn", site) + # conj(ζⱼ) c↓ⱼ† c↓₀ + ℓ += + -conj(outercoup(mc, j)) * + gkslcommutator("Adn", chain_edge_site, jws..., "FAdn†", site) + end + + # Dissipative part + for (j, site) in enumerate(sitenumbers) + # c↑ₖ† ρ c↑ₖ = a↑ₖ† Fₖ₋₁ ⋯ F₁ ρ F₁ ⋯ Fₖ₋₁ a↑ₖ + # Remember that: + # • Fⱼ = (1 - 2 N↑ₖ) (1 - 2 N↓ₖ); + # • Fⱼ and aₛ,ₖ commute only on different sites; + # • {a↓ₖ, a↓ₖ†} = {a↑ₖ, a↑ₖ†} = 1; + # • Fₖ anticommutes with a↓ₖ, a↓ₖ†, a↑ₖ and a↑ₖ†. + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Aup†⋅ * ⋅Aup"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + # c↓ₖ† ρ c↓ₖ = a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ ρ F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Adn†F⋅ * ⋅FAdn"] + ℓ += (damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))...) + + # c↑ₖ c↑ₖ† = F₁ ⋯ Fₖ₋₁ a↑ₖ a↑ₖ† Fₖ₋₁ ⋯ F₁ = + # = F₁² ⋯ Fₖ₋₁² a↑ₖ a↑ₖ† = + # = a↑ₖ a↑ₖ† = + # = 1 - N↑ₖ + ℓ += -damp(mc, j), "Id", site + ℓ += 0.5damp(mc, j), "Nup⋅", site + ℓ += 0.5damp(mc, j), "⋅Nup", site + # c↓ₖ c↓ₖ† = F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ = + # = F₁² ⋯ Fₖ₋₁² Fₖa↓ₖ a↓ₖ†Fₖ = + # = Fₖa↓ₖ a↓ₖ†Fₖ = + # = Fₖ² a↓ₖ a↓ₖ† = + # = 1 - N↓ₖ + ℓ += -damp(mc, j), "Id", site + ℓ += 0.5damp(mc, j), "Ndn⋅", site + ℓ += 0.5damp(mc, j), "⋅Ndn", site + end + return ℓ +end + +############################################################################################ + +filled_closure_op_adjoint(::SiteType, ::Closure, ::Vector{<:Int}, ::Int, ::Int) = nothing + +""" + filled_closure_op_adjoint( + mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int, gradefactor::Int + ) + +Return an OpSum object encoding the Markovian closure operators with parameters given by +`mc`, on sites `sites`, linked to the main TEDOPA/thermofield chain on site +`chain_edge_site`. The integer `gradefactor` is the parity (1: even, -1: odd) of the +operator subject to the time evolution. +This closure replaces a chain starting from a completely filled state. +""" +function filled_closure_op_adjoint( + mc::Closure, sites::Vector{<:Index}, chain_edge_site::Int, gradefactor::Int +) + @assert length(mc) == length(sites) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = filled_closure_op_adjoint( + st, mc, sitenumber.(sites), chain_edge_site, gradefactor + ) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"filled_closure_op\" function not found for " * + "Index tags $(tags(sites[1]))", + ), + ) +end + +function filled_closure_op_adjoint( + st::SiteType"vFermion", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = adjoint_spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + ℓ += -( + outercoup(mc, j) * gkslcommutator("A†", chain_edge_site, jws..., "A", site) + + conj(outercoup(mc, j)) * + gkslcommutator("A", chain_edge_site, jws..., "A†", site) + ) + end + + 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)))...) + # -0.5 (a a† ρ + ρ a a†) = -0.5 (ρ - a† a ρ + ρ a† a) + ℓ += 0.5damp(mc, j), "N⋅", site + ℓ += 0.5damp(mc, j), "⋅N", site + ℓ += -damp(mc, j), "Id", site + end + return ℓ +end + +function filled_closure_op_adjoint( + st::SiteType"vS=1/2", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = adjoint_spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + ℓ += -( + outercoup(mc, j) * gkslcommutator("σ+", chain_edge_site, "σ-", site) + + conj(outercoup(mc, j)) * gkslcommutator("σ-", chain_edge_site, "σ+", site) + ) + end + + for (j, site) in enumerate(sitenumbers) + # a ρ a† + ℓ += damp(mc, j), "σ-⋅ * ⋅σ+", site + # -0.5 (a a† ρ + ρ a a†) = -0.5 (2ρ - a† a ρ + ρ a† a) + ℓ += 0.5damp(mc, j), "N⋅", site + ℓ += 0.5damp(mc, j), "⋅N", site + ℓ += -damp(mc, j), "Id", site + end + return ℓ +end + +function filled_closure_op_adjoint( + st::SiteType"vElectron", + mc::Closure, + sitenumbers::Vector{<:Int}, + chain_edge_site::Int, + gradefactor::Int, +) + ℓ = adjoint_spin_chain(st, freqs(mc), innercoups(mc), sitenumbers) + + for (j, site) in enumerate(sitenumbers) + jws = jwstring(; start=chain_edge_site, stop=site) + # ζⱼ c↑₀† c↑ⱼ (0 = chain edge, j = pseudomode) + ℓ += + -outercoup(mc, j) * + gkslcommutator("Aup†F", chain_edge_site, jws..., "Aup", site) + # conj(ζⱼ) c↑ⱼ† c↑₀ + ℓ += + conj(outercoup(mc, j)) * + gkslcommutator("AupF", chain_edge_site, jws..., "Aup†", site) + # ζⱼ c↓₀† c↓ⱼ + ℓ += + -outercoup(mc, j) * + gkslcommutator("Adn†", chain_edge_site, jws..., "FAdn", site) + # conj(ζⱼ) c↓ⱼ† c↓₀ + ℓ += + conj(outercoup(mc, j)) * + gkslcommutator("Adn", chain_edge_site, jws..., "FAdn†", site) + end + + # Dissipative part + for (j, site) in enumerate(sitenumbers) + # Remember that: + # • Fⱼ = (1 - 2 N↑ₖ) (1 - 2 N↓ₖ); + # • Fⱼ and aₛ,ₖ commute only on different sites; + # • {a↓ₖ, a↓ₖ†} = {a↑ₖ, a↑ₖ†} = 1; + # • Fₖ anticommutes with a↓ₖ, a↓ₖ†, a↑ₖ and a↑ₖ†. + # + # c↑ₖ X c↑ₖ† = F₁ ⋯ Fₖ₋₁ a↑ₖ X a↑ₖ† Fₖ₋₁ ⋯ F₁ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "Aup⋅ * ⋅Aup†"] + ℓ += ( + gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))... + ) + # c↓ₖ X c↓ₖ† = F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ X a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ + opstring = [repeat(["F⋅ * ⋅F"], site - 1); "FAdn⋅ * ⋅Adn†F"] + ℓ += ( + gradefactor * damp(mc, j), collect(Iterators.flatten(zip(opstring, 1:site)))... + ) + + # c↑ₖ c↑ₖ† = F₁ ⋯ Fₖ₋₁ a↑ₖ a↑ₖ† Fₖ₋₁ ⋯ F₁ = + # = F₁² ⋯ Fₖ₋₁² a↑ₖ a↑ₖ† = + # = a↑ₖ a↑ₖ† = + # = 1 - N↑ₖ + ℓ += -damp(mc, j), "Id", site + ℓ += 0.5damp(mc, j), "Nup⋅", site + ℓ += 0.5damp(mc, j), "⋅Nup", site + # c↓ₖ c↓ₖ† = F₁ ⋯ Fₖ₋₁ Fₖa↓ₖ a↓ₖ†Fₖ Fₖ₋₁ ⋯ F₁ = + # = F₁² ⋯ Fₖ₋₁² Fₖa↓ₖ a↓ₖ†Fₖ = + # = Fₖa↓ₖ a↓ₖ†Fₖ = + # = Fₖ² a↓ₖ a↓ₖ† = + # = 1 - N↓ₖ + ℓ += -damp(mc, j), "Id", site + ℓ += 0.5damp(mc, j), "Ndn⋅", site + ℓ += 0.5damp(mc, j), "⋅Ndn", site + end + return ℓ +end + +const filled_closure_op′ = filled_closure_op_adjoint diff --git a/src/current_operators.jl b/src/current_operators.jl new file mode 100644 index 0000000..c1de6e0 --- /dev/null +++ b/src/current_operators.jl @@ -0,0 +1,78 @@ +# 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/exchange_interaction.jl b/src/exchange_interaction.jl new file mode 100644 index 0000000..6afe004 --- /dev/null +++ b/src/exchange_interaction.jl @@ -0,0 +1,279 @@ +export exchange_interaction, exchange_interaction_adjoint, exchange_interaction′ + +exchange_interaction(::SiteType, ::SiteType, ::Int, ::Int; kwargs...) = nothing + +""" + exchange_interaction(s1::Index, s2::Index; kwargs...) + +Return an OpSum object encoding the Hamiltonian part ``-i[H, –]`` of an exchange interaction +between sites `s1` and `s2` term in a GKSL equation. +""" +function exchange_interaction(s1::Index, s2::Index; kwargs...) + for (st1, st2) in zip(sitetypes(s1), sitetypes(s2)) + ℓ = exchange_interaction(st1, st2, sitenumber(s1), sitenumber(s2); kwargs...) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"exchange_interaction\" function not found for " * + "Index tags $(tags(s1)) and $(tags(s2))", + ), + ) +end + +function exchange_interaction( + ::SiteType"Fermion", + ::SiteType"Fermion", + site1::Int, + site2::Int; + coupling_constant::Real=1.0, +) + h = OpSum() + + h += coupling_constant, "c†", site1, "c", site2 + h += coupling_constant, "c†", site2, "c", site1 + + return h +end + +function exchange_interaction( + ::SiteType"vFermion", + ::SiteType"vFermion", + site1::Int, + site2::Int; + coupling_constant::Real=1.0, +) + ℓ = OpSum() + jws = jwstring(; start=site1, stop=site2) + ℓ += ( + coupling_constant * gkslcommutator("a†", site1, jws..., "a", site2) + + coupling_constant * gkslcommutator("a", site1, jws..., "a†", site2) + ) + return ℓ +end + +function exchange_interaction( + ::SiteType"vS=1/2", + ::SiteType"vS=1/2", + site1::Int, + site2::Int; + coupling_constant::Real=1.0, +) + ℓ = OpSum() + ℓ += ( + coupling_constant * gkslcommutator("σ+", site1, "σ-", site2) + + coupling_constant * gkslcommutator("σ-", site1, "σ+", site2) + ) + return ℓ +end + +function exchange_interaction( + ::SiteType"Electron", + ::SiteType"Electron", + site1::Int, + site2::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + h = OpSum() + + h += coupling_constant_up, "c†↑", site1, "c↑", site2 + h += coupling_constant_up, "c†↑", site2, "c↑", site1 + + h += coupling_constant_dn, "c†↓", site1, "c↓", site2 + h += coupling_constant_dn, "c†↓", site2, "c↓", site1 + + return h +end + +function exchange_interaction( + ::SiteType"vElectron", + ::SiteType"vElectron", + site1::Int, + site2::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + ℓ = OpSum() + jws = jwstring(; start=site1, stop=site2) + + ℓ += coupling_constant_up * gkslcommutator("Aup†F", site1, jws..., "Aup", site2) + ℓ += -coupling_constant_up * gkslcommutator("AupF", site1, jws..., "Aup†", site2) + + ℓ += coupling_constant_dn * gkslcommutator("Adn†", site1, jws..., "FAdn", site2) + ℓ += -coupling_constant_dn * gkslcommutator("Adn", site1, jws..., "FAdn†", site2) + # Look out for minus signs: + # (F * a↓)† = a†↓ * F† = -F * a†↓ + # (a†↑ * F)† = F† * a↑ = F * a↑ = -a↑ * F + + return ℓ +end + +function exchange_interaction( + ::SiteType"Electron", + ::SiteType"Fermion", + electron_site::Int, + fermion_site::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + h = OpSum() + jws = jwstring(; start=electron_site, stop=fermion_site) + + h += coupling_constant_up, "a†↑ * F", electron_site, jws..., "a", fermion_site + h += -coupling_constant_up, "a↑ * F", electron_site, jws..., "a†", fermion_site + # Look out for the minus sign: (a†↑ * F)† = F† * a↑ = F * a↑ = -a↑ * F + + h += coupling_constant_dn, "a†↓", electron_site, jws..., "a", fermion_site + h += coupling_constant_dn, "a↓", electron_site, jws..., "a†", fermion_site + + return h +end + +function exchange_interaction( + ::SiteType"vElectron", + ::SiteType"vFermion", + electron_site::Int, + fermion_site::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + ℓ = OpSum() + jws = jwstring(; start=electron_site, stop=fermion_site) + + ℓ += + coupling_constant_up * + gkslcommutator("Aup†F", electron_site, jws..., "a", fermion_site) + ℓ += + -coupling_constant_up * + gkslcommutator("AupF", electron_site, jws..., "a†", fermion_site) + # Look out for the minus sign: (a†↑ * F)† = F† * a↑ = F * a↑ = -a↑ * F + + ℓ += + coupling_constant_dn * + gkslcommutator("Adn†", electron_site, jws..., "a", fermion_site) + ℓ += + coupling_constant_dn * + gkslcommutator("Adn", electron_site, jws..., "a†", fermion_site) + + return ℓ +end + +function exchange_interaction( + ::SiteType"Fermion", + ::SiteType"Electron", + fermion_site::Int, + electron_site::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + h = OpSum() + jws = jwstring(; start=electron_site, stop=fermion_site) + + h += coupling_constant_up, "a†", fermion_site, jws..., "a↑", electron_site + h += coupling_constant_up, "a", fermion_site, jws..., "a†↑", electron_site + + h += coupling_constant_dn, "a†", fermion_site, jws..., "F * a↓", electron_site + h += -coupling_constant_dn, "a", fermion_site, jws..., "F * a†↓", electron_site + # Look out for the minus sign: (F * a↓)† = a†↓ * F† = -F * a†↓ + + return h +end + +function exchange_interaction( + ::SiteType"vFermion", + ::SiteType"vElectron", + fermion_site::Int, + electron_site::Int; + coupling_constant_up::Real=1.0, + coupling_constant_dn::Real=1.0, +) + ℓ = OpSum() + jws = jwstring(; start=electron_site, stop=fermion_site) + + ℓ += + coupling_constant_up * + gkslcommutator("a†", fermion_site, jws..., "Aup", electron_site) + ℓ += + coupling_constant_up * + gkslcommutator("a", fermion_site, jws..., "Aup†", electron_site) + + ℓ += + coupling_constant_dn * + gkslcommutator("a†", fermion_site, jws..., "FAdn", electron_site) + ℓ += + -coupling_constant_dn * + gkslcommutator("a", fermion_site, jws..., "FAdn†", electron_site) + # Look out for the minus sign: (F * a↓)† = a†↓ * F† = -F * a†↓ + + return ℓ +end + +############################################################################################ + +exchange_interaction_adjoint(::SiteType, ::SiteType, ::Int, ::Int; kwargs...) = nothing + +""" + exchange_interaction_adjoint(s1::Index, s2::Index; kwargs...) + +Return an OpSum object encoding the adjoint of the Hamiltonian part ``-i[H, –]`` of an +exchange interaction between sites `s1` and `s2` term in a GKSL equation. +""" +function exchange_interaction_adjoint(s1::Index, s2::Index; kwargs...) + for (st1, st2) in zip(sitetypes(s1), sitetypes(s2)) + ℓ = exchange_interaction_adjoint( + st1, st2, sitenumber(s1), sitenumber(s2); kwargs... + ) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"exchange_interaction_adjoint\" function not found for " * + "Index tags $(tags(s1)) and $(tags(s2))", + ), + ) +end + +# In the GKSL equation, the adjoint of -i[H,-] is simply i[H,-]. + +function exchange_interaction_adjoint( + st1::SiteType"vFermion", st2::SiteType"vFermion", site1::Int, site2::Int; kwargs... +) + return -exchange_interaction(st1, st2, site1, site2; kwargs...) +end + +function exchange_interaction_adjoint( + st1::SiteType"vS=1/2", st2::SiteType"vS=1/2", site1::Int, site2::Int; kwargs... +) + return -exchange_interaction(st1, st2, site1, site2; kwargs...) +end + +function exchange_interaction_adjoint( + st1::SiteType"vElectron", st2::SiteType"vElectron", site1::Int, site2::Int; kwargs... +) + return -exchange_interaction(st1, st2, site1, site2; kwargs...) +end + +function exchange_interaction_adjoint( + st1::SiteType"vElectron", st2::SiteType"vFermion", site1::Int, site2::Int; kwargs... +) + return -exchange_interaction(st1, st2, site1, site2; kwargs...) +end + +function exchange_interaction_adjoint( + st1::SiteType"vFermion", st2::SiteType"vElectron", site1::Int, site2::Int; kwargs... +) + return -exchange_interaction(st1, st2, site1, site2; kwargs...) +end + +const exchange_interaction′ = exchange_interaction_adjoint diff --git a/src/operators.jl b/src/operators.jl new file mode 100644 index 0000000..ec77769 --- /dev/null +++ b/src/operators.jl @@ -0,0 +1,535 @@ +@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...) + isodd(length(args)) && + throw(error("Cannot make pairs out of an odd number of arguments.")) + return [(args[i], args[i + 1]) for i in eachindex(args)[1:2:(end - 1)]] +end + +""" + gkslcommutator(operators::Tuple{String,Int}...) + +Given one or more tuples `(a1, n1), (a2, n2), …`, return an OpSum representing the operation +``x ↦ -i[A_1 A_2 …, x]`` where ``A_i`` is an operator which consists of `ai` at site `ni` +and the identity elsewhere. The string `a1` must be an existing ITensor OpName whose +variants `a1⋅` and `⋅a1` are defined (this function, however, doesn't perform any checks). + +# Examples +```julia-repl +julia> gkslcommutator(("σ+", 1), ("σ-", 2)) +sum( + 0.0 - 1.0im σ+⋅(1,) σ-⋅(2,) + 0.0 + 1.0im ⋅σ+(1,) ⋅σ-(2,) +) +``` +""" +function gkslcommutator(operators::Tuple{String,Int}...) + l = OpSum() + opnames = first.(operators) + sites = last.(operators) + l += (-im, collect(Iterators.flatten(zip(opnames .* "⋅", sites)))...) + l += (+im, collect(Iterators.flatten(zip("⋅" .* opnames, sites)))...) + return l +end + +""" + gkslcommutator(args...) + +Given a sequence `a1, n1, a2, n2, …` where each `ai` is a `String` and each `ni` is an +`Int`, return an OpSum representing the operation ``x ↦ -i[A_1 A_2 …, x]`` where ``A_i`` +is an operator which consists of `ai` at site `ni` and the identity elsewhere. The string +`a1` must be an existing ITensor OpName whose variants `a1⋅` and `⋅a1` are defined +(this function, however, doesn't perform any checks). + +# Examples +```julia-repl +julia> gkslcommutator("σ+", 1, "σ-", 2) +sum( + 0.0 - 1.0im σ+⋅(1,) σ-⋅(2,) + 0.0 + 1.0im ⋅σ+(1,) ⋅σ-(2,) +) +``` +""" +gkslcommutator(args...) = gkslcommutator(makeopsumpairs(args...)...) + +""" + gkslcommutator_itensor(sites::Vector{<:Index}, operators::Tuple{String,Int}...) + +Given a vector of ITensors site indices `sites` and a sequence `a1, n1, a2, n2, …` where +each `ai` is a `String` and each `ni` is an `Int`, return an ITensor with the site indices +`s[n1]`, `s[n1]'`, `s[n2]`, `s[n2]'` and so on which represents the operation +``x ↦ -i[A_1 A_2 …, x]`` where ``A_i`` is an operator consisting of `ai` at site +`ni` and the identity elsewhere. Each `ai` string must be an existing ITensors OpName whose +variants `a1⋅` and `⋅a1` are defined (this function, however, doesn't perform any checks). + + +# Examples + +We start from a system of three 1/2-spins: +```julia +sites = siteinds("vS=1/2", 3) +``` + +Take an operator ``U`` which is ``Sˣ`` on the second site and the identity on the +others. The commutator ``ρ ↦ -i[U,ρ]`` is given by +```julia +gkslcommutator_itensor(sites, "Sx", 2) +``` + +An operator ``V`` which is ``Sy`` on the first site, the identity on the second one, and +``Sz`` on the third one. The commutator ``ρ ↦ -i[V,ρ]`` is given by +```julia +gkslcommutator_itensor(sites, "Sy", 1, "Sz", 3) +``` +""" +function gkslcommutator_itensor(sites::Vector{<:Index}, operators::Tuple{String,Int}...) + operator_names = first.(operators) + site_numbers = last.(operators) + if !allunique(sites) + error("Some sites are repeated in the list. Please use unique site indices.") + # This possibility is not allowed for now, since it would create issues in the + # multiplication loop below. Basically, if two ITensors with the same indices + # are multiplied together (with a simple `*`) then both indices get contracted + # and we end up with a scalar. We should use `apply` instead, but it does not + # work with OneITensor objects... + end + lmult = ITensors.OneITensor() + rmult = ITensors.OneITensor() + for (on, j) in zip(operator_names, site_numbers) + lmult *= op("$on⋅", sites, j) + rmult *= op("⋅$on", sites, j) + end + return -im * (lmult - rmult) +end + +function gkslcommutator_itensor(sites::Vector{<:Index}, args...) + return gkslcommutator_itensor(sites, makeopsumpairs(args...)...) +end diff --git a/src/site_types/electron.jl b/src/site_types/electron.jl new file mode 100644 index 0000000..807ea97 --- /dev/null +++ b/src/site_types/electron.jl @@ -0,0 +1,41 @@ +# Additional ITensor states and operators for the electron SiteType. +function ITensors.op(::OpName"Id", ::SiteType"Electron") # ITensors doesn't define this + return Matrix(1.0I, 4, 4) +end + +ITensors.op(::OpName"Ntot^2", st::SiteType"Electron") = ITensors.op(OpName("Ntot"), st)^2 +ITensors.op(::OpName"ntot^2", st::SiteType"Electron") = ITensors.op(OpName("Ntot^2"), st) + +# These aliases aren't defined in ITensors, but we use them a lot in our scripts, +# so we keep them for compatibility reasons. +ITensors.op(::OpName"NupNdn", st::SiteType"Electron") = ITensors.op(OpName("Nupdn"), st) +ITensors.op(::OpName"Aup†", st::SiteType"Electron") = ITensors.op(OpName("Adagup"), st) +ITensors.op(::OpName"Adn†", st::SiteType"Electron") = ITensors.op(OpName("Adagdn"), st) + +# Here are some compositions of operators that are commonly used in defining exchange +# interactions with Electron site types: we need them, actually, for `vElectron` types, +# where multiplication between operators is trickier to use. +function ITensors.op(::OpName"Aup†F", st::SiteType"Electron") + return ITensors.op(OpName("Adagup"), st) * ITensors.op(OpName("F"), st) +end +ITensors.op(::OpName"a†↑F", st::SiteType"Electron") = ITensors.op(OpName("Aup†F"), st) + +function ITensors.op(::OpName"AupF", st::SiteType"Electron") + return ITensors.op(OpName("Aup"), st) * ITensors.op(OpName("F"), st) +end +ITensors.op(::OpName"a↑F", st::SiteType"Electron") = ITensors.op(OpName("AupF"), st) + +function ITensors.op(::OpName"FAdn", st::SiteType"Electron") + return ITensors.op(OpName("F"), st) * ITensors.op(OpName("Adn"), st) +end +ITensors.op(::OpName"Fa↓", st::SiteType"Electron") = ITensors.op(OpName("FAdn"), st) + +function ITensors.op(::OpName"FAdn†", st::SiteType"Electron") + return ITensors.op(OpName("F"), st) * ITensors.op(OpName("Adagdn"), st) +end +ITensors.op(::OpName"Fa†↓", st::SiteType"Electron") = ITensors.op(OpName("FAdn†"), st) + +function ITensors.op(::OpName"Adn†F", st::SiteType"Electron") + return ITensors.op(OpName("Adagdn"), st) * ITensors.op(OpName("F"), st) +end +ITensors.op(::OpName"a†↓F", st::SiteType"Electron") = ITensors.op(OpName("Adn†F"), st) diff --git a/src/site_types/fermion.jl b/src/site_types/fermion.jl new file mode 100644 index 0000000..3a457fa --- /dev/null +++ b/src/site_types/fermion.jl @@ -0,0 +1,23 @@ +function ITensors.op(::OpName"Id", ::SiteType"Fermion") + return Matrix(1.0I, 2, 2) +end + +function ITensors.op!(Op::ITensor, ::OpName"A", ::SiteType"Fermion", s::Index) + 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 +end +function ITensors.op!(Op::ITensor, on::OpName"adag", st::SiteType"Fermion", s::Index) + return ITensors.op!(Op, OpName("Adag"), st, s) +end +function ITensors.op!(Op::ITensor, on::OpName"A†", st::SiteType"Fermion", s::Index) + return ITensors.op!(Op, OpName("Adag"), st, s) +end +function ITensors.op!(Op::ITensor, on::OpName"a†", st::SiteType"Fermion", s::Index) + return ITensors.op!(Op, OpName("Adag"), st, s) +end diff --git a/src/site_types/fermionic_dot3.jl b/src/site_types/fermionic_dot3.jl new file mode 100644 index 0000000..9105d5b --- /dev/null +++ b/src/site_types/fermionic_dot3.jl @@ -0,0 +1,152 @@ +""" + space(::SiteType"FDot3") + +Create the Hilbert space for a site of type "FDot3". + +No conserved symmetries and quantum number labels are provided for this space. +""" +function ITensors.space(::SiteType"FDot3") + return 2^3 +end + +occ = [ + 0 + 1 +] +emp = [ + 1 + 0 +] + +up = [ + 0 0 + 1 0 +] +dn = [ + 0 1 + 0 0 +] +id = [ + 1 0 + 0 1 +] +fs = [ + 1 0 + 0 -1 +] # fermion string operator + +# basis order: +# e_1 -> |∅⟩ +# e_2 -> c₁†|∅⟩ +# e_3 -> c₂†|∅⟩ +# e_4 -> c₁† c₂†|∅⟩ +# e_5 -> c₃†|∅⟩ +# e_6 -> c₁† c₃†|∅⟩ +# e_7 -> c₂† c₃†|∅⟩ +# e_8 -> c₁† c₂† c₃†|∅⟩ +ITensors.state(::StateName"Emp", ::SiteType"FDot3") = kron(emp, emp, emp) +ITensors.state(::StateName"0", st::SiteType"FDot3") = state(StateName("Emp"), st) +ITensors.state(::StateName"Vac", st::SiteType"FDot3") = state(StateName("Emp"), st) +ITensors.state(::StateName"Vacuum", st::SiteType"FDot3") = state(StateName("Emp"), st) + +# Modes: +# (3) (2) (1) +ITensors.state(::StateName"1", ::SiteType"FDot3") = kron(emp, emp, occ) +ITensors.state(::StateName"2", ::SiteType"FDot3") = kron(emp, occ, emp) +ITensors.state(::StateName"12", ::SiteType"FDot3") = kron(emp, occ, occ) +ITensors.state(::StateName"3", ::SiteType"FDot3") = kron(emp, emp, occ) +ITensors.state(::StateName"13", ::SiteType"FDot3") = kron(occ, emp, occ) +ITensors.state(::StateName"23", ::SiteType"FDot3") = kron(occ, occ, emp) +ITensors.state(::StateName"123", ::SiteType"FDot3") = kron(occ, occ, occ) + +function ITensors.op(::OpName"Id", ::SiteType"FDot3") + return Matrix(1.0I, 2^3, 2^3) +end + +function ITensors.op(::OpName"c1", ::SiteType"FDot3") + return kron(id, id, dn) +end +function ITensors.op(::OpName"c2", ::SiteType"FDot3") + return kron(id, dn, fs) +end +function ITensors.op(::OpName"c3", ::SiteType"FDot3") + return kron(dn, fs, fs) +end + +function ITensors.op(::OpName"c†1", ::SiteType"FDot3") + return kron(id, id, up) +end +function ITensors.op(::OpName"c†2", ::SiteType"FDot3") + return kron(id, up, fs) +end +function ITensors.op(::OpName"c†3", ::SiteType"FDot3") + return kron(up, fs, fs) +end + +function ITensors.op(::OpName"n1", st::SiteType"FDot3") + return ITensors.op(OpName("c†1"), st) * ITensors.op(OpName("c1"), st) +end +function ITensors.op(::OpName"n2", st::SiteType"FDot3") + return ITensors.op(OpName("c†2"), st) * ITensors.op(OpName("c2"), st) +end +function ITensors.op(::OpName"n3", st::SiteType"FDot3") + return ITensors.op(OpName("c†3"), st) * ITensors.op(OpName("c3"), st) +end +function ITensors.op(::OpName"ntot", st::SiteType"FDot3") + return sum([ITensors.op(OpName("n$k"), st) for k in 1:3]) +end +function ITensors.op(::OpName"ntot^2", st::SiteType"FDot3") + return (ITensors.op(OpName("ntot"), st))^2 +end + +function ITensors.op(::OpName"F1", st::SiteType"FDot3") + return kron(id, id, fs) +end +function ITensors.op(::OpName"F2", st::SiteType"FDot3") + return kron(id, fs, id) +end +function ITensors.op(::OpName"F3", st::SiteType"FDot3") + return kron(fs, id, id) +end +function ITensors.op(::OpName"F", st::SiteType"FDot3") + return kron(fs, fs, fs) +end + +function dot_hamiltonian(::SiteType"FDot3", energies, coulomb_repulsion, sitenumber::Int) + E = OpSum() + for k in 1:3 + E += (energies[k], "n$k", sitenumber) + end + + N = OpSum() + N += ("ntot", sitenumber) + + return E + 0.5coulomb_repulsion * (Ops.expand(N^2) - N) +end + +function exchange_interaction(::SiteType"FDot3", s1::Index, s2::Index) + stypes1 = sitetypes(s1) + stypes2 = sitetypes(s2) + if (SiteType("FDot3") in stypes1) && !(SiteType("FDot3") in stypes2) + return exchange_interaction(st, sitenumber(s1), sitenumber(s2)) + elseif (SiteType("FDot3") in stypes2) && !(SiteType("FDot3") in stypes1) + return exchange_interaction(st, sitenumber(s2), sitenumber(s1)) + else + # Return an error if no implementation is found for any type. + return throw( + ArgumentError("No FDot3 site type found in either $(tags(s1)) or $(tags(s2)).") + ) + end +end + +function exchange_interaction(::SiteType"FDot3", dot_site::Int, other_site::Int) + h = OpSum() + + for k in 1:3 + h += "c†$k", dot_site, "c", other_site + h += "c†", other_site, "c$k", dot_site + end + # Doesn't work... we can't mix "c" operators from the Fermion SiteType with other + # SiteType operators. + return h +end diff --git a/src/site_types/oscillator.jl b/src/site_types/oscillator.jl new file mode 100644 index 0000000..45dbb57 --- /dev/null +++ b/src/site_types/oscillator.jl @@ -0,0 +1,68 @@ +# Space of harmonic oscillators +# ============================= + +ITensors.alias(::SiteType"Osc") = SiteType"Qudit"() + +""" + ITensors.space(st::SiteType"Osc"; + dim = 2, + conserve_qns = false, + conserve_number = false, + qnname_number = "Number") + +Create the Hilbert space for a site of type "Osc". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +ITensors.space(st::SiteType"Osc"; kwargs...) = space(ITensors.alias(st); kwargs...) + +# Forward "Osc" definitions to "Qudit" states. +# This way each state, operator and so on defined for the "Qudit" SiteType is already +# available for the "Osc" type. +ITensors.val(vn::ValName; st::SiteType"Osc") = val(vn, ITensors.alias(st)) + +function ITensors.state(sn::StateName, st::SiteType"Osc", s::Index; kwargs...) + return state(sn, ITensors.alias(st), s; kwargs...) +end + +function ITensors.op(on::OpName, st::SiteType"Osc", dims::Int...; kwargs...) + return op(on, ITensors.alias(st), dims...; kwargs...) +end + +# ITensors.op functions all require an additional parameter, the dimension: for this reason, +# we need to compute ITensors.dim.(rs), i.e. the dimensions of the Indices of the operator, +# and append them to the function arguments. +function ITensors.op(on::OpName, st::SiteType"Osc", s1::Index, s_tail::Index...; kwargs...) + rs = reverse((s1, s_tail...)) + dims = ITensors.dim.(rs) + opmat = ITensors.op(on, st, dims...; kwargs...) + return ITensors.itensor(opmat, prime.(rs)..., dag.(rs)...) +end + +# Operators +# --------- +function ITensors.op(::OpName"Asum", st::SiteType"Osc", d::Int) + return ITensors.op(OpName("A"), st, d) + ITensors.op(OpName("Adag"), st, d) +end +function ITensors.op(on::OpName"X", st::SiteType"Osc", d::Int) + return ITensors.op(OpName("Asum"), st, d) / sqrt(2) +end +function ITensors.op(on::OpName"Y", st::SiteType"Osc", d::Int) + return im / sqrt(2) * + (ITensors.op(OpName("Adag"), st, d) - ITensors.op(OpName("A"), st, d)) +end + +# Aliases (for backwards compatibility) +ITensors.alias(::OpName"a-") = OpName"A"() +ITensors.alias(::OpName"a+") = OpName"Adag"() +ITensors.alias(::OpName"asum") = OpName"Asum"() + +function ITensors.op(on::OpName"a+", st::SiteType"Osc", d::Int) + return ITensors.op(ITensors.alias(on), st, d) +end +function ITensors.op(on::OpName"a-", st::SiteType"Osc", d::Int) + return ITensors.op(ITensors.alias(on), st, d) +end +function ITensors.op(on::OpName"asum", st::SiteType"Osc", d::Int) + return ITensors.op(ITensors.alias(on), st, d) +end diff --git a/src/site_types/spinhalf.jl b/src/site_types/spinhalf.jl new file mode 100644 index 0000000..9dc9792 --- /dev/null +++ b/src/site_types/spinhalf.jl @@ -0,0 +1,17 @@ +# Additional ITensor states and operators for the spin-half SiteType. +# +# Note that while the following definitions use "S=1/2" to denote the site type, they hold +# for all other aliases of the spin-half type, such as "S=½", automatically. + +# Number operator, aka |↑⟩⟨↑| +ITensors.op(::OpName"N", st::SiteType"S=1/2") = ITensors.op(OpName("ProjUp"), st) +# Identity matrices +ITensors.op(::OpName"Id", ::SiteType"S=1/2") = Matrix(1.0I, 2, 2) +# Jordan-Wigner string operator (aka -σᶻ) +ITensors.op(::OpName"F", st::SiteType"S=1/2") = -ITensors.op(OpName("σz"), st) + +# Ladder operators (aliases) +ITensors.op(::OpName"σ+", st::SiteType"S=1/2") = op(OpName("S+"), st) +ITensors.op(::OpName"σ-", st::SiteType"S=1/2") = op(OpName("S-"), st) +ITensors.op(::OpName"plus", st::SiteType"S=1/2") = op(OpName("S+"), st) +ITensors.op(::OpName"minus", st::SiteType"S=1/2") = op(OpName("S-"), st) diff --git a/src/site_types/vectorized_electron.jl b/src/site_types/vectorized_electron.jl new file mode 100644 index 0000000..f46fd53 --- /dev/null +++ b/src/site_types/vectorized_electron.jl @@ -0,0 +1,96 @@ +# Space of electrons (vectorised) +# ======================================== +""" + ITensors.space(st::SiteType"vElectron") + +Create the Hilbert space for a site of type "vElectron", i.e. a mixed state describing a +site with a 1/2-spin degree of freedom. +The density matrix is represented in the generalised Gell-Mann basis of `Mat(ℂ⁴)`, composed +of Hermitian traceless matrices together with the identity matrix. +""" +ITensors.space(::SiteType"vElectron") = 16 + +# Shorthand notation: +function vstate(sn::StateName, ::SiteType"vElectron") + v = ITensors.state(sn, SiteType("Electron")) + return vec(kron(v, v'), gellmannbasis(4)) +end +function vop(sn::StateName, ::SiteType"vElectron") + sn = statenamestring(sn) + on = sn[1] == 'v' ? sn[2:end] : sn + return vec(ITensors.op(OpName(on), SiteType("Electron")), gellmannbasis(4)) +end + +# States (actual ones) +# -------------------- +ITensors.state(sn::StateName"Emp", st::SiteType"vElectron") = vstate(sn, st) +ITensors.state(sn::StateName"Up", st::SiteType"vElectron") = vstate(sn, st) +ITensors.state(sn::StateName"Dn", st::SiteType"vElectron") = vstate(sn, st) +ITensors.state(sn::StateName"UpDn", st::SiteType"vElectron") = vstate(sn, st) + +# States representing vectorised operators +# ---------------------------------------- +ITensors.state(sn::StateName"vId", st::SiteType"vElectron") = vop(sn, st) +ITensors.state(sn::StateName"vNup", st::SiteType"vElectron") = vop(sn, st) +ITensors.state(sn::StateName"vNdn", st::SiteType"vElectron") = vop(sn, st) +ITensors.state(sn::StateName"vNtot", st::SiteType"vElectron") = vop(sn, st) +ITensors.state(sn::StateName"vNupNdn", st::SiteType"vElectron") = vop(sn, st) + +function ITensors.state(::StateName"vecId", ::SiteType"vElectron") + return ITensors.state(StateName("vId"), SiteType("vElectron")) +end + +# Operators acting on vectorised spins +# ------------------------------------ +function premultiply(mat, ::SiteType"vElectron") + return PseudomodesTTEDOPA.vec(x -> mat * x, gellmannbasis(4)) +end +function postmultiply(mat, ::SiteType"vElectron") + return PseudomodesTTEDOPA.vec(x -> x * mat, gellmannbasis(4)) +end + +# The goal here is to define operators "A⋅" and "⋅A" in an automatic way whenever the +# OpName "A" is defined for the Electron site type. +# This is handy, but unless we find a better way to define this function this means that +# _every_ operator has to be written this way; we cannot just return op(on, st) at the end +# if no "⋅" is found, otherwise an infinite loop would be entered. +# We make an exception, though, for "Id" since it is an essential operator, and something +# would probably break if it weren't defined. +function ITensors.op(on::OpName, st::SiteType"vElectron"; kwargs...) + name = strip(String(ITensors.name(on))) # Remove extra whitespace + if name == "Id" + return Matrix(1.0I, 16, 16) + end + dotloc = findfirst("⋅", name) + # This returns the position of the cdot in the operator name String. + # It is `nothing` if no cdot is found. + if !isnothing(dotloc) + on1, on2 = nothing, nothing + on1 = name[1:prevind(name, dotloc.start)] + on2 = name[nextind(name, dotloc.start):end] + # If the OpName `on` is written correctly, i.e. it is either "A⋅" or "⋅A" for some + # A, then either `on1` or `on2` has to be empty (not both, not neither of them). + if (on1 == "" && on2 == "") || (on1 != "" && on2 != "") + throw( + ArgumentError( + "Invalid operator name: $name. Operator name is not \"Id\" or of the " * + "form \"A⋅\" or \"⋅A\"", + ), + ) + end + # name == "⋅A" -> on1 is an empty string + # name == "A⋅" -> on2 is an empty string + if on1 == "" + mat = try_op(OpName(on2), SiteType("Electron"); kwargs...) + return postmultiply(mat, st) + elseif on2 == "" + mat = try_op(OpName(on1), SiteType("Electron"); kwargs...) + return premultiply(mat, st) + else + # This should logically never happen but, just in case, we throw an error. + error("Unknown error with operator name $name") + end + else + error("Operator name $name is not \"Id\" or of the form \"A⋅\" or \"⋅A\"") + end +end diff --git a/src/site_types/vectorized_fermion.jl b/src/site_types/vectorized_fermion.jl new file mode 100644 index 0000000..7b73e7b --- /dev/null +++ b/src/site_types/vectorized_fermion.jl @@ -0,0 +1,113 @@ +# Space of spin-1/2 particles (vectorised) +# ======================================== +""" + ITensors.space(st::SiteType"vFermion") + +Create the Hilbert space for a site of type "vFermion", i.e. a vectorised +spin-1/2 particle, where the vectorisation is performed wrt the generalised +Gell-Mann basis of `Mat(ℂ²)`, composed of Hermitian traceless matrices +together with the identity matrix. +""" +ITensors.space(::SiteType"vFermion") = 4 + +# Elements of and operators on Mat(ℂ²) are expanded wrt the basis {Λᵢ}ᵢ₌₁⁴ of +# generalised Gell-Mann matrices (plus a multiple of the identity). +# An element A ∈ Mat(ℂ²) is representeb by the a vector v such that +# vᵢ = tr(Λᵢ A), +# while a linear map L : Mat(ℂ²) → Mat(ℂ²) by the matrix ℓ such that +# ℓᵢⱼ = tr(Λᵢ L(Λⱼ)). + +# Shorthand notation: +function vstate(sn::StateName, ::SiteType"vFermion") + v = ITensors.state(sn, SiteType("Fermion")) + return PseudomodesTTEDOPA.vec(kron(v, v'), gellmannbasis(2)) +end +function vop(sn::StateName, ::SiteType"vFermion") + sn = statenamestring(sn) + on = sn[1] == 'v' ? sn[2:end] : sn + return PseudomodesTTEDOPA.vec(try_op(OpName(on), SiteType("Fermion")), gellmannbasis(2)) +end + +# States (actual ones) +# -------------------- +ITensors.state(sn::StateName"Emp", st::SiteType"vFermion") = vstate(sn, st) +ITensors.state(sn::StateName"Occ", st::SiteType"vFermion") = vstate(sn, st) + +function ITensors.state(::StateName"Up", st::SiteType"vFermion") + return ITensors.state(StateName("Occ"), st) +end +function ITensors.state(::StateName"Dn", st::SiteType"vFermion") + return ITensors.state(StateName("Emp"), st) +end + +# States representing vectorised operators +# ---------------------------------------- +ITensors.state(sn::StateName"vId", st::SiteType"vFermion") = vop(sn, st) +ITensors.state(sn::StateName"vN", st::SiteType"vFermion") = vop(sn, st) + +ITensors.state(sn::StateName"vA", st::SiteType"vFermion") = vop(sn, st) +ITensors.state(sn::StateName"va", st::SiteType"vFermion") = vop(sn, st) + +ITensors.state(sn::StateName"vAdag", st::SiteType"vFermion") = vop(sn, st) +ITensors.state(sn::StateName"vadag", st::SiteType"vFermion") = vop(sn, st) +ITensors.state(sn::StateName"vA†", st::SiteType"vFermion") = vop(sn, st) +ITensors.state(sn::StateName"va†", st::SiteType"vFermion") = vop(sn, st) + +function ITensors.state(::StateName"vecId", ::SiteType"vFermion") + return ITensors.state(StateName("vId"), SiteType("vFermion")) +end + +# Operator dispatch +# ================= +function premultiply(mat, ::SiteType"vFermion") + return PseudomodesTTEDOPA.vec(x -> mat * x, gellmannbasis(2)) +end +function postmultiply(mat, ::SiteType"vFermion") + return PseudomodesTTEDOPA.vec(x -> x * mat, gellmannbasis(2)) +end + +# The goal here is to define operators "A⋅" and "⋅A" in an automatic way whenever the +# OpName "A" is defined for the Fermion site type. +# This is handy, but unless we find a better way to define this function this means that +# _every_ operator has to be written this way; we cannot just return op(on, st) at the end +# if no "⋅" is found, otherwise an infinite loop would be entered. +# We make an exception, though, for "Id" since it is an essential operator, and something +# would probably break if it weren't defined. +function ITensors.op(on::OpName, st::SiteType"vFermion"; kwargs...) + name = strip(String(ITensors.name(on))) # Remove extra whitespace + if name == "Id" + return Matrix(1.0I, 4, 4) + end + dotloc = findfirst("⋅", name) + # This returns the position of the cdot in the operator name String. + # It is `nothing` if no cdot is found. + if !isnothing(dotloc) + on1, on2 = nothing, nothing + on1 = name[1:prevind(name, dotloc.start)] + on2 = name[nextind(name, dotloc.start):end] + # If the OpName `on` is written correctly, i.e. it is either "A⋅" or "⋅A" for some + # A, then either `on1` or `on2` has to be empty (not both, not neither of them). + if (on1 == "" && on2 == "") || (on1 != "" && on2 != "") + throw( + ArgumentError( + "Invalid operator name: $name. Operator name is not \"Id\" or of the " * + "form \"A⋅\" or \"⋅A\"", + ), + ) + end + # name == "⋅A" -> on1 is an empty string + # name == "A⋅" -> on2 is an empty string + if on1 == "" + mat = try_op(OpName(on2), SiteType("Fermion"); kwargs...) + return postmultiply(mat, st) + elseif on2 == "" + mat = try_op(OpName(on1), SiteType("Fermion"); kwargs...) + return premultiply(mat, st) + else + # This should logically never happen but, just in case, we throw an error. + error("Unknown error with operator name $name") + end + else + error("Operator name $name is not \"Id\" or of the form \"A⋅\" or \"⋅A\"") + end +end diff --git a/src/site_types/vectorized_fermionic_dot3.jl b/src/site_types/vectorized_fermionic_dot3.jl new file mode 100644 index 0000000..90f58b2 --- /dev/null +++ b/src/site_types/vectorized_fermionic_dot3.jl @@ -0,0 +1,147 @@ +""" + space(::SiteType"vFDot3") + +Create the Hilbert space for a site of type "vFDot3", i.e. a mixed state describing a +site with a fermionic three-level quantum dot. +The density matrix is represented in the generalised Gell-Mann basis, composed +of Hermitian traceless matrices together with the identity matrix. + +No conserved symmetries and quantum number labels are provided for this space. +""" +function ITensors.space(::SiteType"vFDot3") + return (2^3)^2 +end + +# An element A ∈ Mat(ℂ⁴) is representeb by the a vector v such that +# vᵢ = tr(Λᵢ A), +# while a linear map L : Mat(ℂ⁴) → Mat(ℂ⁴) by the matrix ℓ such that +# ℓᵢⱼ = tr(Λᵢ L(Λⱼ)). + +# Shorthand notation: +function vstate(sn::AbstractString, ::SiteType"vFDot3") + v = ITensors.state(StateName(sn), SiteType("FDot3")) + return PseudomodesTTEDOPA.vec(kron(v, v'), gellmannbasis(2^3)) +end +function vop(on::AbstractString, ::SiteType"vFDot3") + return PseudomodesTTEDOPA.vec( + ITensors.op(OpName(on), SiteType("FDot3")), gellmannbasis(2^3) + ) +end + +# basis order: +# e_1 -> |∅⟩ +# e_2 -> c₁†|∅⟩ +# e_3 -> c₂†|∅⟩ +# e_4 -> c₁† c₂†|∅⟩ +# e_5 -> c₃†|∅⟩ +# e_6 -> c₁† c₃†|∅⟩ +# e_7 -> c₂† c₃†|∅⟩ +# e_8 -> c₁† c₂† c₃†|∅⟩ +ITensors.state(::StateName"Emp", st::SiteType"vFDot3") = vstate("Emp", st) +ITensors.state(::StateName"0", st::SiteType"vFDot3") = state(StateName("Emp"), st) +ITensors.state(::StateName"Vac", st::SiteType"vFDot3") = state(StateName("Emp"), st) +ITensors.state(::StateName"Vacuum", st::SiteType"vFDot3") = state(StateName("Emp"), st) + +ITensors.state(::StateName"1", st::SiteType"vFDot3") = vstate("1", st) +ITensors.state(::StateName"2", st::SiteType"vFDot3") = vstate("2", st) +ITensors.state(::StateName"12", st::SiteType"vFDot3") = vstate("12", st) +ITensors.state(::StateName"3", st::SiteType"vFDot3") = vstate("3", st) +ITensors.state(::StateName"13", st::SiteType"vFDot3") = vstate("13", st) +ITensors.state(::StateName"23", st::SiteType"vFDot3") = vstate("23", st) +ITensors.state(::StateName"123", st::SiteType"vFDot3") = vstate("123", st) + +ITensors.state(::StateName"vId", st::SiteType"vFDot3") = vop("Id", st) +ITensors.state(::StateName"vecId", st::SiteType"vFDot3") = vop("Id", st) +ITensors.state(::StateName"vn1", st::SiteType"vFDot3") = vop("n1", st) +ITensors.state(::StateName"vn2", st::SiteType"vFDot3") = vop("n2", st) +ITensors.state(::StateName"vn3", st::SiteType"vFDot3") = vop("n3", st) +ITensors.state(::StateName"vntot", st::SiteType"vFDot3") = vop("ntot", st) + +# Operator dispatch +# ================= +function premultiply(mat, ::SiteType"vFDot3") + return PseudomodesTTEDOPA.vec(x -> mat * x, gellmannbasis(2^3)) +end +function postmultiply(mat, ::SiteType"vFDot3") + return PseudomodesTTEDOPA.vec(x -> x * mat, gellmannbasis(2^3)) +end + +# The goal here is to define operators "A⋅" and "⋅A" in an automatic way whenever the +# OpName "A" is defined for the FDot3 site type. +# This is handy, but unless we find a better way to define this function this means that +# _every_ operator has to be written this way; we cannot just return op(on, st) at the end +# if no "⋅" is found, otherwise an infinite loop would be entered. +# We make an exception, though, for "Id" since it is an essential operator, and something +# would probably break if it weren't defined. +function ITensors.op(on::OpName, st::SiteType"vFDot3"; kwargs...) + name = strip(String(ITensors.name(on))) # Remove extra whitespace + if name == "Id" + return Matrix(1.0I, 2^6, 2^6) + end + dotloc = findfirst("⋅", name) + # This returns the position of the cdot in the operator name String. + # It is `nothing` if no cdot is found. + if !isnothing(dotloc) + on1, on2 = nothing, nothing + on1 = name[1:prevind(name, dotloc.start)] + on2 = name[nextind(name, dotloc.start):end] + # If the OpName `on` is written correctly, i.e. it is either "A⋅" or "⋅A" for some + # A, then either `on1` or `on2` has to be empty (not both, not neither of them). + if (on1 == "" && on2 == "") || (on1 != "" && on2 != "") + throw( + ArgumentError( + "Invalid operator name: $name. Operator name is not \"Id\" or of the " * + "form \"A⋅\" or \"⋅A\"", + ), + ) + end + # name == "⋅A" -> on1 is an empty string + # name == "A⋅" -> on2 is an empty string + if on1 == "" + mat = try_op(OpName(on2), SiteType("FDot3"); kwargs...) + return postmultiply(mat, st) + elseif on2 == "" + mat = try_op(OpName(on1), SiteType("FDot3"); kwargs...) + return premultiply(mat, st) + else + # This should logically never happen but, just in case, we throw an error. + error("Unknown error with operator name $name") + end + else + error("Operator name $name is not \"Id\" or of the form \"A⋅\" or \"⋅A\"") + end +end + +function dot_hamiltonian(::SiteType"vFDot3", energies, coulomb_repulsion, sitenumber::Int) + E = OpSum() + for k in 1:3 + E += energies[k] * gkslcommutator("n$k", sitenumber) + end + + N² = gkslcommutator("ntot^2", sitenumber) + N = gkslcommutator("ntot", sitenumber) + + return E + 0.5coulomb_repulsion * (N² - N) +end + +function exchange_interaction(::SiteType"vFDot3", s1::Index, s2::Index) + stypes1 = sitetypes(s1) + stypes2 = sitetypes(s2) + if (SiteType("vFDot3") in stypes1) && !(SiteType("vFDot3") in stypes2) + return exchange_interaction(st, sitenumber(s1), sitenumber(s2)) + elseif (SiteType("vFDot3") in stypes2) && !(SiteType("vFDot3") in stypes1) + return exchange_interaction(st, sitenumber(s2), sitenumber(s1)) + else + # Return an error if no implementation is found for any type. + throw( + ArgumentError("No vFDot3 site type found in either $(tags(s1)) or $(tags(s2)).") + ) + end +end + +function exchange_interaction(::SiteType"vFDot3", dot_site::Int, other_site::Int) + return sum([ + gkslcommutator("c†$k", dot_site, "σ-", other_site) + + gkslcommutator("c$k", dot_site, "σ+", other_site) for k in 1:3 + ]) +end diff --git a/src/site_types/vectorized_oscillator.jl b/src/site_types/vectorized_oscillator.jl new file mode 100644 index 0000000..662b893 --- /dev/null +++ b/src/site_types/vectorized_oscillator.jl @@ -0,0 +1,294 @@ +# Space of harmonic oscillators (vectorised) +# ========================================== +""" + ITensors.space(st::SiteType"vOsc"; dim = 2) + +Create the Hilbert space for a site of type "vOsc", i.e. a vectorised +harmonic oscillator of dimension `dim` (this means that the space has dimension +`dim^2`), where the vectorisation is performed wrt the (Hermitian) Gell-Mann +basis of `Mat(ℂᵈⁱᵐ)`. +""" +function ITensors.space(::SiteType"vOsc"; dim=2) + return dim^2 +end + +# Just as with the "Osc" site type, "vOsc" operator and states require that we specify the +# dimension of the space, so we need to compute ITensors.dim.(rs), i.e. the dimensions of +# the Indices of the operator, and append them to the function arguments. + +function ITensors.state(sn::StateName, st::SiteType"vOsc", s::Index; kwargs...) + d = isqrt(ITensors.dim(s)) + stvec = ITensors.state(sn, st, d; kwargs...) + return ITensors.itensor(stvec, s) +end + +function ITensors.op(on::OpName, st::SiteType"vOsc", s1::Index, s_tail::Index...; kwargs...) + rs = reverse((s1, s_tail...)) + dims = isqrt.(ITensors.dim.(rs)) + # Use `isqrt` which takes an Int and returns an Int; the decimal part is + # truncated, but it's not a problem since by construction dim.(rs) are + # perfect squares, so the decimal part should be zero. + opmat = ITensors.op(on, st, dims...; kwargs...) + return ITensors.itensor(opmat, prime.(rs)..., dag.(rs)...) +end + +# Aliases (for backwards compatibility) +ITensors.alias(::SiteType"HvOsc") = SiteType"vOsc"() +ITensors.alias(::SiteType"vecOsc") = SiteType"vOsc"() + +function ITensors.space(st::SiteType"HvOsc"; kwargs...) + return ITensors.space(ITensors.alias(st); kwargs...) +end +function ITensors.space(st::SiteType"vecOsc"; kwargs...) + return ITensors.space(ITensors.alias(st); kwargs...) +end + +ITensors.val(vn::ValName, st::SiteType"HvOsc") = ITensors.val(vn, ITensors.alias(st)) +ITensors.val(vn::ValName, st::SiteType"vecOsc") = ITensors.val(vn, ITensors.alias(st)) + +function ITensors.state(sn::StateName, st::SiteType"vecOsc", s::Index; kwargs...) + return ITensors.state(sn, ITensors.alias(st), s; kwargs...) +end +function ITensors.state(sn::StateName, st::SiteType"HvOsc", s::Index; kwargs...) + return ITensors.state(sn, ITensors.alias(st), s; kwargs...) +end + +function ITensors.op( + on::OpName, st::SiteType"vecOsc", s1::Index, s_tail::Index...; kwargs... +) + return ITensors.op(on, ITensors.alias(st), s1, s_tail...; kwargs...) +end +function ITensors.op( + on::OpName, st::SiteType"HvOsc", s1::Index, s_tail::Index...; kwargs... +) + return ITensors.op(on, ITensors.alias(st), s1, s_tail...; kwargs...) +end + +# Shorthand notation: +function vstate(sn::StateName, ::SiteType"vOsc", d::Int) + v = ITensors.state(sn, SiteType("Osc")) + return PseudomodesTTEDOPA.vec(kron(v, v'), gellmannbasis(d)) +end +function vop(sn::StateName, ::SiteType"vOsc", d::Int) + sn = statenamestring(sn) + on = sn[1] == 'v' ? sn[2:end] : sn + return PseudomodesTTEDOPA.vec(try_op(OpName(on), SiteType("Osc"), d), gellmannbasis(d)) +end + +# States +# ------ +function ITensors.state(::StateName{N}, ::SiteType"vOsc", d::Int) where {N} + # Eigenstates êₙ ⊗ êₙ of the number operator, wrt the Hermitian basis. + n = parse(Int, String(N)) + v = zeros(d) + v[n + 1] = 1.0 + return vec(kron(v, v'), gellmannbasis(d)) +end + +function ITensors.state( + ::StateName"ThermEq", st::SiteType"vOsc", d::Int; frequency::Real, temperature::Real +) + if temperature == 0 + return ITensors.state(StateName("0"), st, d) + else + numop = ITensors.op(OpName("N"), SiteType("Osc"), d) + # We don't need to define our own matrix for the number operator when + # we can call this one instead. + ρ_eq = exp(-frequency / temperature * numop) + ρ_eq /= tr(ρ_eq) + return vec(ρ_eq, gellmannbasis(d)) + end +end + +# Product of X = (a+a†)/√2 and of the thermal equilibrium state Z⁻¹vec(exp(-βH)). +# It is used in the computation of the correlation function of the bath. +function ITensors.state( + ::StateName"X⋅Therm", st::SiteType"vOsc", d::Int; frequency::Real, temperature::Real +) + xop = ITensors.op(OpName("X"), SiteType("Osc"), d) + if temperature == 0 + ρ_eq = zeros(Float64, d, d) + ρ_eq[1, 1] = 1.0 + else + numop = ITensors.op(OpName("N"), SiteType("Osc"), d) + ρ_eq = exp(-frequency / temperature * numop) + ρ_eq /= tr(ρ_eq) + end + return vec(xop * ρ_eq, gellmannbasis(d)) +end + +# States representing vectorised operators +# ---------------------------------------- +ITensors.state(sn::StateName"vAdag", st::SiteType"vOsc", d::Int) = vop(sn, st, d) +ITensors.state(sn::StateName"vA", st::SiteType"vOsc", d::Int) = vop(sn, st, d) +ITensors.state(sn::StateName"vN", st::SiteType"vOsc", d::Int) = vop(sn, st, d) +ITensors.state(sn::StateName"vId", st::SiteType"vOsc", d::Int) = vop(sn, st, d) +ITensors.state(sn::StateName"vX", st::SiteType"vOsc", d::Int) = vop(sn, st, d) +ITensors.state(sn::StateName"vY", st::SiteType"vOsc", d::Int) = vop(sn, st, d) + +# Aliases (for backwards compatibility) +function ITensors.state(::StateName"veca+", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vAdag"), st, d) +end +function ITensors.state(::StateName"veca-", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vA"), st, d) +end +function ITensors.state(::StateName"vecplus", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vAdag"), st, d) +end +function ITensors.state(::StateName"vecminus", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vA"), st, d) +end +function ITensors.state(::StateName"vecN", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vN"), st, d) +end +function ITensors.state(::StateName"vecId", st::SiteType"vOsc", d::Int) + return ITensors.state(StateName("vId"), st, d) +end + +# Operator dispatch +# ================= +function premultiply(mat, ::SiteType"vOsc", d::Int) + return PseudomodesTTEDOPA.vec(x -> mat * x, gellmannbasis(d)) +end +function postmultiply(mat, ::SiteType"vOsc", d::Int) + return PseudomodesTTEDOPA.vec(x -> x * mat, gellmannbasis(d)) +end + +# The goal here is to define operators "A⋅" and "⋅A" in an automatic way whenever the +# OpName "A" is defined for the Osc site type. +# This is handy, but unless we find a better way to define this function this means that +# _every_ operator has to be written this way; we cannot just return op(on, st) at the end +# if no "⋅" is found, otherwise an infinite loop would be entered. +# We make an exception, though, for "Id" since it is an essential operator, and something +# would probably break if it weren't defined. +function ITensors.op(on::OpName, st::SiteType"vOsc", d::Int; kwargs...) + name = strip(String(ITensors.name(on))) # Remove extra whitespace + if name == "Id" + return Matrix(1.0I, d^2, d^2) + end + dotloc = findfirst("⋅", name) + # This returns the position of the cdot in the operator name String. + # It is `nothing` if no cdot is found. + if !isnothing(dotloc) + on1, on2 = nothing, nothing + on1 = name[1:prevind(name, dotloc.start)] + on2 = name[nextind(name, dotloc.start):end] + # If the OpName `on` is written correctly, i.e. it is either "A⋅" or "⋅A" for some + # A, then either `on1` or `on2` has to be empty (not both, not neither of them). + if (on1 == "" && on2 == "") || (on1 != "" && on2 != "") + throw( + ArgumentError( + "Invalid operator name: $name. Operator name is not \"Id\" or of the " * + "form \"A⋅\" or \"⋅A\"", + ), + ) + end + # name == "⋅A" -> on1 is an empty string + # name == "A⋅" -> on2 is an empty string + if on1 == "" + mat = try_op(OpName(on2), SiteType("Osc"), d; kwargs...) + return postmultiply(mat, st, d) + elseif on2 == "" + mat = try_op(OpName(on1), SiteType("Osc"), d; kwargs...) + return premultiply(mat, st, d) + else + # This should logically never happen but, just in case, we throw an error. + error("Unknown error with operator name $name") + end + else + error("Operator name $name is not \"Id\" or of the form \"A⋅\" or \"⋅A\"") + end +end + +# GKSL equation terms +# ------------------- +# For reference: given an Index `i`, +# op(i, "B * A") == replaceprime(op(i', "B") * op(i, "A"), 2, 1) +# so the composition is performed from right to left. +# Look out for the correct order of operations: +# vec(ABx) == op(A⋅) * vec(Bx) == op(A⋅) * (op(B⋅) * vec(x)) +# vec(xAB) == op(⋅B) * vec(xA) == op(⋅B) * (op(⋅A) * vec(x)) + +# Separate absorption and dissipation terms in GKSL equation +""" + dissipator_gain(n::Int) + +Return an OpSum object representing a dissipator operator (as in a GKSL equation) on site +`n`, associated to the jump operator `Adag` (the creation operator). +The OpNames `A` and `Adag` must be defined for the site type in use. +""" +function dissipator_gain(n::Int) + dsp = OpSum() + # a† ρ a - ½ a a† ρ - ½ ρ a a† + dsp += "Adag⋅ * ⋅A", n + dsp += -0.5, "A⋅ * Adag⋅", n + dsp += -0.5, "⋅Adag * ⋅A", n + return dsp +end + +""" + dissipator_loss(n::Int) + +Return an OpSum object representing a dissipator operator (as in a GKSL equation) on site +`n`, associated to the jump operator `A` (the annihilation operator). +The OpNames `A` and `Adag` must be defined for the site type in use. +""" +function dissipator_loss(n::Int) + dsp = OpSum() + # a ρ a† - ½ a† a ρ - ½ ρ a† a + dsp += "A⋅ * ⋅Adag", n + dsp += -0.5, "N⋅", n + dsp += -0.5, "⋅N", n + return dsp +end + +""" + dissipator(n::Int, frequency::Real, temperature::Real) + +Return an OpSum object which represents the dissipator in the GKSL equation for given +`frequency` and `temperature`, i.e. with dissipation coefficient `1/(exp(freq/temp) - 1)`. +""" +function dissipator(n::Int, frequency::Real, temperature::Real) + if temperature == 0 + return dissipator_loss(n) + else + # Use `expm1(x)` instead of `e^x-1` for better precision. + avgn = 1 / expm1(frequency / temperature) + return (avgn + 1) * dissipator_loss(n) + avgn * dissipator_gain(n) + end +end + +""" + mixedlindbladplus(n1::Int, n2::Int) + +Return on OpSum expression which represents a mixed dissipation operator, appearing in the +equation for two pseudomodes, on sites at positions `n1` and `n2` in the system. +""" +function mixedlindbladplus(n1::Int, n2::Int) + x = OpSum() + x += "A⋅", n1, "⋅Adag", n2 + x += "⋅Adag", n1, "A⋅", n2 + x += -0.5, "Adag⋅", n1, "A⋅", n2 + x += -0.5, "A⋅", n1, "Adag⋅", n2 + x += -0.5, "⋅Adag", n1, "⋅A", n2 + x += -0.5, "⋅A", n1, "⋅Adag", n2 + return x +end + +""" + mixedlindbladminus(n1::Int, n2::Int) + +Return on OpSum expression which represents a mixed dissipation operator, appearing in the +equation for two pseudomodes, on sites at positions `n1` and `n2` in the system. +""" +function mixedlindbladminus(n1::Int, n2::Int) + x = OpSum() + x += "Adag⋅", n1, "⋅A", n2 + x += "Adag⋅", n2, "⋅A", n1 + x += -0.5, "A⋅", n1, "Adag⋅", n2 + x += -0.5, "Adag⋅", n1, "A⋅", n2 + x += -0.5, "⋅A", n1, "⋅Adag", n2 + x += -0.5, "⋅Adag", n1, "⋅A", n2 + return x +end diff --git a/src/site_types/vectorized_spinhalf.jl b/src/site_types/vectorized_spinhalf.jl new file mode 100644 index 0000000..8e3ba5f --- /dev/null +++ b/src/site_types/vectorized_spinhalf.jl @@ -0,0 +1,158 @@ +# Space of spin-1/2 particles (vectorised) +# ======================================== +""" + ITensors.space(st::SiteType"vS=1/2"; dim = 2) + +Create the Hilbert space for a site of type "vS=1/2", i.e. a vectorised +spin-1/2 particle, where the vectorisation is performed wrt the generalised +Gell-Mann basis of `Mat(ℂ²)`, composed of Hermitian traceless matrices +together with the identity matrix. +""" +ITensors.space(::SiteType"vS=1/2") = 4 + +# Elements of and operators on Mat(ℂ²) are expanded wrt the basis {Λᵢ}ᵢ₌₁⁴ of +# generalised Gell-Mann matrices (plus a multiple of the identity). +# An element A ∈ Mat(ℂ²) is representeb by the a vector v such that +# vᵢ = tr(Λᵢ A), +# while a linear map L : Mat(ℂ²) → Mat(ℂ²) by the matrix ℓ such that +# ℓᵢⱼ = tr(Λᵢ L(Λⱼ)). + +# Aliases (for backwards compatibility) +ITensors.alias(::SiteType"HvS=1/2") = SiteType"vS=1/2"() +ITensors.alias(::SiteType"vecS=1/2") = SiteType"vS=1/2"() + +ITensors.space(st::SiteType"HvS=1/2") = ITensors.space(ITensors.alias(st)) +ITensors.space(st::SiteType"vecS=1/2") = ITensors.space(ITensors.alias(st)) + +ITensors.val(vn::ValName, st::SiteType"HvS=1/2") = ITensors.val(vn, ITensors.alias(st)) +ITensors.val(vn::ValName, st::SiteType"vecS=1/2") = ITensors.val(vn, ITensors.alias(st)) + +# ! +# We need to replicate the signatures of the functions below: since the states are defined +# using ITensors.state(::StateName, ::SiteType"vS=1/2"), the aliasing function must follow +# the same structure, otherwise it doesn't pick up the definitions for "vS=1/2". +# The same goes for the operators. +function ITensors.state(sn::StateName, st::SiteType"vecS=1/2"; kwargs...) + return ITensors.state(sn, ITensors.alias(st); kwargs...) +end +function ITensors.state(sn::StateName, st::SiteType"HvS=1/2"; kwargs...) + return ITensors.state(sn, ITensors.alias(st); kwargs...) +end + +function ITensors.op(on::OpName, st::SiteType"vecS=1/2"; kwargs...) + return ITensors.op(on, ITensors.alias(st); kwargs...) +end +function ITensors.op(on::OpName, st::SiteType"HvS=1/2"; kwargs...) + return ITensors.op(on, ITensors.alias(st); kwargs...) +end + +# Shorthand notation: +function vstate(sn::StateName, ::SiteType"vS=1/2") + v = ITensors.state(sn, SiteType("S=1/2")) + return PseudomodesTTEDOPA.vec(kron(v, v'), gellmannbasis(2)) +end +function vop(sn::StateName, ::SiteType"vS=1/2") + sn = statenamestring(sn) + on = sn[1] == 'v' ? sn[2:end] : sn + return PseudomodesTTEDOPA.vec(try_op(OpName(on), SiteType("S=1/2")), gellmannbasis(2)) +end + +# States (actual ones) +# -------------------- +ITensors.state(sn::StateName"Up", st::SiteType"vS=1/2") = vstate(sn, st) +ITensors.state(sn::StateName"Dn", st::SiteType"vS=1/2") = vstate(sn, st) + +# States representing vectorised operators +# ---------------------------------------- +ITensors.state(sn::StateName"vSx", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vSy", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vSz", st::SiteType"vS=1/2") = vop(sn, st) + +ITensors.state(sn::StateName"vX", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vY", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vZ", st::SiteType"vS=1/2") = vop(sn, st) + +ITensors.state(sn::StateName"vσx", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vσy", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vσz", st::SiteType"vS=1/2") = vop(sn, st) + +ITensors.state(sn::StateName"vId", st::SiteType"vS=1/2") = vop(sn, st) +ITensors.state(sn::StateName"vN", st::SiteType"vS=1/2") = vop(sn, st) + +# Aliases (for backwards compatibility) +function ITensors.state(::StateName"vecσx", st::SiteType"vS=1/2") + return ITensors.state(StateName("vσx"), st) +end +function ITensors.state(::StateName"vecσy", st::SiteType"vS=1/2") + return ITensors.state(StateName("vσy"), st) +end +function ITensors.state(::StateName"vecσz", st::SiteType"vS=1/2") + return ITensors.state(StateName("vσz"), st) +end +function ITensors.state(::StateName"vecplus", ::SiteType"vS=1/2") + return vec(ITensors.op(OpName("S+"), SiteType("S=1/2")), gellmannbasis(2)) +end +function ITensors.state(::StateName"vecminus", ::SiteType"vS=1/2") + return vec(ITensors.op(OpName("S-"), SiteType("S=1/2")), gellmannbasis(2)) +end +function ITensors.state(::StateName"vecN", st::SiteType"vS=1/2") + return ITensors.state(StateName("vN"), st) +end +function ITensors.state(::StateName"vecId", st::SiteType"vS=1/2") + return ITensors.state(StateName("vId"), st) +end + +# Operator dispatch +# ================= +function premultiply(mat, ::SiteType"vS=1/2") + return PseudomodesTTEDOPA.vec(x -> mat * x, gellmannbasis(2)) +end +function postmultiply(mat, ::SiteType"vS=1/2") + return PseudomodesTTEDOPA.vec(x -> x * mat, gellmannbasis(2)) +end + +# The goal here is to define operators "A⋅" and "⋅A" in an automatic way whenever the +# OpName "A" is defined for the S=1/2 site type. +# This is handy, but unless we find a better way to define this function this means that +# _every_ operator has to be written this way; we cannot just return op(on, st) at the end +# if no "⋅" is found, otherwise an infinite loop would be entered. +# We make an exception, though, for "Id" since it is an essential operator, and something +# would probably break if it weren't defined. +function ITensors.op(on::OpName, st::SiteType"vS=1/2"; kwargs...) + name = strip(String(ITensors.name(on))) # Remove extra whitespace + if name == "Id" + return Matrix(1.0I, 4, 4) + end + dotloc = findfirst("⋅", name) + # This returns the position of the cdot in the operator name String. + # It is `nothing` if no cdot is found. + if !isnothing(dotloc) + on1, on2 = nothing, nothing + on1 = name[1:prevind(name, dotloc.start)] + on2 = name[nextind(name, dotloc.start):end] + # If the OpName `on` is written correctly, i.e. it is either "A⋅" or "⋅A" for some + # A, then either `on1` or `on2` has to be empty (not both, not neither of them). + if (on1 == "" && on2 == "") || (on1 != "" && on2 != "") + throw( + ArgumentError( + "Invalid operator name: $name. Operator name is not \"Id\" or of the " * + "form \"A⋅\" or \"⋅A\"", + ), + ) + end + # name == "⋅A" -> on1 is an empty string + # name == "A⋅" -> on2 is an empty string + if on1 == "" + mat = try_op(OpName(on2), SiteType("S=1/2"); kwargs...) + return postmultiply(mat, st) + elseif on2 == "" + mat = try_op(OpName(on1), SiteType("S=1/2"); kwargs...) + return premultiply(mat, st) + else + # This should logically never happen but, just in case, we throw an error. + error("Unknown error with operator name $name") + end + else + error("Operator name $name is not \"Id\" or of the form \"A⋅\" or \"⋅A\"") + end +end diff --git a/src/spin_chain.jl b/src/spin_chain.jl new file mode 100644 index 0000000..5e26477 --- /dev/null +++ b/src/spin_chain.jl @@ -0,0 +1,220 @@ +export spin_chain, spin_chain_adjoint, spin_chain′ + +spin_chain(::SiteType, ::Vector{<:Real}, ::Vector{<:Real}, ::Vector{<:Int}) = nothing + +""" + spin_chain( + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sites::Vector{<:Index} + ) + +Return an OpSum object encoding the Hamiltonian part ``-i[H, –]`` of the GKSL equation +for a spin chain of frequencies `freqs` and coupling constants `coups`, on `sites`. +""" +function spin_chain(freqs::Vector{<:Real}, coups::Vector{<:Real}, sites::Vector{<:Index}) + @assert length(freqs) == length(coups) + 1 + @assert length(sites) == length(freqs) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = spin_chain(st, freqs, coups, sitenumber.(sites)) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"spin_chain\" function not found for Index tags $(tags(sites[1]))" + ), + ) +end + +function spin_chain( + ::SiteType"Fermion", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + h = OpSum() + + for (j, site) in enumerate(sitenumbers) + h += freqs[j], "n", site + end + + for (j, (site1, site2)) in enumerate(partition(sitenumbers, 2, 1)) + h += coups[j], "c†", site1, "c", site2 + h += coups[j], "c†", site2, "c", site1 + end + + return h +end + +function spin_chain( + ::SiteType"Electron", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + h = OpSum() + + for (j, site) in enumerate(sitenumbers) + h += freqs[j], "n↓", site + h += freqs[j], "n↑", site + end + + for (j, (site1, site2)) in enumerate(partition(sitenumbers, 2, 1)) + h += coups[j], "c†↓", site1, "c↓", site2 + h += coups[j], "c†↓", site2, "c↓", site1 + + h += coups[j], "c†↑", site1, "c↑", site2 + h += coups[j], "c†↑", site2, "c↑", site1 + end + + return h +end + +function spin_chain( + ::SiteType"vFermion", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + ℓ = OpSum() + for (j, site) in enumerate(sitenumbers) + ℓ += freqs[j] * gkslcommutator("N", site) + end + for (j, (site1, site2)) in enumerate(partition(sitenumbers, 2, 1)) + jws = jwstring(; start=site1, stop=site2) + ℓ += + coups[j] * ( + gkslcommutator("A†", site1, jws..., "A", site2) + + gkslcommutator("A", site1, jws..., "A†", site2) + ) + end + return ℓ +end + +function spin_chain( + ::SiteType"vS=1/2", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + ℓ = OpSum() + for (j, site) in enumerate(sitenumbers) + ℓ += freqs[j] * gkslcommutator("N", site) + end + for (j, (site1, site2)) in enumerate(partition(sitenumbers, 2, 1)) + ℓ += + coups[j] * ( + gkslcommutator("σ+", site1, "σ-", site2) + + gkslcommutator("σ-", site1, "σ+", site2) + ) + end + return ℓ +end + +function spin_chain( + ::SiteType"vElectron", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + ℓ = OpSum() + for (j, site) in enumerate(sitenumbers) + ℓ += freqs[j] * gkslcommutator("Ntot", site) + end + for (j, (site1, site2)) in enumerate(partition(sitenumbers, 2, 1)) + jws = jwstring(; start=site1, stop=site2) + ℓ += + coups[j] * ( + gkslcommutator("Aup†F", site1, jws..., "Aup", site2) - + gkslcommutator("AupF", site1, jws..., "Aup†", site2) + + gkslcommutator("Adn†", site1, jws..., "FAdn", site2) - + gkslcommutator("Adn", site1, jws..., "FAdn†", site2) + ) + end + return ℓ +end + +############################################################################################ + +function spin_chain_adjoint(::SiteType, ::Vector{<:Real}, ::Vector{<:Real}, ::Vector{<:Int}) + return nothing +end + +""" + spin_chain_adjoint( + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sites::Vector{<:Index} + ) + +Return an OpSum object encoding the adjoint of the Hamiltonian part ``-i[H, –]`` of the +GKSL equation for a spin chain of frequencies `freqs` and coupling constants `coups`, on +`sites`. +""" +function spin_chain_adjoint( + freqs::Vector{<:Real}, coups::Vector{<:Real}, sites::Vector{<:Index} +) + @assert length(freqs) == length(coups) + 1 + @assert length(sites) == length(freqs) + stypes = sitetypes(first(sites)) + for st in stypes + # Check if all sites have this type (otherwise skip this tag). + if all(i -> st in sitetypes(i), sites) + # If the type is shared, then try calling the function with it. + ℓ = spin_chain_adjoint(st, freqs, coups, sitenumber.(sites)) + # If the result is something, return that result. + if !isnothing(ℓ) + return ℓ + end + end + # Otherwise, try again with another type from the initial ones. + end + # Return an error if no implementation is found for any type. + return throw( + ArgumentError( + "Overload of \"spin_chain_adjoint\" function not found for " * + "Index tags $(tags(sites[1]))", + ), + ) +end + +# In the GKSL equation, the adjoint of -i[H,-] is simply i[H,-]. +# +function spin_chain_adjoint( + st::SiteType"vFermion", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + return -spin_chain(st, freqs, coups, sitenumbers) +end + +function spin_chain_adjoint( + st::SiteType"vS=1/2", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + return -spin_chain(st, freqs, coups, sitenumbers) +end + +function spin_chain_adjoint( + st::SiteType"vElectron", + freqs::Vector{<:Real}, + coups::Vector{<:Real}, + sitenumbers::Vector{<:Int}, +) + return -spin_chain(st, freqs, coups, sitenumbers) +end + +const spin_chain′ = spin_chain_adjoint diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..8e8eab0 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,886 @@ +""" + statenamestring(sn::StateName{T}) where T + +Return the name of the state "sn" as a string. +""" +function statenamestring(sn::StateName{T}) where T + return string(T) +end + +""" + sitenumber(i::Index) + +Return the site number of the given Index, i.e. the number in the "n=N" tag. +""" +function sitenumber(i::Index) + st = split(replace(string.(tags.(i)), "\"" => ""), ",") + sitetag = first(filter!(s -> occursin("n=", s), st)) + # TODO what if sitetag is empty? + siten = replace(sitetag, "n=" => "") + return parse(Int, siten) +end + +""" + jwstring(; start, stop, op::AbstractString="F") + +Return the vector ``(op, start + 1, op, start + 2, ..., op, stop - 1)``. +""" +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 +# of course be orthonormal wrt this inner product. +# The canonical basis or the Gell-Mann one are okay. + +""" + vec(A::Matrix, basis::Vector) + +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) + +Compute the matrix of coefficients of the linear map `L` wrt the basis `basis`. +The linearity of the map is not checked, so using this function with non-linear +functions leads to undefined results. +""" +function vec(L::Function, basis::Vector) + return [tr(bi' * L(bj)) for (bi, bj) in Base.product(basis, basis)] +end + +""" + partialtrace(sites::Vector{Index{Int64}}, v::MPS, j::Int) + +Compute the partial trace, on the `j`-th site of `sites`, of the matrix +represented, as vectorised, by the MPS `v`. +The result is a `Vector` containing the coordinates of the partial trace +(i.e. the result is still in vectorised form). +""" +function partialtrace(sites::Vector{Index{Int64}}, v::MPS, j::Int) + # Custom contraction sequence: we start from v[end] and we contract it with + # a vecId state; we contract the result with v[end-1] and so on, until we + # get to v[j]. # Then we do the same starting from v[1]. + x = ITensor(1.0) + for i in length(sites):-1:(j + 1) + x = v[i] * state(sites[i], "vecId") * x + end + y = ITensor(1.0) + for i in 1:(j - 1) + y = v[i] * state(sites[i], "vecId") * y + end + z = y * v[j] * x + # Now `vector(z)` is the coordinate vector of the partial trace. + return vector(z) +end + +# Generalised Gell-Mann matrices +# ============================== +""" + gellmannmatrix(j, k, dim) + +Return the `(j,k)` generalised Gell-Mann matrix of dimension `dim`, normalised +wrt the Hilbert-Schmidt inner product ``(A,B) = tr(A†B)``. +The matrices are indexed as follows: + + * if ``j > k`` the matrix is symmetric and traceless; + * if ``j < k`` the matrix is antisymmetric; + * if ``j = k`` the matrix is diagonal. + +In particular, ``j = k = dim`` gives a matrix proportional to the identity. +The two indices `j` and `k` determine the non-zero coefficients of the matrix. +The whole set of (different) Gell-Mann matrices that can be generated with this +function is a basis of ``Mat(ℂᵈⁱᵐ)``. +""" +function gellmannmatrix(j, k, dim) + if j > dim || k > dim || j < 0 || k < 0 + throw(DomainError) + end + m = zeros(ComplexF64, dim, dim) + if j > k + m[j, k] = 1 / sqrt(2) + m[k, j] = 1 / sqrt(2) + elseif k > j + m[j, k] = -im / sqrt(2) + m[k, j] = im / sqrt(2) + elseif j == k && j < dim + for i in 1:j + m[i, i] = 1 + end + m[j + 1, j + 1] = -j + m .*= sqrt(1 / (j * (j + 1))) + else + for i in 1:dim + m[i, i] = 1 / sqrt(dim) + end + end + return m +end + +""" + gellmannbasis(dim) + +Return a list containing the "Hermitian basis" of ``Mat(ℂᵈⁱᵐ)``, i.e. composed +of the ``dim²`` generalised Gell-Mann matrices. +""" +function gellmannbasis(dim) + return [gellmannmatrix(j, k, dim) for (j, k) in [Base.product(1:dim, 1:dim)...]] + # We need to splat the result from `product` so that the result is a list + # of matrices (a Vector) and not a Matrix. +end + +""" + canonicalmatrix(i, j, dim) + +Return the (`i`,`j`) element of the canonical basis of ``Mat(ℂᵈⁱᵐ)``, i.e. a +`dim`×`dim` matrix whose element on the `i`-th row and `j`-th column is ``1``, +and zero elsewhere. +""" +function canonicalmatrix(i, j, dim) + m = zeros(ComplexF64, dim, dim) + m[i, j] = 1 + return m +end + +""" + canonicalbasis(dim) + +Return a list of the matrices in the canonical basis of ``Mat(ℂᵈⁱᵐ)``. +The list is ordered corresponding to column-based vectorisation, i.e. + + canonicalbasis(dim)[j] = canonicalmatrix((j-1)%dim + 1, (j-1)÷dim + 1, dim) + +with ``j ∈ {1,…,dim²}``. With this ordering, +``vec(A)ⱼ = tr(canonicalbasis(dim)[j]' * A)``. +""" +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 +# ===================== +""" + vonneumannentropy(ψ::MPS, sites::Vector{Index{Int64}}, n::Int) + +Compute the entanglement entropy of the biparition ``(1,…,n)|(n+1,…,N)`` of +the system in state described by the MPS `ψ` (defined on the sites `sites`), +using its Schmidt decomposition. +""" +function vonneumannentropy(ψ₀::MPS, sites::Vector{Index{Int64}}, n::Int) + ψ = orthogonalize(ψ₀, n) + # Decompose ψ[n] in singular values, treating the Link between sites n-1 and n + # and the physical index as "row index"; the remaining index, the Link + # between n and n+1, is the "column index". + _, S, _ = svd(ψ[n], (linkind(ψ, n - 1), sites[n])) + # Compute the square of the singular values (aka the Schmidt coefficients + # of the bipartition), and from them the entropy. + sqdiagS = [S[j, j]^2 for j in ITensors.dim(S, 1)] + return -sum(p -> p * log(p), sqdiagS; init=0.0) +end + +# Chop (from Mathematica) +# ======================= +# Imitation of Mathematica's "Chop" function + +""" + chop(x::Real; tolerance=1e-10) + +Truncates `x` to zero if it is less than `tolerance`. +""" +function chop(x::Real; tolerance=1e-10) + return abs(x) > tolerance ? x : zero(x) +end + +""" + chop(x::Complex; tolerance=1e-10) + +Truncates the real and/or the imaginary part of `x` to zero if they are less +than `tolerance`. +""" +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) + +Concatenate `left` and `right`, returning `left` ``⊗`` `right`. +""" +function chain(left::MPS, right::MPS) + # This function is like mpnum's `chain`: it takes two MPSs ans concatenates + # them. The code is "inspired" from `ITensors.MPS` in mps.jl:308. + + midN = length(left) # The site with the missing link between the two MPSs. + # First of all we shift the Link tags of the given MPSs, so that the final + # enumeration of the tags is correct. + # Note that in each MPS the numbers in the Link tags do not follow the + # numbering of the Sites on which it is based, they always start from 1. + for j in eachindex(left) + replacetags!(left[j], "l=$j", "l=$j"; tags="Link") + replacetags!(left[j], "l=$(j-1)", "l=$(j-1)"; tags="Link") + end + for j in eachindex(right) + replacetags!(right[j], "l=$j", "l=$(midN+j)"; tags="Link") + replacetags!(right[j], "l=$(j-1)", "l=$(midN+j-1)"; tags="Link") + end + # "Shallow" concatenation of the MPSs (there's still a missing link). + M = MPS([left..., right...]) + # We create a "trivial" Index of dimension 1 and add it to the two sites + # which are not yet connected. + # The Index has dimension 1 because this is a tensor product between the + # states so there's no correlation between them. + missing_link = Index(1; tags="Link,l=$midN") + M[midN] = M[midN] * state(missing_link, 1) + M[midN + 1] = state(dag(missing_link), 1) * M[midN + 1] + + return M +end + +""" + chain(left::MPO, right::MPO) + +Concatenate `left` and `right`, returning `left` ``⊗`` `right`. +""" +function chain(left::MPO, right::MPO) + # Like the previous `chain`, but for MPOs. + midN = length(left)# The site with the missing link between the two MPOs. + for j in eachindex(right) + replacetags!(right[j], "l=$j", "l=$(midN+j)"; tags="Link") + replacetags!(right[j], "l=$(j-1)", "l=$(midN+j-1)"; tags="Link") + end + M = MPO([left..., right...]) + missing_link = Index(1; tags="Link,l=$midN") + M[midN] = M[midN] * state(missing_link, 1) + M[midN + 1] = M[midN + 1] * state(dag(missing_link), 1) + # The order of the Indexes in M[midN] and M[midN+1] ns not what we would get + # if we built an MPO on the whole system as usual; namely, the two "Link" + # Indexes are swapped. + # This however should not matter since ITensor does not care about the order. + return M +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...) + +Concatenate the given MPOs into a longer MPO, returning their tensor product. +""" +chain(a::MPO, b...) = chain(a, chain(b...)) + +""" + embed_slice(sites::Array{Index{Int64}}, range::UnitRange{Int}, slice::MPO) + +Embed `slice`, defined on a subset `range` of `sites`, into a MPO which covers +the whole `sites`. + +The MPO is extended by filling the empty spots with an "Id" operator, therefore +an operator with OpName "Id" is required to be defined for the SiteTypes of the +remaining sites. + +# Arguments +- `sites::Array{Index{Int64}}`: the sites of the whole system. +- `range::UnitRange{Int}`: the range spanned by `slice`. +- `slice::MPO`: the MPO to be extended. +""" +function embed_slice(sites::Array{Index{Int64}}, range::UnitRange{Int}, slice::MPO) + # TODO: compute automatically on which sites the MPO is defined, without + # having to supply the range explicitly as an argument. + if length(slice) != length(range) + throw(DimensionMismatch("slice and range must have the same size.")) + end + if !issubset(range, eachindex(sites)) + throw(BoundsError(range, sites)) + end + + if range[begin] == 1 && range[end] == length(sites) + mpo = slice + elseif range[begin] == 1 + mpo = chain(slice, MPO(sites[(range[end] + 1):end], "Id")) + elseif range[end] == length(sites) + mpo = chain(MPO(sites[1:(range[begin] - 1)], "Id"), slice) + else + mpo = chain( + MPO(sites[1:(range[begin] - 1)], "Id"), + slice, + MPO(sites[(range[end] + 1):end], "Id"), + ) + end + return mpo +end + +""" + embed_slice(sites::Array{Index{Int64}}, range::UnitRange{Int}, slice::MPS) + +Embed `slice`, defined on a subset `range` of `sites`, into a MPS which covers +the whole `sites` (to be interpreted as a vectorised operator). + +The MPS is extended by filling the empty spots with a "vecId" operator, +therefore an operator with OpName "vecId" is required to be defined for the +SiteTypes of the remaining sites. + +# Arguments +- `sites::Array{Index{Int64}}`: the sites of the whole system. +- `range::UnitRange{Int}`: the range spanned by `slice`. +- `slice::MPS`: the MPS to be extended. +""" +function embed_slice(sites::Array{Index{Int64}}, range::UnitRange{Int}, slice::MPS) + if length(slice) != length(range) + throw(DimensionMismatch("slice e range must have the same size.")) + end + if !issubset(range, eachindex(sites)) + throw(BoundsError(range, sites)) + end + + if range[begin] == 1 && range[end] == length(sites) + mpo = slice + elseif range[begin] == 1 + mpo = chain(slice, MPS(sites[(range[end] + 1):end], "vecId")) + elseif range[end] == length(sites) + mpo = chain(MPS(sites[1:(range[begin] - 1)], "vecId"), slice) + else + mpo = chain( + MPS(sites[1:(range[begin] - 1)], "vecId"), + slice, + MPS(sites[(range[end] + 1):end], "vecId"), + ) + end + return mpo +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) + +Return a list of Strings "(a,b)" formed by all adjacent items in `v`. +""" +function consecutivepairs(v::AbstractVector) + return string.("(", v[1:(end - 1)], ",", v[2:end], ")") +end + +""" + try_op(on::OpName, st::SiteType; kwargs...) + +Return the matrix of an ITensor operator, if it exists, trying first the +```julia +op(::OpName, ::SiteType; kwargs...) +``` +syntax, and then +```julia +op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) +``` +if the former returns nothing. +""" +function try_op(on::OpName, st::SiteType; kwargs...) + stname = String(ITensors.tag(st)) + opstring = String(ITensors.name(on)) + # Try calling a function of the form: + # op(::OpName, ::SiteType; kwargs...) + # which returns a Julia matrix + mat = ITensors.op(on, st; kwargs...) + if isnothing(mat) + # Otherwise try calling a function of the form + # op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) + dummy = siteind(stname) + Op = ITensor(prime(dummy), ITensors.dag(dummy)) + r = ITensors.op!(Op, on, st, dummy; kwargs...) + if isnothing(r) + throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for " * + "operator name \"$opstring\" and Index tag $(tags(dummy)).", + ), + ) + end + mat = matrix(Op) + end + if isnothing(mat) + throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for operator " * + "name \"$opstring\" and Index tag \"$stname\".", + ), + ) + end + return mat +end + +""" + try_op(on::OpName, st::SiteType, d::Int; kwargs...) + +Like try_op(on::OpName, st::SiteType; kwargs...) (see [`try_op`](@ref)), but with an +additional Int argument so that it can be used by SiteTypes without a fixed dimension. +""" +function try_op(on::OpName, st::SiteType, d::Int; kwargs...) + stname = String(ITensors.tag(st)) + opstring = String(ITensors.name(on)) + # Try calling a function of the form: + # op(::OpName, ::SiteType, ::Int; kwargs...) + # which returns a Julia matrix + mat = ITensors.op(on, st, d; kwargs...) + if isnothing(mat) + # Otherwise try calling a function of the form + # op!(::ITensor, ::OpName, ::SiteType, ::Index..., ::Int; kwargs...) + dummy = siteind(stname) + Op = ITensor(prime(dummy), ITensors.dag(dummy)) + r = ITensors.op!(Op, on, st, dummy, d; kwargs...) + if isnothing(r) + throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for " * + "operator name \"$opstring\" and Index tag $(tags(dummy)).", + ), + ) + end + mat = matrix(Op) + end + if isnothing(mat) + throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for operator " * + "name \"$opstring\" and Index tag \"$stname\".", + ), + ) + end + return mat +end