Updated doc with convergence check page

docs/make.jl

@@ -13,6 +13,7 @@ makedocs(
         "user-guide.md",
         "nutshell.md",
         "Examples" => ["./examples/sbm.md", "./examples/puredephasing.md", "./examples/timedep.md", "./examples/anderson-model.md", "./examples/bath-observables.md", "./examples/protontransfer.md"],
+        "convergence.md",
         "theory.md",
         "Methods" => "methods.md",
         "dev.md"

docs/src/convergence.md

@@ -0,0 +1,65 @@
+# Convergence checks
+
+In this section, we analyze the numerical accuracy of our simulations by outlining the necessary approximations and providing practical guidelines to ensure convergence of results.
+
+## Approximations
+
+### Chain length $N$
+In numerical simulations, a truncation on the number of chain modes (and therefore the chain length) is introduced to handle a finite chain instead of a semi-infinite one. This truncation, denoted as $N$, corresponds to a sampling of the modes in the original environment.
+
+The chain length $N$ is directly connected to the simulation time as follows: excitations injected into the chain propagate as a wavefront traveling along the chain. Perturbations to the initial state outside this wavefront are exponentially suppressed. Therefore, truncating the chain basis beyond the expanding wavefront ensures that the resulting sampling error is also exponentially small  [^woods_simulating_2015][^DeVega_howto_2015].
+
+To optimize the chain length for a given simulation time and set of chain coefficients, the following procedure can be used:
+
+1. Set the total simulation time `T` and compute the chain coefficients for an intentionally oversized chain `N_huge` (longer than needed):
+```
+T = 100 
+N_huge = 1000
+cpars = chaincoeffs_flat(N_huge, αchain; ωc = 1.0)
+```
+2. Use the built-in function [`findchainlength`](@ref) to determine the optimal chain length `N_opt` based on the propagation speed of the wavefront (given by the hopping coefficients $t_n$):
+```
+N_opt = findchainlength(T, cpars; eps=10^-4, verbose=false)
+```
+2bis. A rule-of-thumb estimate `N_est` can also be obtained with [`findchainlength`](@ref) using universal asymptotic properties of chain mapped environment:
+```
+N_est = findchainlength(T, ωc, β)
+```
+
+### Local dimension $d$
+Another critical truncation is imposed on the local dimension $d$ of each tensor in the MPS that undergoes dynamic evolution. The local dimension $d$ corresponds to the number of Fock states retained in the Hilbert space of each chain mode. Since harmonic oscillators are, in principle, infinite-dimensional systems, truncating their Hilbert space to a finite $d$ is necessary for numerical computations. The choice of $d$ determines the maximum number of excitations per site in the MPS and must be carefully tuned to capture the relevant physics while ensuring numerical convergence.
+For a bath that initially contains only a finite number of particles, i.e. any bath in practice, the error originating from the local Hilbert space truncation vanishes exponentially as $d$ increases[^woods_simulating_2015]. 
+
+A practical guideline for choosing $d$ can be derived by considering the physical states being simulated. For example, for coherent states, the occupation number follows a Poisson distribution. In such cases, the average occupation number is given by the mean $\langle n \rangle$, and the probability of observing a state with occupation number $n$ decreases exponentially for $n \gg \langle n \rangle$. Thus, the truncation $d$ should be chosen such that the cumulative probability of truncation error is negligible.
+
+### Bond dimensions and time evolution
+Different time-evolution algorithms are implemented in MPSDynamics, and can be selected as options in the function [`runsim`](@ref). Detailed descriptions of these algorithms are provided in the [Theoretical Background](@ref) section of this documentation. The key parameter for ensuring convergence in these methods is controlled via the `convparams` option, which accepts arrays of multiple parameter values to test until convergence is achieved.
+
+Below are the main time-evolution algorithms and their corresponding convergence parameters:
+- One-Site Time-Dependent Variational Principle (1TDVP):
+  - Selected using: `method = :TDVP1`.
+  - The convergence parameter is the bond dimension of the MPS, which defines the manifold on which the time evolution is constrained. To ensure convergence, multiple bond dimensions can be tested: `convparams = [bond1, bond2, bond3]`.
+- Two-Site Time-Dependent Variational Principle (2TDVP):
+  - Selected using: `method = :TDVP2`.
+  - The convergence parameter is the SVD singular value truncation threshold, which determines the precision of the truncation during simulations. Multiple truncation thresholds can be tested: `convparams = [trunc1, trunc2, trunc3]`.
+  - The bond dimensions at each chain-site and at each time-step can be saved by specifying `savebonddims = true`. The maximum allowed value of bond dimension can be chosen by setting the `Dlim` option.
+- Adaptive Time-Dependent Variational Principle (DTDVP):
+  - Selected using: `method = :DTDVP`.
+  - The convergence parameter is the target precision for the adaptive algorithm. Multiple precision values can be provided to test convergence: `convparams = [prec1, prec2, prec3]`.
+  - The bond dimensions at each chain-site and at each time-step can be saved by specifying `savebonddims = true`. The maximum allowed value of bond dimension can be chosen with the `Dlim` option.
+  - In some cases, it might be better to embed the initial time MPS in a manifold of higher bond-dimension, which can be done with the function `mpsembed!`, to avoid for the dynamics to get stuck.
+
+### Time-step 
+The time step used in simulations must be small enough to accurately capture the dynamics of the system. To ensure convergence, gradually reduce the time step until key observables (e.g., energy, population dynamics, or correlation functions) stabilize. It is often practical to test a range of time steps and assess their impact on results to determine an optimal balance between accuracy and computational cost.
+
+## Common pitfalls
+
+### Hard cut-off in the spectral density function
+The cutoff frequency $\omega_c$ in the spectral density function has to be selected in order to incorporate all of the relevant frequencies. It is to be noted that, for spectral densities belonging to the Szegö class[^chin_exact_2010], the cutoff frenquency $\omega_c$ also determines the asymptotic values to which the chain coefficients converge [^woods_mappings_2014]. Therefore, if we want to reach translational invariance (converged values) in the chain coefficients, then we need to specify a cutoff. To avoid ringing effects, it is a good choice to select a spectral density function that approaches zero before the cutoff frequency.
+
+
+## References
+[^woods_simulating_2015]: Woods, M. P.; Cramer, M.; Plenio, M. B. Simulating Bosonic Baths with Error Bars. Phys. Rev. Lett. 2015, 115 (13), 130401. https://doi.org/10.1103/PhysRevLett.115.130401.
+[^DeVega_howto_2015]: De Vega, I.; Schollwöck, U.; Wolf, F. A. How to Discretize a Quantum Bath for Real-Time Evolution. Phys. Rev. B 2015, 92 (15), 155126. https://doi.org/10.1103/PhysRevB.92.155126.
+[^woods_mappings_2014]: Woods,  M. P. et al. “Mappings of open quantum systems onto chain representations and Markovian embeddings”. In: Journal of Mathematical Physics 55.3 (Mar. 2014), p. 032101. issn: 0022-2488, 1089-7658. doi: 10.1063/1.4866769
+[^chin_exact_2010]: Chin, A. W.  et al. “Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials”. In: Journal of Mathematical Physics 51.9 (Sept. 2010), p. 092109. issn: 0022-2488, 1089-7658. doi: 10.1063/1.3490188. arXiv: 1006.4507

docs/src/examples/anderson-model.md

@@ -275,6 +275,8 @@ plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
 
 ```
 
+![Output of this program](double_output.svg)
+
 
 #### Interleaved chain geometry
 The Hamiltonian is defined using the [`MPSDynamics.interleaved_tightbinding_mpo`](@ref) function, which takes as an input the number of modes of each chain `N`, the defect's energy `ϵd`, and the chain coefficients of the first `chainparams1` and second `chainparams2` chain. The MPS for the initial state is a factorized state (bond dimension 1) made of: a filled impurity, and 2N alternate filled-empty states. 
@@ -406,6 +408,8 @@ end
 plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
 ```
 
+![Output of this program](interleaved_output.svg)
+
 ________________
 ## Bibliography
 

docs/src/examples/bath-observables.md

@@ -21,8 +21,8 @@ s = 1       # ohmicity
 ```
 We set the specifics of the simulation:
 ```julia
-method = :TDVP1         # time-evolution method
-conv = 3                # bond dimension for the TDVP1
+method = :TDVP2         # time-evolution method
+conv = 0.001            # Allowed SVD singular value truncation
 dt = 0.5                # time step
 tfinal = 60.0           # simulation time
 ```
@@ -62,6 +62,8 @@ A, dat = runsim(dt, tfinal, A, H, prec=1E-4;
                 convparams = conv,
                 reduceddensity = true,
                 verbose = false,
+                savebonddims = true,
+                Dlim = 100,
                 save = false,
                 plot = true,
                 );
@@ -193,6 +195,8 @@ p4 = plot(omeg[Mhalf:M], bath_occup_phys, lw=4,
 plot(p1, p2, p3, p4, layout = (2, 2), size = (1400, 1200))
 ```
 
+![Output of this program](bath_obs_output.svg)
+
 ___________________
 ## Bibliography
 

docs/src/index.md

@@ -17,7 +17,7 @@ The package may be installed by typing the following into a Julia REPL
 ## Table of Contents
 
 ```@contents
-Pages = ["index.md", "user-guide.md", "examples/sbm.md", "examples/puredephasing.md", "examples/timedep.md", "examples/anderson-model.md", "examples/bath-observables.md", "examples/protontransfer.md", "theory.md", "methods.md", "dev.md"]
+Pages = ["index.md", "user-guide.md", "examples/sbm.md", "examples/puredephasing.md", "examples/timedep.md", "examples/anderson-model.md", "examples/bath-observables.md", "examples/protontransfer.md", "convergence.md", "theory.md", "methods.md", "dev.md"]
 Depth = 3
 ```
 

examples/anderson_model_double.jl

@@ -19,12 +19,15 @@ function ϵ(x)
     return x
 end
 
-function V(x)
+function J(x)
     return sqrt(1 - x^2) # semi-circular density of states
 end
 
-chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, V, save=false) # empty
-chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, V, save=false) # filled
+chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, J, save=false) # empty
+chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, J, save=false) # filled
+
+c1 = chainparams1[3]
+c2 = chainparams2[3]  
 
 #=
 curdir = @__DIR__
