Update doc, correct typos

docs/src/nutshell.md

@@ -18,7 +18,7 @@ This enables efficient representation of the full {System + Environment} state a
 
 To go further, you will find in this documentation:
 * a [User Guide](@ref) explaining how to use the package
-* several examples illustrating different features of the package: adaptive time-evolution methods, [Time-dependent Hamiltonian](@ref), or fermionic environments with [The Anderson Impurity Model](@ref)
+* several examples illustrating different features of the package: adaptive time-evolution methods in [The Spin-Boson Model](@ref), [Time-dependent Hamiltonian](@ref), or fermionic environments with [The Anderson Impurity Model](@ref)
 * a [Theoretical Background](@ref) covering the chain mapping procedure, flying over the basics of tensor networks, and summarizing the TDVP methods
 
 If this package was useful in your work, do not forget [Citation](@ref).

docs/src/theory.md

@@ -32,7 +32,7 @@ where $P_n(\omega)$ are orthonormal polynomials such that
 and the inverse transformation is
 
 ```math
-    \hat{b}_n = \int_0^{+\infty} U_n(\omega)\hat{a}_\omega\mathrm{d}\omega\ .
+    \hat{b}_n = \int_0^{+\infty} U^*_n(\omega)\hat{a}_\omega\mathrm{d}\omega\ .
 ```
 
 Note that the orthonormality of the polynomials ensures the unitarity of the transformation defined in Eq.(2).
@@ -41,12 +41,12 @@ The mapping from a continuous set of modes to a (still infinite) discrete set mi
 Under this transformation, the Hamiltonian in Eq.(1) becomes
 
 ```math
-    \hat{H}= \hat{H}_S + \sum_{n=0}^{+\infty}\varepsilon_n\hat{b}_n^\dagger\hat{b}_n + t_n(\hat{b}_{n+1}^\dagger\hat{b}_n + \mathrm{h.c.}) + \kappa\hat{A}_S(\hat{b}_0 + \hat{b}_0^\dagger)\ .
+    \hat{H}= \hat{H}_S + \sum_{n=0}^{+\infty}\varepsilon_n\hat{b}_n^\dagger\hat{b}_n + t_n(\hat{b}_{n+1}^\dagger\hat{b}_n + \mathrm{h.c.}) + c_0 \hat{A}_S(\hat{b}_0 + \hat{b}_0^\dagger)\ .
 ```
 
 Hence, this mapping transforms the normal bath Hamiltonian into a tight-binding Hamiltonian with on-site energies $\varepsilon_n$ and hopping energies $t_n$.
 Another important consequence of this mapping is that now the system only interacts with the first mode $n = 0$ of the chain-mapped environment.
-The chain coefficients $\varepsilon_n$, $t_n$, and the coupling $\kappa$ depend solely on the SD.
+The chain coefficients $\varepsilon_n$, $t_n$, and the coupling $c_0$ depend solely on the SD.
 
 This makes chain mapping a tool of choice for describing systems coupled to environment with highly structured SD (e.g. experimentally measured or calculated *ab initio*)[^chin_role_2013][^alvertis_nonequilibrium_2019][^dunnett_influence_2021][^caycedosoler_exact_2022].
 In this new representation, the Hamiltonian in Eq.(5) has naturally a 1D chain topology.
@@ -128,7 +128,7 @@ and it can be now used as pure an initial state for both of the environments. Mo
 Therefore, solving the dynamics given by the original Hamiltonian $\hat H$, starting from the initial condition $\hat \rho_S(0) \otimes \hat \rho_E(\beta)$, is equivalent to solving the dynamics given by the following Hamiltonian:
 ```math
        \hat H = \hat H_S +\hat H_E +\hat H_I  = \\
-       = \overbrace{\hat A_S}^{\hat H_S} + \overbrace{\sum_k \omega_k \big(\hat a_{1k}^\dagger \hat a_{1k} - \hat a_{2k}^\dagger \hat a_{2k} \big)}^{H_E} + \overbrace{\hat L_S \otimes \sum_k g_{1k}(\hat a_{1k}^\dagger + \hat a_{1k})+\hat L_S \otimes \sum_k g_{2k}(\hat a_{2k}^\dagger + \hat a_{2k})}^{H_I}.
+       = \overbrace{\hat A_S}^{\hat H_S} + \overbrace{\sum_k \omega_k \big(\hat a_{1k}^\dagger \hat a_{1k} - \hat a_{2k}^\dagger \hat a_{2k} \big)}^{\hat H_E} + \overbrace{\hat L_S \otimes \sum_k g_{1k}(\hat a_{1k}^\dagger + \hat a_{1k})+\hat L_S \otimes \sum_k g_{2k}(\hat a_{2k}^\dagger + \hat a_{2k})}^{\hat H_I}.
 ```
 where $\hat L_S = \hat L_S^\dagger$, considering $\hat \rho_S(0) \otimes \ket{\Omega}\bra{\Omega}$ as the initial state of system and environment. It is this Hamiltonian that is mapped on two chains with TEDOPA, to be able to perform the time evolution using tensor network techniques. 
 
