Doc writing

.github/workflows/docs.yml

@@ -3,7 +3,7 @@ name: Documentation
 on:
   push:
     branches:
-      - master
+      - doc-writing
     tags: '*'
   pull_request:
 

docs/make.jl

@@ -7,11 +7,21 @@ makedocs(
     modules = [MPSDynamics],
     highlightsig = true,
     sitename = "MPSDynamics.jl",
-    authors = "Angus Dunnett, Thibaut Lacroix, Brieuc Le Dé",
-    expandfirst = []
+    authors = "Angus Dunnett, Thibaut Lacroix, Brieuc Le Dé, Angela Riva",
+    pages = [
+        "index.md",
+        "user-guide.md",
+        "Examples" => ["./examples/sbm.md", "./examples/puredephasing.md"],
+        "theory.md",
+        "Methods" => "methods.md",
+        "dev.md"
+    ],
+    expandfirst = [],
+    remotes = nothing
 )
 
 deploydocs(
     repo = "github.com/shareloqs/MPSDynamics.git",
-    devurl = "docs"
+    devurl = "docs",
+    devbranch = "doc-writing"
 )

docs/src/dev.md

@@ -0,0 +1,15 @@
+# Developpers
+
+## Simulation Workflow
+
+The flow chart will go here
+
+## How to Contribute
+
+Contributions are welcome! Don't hesitate to contact us if you
+
+* found a bug;
+* have a suggestion on how to improve the code and/or documentation;
+* would like to get involved in writing code and/or documentation.
+
+You can contact us by raising an [issue on Github](https://github.com/shareloqs/MPSDynamics/issues), or by writing to one of the developpers.

docs/src/examples/puredephasing.md

@@ -0,0 +1,5 @@
+# Pure-Dephasing
+
+## Zero Temperature
+
+## Finite Temperature

docs/src/examples/sbm.md

@@ -0,0 +1 @@
+# The Spin-Boson Model

docs/src/index.md

@@ -6,10 +6,8 @@ The package also provides strong support for the measurement of observables, as
 The package has been developed with the aim of studying non-Markovian open system dynamics at finite temperature using the state-of-the-art numerically exact Thermalized-Time Evolving Density operator with Orthonormal Polynomials Algorithm (T-TEDOPA) based on environment chain mapping.
 However the methods implemented can equally be applied to other areas of physics.
 
-## Table of Contents
-
-```@contents
-```
+!!! warning
+	The documentation is currently undergoing massive restructurations/improvement. It's a work in progress until the next release scheduled for May, 2024.
 
 ## Installation
 
@@ -19,8 +17,26 @@ The package may be installed by typing the following into a Julia REPL
     ] add https://github.com/shareloqs/MPSDynamics.git
 ```
 
-## Functions
+## Table of Contents
+
+```@contents
+Pages = ["index.md", "user-guide.md", "./examples/sbm.md", "./examples/puredephasing.md", "theory.md", "methods.md", "dev.md"]
+Depth = 3
+```
 
-```@autodocs
-Modules = [MPSDynamics]
+## Citation
+If you use the package in your research, please consider citing it.
+You can add the Zenodo record to your BibTex file:
+
+```tex
+@misc{mpsdynamics_zenodo2021,
+	title = {shareloqs/{MPSDynamics}},
+	shorttitle = {{MPSDynamics.jl}},
+	url = {https://zenodo.org/record/5106435},
+	abstract = {Tensor network simulations for finite temperature, open quantum system dynamics},
+	publisher = {Zenodo},
+	author = {Dunnett, Angus and Lacroix, Thibaut and Le Dé, Brieuc and Riva, Angela},
+	year = {2021},
+	doi = {10.5281/zenodo.5106435},
+}
 ```

docs/src/methods.md

@@ -0,0 +1,5 @@
+# List of all methods
+
+```@autodocs
+Modules = [MPSDynamics]
+```

docs/src/user-guide.md

