Update anderson-model.md

added the commented code to make the example more tutorial like

docs/src/examples/anderson-model.md

@@ -1,12 +1,24 @@
 # 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 perform tensor network simulations of the Single Impurity Anderson Model (SIAM). The SIAM Hamiltonian is defined as:
+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. 
+
+
+## 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}}.
 ```
 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, the chain labelled $2$ for the filled modes.
+[^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 \\
+\hat f_{2k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k)  \hat c_{1k} +\sinh(\theta_k)  \hat c_{2k}^\dagger.
+```
+We remark that this is the same Bogoliubov transformation used in the thermofield[^devega_thermo_2015] for the bosonic case: the only thing that changes is a plus sign, that takes into account the fermionic anti-commutation relations. With these new fermionic operators we obtain the thermofield-transformed Hamiltonian, where the system interacts with two environments at zero temperature, allowing for pure state simulations (and thus the employement of MPS).
+
+The thermofield-transformed Hamiltonian is then mapped on two chains, defined and constructed using the TEDOPA chain mapping: the chain labelled $1$ is for the empty modes, the chain labelled $2$ for the filled modes.
 The following relations are used to define the functions equivalent to the spectral density of the bosonic case, one for each chain:
 ```math
 \begin{split}
@@ -27,11 +39,11 @@ where the $J_{i,n}$ coefficients are the couplings between the chain sites and t
 ```
 
 
-## Double chain mapping
+### Double chain MPO
 
 ![image](doublechain.png)
 
-In the double chain geometry of Fig. \ref{subfig:double_ferm}, the MPOs bond dimension is: $\chi = 4$. The MPO has the structure defined in Eq. \ref{eq:mpo} and therefore can be seen as the product of the following matrices:
+We can represent the Hamiltonian $\hat H^\text{chain}$ as a MPOs bond dimension is: $\chi = 4$. The MPO has the structure of a double chain, with an impurty site at the center, and therefore can be seen as the product of the following matrices:
 ```math
     H = W_{1N} \cdot...\cdot W_{1 0} \cdot W_d \cdot W_{20} \cdot ... \cdot W_{2 N},  
 ```
@@ -68,9 +80,8 @@ where the matrices are defined as:
 \end{bmatrix}.
 \end{split}
 ```
-The system starts from a filled state, the chain starts as in the Fermi sea.
 
-## Interleaved chain mapping
+### Interleaved chain MPO
 
 ![image](foldedchain.png)
 
@@ -122,6 +133,278 @@ where the matrices are defined as:
 \end{bmatrix} .
 \end{split}
 ```