@@ -145,7 +145,7 @@ This recurrence relation can be used to construct the polynomials with the condi
 
 The energy of the chain mode $n$ is given by $\varepsilon_n = A_n C_n^{-1}$ and $t_n=C_n^{-1}$ is the coupling between mode $n$ and $n+1$[^chin_exact_2010].
 
-The system couples _only_ to the first mode with the coupling strength $\kappa = ||p_0||$.
+The system couples _only_ to the first mode with the coupling strength $c_0 = ||p_0||$.
 
 Explain that for some weight function/SD they are known analytically and that for others we can use the build-in routines inspired by Gautschi or the PolyChaos.jl package.
 

docs/src/user-guide.md

@@ -2,7 +2,7 @@
 
 Here we explain the different steps to perform a simulation.
 
-Examples with detailed explanations can be found in [Examples](@ref).
+Examples with detailed explanations can be found in the `Examples` section.
 
 ## Initial State
 
@@ -49,7 +49,7 @@ In the context of Open Quantum Systems, custom chain coefficients for the enviro
 
 ### Build-in Hamiltonian
 
-MPSDynamics provides several topical Hamiltonians directly in the form of MPO or Tree Tensor Networks such as the Ising model [`MPSDynamics.isingmpo`](@ref), the XYZ Hamiltonian [`MPSDynamics.xyzmpo`](@ref), the Spin Boson Model [`MPSDynamics.spinbosonmpo`](@ref), a spin coupled to two bosonic baths `MPSDynamics.twobathspinmpo`, nearest neighbour interactions Hamiltonian `MPSDynamics.nearestneighbourmpo`, the independent boson model `MPSDynamics.ibmmpo`, (non-)uniform tight-binding chain Hamiltonian [`MPSDynamics.tightbindingmpo`](@ref).
+MPSDynamics provides several topical Hamiltonians directly in the form of MPO or Tree Tensor Networks such as the Ising model [`MPSDynamics.isingmpo`](@ref), the XYZ Hamiltonian [`MPSDynamics.xyzmpo`](@ref), the Spin Boson Model [`MPSDynamics.spinbosonmpo`](@ref), a spin coupled to two bosonic baths [`MPSDynamics.twobathspinmpo`](@ref), nearest neighbour interactions Hamiltonian [`MPSDynamics.nearestneighbourmpo`](@ref), the independent boson model [`MPSDynamics.ibmmpo`](@ref), (non-)uniform tight-binding chain Hamiltonian [`MPSDynamics.tightbindingmpo`](@ref).
 
 ### Convert a MPO from ITensor
 
@@ -64,14 +64,14 @@ MPOs are fundamentally a lists of rank-4 tensors such that the right bond dimens
 
 In the T-TEDOPA framework (see [Theoretical Background](@ref) for more details) finite temperature simulations are done with an effective pure state description of the system and the environment where the coupling coefficients (or the SD) is temperature-dependent.
 
-The corresponding chain coefficients for an Ohmic or a user provided spectral density (that can thus in pratice be either at zero or finite temperature) are computed with the `[`chaincoeffs_finiteT`](@ref)`.
+The corresponding chain coefficients for an Ohmic or a user provided spectral density (that can thus in pratice be either at zero or finite temperature) are computed with the [`chaincoeffs_finiteT`](@ref).
 This method is based on the `ORTHOPOL` routines[^Gautschi]
 
 ## Observables
 
 System and environment observables can be computed, as well as system-and-environment 'non-local' observables.
 
