Move files over from PseudomodesTTEDOPA

Project.toml

@@ -4,8 +4,19 @@ authors = ["Davide Ferracin <davide.ferracin@protonmail.com> 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"
 

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"

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 <davide.ferracin@protonmail.com> and contributors",
     sitename="LindbladVectorizedTensors.jl",
     format=Documenter.HTML(;

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, ρ)
+```

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↓`
+

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′
 ```