@@ -71,7 +74,7 @@ A, dat = runsim(dt, T, A, H;
                 method = method,
                 obs = [ob1, ob2, ob3], 
                 convobs = [ob1],
-                params = @LogParams(N, ϵd, β),
+                params = @LogParams(N, ϵd, β, c1, c2),
                 convparams = [prec],   
                 Dlim = Dmax,          
                 savebonddims = true,   # we want to save the bond dimension

examples/anderson_model_interleaved.jl

@@ -19,12 +19,15 @@ function ϵ(x)
     return x
 end
 
-function V(x)
+function J(x)
     return sqrt(1 - x^2) # semi-circular density of states
 end
 
-chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, V, save=false) # empty
-chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, V, save=false) # filled
+chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, J, save=false) # empty
+chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, J, save=false) # filled
+
+c1 = chainparams1[3]
+c2 = chainparams2[3]  
 
 #=
 curdir = @__DIR__
@@ -77,7 +80,7 @@ A, dat = runsim(dt, T, A, H;
                 method = method,
                 obs = [ob1, ob2], 
                 convobs = [ob1],
-                params = @LogParams(N, ϵd, β),
+                params = @LogParams(N, ϵd, β, c1, c2),
                 convparams = [prec],   
                 Dlim = Dmax,          
                 savebonddims = true,   # we want to save the bond dimension

src/finitetemperature.jl

@@ -118,7 +118,7 @@ function chaincoeffs_finiteT(nummodes, β, ohmic=true; α=1, s=1, J=nothing, ωc
 end
 
 """
-    chaincoeffs_fermionic(nummodes, β, chain; ϵ=x, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true)
+ 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 β.