-Observables that will be passed to `MPSDynamics.runsim`(@ref) to have their expectation value computated at each time step are defined with the [`OneSiteObservable`](@ref) and `[`TwoSiteObservable`](@ref)`.
+Observables that will be passed to [`MPSDynamics.runsim`](@ref) to have their expectation value computated at each time step are defined with the [`OneSiteObservable`](@ref) and [`TwoSiteObservable`](@ref).
 
 One-site and two-site obsevables work similarly, they need to be given a name, an (pair of) operator(s) and the (list of) site(s) on which they are evaluated.
 
@@ -141,7 +141,7 @@ Ahead of any time evolution, one computes new bond dimensions of the MPS if the
 Then, TDVP1 is performed, using projectors with sub-spaces expanded accordingly, to produce an MPS evolved by one time step with the new,increased, bond dimensions.
 This version is faster than TDVP1 due to the acceleration gained from having more optimised bond dimensions; the bond update step is cheap and so its cost should not normally outweigh this advantage.
 However, two things need to be noticed:
-* DTDVP can get stuck into local minima from the initial time. If the bond dimension does not change during the time-evolution, consider embeding your MPS in a larger manifold with [`mpsembed!`](@ref) before time-evolving it.
+* DTDVP can get stuck into local minima from the initial time. If the bond dimension does not change during the time-evolution, consider embeding your MPS in a larger manifold with [`MPSDynamics.mpsembed!`](@ref) before time-evolving it.
 * For models with analytical solutions that are MPS of a given bond dimension (such as the independent boson model), DTDVP can overshoot the analytical bond dimension because the relative rate of change of the projection error becomes dominated by random numerical fluctuations.
 
 The convergence parameter is a threshold value `p` for the rate of change of the projection error with respect to the bond dimension.

src/MPSDynamics.jl

@@ -25,7 +25,6 @@ include("run_all.jl")
 include("run_1TDVP.jl")
 include("run_2TDVP.jl")
 include("run_DTDVP.jl")
-include("run_A1TDVP.jl")
 include("chainA1TDVP.jl")
 include("switchmpo.jl")
 include("finitetemperature.jl")

src/finitetemperature.jl

@@ -138,8 +138,6 @@ Generate chain coefficients ``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]`` for a fermio
 * Mmax: maximum number of integration points
 * save: if true the coefficients are saved
 """
-
-
 function 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)
 
     N = nummodes # Number of bath modes
@@ -209,4 +207,4 @@ function chaincoeffs_fermionic(nummodes, β, chain; ϵ=nothing, J=nothing, ωc=1
     end
 
     return [jacerg[:,1], jacerg[1:N-1,2],jacerg[N,2]]
-end
+end

src/fundamentals.jl

@@ -3,7 +3,15 @@ default(size = (800,600), reuse = true)
 crea(d) = diagm(-1 => [sqrt(i) for i=1:d-1])
 anih(d) = Matrix(crea(d)')
 numb(d) = crea(d)*anih(d)
+"""
+    function disp(d)
+
+"""
 disp(d) = (1/sqrt(2))*(crea(d)+anih(d))
+"""
+    function disp(d,ωvib,m)
+
+"""
 disp(d,ωvib,m) = (1/(2*sqrt(m*ωvib/2)))*(crea(d)+anih(d))
 mome(d) = (1/sqrt(2))*(im*(crea(d)-anih(d)))
 sx = [0. 1.; 1. 0.]
@@ -75,11 +83,10 @@ function choose(x, n)
 end
 
 """
-function therHam(psi, site1, site2)
+     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
+Calculates Hβ such that ρ = e^(-βH) for some density matrix ρ obatined from tracing out everything outside the range [site1,site2] in the MPS psi
 """
-
 function therHam(psi, site1, site2)
     pmat = ptracemps(psi, site1, site2)
     pmat = 0.5 * (pmat + pmat')
@@ -199,7 +206,6 @@ end
     Hamiltonian into the chain representation. 
 
 """
-
 function measurecorrs(oper, e::Vector, t::Vector)
     N = size(oper)[1]
     hmat = diagm(0=>e[1:N], 1=>t[1:N-1], -1=>t[1:N-1])
