List of all methods

apply1siteoperator!(A, O, sites::Int)

Apply an operator O on the MPS A. O is acting on only one site ::Int. The resulting MPS A is the MPS modified by the operator O.

apply1siteoperator!(A, O, sites::Vector{Int})

Apply an operator O on the MPS A. O is acting on several sites ::Vector{Int}. The resulting MPS A is the MPS modified by the operator O.

applympo!(A, H)

Apply an MPO H on the MPS A. H must have the same number of site than A. The resulting MPS A is the MPS modified by the MPO H.

chainmps(N::Int, site::Int, numex::Int)

Generate an MPS with numex excitations on site

The returned MPS will have bond-dimensions and physical dimensions numex+1

chainmps(N::Int, sites::Vector{Int}, numex::Int)

Generate an MPS with numex excitations of an equal super-position over sites

The returned MPS will have bond-dimensions and physical dimensions numex+1

displacedchainmps(A::Vector{Any}, N::Int, Nm::Int, γ::Any)

Given a MPS A, return a MPS B where the Nm-long chain is displaced by γ without displacing the N-long system.

electronkmps(N::Int, k::Vector{Int}, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])

Generate an MPS with 2 electrons in k-states k1 and k2.

electronkmps(N::Int, k::Int, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])

Generate an MPS for an electron with momentum k.

elementmpo(M, el...)

Return the element of the MPO M for the set of physical states el...

elementmps(A, el...)

Return the element of the MPS A for the set of physical states el...

Examples

julia> A = chainmps(6, [2,4], 1);

julia> elementmps(A, 1, 2, 1, 1, 1, 1)
0.7071067811865475

julia> elementmps(A, 1, 1, 1, 2, 1, 1)
0.7071067811865475

julia> elementmps(A, 1, 2, 1, 2, 1, 1)
0.0

julia> elementmps(A, 1, 1, 1, 1, 1, 1)
0.0

entanglemententropy(A)

For a list of tensors A representing a right orthonormalized MPS, compute the entanglement entropy for a bipartite cut for every bond.

modemps(N::Int, k::Vector{Int}, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])

Generate an MPS with numex excitations of an equal superposition of modes k of a bosonic tight-binding chain.

modemps(N::Int, k::Int, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])

Generate an MPS with numex excitations of mode k of a bosonic tight-binding chain.

chainparams takes the form [e::Vector, t::Vector] where e are the on-site energies and t are the hoppping parameters.

The returned MPS will have bond-dimensions and physical dimensions numex+1

mpsembed!(A::Vector, Dmax::Int)

Embed MPS A in manifold of max bond-dimension Dmax

mpsleftnorm!(A::Vector, jq::Int=length(A))

Left orthoganalise MPS A up to site jq.

mpsmixednorm!(A::Vector, OC::Int)

Put MPS A into mixed canonical form with orthogonality centre on site OC.

mpsrightnorm!(A::Vector, jq::Int=1)

Right orthoganalise MPS A up to site jq.

multiply(M1::Vector, M2::Vector)

Calculates M1*M2 where M1 and M2 are MPOs

normmps(A::Vector; mpsorthog=:None)

Calculate norm of MPS A.

Setting mpsorthog=:Right/:Left will calculate the norm assuming right/left canonical form. Setting mpsorthog=OC::Int will cause the norm to be calculated assuming the orthoganility center is on site OC. If mpsorthog is :None the norm will be calculated as an MPS-MPS product.

orthcentersmps(A::Vector)

Compute the orthoganality centres of MPS A.

Return value is a list in which each element is the corresponding site tensor of A with the orthogonality centre on that site. Assumes A is right normalised.

physdims(M)

Return the physical dimensions of an MPS or MPO M.

productstatemps(physdims::Dims, Dmax=1; state=:Vacuum, mpsorthog=:Right)

Return an MPS representing a product state with local Hilbert space dimensions given by physdims.

By default all bond-dimensions will be 1 since the state is a product state. However, to embed the product state in a manifold of greater bond-dimension, Dmax can be set accordingly.

The indvidual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state (where the state of each site is represented by a column vector with a 1 in the first row and zeros elsewhere). Setting state=:FullOccupy will produce an MPS in which each site is fully occupied (ie. a column vector with a 1 in the last row and zeros elsewhere).