+
+## Code implementation
+
+The fermionic mapping can be strightforwardly implemented with the methods of `MPSDyanmics.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.
+
+Both of the examples start with the definition of the physical parameters,
+```
+N = 40      # number of chain sites
+β = 10.0     # inverse temperature
+μ = 0.      # chemical potential
+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 `\beta`, 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$.
+```
+function ϵ(x)
+    return x
+end
+
+function J(x)
+    return sqrt(1 - x^2) # semi-circular density of states
+end
+
+chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, J, save=false) # empty
+chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, J, save=false) # filled
+```
+We then specify the simulation parameters
+```
+dt = 0.25           # time step
+T = 15.0            # simulation time
+method = :DTDVP     # time-evolution method
+Dmax = 150          # MPS max bond dimension
+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 `\epsilon_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. 
+```
+H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
+
+ψ =  unitcol(2,2) # (0,1) filled impurity state
+A = productstatemps(physdims(H), state=[fill(unitcol(2,2), N)..., ψ, fill(unitcol(1,2), N)...]) # MPS
+```
+To avoid the `DTDVP` algorithm from getting stuck in a local minimum, it is better to embed the MPS in a manifold of bond dimension 2 (or more):
+```
+mpsembed!(A, 2) # to embed the MPS in a manifold of bond dimension 2
+```
+We can now define the observables for the two chains and for the impurity
+```
+ob1 = OneSiteObservable("chain1_filled_occup", numb(2), (1,N))
+ob2 = OneSiteObservable("chain2_empty_occup", numb(2), (N+2, 2N+1))
+ob3 = OneSiteObservable("system_occup", numb(2), N+1)
+```
+and run the simulation
+```
+A, dat = runsim(dt, T, A, H;
+                name = "Anderson impurity problem (folded chain)",
+                method = method,
+                obs = [ob1, ob2, ob3], 
+                convobs = [ob1],
+                params = @LogParams(N, ϵd, β, c1, c2),
+                convparams = [prec],   
+                Dlim = Dmax,          
+                savebonddims = true,   # we want to save the bond dimension
+                verbose = false,
+                save = false,
+                plot = true,
+                );
+```
+With very minimal post-processing of the data
+```
+# Reshaping the vector to a column matrix and horizontal concatenation
+system_occup_col = reshape(dat["data/system_occup"], :, 1)
+occ = hcat(dat["data/chain1_filled_occup"]', system_occup_col)
+occ = vcat(occ', dat["data/chain2_empty_occup"])
+```
+we plot the results:
+```
+# Plot the system occupation    
+p1 = plot(
+    dat["data/times"],
+    dat["data/system_occup"],
+    xlabel = L"$t$",
+    ylabel = L"$n_d$",
+    title = "System Occupation",
+    size = (700, 500)
+)
+
+# Plot the occupation of the chain sites
+p2 = heatmap(
+    collect(1:2*N+1),
+    dat["data/times"],
+    transpose(occ),  # Use the matrix form
+    cmap = :coolwarm,
+    aspect_ratio = :auto,
+    xlabel = L"$N_{i,j}$ chain sites",
+    ylabel = L"$t$",
+    title = "Chain Occupation",
+    colorbar = true,
+    size = (700, 500)
+)
+
+# Plot the bond dimensions
+p3 = heatmap(
+    collect(1:2*N+2),
+    dat["data/times"],
+    transpose(dat["data/bonddims"]),
+    cmap = :magma,
+    aspect_ratio = :auto,
+    xlabel = L"$N_{i,j}$ chain sites",
+    ylabel = L"$t$",
+    title = "Bond Dimensions",
+    colorbar = true,
+    size = (700, 500)
+)
+
+# Define indices for columns to be plotted
+columns_to_plot = [1, 5, 10, 15, 20]
+
+# Plot vertical slices for occupancy
+p4 = plot(title = "Chain occupation")
+for col in columns_to_plot
+    plot!(p4, occ[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = "chain occupation")
+end
+
+# Plot vertical slices for bond dimensions
+p5 = plot(title = "Bond Dimensions")
+for col in columns_to_plot
+    plot!(p5, dat["data/bonddims"][:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = L"$\chi$")
+end
+
+# Display the plots
+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 `\epsilon_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. 
+```
+H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
+
+ψ =  unitcol(2,2) # (0,1) filled impurity state
+Tot = Any[]
+push!(Tot, ψ)
+for i in 1:(N)
+    push!(Tot, unitcol(2,2))
+    push!(Tot, unitcol(1,2))
+end
+
+A = productstatemps(physdims(H), state=Tot) # MPS
+```
+To avoid the `DTDVP` algorithm from getting stuck in a local minimum, it is better to embed the MPS in a manifold of bond dimension 2 (or more):
+```
+mpsembed!(A, 2) # to embed the MPS in a manifold of bond dimension 2
+```
+We can finally define the observables for the interleaved chain
+```
+ob1 = OneSiteObservable("system_occup", numb(2), 1)
+ob2 = OneSiteObservable("folded_chain_occup", numb(2), (2,2N+1))
+```
+and run the simulation
+```
+A, dat = runsim(dt, T, A, H;
+                name = "Anderson impurity problem (folded chain)",
+                method = method,
+                obs = [ob1, ob2], 
+                convobs = [ob1],
+                params = @LogParams(N, ϵd, β, c1, c2),
+                convparams = [prec],   
+                Dlim = Dmax,          
+                savebonddims = true,   # we want to save the bond dimension
+                verbose = false,
+                save = false,
+                plot = true,
+                );
+```
+To show the data in a clear way, we do a bit of post-processing of the data to _unfold_ the chain, and move back to the double chain representation:
+```
+unfolded_occ = Vector{Vector{Float64}}()  # Assuming the elements are of type Float64
+unfolded_bonds = Vector{Vector{Float64}}()  # Adjust the type based on actual data
+
+# Populate unfolded_occ by iterating in the specific order mentioned
+for i in 1:N  # Adjusted for 1-based indexing
+    push!(unfolded_occ, dat[ "data/folded_chain_occup"][2N + 1 - 2i, :])
+end
+
+push!(unfolded_occ, dat["data/folded_chain_occup"][1,:])
+
+for i in 2:N
+    push!(unfolded_occ, dat["data/folded_chain_occup"][2i,:])
+end
+
+# Populate unfolded_bonds similarly
+for i in 1:(N+1)  # Adjusted for 1-based indexing
+    push!(unfolded_bonds, dat["data/bonddims"][2N + 3 - 2i,:])  # Assuming bonddims is directly accessible
+end
+
+push!(unfolded_bonds, dat["data/bonddims"][1,:])
+
+for i in 2:(N+1)
+    push!(unfolded_bonds, dat["data/bonddims"][2i,:])
+end
+
+unfolded_bonds_matrix = hcat(unfolded_bonds...)'
+unfolded_occ_matrix = hcat(unfolded_occ...)'
+```
+We conclude by plotting the data
+```
+# Plot the system occupation    
+p1 = plot(
+    dat["data/times"],
+    dat["data/system_occup"],
+    xlabel = L"$t$",
+    ylabel = L"$n_d$",
+    title = "System Occupation",
+    size = (700, 500)
+)
+
+# Plot the occupation of the chain sites
+p2 = heatmap(
+    collect(1:2*N),
+    dat["data/times"],
+    transpose(unfolded_occ_matrix),  # Use the matrix form
+    cmap = :coolwarm,
+    aspect_ratio = :auto,
+    xlabel = L"$N_{i,j}$ chain sites",
+    ylabel = L"$t$",
+    title = "Chain Occupation",
+    colorbar = true,
+    size = (700, 500)
+)
+
+# Plot the bond dimensions
+p3 = heatmap(
+    collect(1:2*N+2),
+    dat["data/times"],
+    transpose(unfolded_bonds_matrix),
+    cmap = :magma,
+    aspect_ratio = :auto,
+    xlabel = L"$N_{i,j}$ chain sites",
+    ylabel = L"$t$",
+    title = "Bond Dimensions",
+    colorbar = true,
+    size = (700, 500)
+)
+
+# Define indices for columns to be plotted
+columns_to_plot = [1, 5, 10, 15, 20]
+
+
+# Plot vertical slices for occupancy
+p4 = plot(title = "Chain occupation")
+for col in columns_to_plot
+    plot!(p4, unfolded_occ_matrix[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = "chain occupation")
+end
+
+# Plot vertical slices for bond dimensions
+p5 = plot(title = "Bond Dimensions")
+for col in columns_to_plot
+    plot!(p5, unfolded_bonds_matrix[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = L"$\chi$")
+end
+
+# Display the plots
+plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
+```
+
 ________________
 ### References