@@ -222,7 +228,6 @@ end
     # Returns
     - `Float64`: The result of the modified cosine function.
 """
-
 function cosineh(omega, bet)
     return 1/sqrt(1 - exp(-omega * (bet)))
 end
@@ -240,7 +245,6 @@ end
     # Returns
     - `Float64`: The result of the modified cosine function.
 """
-
 function sineh(omega, bet)
     return 1/sqrt(-1 + exp(omega * float(bet)))
 end
@@ -263,7 +267,6 @@ end
     # Returns
     - `Vector{Float64}`: The physical occupation values.
 """
-
 function physical_occup(corr_constr, corr_destr, omega, occup, b, M)
     x = range(-1, stop=1, length=M)
 

src/measure.jl

@@ -24,7 +24,6 @@ Computes the local expectation value of the one-site operator `op` on the every
 one-site observables that are obs and convobs parameters for the `runsim` function.
 
 """
-
 OneSiteObservable(name, op) = OneSiteObservable(name, op, nothing, ishermitian(op), true)
 
 struct TwoSiteObservable <: Observable
@@ -57,7 +56,6 @@ Computes the local expectation value of operators `op1` and `op2` where `op1` ac
 several-site observables that are obs and convobs parameters for the `runsim` function.
 
 """
-
 function TwoSiteObservable(name, op1, op2, sites1=nothing, sites2=nothing)
     return TwoSiteObservable(name, op1, op2, sites1, sites2, sites1==nothing && sites2==nothing)
 end
@@ -100,7 +98,6 @@ For a list of tensors `A` representing a right orthonormalized MPS, compute the
 value of the MPO M on every site.
 
 """
-
 function measurempo(A::Vector, M::Vector)
     N = length(M)
     N == length(A) || throw(ArgumentError("MPO has $N site while MPS has $(length(A)) sites"))
@@ -117,7 +114,6 @@ For a list of tensors `A` representing a right orthonormalized MPS, compute the
 value of the MPO M on specified sites.
 
 """
-
 function measurempo(A::Vector, M::Vector, sites::Tuple{Int, Int})
     N = sites[2] - sites[1] + 1
     F = fill!(similar(M[1], (1,1,1)), 1)
@@ -219,7 +215,6 @@ For a list of tensors `A` representing a right orthonormalized MPS, compute the
 value of a one-site operator O for a chainsection.
 
 """
-
 function measure1siteoperator(A::Vector, O, chainsection::Tuple{Int64,Int64})
     ρ = ones(ComplexF64, 1, 1)
 
@@ -276,7 +271,6 @@ For a list of tensors `A` representing a right orthonormalized MPS, compute the
 value of a one-site operator O for a single site.
 
 """
-
 function measure1siteoperator(A::Vector, O, site::Int)
     ρ = ones(ComplexF64, 1, 1)
     T = ishermitian(O) ? Float64 : ComplexF64
@@ -473,7 +467,6 @@ end
 Caculate expectation of M1*M2 where M1 acts on sites1 and M2 acts on sites2, assumes A is right normalised.
 
 """
-
 function measure2siteoperator(A::Vector, M1, M2, sites1::Vector{Int}, sites2::Vector{Int})
     if size(M1) == size(M2)
         herm_cis = ishermitian(M1*M2)
@@ -993,8 +986,6 @@ end
 Caculate expectation of Os on MPS A.
 
 """
-
-
 function measure(A::Vector, Os::Vector; kwargs...)
     numobs = length(Os)
     numobs==0 && return Any[]
@@ -1034,7 +1025,6 @@ end
 Caculate the reduced density matrix of the MPS A at the specified site.
 
 """
-
 function rhoreduced_1site(A::Vector, site::Int=1)
     N = length(A)
     ρR = Vector{Any}(undef, N-site+1)
@@ -1060,7 +1050,6 @@ Caculate the reduced density matrix of the MPS A of two neigbour sites. The resu
 corresponding to the dimensions  of the two sites
 
 """
-
 function rhoreduced_2sites(A::Vector, sites::Tuple{Int, Int})
     N = length(A)
     site1, site2=sites

src/models.jl

@@ -871,7 +871,6 @@ end
 
     The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``.
 """
-
 function interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
 
     e1 = chainparams1[1]

src/mpsBasics.jl