The argument mpsorthog can be used to set the gauge of the resulting MPS.

Example

julia> ψ = unitcol(1,2); d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1

julia> cpars = chaincoeffs_ohmic(N, α, s)

julia> H = spinbosonmpo(ω0, Δ, d, N, cpars)

julia> A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>

productstatemps(N::Int, d::Int, Dmax=1; state=:Vacuum, mpsorthog=:Right)

Return an N-site MPS with all local Hilbert space dimensions given by d.

randmps(N::Int, d::Int, Dmax::Int, T=Float64)

Construct a random, N-site, right-normalised MPS with all local Hilbert space dimensions given by d.

randmps(physdims::Dims{N}, Dmax::Int, T::Type{<:Number} = Float64) where {N}

Construct a random, right-normalised MPS with local Hilbert space dimensions given by physdims and max bond-dimension given by Dmax.

T specifies the element type, eg. use T=ComplexF64 for a complex valued MPS.

reversempo!(M)

Reverse the left and right dimensions of the MPO M. The resulting MPO M is the reversed MPO.

reversempo(M)

Reverse the left and right dimensions of the MPO M.

reversemps!(A)

Reverse the left and right dimensions of the MPS A. The resulting MPS A is the reversed MPS.

reversemps(A)

Reverse the left and right bond-dimensions of the MPS A.

svdmps(A)

For a right normalised mps A compute the full svd spectrum for a bipartition at every bond.

Base.reshape(x::Number, dims...)

Reshape any matrix with dimensions "dims"

Example

julia> A = rand(2,3,2)

julia> B = reshape(A,4,3)

julia> size(B) == (4,3) 

U, S, Vd = svdtrunc(A; truncdim = max(size(A)...), truncerr = 0.)

Perform a truncated SVD, with maximum number of singular values to keep equal to truncdim or truncating any singular values smaller than truncerr. If both options are provided, the smallest number of singular values will be kept. Unlike the SVD in Julia, this returns matrix U, a diagonal matrix (not a vector) S, and Vt such that A ≈ U * S * Vt

MPOtoVector(mpo::MPO)

Convert an ITensors chain MPO into a form compatible with MPSDynamics

chainprop(t, cparams)

Propagate an excitation placed initially on the first site of a tight-binding chain with parameters given by cparams for a time t and return occupation expectation for each site.

cosineh(omega, bet)

Calculates the hyperbolic cosine function function based on the input parameters, for the Bogoliubov transformation necessary for the thermofield transformation.

Arguments

• omega::Float64: The frequency parameter.
• bet::Float64: The beta parameter.

Returns

• Float64: The result of the modified cosine function.

disp(d,ωvib,m)

Displacement operator $X = \frac{\sqrt{2}}{2\sqrt{m \omega_{vib}}}(a + a^{\dagger})$

disp(d)

Mass and frequency-weighted displacement operator $X = \frac{1}{2}(a + a^{\dagger})$

dynamap(ps1,ps2,ps3,ps4)

Calculate complete dynamical map to time step at which ps1, ps2, ps3 and ps4 are specified.

Arguments

• ps1 : time evolved system density matrix starting from initial state up
• ps2 : time evolved system density matrix starting from initial state down
• ps3 : time evolved system density matrix starting from initial state (up + down)/sqrt(2)
• ps4 : time evolved system density matrix starting from initial state (up - i*down)/sqrt(2)

eigenchain(cparams; nummodes=nothing)

endsiteocc(t, cparams)

entropy(rho)

findchainlength(T, ωc::Float64, β=Inf)

Estimate length of chain using universal asymptotic properties of chain mapped environments given the simulation time T, the bath cut-off frequency ωc, and the inverse temperature β.

findchainlength(T, cparams::Vector; eps=10^-6, verbose=false)

Estimate length of chain required for a particular set of chain parameters by calculating how long an excitation on the first site takes to reach the end. The chain length is given as the length required for the excitation to have just reached the last site after time T. The initial number of sites in cparams has to be larger than the findchainlength result.

measurecorrs(oper, , e::Vector, t::Vector)

Parameters

oper: Square matrix (Matrix{Float64}) representing the operator to be transformed. e: Vector (Vector{Float64}) of diagonal (on-site energy) chain coefficients. t: Vector (Vector{Float64}) of off-diagonal (hopping terms) chain coefficients.

