Fix typos and add docs

docs/src/examples/anderson-model.md

@@ -1,16 +1,17 @@
 # The Anderson Impurity Model 
 
-Here we give some context on the examples provided in `MPSDynamics/examples/anderson_model_double.jl` and `MPSDynamics/examples/anderson_model_interleaved.jl`. In these two examples, we use the fermionic chain mapping proposed in [^khon_efficient_2021] to show how to perform tensor network simulations of the Single Impurity Anderson Model (SIAM) with the `MPSDynamics.jl` library. 
+Here we give some context on the examples provided in `MPSDynamics/examples/anderson_model_double.jl` and `MPSDynamics/examples/anderson_model_interleaved.jl`. In these two examples, we use the fermionic chain mapping proposed in [^khon_efficient_2021] to show how to perform tensor network simulations of the [Single Impurity Anderson Model](https://en.wikipedia.org/wiki/Anderson_impurity_model) (SIAM) with the `MPSDynamics.jl` library. 
 
 
 ## Basics of the fermionic chain mapping
 
 Before giving the code implementation, a short recap on the problem statement and on the fermionic mapping used to solved it. The SIAM Hamiltonian is defined as:
 
 ```math
-    \hat H^\text{SIAM}  = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{H_I^\text{chain}}.
+    \hat H^\text{SIAM}  = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{\hat H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{\hat H_\text{cond}}.
 ```
-All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation 
+
+This Hamiltonian represents a located impurity ($\hat H_\text{loc}$), conduction electrons ($\hat H_\text{cond}$) and a hybridization term between the impurity and the conduction electrons ($\hat H_\text{hyb}$). All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation 
 [^devega_thermo_2015], performed with fermions: ancillary fermionic operators $\hat c_{2k}$ are defined, one for each of the original fermionic modes $\hat c_{1k}$. A Bogoliubov transformation is then applied, so that two new fermionic modes $\hat f_{1k}$ and $\hat f_{2k}$ are defined as a linear combination of $\hat c_{1k}$ and $\hat c_{2k}$. Two chains are defined: the chain labelled $1$ for the empty modes:
 ```math
 \hat f_{1k}=e^{-iG}  \hat c_k e^{iG}= \cosh(\theta_k) \hat c_{1k} -\sinh(\theta_k)  \hat c_{2k}^\dagger \\
@@ -136,7 +137,7 @@ where the matrices are defined as:
 
 ## Code implementation
 
-The fermionic mapping can be strightforwardly implemented with the methods of `MPSDyanmics.jl`. We provide two examples: 
+The fermionic mapping can be strightforwardly implemented with the methods of `MPSDynamics.jl`. We provide two examples: 
 - `examples/anderson_model_double`, simulating the SIAM with the double-chain geometry
 - `examples/anderson_model_interleaved`, simulating the SIAM with the interleaved-chain geometry
 They only differ on the way the Hamiltonian MPO is defined. We now briefly review the important bits of the code.
@@ -149,7 +150,7 @@ N = 40      # number of chain sites
 Ed = 0.3    # energy of the impurity
 ϵd = Ed - μ # energy of the impurity minus the chemical potential
 ```
-The `chaincoeffs_fermionic` function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain `N`, the inverse temperature `β`, a label to specify if the chain modes are empty (label is `1.0`) or filled (label is `2.0`), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$.
+The [`MPSDynamics.chaincoeffs_fermionic`](@ref) function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain `N`, the inverse temperature `β`, a label to specify if the chain modes are empty (label is `1.0`) or filled (label is `2.0`), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$.
 ```julia
 function ϵ(x)
     return x
@@ -174,7 +175,7 @@ prec = 0.0001       # precision for the adaptive TDVP
 and with this we are ready to construct the Hamiltonian MPO and specify the initial state, which will obviously differ depending on the chosen geometry.
 
 #### Double chain geometry
-The Hamiltonian is defined using the `tightbinding_mpo` 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 made of: N filled states, a filled impurity, and N empty states. 
+The Hamiltonian is defined using the [`MPSDynamics.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 made of: N filled states, a filled impurity, and N empty states. 
 ```julia
 H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
 
@@ -198,7 +199,7 @@ A, dat = runsim(dt, T, A, H;
                 method = method,
                 obs = [ob1, ob2, ob3], 
                 convobs = [ob1],
-                params = @LogParams(N, ϵd, β, c1, c2),
+                params = @LogParams(N, ϵd, β),
                 convparams = [prec],   
                 Dlim = Dmax,          
                 savebonddims = true,   # we want to save the bond dimension
@@ -276,7 +277,7 @@ plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
 
 
 #### Interleaved chain geometry
-The Hamiltonian is defined using the `interleaved_tightbinding_mpo` 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. 
+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. 
 ```julia
 H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
 
@@ -306,7 +307,7 @@ A, dat = runsim(dt, T, A, H;
                 method = method,
                 obs = [ob1, ob2], 
                 convobs = [ob1],
-                params = @LogParams(N, ϵd, β, c1, c2),
+                params = @LogParams(N, ϵd, β),
                 convparams = [prec],   
                 Dlim = Dmax,          
                 savebonddims = true,   # we want to save the bond dimension

docs/src/examples/bath-observables.md

@@ -79,7 +79,7 @@ for t in 1:T
 end
 ```
 
-During a numerical simulation we work with chain of finite length $N$. This truncation on chain modes, introduces a sampling on the modes in the original star-like environment. To recover the frequency modes that are implicitly sampled, one has to diagonalize the tri-diagonal $N\times N$ matrix $H^\text{chain}$, where the diagonal is formed by the $e_n$ coefficients, that is the chain's frequencies, and the upper and lower diagonals by the $N-1$ hopping coefficients $t_n$. The unitary matrix that defines the change of basis from the star-like to the chain-like environment is $U_n$, constituted by the eigenvectors of $H^\text{chain}$. In the code, we use the `eigenchain` function:
+During a numerical simulation we work with chain of finite length $N$. This truncation on chain modes, introduces a sampling on the modes in the original star-like environment. To recover the frequency modes that are implicitly sampled, one has to diagonalize the tri-diagonal $N\times N$ matrix $H^\text{chain}$, where the diagonal is formed by the $e_n$ coefficients, that is the chain's frequencies, and the upper and lower diagonals by the $N-1$ hopping coefficients $t_n$. The unitary matrix that defines the change of basis from the star-like to the chain-like environment is $U_n$, constituted by the eigenvectors of $H^\text{chain}$. In the code, we use the [`MPSDynamics.eigenchain`](@ref) function:
 
 ```julia
 omeg = eigenchain(cpars, nummodes=N).values
@@ -94,7 +94,7 @@ Therefore to calculate the number of modes of the environment associated to a sp
 \end{aligned}
 ```
 
-This is done in the code using the `measuremodes(X, cpars[1], cpars[2])` function, which outputs the vector of the diagonal elements of the operators, in the following way:
+This is done in the code using the [`MPSDynamics.measuremodes`](@ref) function, which outputs the vector of the diagonal elements of the operators, in the following way:
 
 ```julia
 bath_occup = mapslices(X -> measuremodes(X, cpars[1], cpars[2]), dat["data/cdagc"], dims = [1,2])
@@ -126,7 +126,7 @@ It is possible to invert the thermofield transformation (details in [^riva_therm
 \end{aligned}
 ```
 
-We remark that in the thermofield case, a negative frequency $\omega_{2k}$ is associated to each positive frequency $\omega_{1k}$. The sampling is therefore symmetric around zero. This marks a difference with T-TEDOPA, where the sampling of frequencies was obtained through the thermalized measure $d\mu(\beta) = \sqrt{J(\omega, \beta)}d\omega$, and was not symmetric. To recover the results for the physical bath of frequencies starting from the results of our simulations, that were conducted using the T-TEDOPA chain mapping, we need to do an extrapolation for all of the mean values, in order to have their values for each $\omega$ at $-\omega$ as well. This is done in the code with the `physical_occup` function:
+We remark that in the thermofield case, a negative frequency $\omega_{2k}$ is associated to each positive frequency $\omega_{1k}$. The sampling is therefore symmetric around zero. This marks a difference with T-TEDOPA, where the sampling of frequencies was obtained through the thermalized measure $d\mu(\beta) = \sqrt{J(\omega, \beta)}d\omega$, and was not symmetric. To recover the results for the physical bath of frequencies starting from the results of our simulations, that were conducted using the T-TEDOPA chain mapping, we need to do an extrapolation for all of the mean values, in order to have their values for each $\omega$ at $-\omega$ as well. This is done in the code with the [`MPSDynamics.physical_occup`](@ref) function:
 
 ```julia
 bath_occup_phys = physical_occup(correlations_cdag[:,:,T], correlations_c[:,:,T], omeg, bath_occup[:,:,T], β, N)

src/deprecated/fundamentals_deprecated.jl

@@ -0,0 +1,16 @@
+# Doesn't really work because indices change order in ITensors
+function MPStoVector(mps::MPS)
+    N = length(mps)
+    A = [Array(mps[i], mps[i].inds...) for i=1:N]
+    dims = size(A[1])
+    A[1] = reshape(A[1], 1, dims...)
+    A[1] = permutedims(A[1], [1,3,2])
+    for i=2:N-1
+        A[i] = permutedims(A[i], [3,1,2])
+    end
+    dims = size(A[N])
+    A[N] = reshape(A[N], dims..., 1)
+    A[N] = permutedims(A[N], [1,3,2])
+    return A
+end
+

src/deprecated/models_deprecated.jl

@@ -0,0 +1,18 @@
+#Deprecated#
+function getchaincoeffs(nummodes, α, s, beta, ωc=1)
+    matlabdir = ENV["MATDIR"]
+    astr = paramstring(α, 2)
+    bstr = paramstring(beta, 3)
+    datfname = string("jacerg_","a",astr,"s$s","beta",bstr,".csv")
+    chaincoeffs = readdlm(string(matlabdir,datfname),',',Float64,'\n')
+    es = ωc .* chaincoeffs[:,1]
+    ts = ωc .* chaincoeffs[1:end-1,2]
+    c0 = ωc * sqrt(chaincoeffs[end,2]/pi)
+    Nmax = length(es)
+    if Nmax < nummodes
+        throw(ErrorException("no data for nummodes > $Nmax"))
+    end
+    return [es[1:nummodes], ts[1:nummodes-1], c0]
+end
+##
+