@@ -562,7 +562,6 @@ end
 Reverse the left and right dimensions of the MPS A. The resulting MPS A is the reversed MPS.
 
 """
-
 function reversemps!(A)
     N = length(A)
     reverse!(A)
@@ -577,7 +576,6 @@ end
 Reverse the left and right bond-dimensions of the MPS A.
 
 """
-
 function reversemps(A)
     N = length(A)
     Ar = Vector{Any}(undef, N)
@@ -593,7 +591,6 @@ end
 Reverse the left and right dimensions of the MPO M. The resulting MPO M is the reversed MPO.
 
 """
-
 function reversempo!(M)
     N = length(M)
     reverse!(M)
@@ -608,7 +605,6 @@ end
 Reverse the left and right dimensions of the MPO M.
 
 """
-
 function reversempo(M)
     N = length(M)
     Mr = Vector{Any}(undef, N)
@@ -697,10 +693,9 @@ function multiply(M1::Vector, M2::Vector)
 end
 
 """
-    mpsembed(A::Vector, Dmax::Int)
+    mpsembed!(A::Vector, Dmax::Int)
 
 Embed MPS `A` in manifold of max bond-dimension `Dmax`
-
 """
 function mpsembed!(A::Vector, Dmax::Int)
 

src/observables.jl

@@ -23,7 +23,6 @@ Base.ndims(::FockError) = 0
 Return the measure of the observable obs on the MPS A. 
 
 """
-
 function measure(A, obs::FockError; t=0, kwargs...)
     d = physdims(A)[obs.sites[1]:obs.sites[2]][1]
     all(x->x==d, physdims(A)[obs.sites[1]:obs.sites[2]]) || error("MPS has non-uniform local Hilbert space dimensions")

src/reshape.jl

@@ -14,7 +14,6 @@ julia> size(B) == (4,3)
 
 ```
 """
-
 function Base.reshape(x::Number, dims...)
     prod(dims) == 1 || throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size 1"))
     return fill(x, dims...)

src/treeMeasure.jl

@@ -4,7 +4,6 @@
 For a Tree, compute the local expectation value of a one-site operator O for site id.
 
 """
-
 function measure1siteoperator(net::TreeNetwork, O, id::Int)
     ρ = ones(ComplexF64, 1, 1)
     T = ishermitian(O) ? Float64 : ComplexF64
@@ -22,7 +21,6 @@ end
 For a Tree, compute the local expectation value of a one-site operator O for the specified site range.
 
 """
-
 function measure1siteoperator(net::TreeNetwork, O, sites::Tuple{Int,Int})
     ρ = ones(ComplexF64, 1, 1)
     firstsite = sites[1]
@@ -65,7 +63,6 @@ end
 For a Tree, compute the local expectation value of two one-site operators O1 and O2 for the specified site range.
 
 """
-
 function measure2siteoperator(net::TreeNetwork, O1, O2, sites::Tuple{Int,Int})
 
     herm_cis = ishermitian(O1*O2)

src/treeTDVP.jl

@@ -27,7 +27,6 @@ end
 
 Return the physical dimensions of a tree-MPS or tree-MPO `M`.
 """
-
 function physdims(M::TreeNetwork)
     N = length(M)
     res = Vector{Int}(undef, N)
@@ -275,7 +274,6 @@ tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F=nothing; 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.  
 """
-
 function tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F::Vector, id::Int; verbose=false, kwargs...)
 
     children = A.tree[id].children
@@ -565,7 +563,6 @@ 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>
 ```
 """
-
 function productstatemps(tree_::Tree, physdims::Dims, Dmax::Int=1; state=:Vacuum)
     tree = deepcopy(tree_)
     hn = findheadnode(tree)
@@ -613,7 +610,6 @@ productstatemps(tree::Tree, physdims::Int, Dmax::Int; state=:Vacuum) =
 Embed tree-MPS `A` in manifold of max bond-dimension `Dmax`.
 
 """
-
 function mpsembed!(A::TreeNetwork, Dmax::Int)
     tree = deepcopy(A.tree)
     pdims = physdims(A)
@@ -645,7 +641,6 @@ end
 Return the bon-dimension of a tree-MPS `A`.
 
 """
-
 function bonddims(A::TreeNetwork)
     N = length(A)
     mat = zeros(Int, N, N)