Returns

Matrix{Float64}: This matrix is the operator oper transformed back from the chain representation to the representation corresponding to the extended bath. The resulting operator represents quantities like mode occupations or other properties in the basis of environmental modes associated with specific frequencies $\omega_i$.

Description

This function performs a basis transformation of the operator oper. Specifically, this transformation reverses the unitary transformation that maps the extended bath Hamiltonian into the chain representation.

measuremodes(adaga, e=1.0, t=1.0)

measuremodes(A, chainsection::Tuple{Int64,Int64}, U::AbstractArray)

for longer chains it can be worth calculating U in advance

measuremodes(A, chainsection::Tuple{Int64,Int64}, e::Array{Float64,1}, t::Array{Float64,1})

measuremodes(adaga, e::Vector, t::Vector)

physical_occup(corr_constr, corr_destr, omega, occup, b, M)

Calculates the physical occupation based on correlation matrices, omega values, and other parameters. The physical occupation in the original frequency environment is computed by reverting the thermofield transformation.

Arguments

• corr_constr::Matrix{ComplexF64}: The correlation construction matrix.
• corr_destr::Matrix{ComplexF64}: The correlation destruction matrix.
• omega::Vector{Float64}: The omega values.
• occup::Matrix{Float64}: The occupation matrix.
• b::Float64: The beta parameter.
• M::Int: The number of points for interpolation.

Returns

• Vector{Float64}: The physical occupation values.

randisometry([T=Float64], dims...)

Construct a random isometry

rmsd(ob1, ob2)

Calculate the root mean squared difference between two measurements of an observable over the same time period.

sineh(omega, bet)

Calculates the hyperbolic sine function function based on the input parameters, for the Bogoliubov transformation necessary for the thermofield transformation.

Arguments

• omega::Float64: The frequency parameter.
• bet::Float64: The beta parameter.

Returns

• Float64: The result of the modified cosine function.

 therHam(psi, site1, site2)

Calculates Hβ such that ρ = e^(-βH) for some density matrix ρ obatined from tracing out everything outside the range [site1,site2] in the MPS psi

thermaloccupations(β, cparams...)

addchild!(tree::Tree, id::Int)

Add child to node id of tree.

addchildren!(tree::Tree, id::Int, n::Int)

Add n children to node id of tree.

findchild(node::TreeNode, id::Int)

Return integer corresponding to the which number child site id is of node.

randtree(numnodes::Int, maxdegree::Int)

Construct a random tree with nummodes modes and max degree maxdegree.

measure1siteoperator(net::TreeNetwork, O, id::Int)

For a Tree, compute the local expectation value of a one-site operator O for site id.

measure1siteoperator(net::TreeNetwork, O, sites::Tuple{Int,Int})

For a Tree, compute the local expectation value of a one-site operator O for the specified site range.

measure2siteoperator(net::TreeNetwork, O1, O2, sites::Tuple{Int,Int})

For a Tree, compute the local expectation value of two one-site operators O1 and O2 for the specified site range.

bonddims(A::TreeNetwork)

Return the bon-dimension of a tree-MPS A.

mpsembed(A::TreeNetwork, Dmax::Int)

Embed tree-MPS A in manifold of max bond-dimension Dmax.

mpsmixednorm!(A::TreeNetwork, id::Int)

Normalise tree-MPS A such that orthogonality centre is on site id.

mpsmoveoc!(A::TreeNetwork, id::Int)

Move the orthogonality centre of right normalised tree-MPS A to site id.

This function will be more efficient than using mpsmixednorm! if the tree-MPS is already right-normalised.

mpsrightnorm!(A::TreeNetwork)

When applied to a tree-MPS, right normalise towards head-node.

mpsshiftoc!(A::TreeNetwork, newhd::Int)

Shift the orthogonality centre by one site, setting new head-node newhd.

normmps(net::TreeNetwork; mpsorthog=:None)

When applied to a tree-MPS mpsorthog=:Left is not defined.

physdims(M::TreeNetwork)

Return the physical dimensions of a tree-MPS or tree-MPO M.

productstatemps(tree_::Tree, physdims::Dims, Dmax::Int=1; state=:Vacuum)