@@ -0,0 +1,188 @@
+# User Guide
+
+Here we explain the different steps to perform a simulation.
+
+Examples with detailed explanations can be found in [Examples](@ref).
+
+## Initial State
+
+The initial many-body state of the {System + Environment} can be described easily as [Matrix Product States](@ref) (MPS) or [Tree Tensor Networks](@ref) (TTN) state (e.g. when a system is coupled to several environments).
+
+### Matrix Product States
+
+A MPS can be initialized with several methods.
+
+The [`productstatemps`](@ref) method enables  to instantiate arbitrary MPS of fixed uniform bond dimension with non-uniform physical dimensions.
+The individual 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. 
+Setting `state=:FullOccupy` will produce an MPS in which each site is fully occupied.
+The gauge of the MPS can also be set using a keyword argument.
+
+```julia
+julia> ψ = unitcol(1,2) # system initial state
+
+julia> d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1 
+
+julia> cpars = chaincoeffs_ohmic(N, α, s) # chain coefficient for an Ohmic spectral density
+
+julia> H = spinbosonmpo(ω0, Δ, d, N, cpars)
+
+julia> A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
+```
+
+Alternatively, a chain with a specified number of excitation localised on one site, or delocalized accross several sites can be generated with [`MPSDynamics.chainmps`](@ref).
+
+Random MPS can also be generated with the [`randmps`](@ref) method.
+
+For the case of fermionic states (which need to be anti-symmetrized), the [`MPSDynamics.electronkmps`](@ref) method generate an MPS for an electron with momentum `k`, and the [`MPSDynamics.electron2kmps`](@ref) generate an MPS with 2 electrons in k-states `k1` and `k2`.
+
+### Tree Tensor Networks
+
+Write a quick explanation of the how trees are structured: parents, child, nodes; and that most methods to initialize a MPS are overloaded for `TreeNetwork`. 
+
+## Hamiltonian
+
+In order to perform time evolution and have access to the dynamics of the many-body state a Hamiltonian needs to be specified in the form of a Matrix Product Operator (MPO) of as a tree tensor network.
+Either way, this can be done by using a [Build-in Hamiltonian](@ref), [Convert a MPO from ITensor](@ref), or creating a [Tailored MPO](@ref).
+
+In the context of Open Quantum Systems, custom chain coefficients for the environment can be generated for finite temperature simulations, and/or user provided spectral densities (SDs).
+
+### 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).
+
+### Convert a MPO from ITensor
+
+The method [`MPSDynamics.MPOtoVector`](@ref) converts an ITensors chain MPO into a form compatible with MPSDynamics.
+
+### Tailored MPO
+
+One can also construct MPO tailored to the problem one is interested in.
+MPOs are fundamentally a lists of rank-4 tensors such that the right bond dimension of the nth tensor must be equal to the left bond dimension of the n+1th tensor; and the dimension of the physical bonds of the nth tensor must be equal to the corresponding physical bond on the MPS.
+
+### Finite Temperature and Custom SD
+
+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)`.
+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)`.
+
+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.
+
+For instance one can calculated the average number of excitation in each of the ``N`` environmental modes
+
+```julia
+    ob = OneSiteObservable("chain mode occupation", numb(d), (2,N+1))
+```
+
+It is also possible to measure composite system/environment observables, for example
+
+```julia
+    ob = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1))
+```
+which measure the correlation between the spin in the x-direction and the displacement of the bath modes.
+
+Purely environmental 'non-local' observables such as 
+
+```julia
+    ob = TwoSiteObservable("bath coherence", crea(d), anih(d), collect(2:N+1), collect(2:N+1))
+```
+that computes all the chain mode coherences ``\langle\hat{a}_n^\dagger\hat{a}_{m}\rangle`` (and the chain mode occupation when ``n = m``).
+
+We note that if one knows the coherences and populations of all the chain modes, it is then possible to reconstruct the populations of the normal mode environment.
+
+## Time-Evolution
+
+Explain that we do TDVP.
+
+A simulation is runned by the `MPSDynamics.runsim`(@ref) function.
+
+Local one-site and two-site observables, as well as non-local two-site observables, can be efficiently computed for each time-step of the time-evolution.
+
+
+### One-site TDVP
+
+Fixed bond dimension, complexity XXX, preserves unitarity.
+
+The convergence parameter is the bond dimension.
+
+Is set in `MPSDynamics.runsim`(@ref) usind the key word argument `method=:TDVP1`
+
+### Two-Site TDVP
+
+Varying bond dimension, complexity higher than 1TDVP, but breaks unitarity because of a SVD.
+
+The convergence parameter is the threshold of the SVD.
+
+Is set in `MPSDynamics.runsim`(@ref) usind the key word argument `method=:TDVP2`
+
+### Adaptive One-Site TDVP
+
+Keeps the good scaling of the 1TDVP, preserves unitarity but is able to increase the bond dimendion.
+
+The convergence parameter is a threshold value for the rate of change of the projection error with respect to the bond dimension.
+
+Is set in `MPSDynamics.runsim`(@ref) usind the key word argument `method=:DTDVP`
+
+## Data Storage
+
+The data (i.e. observables time-series) is stored in the JLD format which is based on HDF5.
+
+The HDF5 format is natively supported across many platforms and languages (e.g. `Python`, or `Mathematica`).
+
+For the data to be saved to a file after a run, the keyword argument `save=true` needs to be used in `runsim`.
+
+The directory where the data should be saved can be chosen by setting a path with the `savedir` keyword argument.
+
+A `plot` keyword argument can also be used to choose whether plots for 1D observables will be automatically generated and saved along with the data. 
+
+Loading the data in Julia using the [`JLD.jl`](https://github.com/JuliaIO/JLD.jl) package will recover the full type information of the Julia variables that were stored.
+
+```julia
+julia> using JLD
+
+julia> dat = load("filename.jld")
+Dict{String, Any} with 2 entries:
+  "parameters" => Dict{String, Any}("tmax"=>0.2, "method"=>:DTDVP, "dt"=>0.0005…
+  "data"       => Dict{String, Any}("bonddims"=>[1 1 … 1 1; 1 2 … 2 2; … ; 1 1 …
+```
+
+If the data is loaded in an variable `dat` the structure is the folowing:
+- the `parameters` entry contains a dictionary where all the simulation parameters are stored
+- the `data` entry contains a dictionary where are stored simulation time, the observables and (whenever relevent) the bond dimension of the state at each time steps. 
+
+```julia
+julia> dat["parameters"]
+Dict{String, Any} with 11 entries:
+  "tmax"       => 0.2
+  "method"     => :DTDVP
+  "dt"         => 0.0005
+  "name"       => "my model"
+  "Δ"          => 0.0
+  "β"          => 0.0186854
+  "N"          => 300.0
+  "d"          => 15.0
+  "unid"       => "Ovzm6"
+  "ω0"         => 0.0
+  "convparams" => 0.0005
+
+julia> dat["data"]
+Dict{String, Any} with 6 entries:
+  "bonddims" => [1 1 … 1 1; 1 2 … 2 2; … ; 1 1 … 7 7; 1 1 … 1 1]
+  "sx"       => [1.0, 0.912818, 0.741759, 0.605797, 0.528792, 0.492497, 0.47976…
+  "sz"       => [0.0, 0.0871825, 0.25824, 0.394201, 0.471207, 0.507503, 0.52023…
+  "nchain"   => [0.0 0.23466 … 1.84319 1.76098; 0.0 0.00231507 … 0.83105 0.9033…
+  "sy"       => [0.0, -0.0133489, -0.0588887, -0.0858181, -0.0759996, -0.048539…
+  "times"    => [0.0, 0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.00…
+
+```
+
+[^Gautschi]:
+	> Gautschi, W. Algorithm 726: ORTHPOL–a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans. Math. Softw. 20, 21–62 (1994).
+

src/fundamentals.jl

@@ -178,9 +178,9 @@ end
 """
     findchainlength(T, cparams...; eps=10^-6)
 
-Estimate length of chain required for a particular set of chain parameters by calulating how long an excitation on the
+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.
+reached the last site after time T. The initial number of sites in cparams has to be larger than the findchainlength result.
 
 """
 function findchainlength(T, cparams; eps=10^-4, verbose=false)

src/measure.jl

@@ -987,6 +987,14 @@ function measure(A::Vector, O::CdagCdn, ρ::Vector)
     return expval
 end
 
+"""
+     measure(A::Vector, Os::Vector; kwargs...)
+
+Caculate expectation of Os on MPS A.
+
+"""
+
+
 function measure(A::Vector, Os::Vector; kwargs...)
     numobs = length(Os)
     numobs==0 && return Any[]