Return a tree-MPS representing a product state with local Hilbert space dimensions given by physdims.

By default all bond-dimensions will be 1 since the state is a product state. However, to embed the product state in a manifold of greater bond-dimension, Dmax can be set accordingly.

The indvidual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state (where the state of each site is represented by a column vector with a 1 in the first row and zeros elsewhere). Setting state=:FullOccupy will produce an MPS in which each site is fully occupied (ie. a column vector with a 1 in the last row and zeros elsewhere).

Example

julia> ψ = unitcol(1,2); d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1

julia> cpars = chaincoeffs_ohmic(N, α, s)

julia> H = spinbosonmpo(ω0, Δ, d, N, cpars, tree=true)

julia> A = productstatemps(H.tree, physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # tree-MPS representation of |ψ>|Vacuum>

randmps(tree::Tree, physdims, Dmax::Int, T::Type{<:Number} = Float64)

Construct a random, right-normalised, tree-MPS, with structure given by tree and max bond-dimension given by Dmax.

The local Hilbert space dimensions are specified by physdims which can either be of type Dims{length(tree)}, specifying the dimension of each site, or of type Int, in which case the same local dimension is used for every site.

tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F::Vector, id::Int; verbose=false, kwargs...)

Propagates the tree-MPS A with the tree-MPO M following the 1-site TDVP method. The sweep is done back and forth with a time step dt/2. F represents the merged left and right parts of the site being propagated.

OneSiteObservable(name,op,sites)

Computes the local expectation value of the one-site operator op on the specified sites. Used to define one-site observables that are obs and convobs parameters for the runsim function.

OneSiteObservable(name,op)

Computes the local expectation value of the one-site operator op on the every site. Used to define one-site observables that are obs and convobs parameters for the runsim function.

TwoSiteObservable(name,op1,op2,sites1=nothing,sites2=nothing)

Computes the local expectation value of operators op1 and op2 where op1 acts on sites1 and op2 acts on sites2. Used to define several-site observables that are obs and convobs parameters for the runsim function.

measure(A, O; kwargs...)

Measure observable O on mps state A

 measure(A::Vector, Os::Vector; kwargs...)

Caculate expectation of Os on MPS A.

measure1siteoperator(A::Vector, O, site::Int)

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for a single site.

measure1siteoperator(A::Vector, O, chainsection::Tuple{Int64,Int64})

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for a chainsection.

measure1siteoperator(A::Vector, O, sites::Vector{Int})

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site or just one if it is specified.

For calculating operators on single sites this will be more efficient if the site is on the left of the mps.

measure1siteoperator(A::Vector, O)

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site.

 measure2siteoperator(A::Vector, M1, M2, j1, j2)

Caculate expectation of M1*M2 where M1 acts on site j1 and M2 acts on site j2, assumes A is right normalised.

 measure2siteoperator(A::Vector, M1, M2, sites1::Vector{Int}, sites2::Vector{Int})

Caculate expectation of M1*M2 where M1 acts on sites1 and M2 acts on sites2, assumes A is right normalised.

measurempo(A::Vector, M::Vector, sites::Tuples{Int,Int})

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of the MPO M on specified sites.

measurempo(A::Vector, M::Vector)

For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of the MPO M on every site.

 rhoreduced_1site(A::Vector, site::Int=1)

Caculate the reduced density matrix of the MPS A at the specified site.

 rhoreduced_2sites(A::Vector, site::Tuple{Int, Int})

Caculate the reduced density matrix of the MPS A of two neigbour sites. The resulting dimensions will be the four physical dimensions in total, corresponding to the dimensions of the two sites

measure(A, obs::FockError; t=0, kwargs...)

Return the measure of the observable obs on the MPS A.

chaincoeffs_ohmic(N, α, s; ωc=1, soft=false)

Generate chain coefficients $[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$ for an Harmonic bath at zero temperature with a power law spectral density given by:

soft cutoff: $J(ω) = 2αω_c (\frac{ω}{ω_c})^s \exp(-ω/ω_c)$

hard cutoff: $J(ω) = 2αω_c (\frac{ω}{ω_c})^s θ(ω-ω_c)$

The coefficients parameterise the chain Hamiltonian

$H = H_S + c_0 A_S⊗B_0+\sum_{i=0}^{N-1}t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$

which is unitarily equivalent (before the truncation to N sites) to

$H = H_S + A_S⊗\int_0^∞dω\sqrt{J(ω)}B_ω + \int_0^∞dωωb_ω^\dagger b_ω$.

correlatedenvironmentmpo(R::Vector, Nm::Int, d::Int; chainparams, fnamecc::String, s=1, α=1, ωc=1, c_phonon=1, β="inf", issoft=false)

Generate a MPO for a one-dimensional bosonic bath spatially correlated to a multi-component system

$H_B + H_{int} = \int_{-∞}^{+∞} dk ω_k b_k^\dagger b_k + ∑_j \int_{-∞}^{+∞}dk \sqrt{J(k)}(A_j b_k e^{i k R_j} + h.c.)$.

The interactions between the system and the chain-mapped bath are long range, i.e. each site interacts with all the chain modes. The spectral density is assumed to be Ohmic $J(ω) = 2αωc(ω/ωc)^s$.

Arguments

• R: List of system's components positions
• Nm: Number of chain modes. The actual number of mode will be doubled to account for the left and right moving excitations.
• d: Local Hilbert space dimension of the bath modes
• chainparams: chain parameters, of the form chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
• fnamecc: Path to a file containing pre-computed long-range coupling coefficient. If not provided, the coupling coefficients will be computed and stored.
• s: Ohmicity
• α: Kondo parameter
• ωc: Bath cut-off frequency
• c_phonon: Speed of sound in the bath
• β: Inverse temperature
• issoft: Is the cut-off of the Ohmic SD soft or hard?

hbathchain(N::Int, d::Int, chainparams, longrangecc...; tree=false, reverse=false, coupletox=false)

Generate MPO representing a tight-binding chain of N oscillators with d Fock states each. Chain parameters are supplied in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$. The output does not itself represent a complete MPO but will possess an end which is open and should be attached to another tensor site, usually representing the system.

Arguments

• reverse: If reverse=true create a chain were the last (i.e. Nth) site is the site which couples to the system
• coupletox: Used to choose the form of the system coupling. coupletox=true gives a non-number conserving coupling of the form $H_{\text{I}}= A_{\text{S}}(b_{0}^\dagger + b_0)$ where $A_{\text{S}}$ is a system operator, while coupletox=false gives the number-converving coupling $H_{\text{I}}=(A_{\text{S}} b_{0}^\dagger + A_{\text{S}}^\dagger b_0)$
• tree: If true the resulting chain will be of type TreeNetwork; useful for construcing tree-MPOs

Example

One can constuct a system site tensor to couple to a chain by using the function up to populate the tensor. For example, to construct a system site with Hamiltonian Hs and coupling operator As, the system tensor M is constructed as follows for a non-number conserving interaction:

u = one(Hs) # system identity
M = zeros(1,3,2,2)
M[1, :, :, :] = up(Hs, As, u)

The full MPO can then be constructed with:

Hmpo = [M, hbathchain(N, d, chainparams, coupletox=true)...]

Similarly for a number conserving interaction the site tensor would look like:

u = one(Hs) # system identity
M = zeros(1,4,2,2)
M[1, :, :, :] = up(Hs, As, As', u)

And the full MPO would be

Hmpo = [M, hbathchain(N, d, chainparams; coupletox=false)...]

heisenbergmpo(N::Int, J=1.0) = xyzmpo(N; Jx=J)

Generate MPO for the N-spin Heisenberg XXX model, defined by the Hamiltonian

$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - J σ_z^{n} σ_z^{n+1}$

with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.

ibmmpo(ω0, d, N, chainparams; tree=false)

Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators with the interacting boson model (IBM), defined by the Hamiltonian

$H = \frac{ω_0}{2}σ_z + c_0σ_z(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.

The spin is on site 1 of the MPS and the bath modes are to the right.

This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by

$H = \frac{ω_0}{2}σ_z + σ_z\int_0^∞ dω\sqrt{J(ω)}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.

The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.

interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)

Generate MPO for a tight-binding chain of N fermionic sites with a single impurity site (fermionic as well) of energy ϵd. The impurity is coupled to two leads, each described by a set of chain parameters. The interactions are next-nearest-neighbour, with the first site corresponding to the impurity, and the two chains organised in an interleaved fashion.

Arguments

• N::Int: number of sites in the chain
• ϵd::Real: energy of the impurity site at the first site, as Ed - μ, where μ is the chemical potential
• chainparams1::Array{Real,1}: chain parameters for the first lead
• chainparams2::Array{Real,1}: chain parameters for the second lead

The chain parameters are given in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$.

isingmpo(N; J=1.0, h=1.0)

Generate MPO for the N-spin 1D Ising model with external field $\vec{h} = (0,0,h)$, defined by the Hamiltonian

$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} + \sum_{n=1}^{N}(- h_z σ_z^{n})$

with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.

longrange_isingmpo(N::Int, α::Float64=0.; J=1.0, h=1.0) = longrange_xyzmpo(N, α; Jx=J, Jy=0., Jz=0., hz=h, hx=0.)

longrange_xyzmpo(N::Int, α::Float64=0.; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)

Gennerate MPO for the N-spin long-range XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian

multisitempo(N::Int, E=[]::Vector, J=[]::Vector, As=[]::Vector)

Generate a MPO for a N-site multi-component system to be coupled with a correlated environment.

Arguments:

• N: number of sites
• E: list of on-site energies
• J: list of sites tunnelling energies
• As: list of system operators coupling to the environment

nearestneighbourmpo(tree_::Tree, h0, A, Ad = A')

nearestneighbourmpo(N::Int, h0, A, Ad = A')

protontransfermpo(ω0e,ω0k,x0e,x0k, Δ, dRC, d, N, chainparams, RCparams, λreorg)

Generate a MPO for a system described in space with a reaction coordinate (RC) tensor. The Hamiltonian of the two-level system and of the reaction coordinate tensor reads

$H_S + H_{RC} + H_{int}^{S-RC} = \omega^0_{e} |e\rangle \langle e| + \omega^0_{k} |k\rangle \langle k| + \Delta (|e\rangle \langle k| + |k\rangle \langle e|) + \omega_{RC} (d^{\dagger}d + \frac{1}{2}) + g_{e} |e\rangle \langle e|( d + d^{\dagger})+ g_{k} |k \rangle \langle k|( d + d^{\dagger})$ The RC tensor is coupled to a bosonic bath, taking into account the induced reorganization energy $H_B + H_{int}^{RC-B} = \int_{-∞}^{+∞} dk ω_k b_k^\dagger b_k - (d + d^{\dagger})\int_0^∞ dω\sqrt{J(ω)}(b_ω^\dagger+b_ω) + \lambda_{reorg}(d + d^{\dagger})^2$ with $\lambda_{reorg} = \int \frac{J(\omega)}{\omega}d\omega.$

Arguments

• ω0e: enol energy at reaction coordinate value x=0
• ω0k: keto energy at reaction coordinate value x=0
• x0e: enol equilibrium displacement
• x0k: keto equilibrium displacement
• Δ: direct coupling between enol and keto
• dRC: fock space of the RC tensor
• d: number of Fock states of the chain modes
• N: length of the bosonic chain
• chainparams: chain parameters, of the form chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
• RCparams: RC tensor parameter, of the form RCparams=$[ω_RC,-g/x]$
• λreorg: reorganization energy

puredephasingmpo(ΔE, dchain, Nchain, chainparams; tree=false)

Generate MPO for a pure dephasing model, defined by the Hamiltonian $H = \frac{ΔE}{2} σ_z + \frac{σ_z}{2} c_0 (b_0^\dagger + b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N-1} ϵ_i b_i^\dagger b_i$

The spin is on site 1 of the MPS and the bath modes are to the right.

Arguments

• ΔE::Real: energy splitting of the spin
• dchain::Int: physical dimension of the chain sites truncated Hilbert spaces
• Nchain::Int: number of sites in the chain
• chainparams::Array{Real,1}: chain parameters for the bath chain. The chain parameters are given in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$.
• tree::Bool: if true, return a TreeNetwork object, otherwise return a vector of MPO tensors

readchaincoeffs(fdir, params...)

spinbosonmpo(ω0, Δ, d, N, chainparams; rwa=false, tree=false)

Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators, defined by the Hamiltonian

$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0σ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.

The spin is on site 1 of the MPS and the bath modes are to the right.

This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by

$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{J(ω)}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.

The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.

The rotating wave approximation can be made by setting rwa=true.

spinchainmpo(N::Int; J=1.0, hz=1.0, hx=0.0, i=div(N,2))

tightbinding_mpo(N, ϵd, chainparams1, chainparams2)

Generate MPO for a tight-binding chain of N fermionic sites with a single impurity site (fermionic as well) of energy ϵd at the center. The impurity is coupled to two leads, each described by a set of chain parameters. The interactions are nearest-neighbour, with the first N/2-1 sites corresponding to the first lead, the Nth site corresponding to the impurity, and the rest of the sites corresponding to the second lead.

Arguments

• N::Int: number of sites in the chain
• ϵd::Real: energy of the impurity site at the center, as Ed - μ, where μ is the chemical potential
• chainparams1::Array{Real,1}: chain parameters for the first lead
• chainparams2::Array{Real,1}: chain parameters for the second lead

The chain parameters are given in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$.

tightbindingmpo(N::Int, d::Int; J=1.0, e=1.0)

tunnelingmpo(ϵ, delta, α, s, β, d::Int, nummodes::Int; tree=false, ωc=1)

twobathspinmpo(ω0, Δ, Nl, Nr, dl, dr, chainparamsl=[fill(1.0,N),fill(1.0,N-1), 1.0], chainparamsr=chainparamsl; tree=false)

Generate MPO for a spin-1/2 coupled to two chains of harmonic oscillators, defined by the Hamiltonian

$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0^rσ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N_r-1} t_i^r (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N_r} ϵ_i^rb_i^\dagger b_i + c_0^lσ_x(d_0^\dagger+d_0) + \sum_{i=0}^{N_l-1} t_i^l (d_{i+1}^\dagger d_i +h.c.) + \sum_{i=0}^{N_l} ϵ_i^l d_i^\dagger d_i$.

The spin is on site $N_l + 1$ of the MPS, surrounded by the left chain modes and the right chain modes.

This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by

$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{J(ω)}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ωi + σ_x\int_0^∞ dω\sqrt{J^l(ω)}(d_ω^\dagger+d_ω) + \int_0^∞ dω ωd_ω^\dagger d_ω$.

The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature. The two chains can have a different spectral density.

xxzmpo(N::Int, Δ = 1.0, J=1.0) = xyzmpo(N; Jx=J, Jy=J, Jz=J*Δ)

Generate MPO for the N-spin XXZ model, defined by the Hamiltonian

$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - \Delta J σ_z^{n} σ_z^{n+1}$

with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.

xyzmpo(N::Int; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)

Generate MPO for the N-spin XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian

$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} - J_y σ_y^{n} σ_y^{n+1} - J_z σ_z^{n} σ_z^{n+1} + \sum_{n=1}^{N}(- h_x σ_x^{n} - h_z σ_z^{n})$

with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.

chaincoeffs_fermionic(nummodes, β, chain; ϵ=nothing, J=nothing, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true)

Generate chain coefficients $[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$ for a fermionic bath at the inverse temperature β.

Arguments

• nummodes: Number of bath modes
• β: inverse temperature
• chain: 1 if the chain modes are empty, 2 if the chain modes are filled
• ϵ: user-provided dispersion relation. Should be a function f(x) where x is the wavenumber
• J: user-provided spectral density. Should be a function f(x) where x is the wavenumber
• ωc: the maximum frequency allowwed in the spectral density
• mc: the number of component intervals
• mp: the number of points in the discrete part of the measure (mp=0 if there is none)
• iq: a parameter to be set equal to 1, if the user provides his or her own quadrature routine, and different from 1 otherwise
• idelta: a parameter whose default value is 1, but is preferably set equal to 2, if iq=1 and the user provides Gauss-type quadrature routines
• procedure: choice between the Stieltjes and the Lanczos procedure
• AB: component intervals
• Mmax: maximum number of integration points
• save: if true the coefficients are saved

chaincoeffs_finiteT(nummodes, β, ohmic=true; α, s, J, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true)

Generate chain coefficients $[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$ for a harmonic bath at the inverse temperature β.

By default a Ohmic spectral density $J(ω) = 2αω_c (\frac{ω}{ω_c})^s θ(ω-ω_c)$ is considered. Users can provide their own spectral density.

Arguments

• nummodes: Number of bath modes
• β: inverse temperature
• ohmic: true if the spectral density is Ohmic, false if the user provides its own spectral density
• α: Kondo parameter of the Ohmic spectral density
• s: ohmicity
• J: user-provided spectral density. Should be a function f(x,i) where x is the frequency and i ∈ {1,...,mc} labels the intervals on which the SD is defined
• ωc: the maximum frequency of the Ohmic spectral density
• mc: the number of component intervals
• mp: the number of points in the discrete part of the measure (mp=0 if there is none)
• iq: a parameter to be set equal to 1, if the user provides his or her own quadrature routine, and different from 1 otherwise
• idelta: a parameter whose default value is 1, but is preferably set equal to 2, if iq=1 and the user provides Gauss-type quadrature routines
• procedure: choice between the Stieltjes and the Lanczos procedure
• AB: component intervals
• Mmax: maximum number of integration points
• save: if true the coefficients are saved

gauss(N,ab)

Gauss quadrature rule for N sites on an interval ab. Given a weight function w encoded by the Nx2 array ab of the first N recurrence coefficients for the associated orthogonal polynomials, the first column of ab containing the N alpha-coefficients and the second column the N beta-coefficients, the call xw = gauss(N,ab) generates the nodes and weights xw of the N-point Gauss quadrature rule for the weight function w. The nodes, in increasing order, are stored in the first column, the n corresponding weights in the second column, of the Nx2 array xw.

lanczos(N,xw)

Given the discrete inner product whose nodes are contained in the first column, and whose weights are contained in the second column, of the Nx2 array xw, the call ab=lanczos(N,xw) generates the first N recurrence coefficients ab of the corresponding discrete orthogonal polynomials. The N alpha-coefficients are stored in the first column, the N beta-coefficients in the second column, of the Nx2 array ab.

The script is adapted from the routine RKPW in W.B. Gragg and W.J. Harrod, ``The numerically stable reconstruction of Jacobi matrices from spectral data'', Numer. Math. 44 (1984), 317-335.

stieltjes(N,xw)

Discretized Stieltjes procedure. Given the discrete inner product whose nodes are contained in the first column, and whose weights are contained in the second column, of the Nx2 array xw, the call ab=stieltjes(N,xw) generates the first N recurrence coefficients ab of the corresponding discrete orthogonal polynomials. The N alpha- coefficients are stored in the first column, the N beta-coefficients in the second column, of the Nx2 array ab.

runsim(dt, tmax, A, H; 
	method=:TDVP1, 
	machine=LocalMachine(), 
	params=[], 
	obs=[], 
	convobs=[],
            convparams=error("Must specify convergence parameters"),
            save=false,
            plot=save,
            savedir=string(homedir(),"/MPSDynamics/"),
            unid=randstring(5),
            name=nothing,
	kwargs...
	)

Propagate the MPS A with the MPO H up to time tmax in time steps of dt. The final MPS is returned to A and the measurement data is returned to dat

Arguments

• method: Several methods are implemented in MPSDynamics. :TDVP1 refers to 1-site TDVP on tree and chain MPS, :TDVP2 refers to 2-site TDVP on chain MPS, :DTDVP refers to a variant of 1-site TDVP with dynamics bond-dimensions on chain MPS
• machine: LocalMachine() points local ressources, RemoteMachine() points distant ressources
• params: list of parameters written in the log.txt file to describe the dynamics. Can be listed with @LogParams().
• obs: list of observables that will be measured at every time step for the most accurate convergence parameter supplied to convparams
• convobs: list of observables that will be measure at every time step for every convergence parameter supplied to convparams
• convparams: list of convergence parameter with which the propagation will be calculated for every parameter. At each parameter, convobs are measured while obs are measured only for the most accurate dynamics
• save: Used to choose whether the data will also be saved to a file
• plot: Used to choose whether plots for 1D observables will be automatically generated and saved along with the data
• savedir: Used to specify the path where resulting files are stored
• unid: Used to specify the name of the directory containing the resulting files
• name: Used to describe the calculation. This name will appear in the log.txt file