update tutorial & clean deprecated examples

Author: qyli

docs/Project.toml

@@ -3,3 +3,4 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteMPS = "6b572690-f317-4dad-a5fe-5fe73a0b7d42"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"

docs/src/tutorial/Heisenberg.md

@@ -1,14 +1,14 @@
 # Heisenberg Chain
 
 ## Ground state
-In this section, we compute the ground-state, finite-T and dynamical properties of a Heisenberg chain, whose Hamiltonian reads
-$$
-H = J\sum_{i}S_{i}\cdot S_{i+1}.
-$$
+In this section, we compute the ground-state, finite-T and dynamical properties of a Heisenberg chain, whose Hamiltonian reads\
+$H = J\sum_{i}S_{i} S_{i+1}$.
 Here we set $J=1$ as energy unit and use spin SU(2) symmetry.
 ```@example Heisenberg
 using FiniteMPS
-using CairoMakie, Statistics, LsqFit # visualization
+using CairoMakie, Statistics  # visualization
+using LsqFit: curve_fit
+using NumericalIntegration: integrate
 
 mkpath("figs_Heisenberg") # save figures
 
@@ -32,7 +32,7 @@ aspace = Rep[SU₂](i => 1 for i in 0:1/2:1)
 ```
 Here `bspace = Rep[SU₂](0 => 1)` is the space of left boundary bond, `0` is the SU(2) quantum number of the total MPS, and `1` is the multiplicity of the trivial representation, therefore this setup indicates the total MPS is a SU(2) scalar.
 
-`aspace` is the space of bulk bonds. Note the fusing of physical space and bond space leads to a constrain due to symmetry. For example, the SU(2) quantum numbers of the bonds  exhibit a integer/half integer oscillation, as the physical space exactly shifts the quantum number by 1/2. Here we use a larger (with redundancy) initial bond space so that contraction of bond indices gives a non-vanished result. 
+`aspace` is the space of bulk bonds. Note the fusion of physical space and bond space leads to a constrain due to symmetry. For example, the SU(2) quantum numbers of the bonds  exhibit a integer/half integer oscillation, as the physical space exactly shifts the bond quantum number by 1/2. Here we use a larger (with redundancy) initial bond space so that contraction of bond indices gives a non-vanished result. 
 
 ```@example Heisenberg
 # DMRG 
@@ -113,7 +113,7 @@ $$
 $$
 
 ## Finite temperature
-Now we move to the finite-temperature properties via [tanTRG](https://doi.org/10.1103/PhysRevLett.130.226502), which belongs to an imaginary-time-evolution method based on MPO TDVP. Note we will use [CBE-TDVP](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.133.026401) to accelerate the computation.
+Now we move to the finite-temperature properties via [tanTRG](https://doi.org/10.1103/PhysRevLett.130.226502), which belongs to an imaginary-time-evolution method based on TDVP of MPO. Note we will use [CBE-TDVP](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.133.026401) to accelerate the computation.
 
 ```@example Heisenberg
 # define a beta list
@@ -150,7 +150,7 @@ Obs = convert(Dict, ObsTree)
 lsSM[1] = calSM(Obs)
 
 ```
-First we define a $\beta$ list, which determines the step length of imaginary-time cooling. [SETTN](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.161104) is adopted to initializes a high-temperature MPO, where `CBEAlg` indicates the algorithm to implement CBE. Currently only `NaiveCBE` is valid, where we directly find the optimal subspace via a svd (use random svd to accelerate), and the expanded bond dimension is set as `D + div(D, 4)`. `lsnoise` set the noise applied in the first several sweeps of variational multiplication. Next we use TDVP to cool down the system.
+First we define a $\beta$ list, which determines the step length of imaginary-time cooling. [SETTN](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.161104) is adopted to initialize a high-temperature MPO, where `CBEAlg` indicates the algorithm to implement CBE. Currently only `NaiveCBE` is valid, where we directly find the optimal subspace via a svd (use random svd to accelerate), and the expanded bond dimension is set as `D + div(D, 4)`. `lsnoise` sets the noise applied in the first several sweeps of variational multiplication. Next we use TDVP to cool down the system.
 
 ```@example Heisenberg
 # TDVP cooling
@@ -169,13 +169,12 @@ for idx in 2:length(lsβ)
 
 	# update data stored in ObsTree 
 	calObs!(ObsTree, ρ)
-     Obs = convert(Dict, ObsTree)
-	lsSM[idx] = calSM(Obs)
+	lsSM[idx] = calSM(convert(Dict, ObsTree))
 end
 ```
-The key function is this part is `TDVPSweep1!` that performs a single left-to-right and right-to-left 1-TDVP sweep. `GCsweep = true` indicates that a manual `GC.gc()` is called per sweep. If you suffer memory problem when using FiniteMPS.jl, the first thing to try is setting `CGsweep = true` and a stronger `GCstep = true` in the main sweeping function (e.g. `DMRGSweep1!` and `TDVPSweep1!`).
+The key function in this part is `TDVPSweep1!` that performs a single left-to-right and right-to-left 1-TDVP sweep. `GCsweep = true` indicates that a manual `GC.gc()` is called per sweep. If you suffer memory problem when using FiniteMPS.jl, the first thing to try is setting `CGsweep = true` and a stronger `GCstep = true` in the main sweeping function (e.g. `DMRGSweep1!` and `TDVPSweep1!`).
 
-In each sweep, the normalization factor after an imaginary-time evolution is extracted to calculate the partition function `lnZ`. In order to calculate spin correlations, regenerating `ObsTree` is not needed, just use `calObs!` again to trigger the in-place update with the new MPO `rho`.
+In each sweep, the normalization factor after an imaginary-time evolution is extracted to calculate the partition function `lnZ`. In order to calculate spin correlations, regenerating `ObsTree` is not necessary, just use `calObs!` again to trigger the in-place update with the new MPO `ρ`.
 
 Below is a simple script for visualization, where the temperature dependence of energy $e$, specific heat $c_V$ and AFM structure factor $S(\pi, \pi)$ are shown. 
 ```@example Heisenberg
@@ -184,7 +183,7 @@ Below is a simple script for visualization, where the temperature dependence of
 # compute C = - ∂S / ∂lnβ
 lsS = lsβ .* lsE .+ lslnZ 
 lslnβ = log.(lsβ)
-lsCe = - diff(lsS) ./ diff(lslnβ)
+lsCv = - diff(lsS) ./ diff(lslnβ)
 lsβ_c = exp.((lslnβ[1:end-1] + lslnβ[2:end])/2)
 
 fig = Figure(size = (480, 400))
@@ -204,7 +203,7 @@ end
 scatterlines!(ax1, 1 ./ lsβ, lsE ./ L)
 lines!(ax1, [0.05, 1], [Eg/L, Eg/L]; color = :red, label = "DMRG")
 axislegend(ax1; position = (0, 0.5))
-scatterlines!(ax2, 1 ./ lsβ_c, lsCe ./ L)
+scatterlines!(ax2, 1 ./ lsβ_c, lsCv ./ L)
 scatterlines!(ax3, 1 ./ lsβ, lsSM)
 lines!(ax3, [0.05, 1], [SM_GS, SM_GS]; color = :red)
 
@@ -214,8 +213,131 @@ save("figs_Heisenberg/FiniteT.png", fig)
 From this example we see that the low-temperature limit of tanTRG does shake hands with the ground state DMRG.
 
 ## Spin dynamics
-TODO
+In this section we will compute the ground-state dynamical spin structure factor (DSF)
+$S(k,\omega) = \sum_{i = 1}^L e^{-ik(r_i-r_j)} \int_{-\infty}^\infty dt e^{i\omega t}\langle S_i^z(t)S_j^z\rangle$,
+where $j$ is chosen as a reference site. More specifically, we will use TDVP to compute $\langle S_i^z(t)S_j^z\rangle = e^{iE_gt}\langle \Psi| S_i^z |\Phi(t)\rangle$ where $|\Phi(t)\rangle = e^{-iHt}S_j^z|\Psi\rangle$.
+```@example Heisenberg
+lst = 0:1.0:10 # time list
+matSijt = zeros(ComplexF64, L, length(lst)) # S_{i,j_ref}(t) 
+
+# obtain |Φ(0)⟩ = S_j|Ψ⟩ 
+j_ref = div(L, 2) 
+
+# fuse the additional bonds due to non-abelian symmetry
+aspace_S = codomain(SU2Spin.SS[2])[1]
+aspace_Ψ = codomain(Ψ[1])[1]
+aspace_Φ = fuse(aspace_S, aspace_Ψ)
+
+# wrap the operator S_j to a MPO 
+Tree = InteractionTree(L)
+addIntr!(Tree, SU2Spin.SS[2], j_ref, 1.0; name = :S)
+S_MPO = AutomataMPO(Tree)
+
+# note Φ should be a complex MPS
+Φ = randMPS(ComplexF64, fill(SU2Spin.pspace, L), vcat(aspace_Φ, fill(aspace, L - 1)))
+
+# variationally find |Φ⟩ = S_j|Ψ⟩ 
+mul!(Φ, S_MPO, Ψ;
+	CBEAlg = NaiveCBE(D + div(D, 4), 1e-8; rsvd = true),
+	trunc = truncdim(D), GCsweep = true,
+     lsnoise = [0.1, 0.01, 0.001], tol = 1e-12,
+)
+```
+Note we cannot directly compute $\langle S_i^z(t)S_j^z\rangle$, as the operator $S^z$ breaks the SU(2) symmetry. Therefore, what we actually compute is $\langle S_i(t) \cdot S_j\rangle / 3$, where $S_j$ operator is a SU(2) spinor with an additional index that labels the representation space. After preparing the initial `Φ` with correct symmetry index, we call `mul!` function to perform variational multiplication. 
+
+```@example Heisenberg
+# define a new ObsTree to calculate each inner product ⟨Ψ|S_i|Φ⟩
+ObsTree = ObservableTree(L)
+for i in 1:L
+	addObs!(ObsTree, SU2Spin.SS[1], i; name = :S)
+end
+
+iso = isometry(aspace_Φ, aspace_S ⊗ aspace_Ψ)
+El = permute(iso', ((2, 1), (3,)))
+calObs!(ObsTree, Ψ, Φ; El = El)
+matSijt[:, 1] = [ObsTree.Refs["S"][(i,)][] / 3 for i in 1:L]
+```
+Emphasize that the isometry to fuse the space of $S_j$ and $|\Psi\rangle$ must be considered carefully. We exactly insert a pair of isometries without changing the contraction result. One is implicitly applied in `mul!` to fuse the two additonal indices to the new one of `Φ`, and the other becomes the left boundary environment tensor `El` sent to `calObs!`. 
+
+```@example Heisenberg
+# time evolution
+
+# construct the trilayer environment for TDVP 
+Env = Environment(Φ', H, Φ)
+# time evolution
+for idx in 2:length(lst)
+	dt = lst[idx] - lst[idx-1]
+
+	TDVPSweep1!(Env, -im * dt;
+		CBEAlg = NaiveCBE(D + div(D, 4), 1e-8; rsvd = true),
+		GCsweep = true, trunc = truncdim(D)
+	)
+
+	calObs!(ObsTree, Ψ, Φ; El = El)
+	# e^{iEgt}⟨Ψ|S_i^+|Φ⟩
+	matSijt[:, idx] = exp(im * Eg * lst[idx]) * map(1:L) do i
+		ObsTree.Refs["S"][(i,)][] / 3
+	end
+end
+```
+In this part we use CBE-TDVP to perform real-time evolution of `Φ` and calculate the time-dependent spin correlations as some inner products. 
+
+```@example Heisenberg
+# spatial FT
+lsk = 0:2/L:2 # unit = π
+matcoef = map([(k, i) for k in lsk, i in 1:L]) do (k, i)
+	exp(-im * k * (i - j_ref) * π)
+end
+matSkt = matcoef * matSijt
+
+# time FT
+lsω = 0:0.1:4 # frequency list
+
+# use a Parzen window to suppress the non-physical oscillation
+function ParzenWindow(x::Float64)
+	if abs(x) < 1 / 2
+		return 1 - 6 * x^2 + 6 * abs(x)^3
+	elseif abs(x) < 1
+		return 2 * (1 - abs(x))^3
+	else
+		return 0.0
+	end
+end
+lsWt = ParzenWindow.(lst ./ maximum(lst))
+
+matSkω = zeros(length(lsk), length(lsω))
+for iω in 1:length(lsω)
+	lscoef = exp.(im * lsω[iω] .* lst) .* lsWt
+	for ik in 1:length(lsk)
+		matSkω[ik, iω] = 2 * integrate(lst, real.(matSkt[ik, :] .* lscoef))
+	end
+end
+```
+After collecting all real-space time-dependent spin correlations, we perform FT to obtain the DSF. Note time-reversal symmetry is used in the numerical integration so that only $t>0$ data is required.
+Moreover, we multiply a parzen window function to the time domain to suppress the non-physical oscillation in frequency domain.
 
+```@example Heisenberg
+# visualization
+fig = Figure(size = (480, 240))
+ax = Axis(fig[1, 1];
+     xlabel = L"k / \pi", 
+     ylabel = L"\omega",
+     limits = ((0, 2), extrema(lsω))
+)
+
+hm = heatmap!(ax, lsk, lsω, matSkω)
+Colorbar(fig[1, 2], hm; label = L"S(k, \omega)")
+
+# upper and lower boundaries
+lsω_upper = π .* sin.(lsk .* (π/2))
+lsω_lower = (π/2) .* abs.(sin.(lsk .* π)) 
+lines!(ax, lsk, lsω_upper; color = :white, linestyle = :dash)
+lines!(ax, lsk, lsω_lower; color = :white, linestyle = :dash)
+
+save("figs_Heisenberg/Skomega.png", fig)
+```
+![](./figs_Heisenberg/Skomega.png)
+The exact (thermodynamical limit) upper and lower boundaries are also shown. Our numerial results align well with them, up to the frequency resolution (FWHM $\approx 8 / t_\textrm{max}$) due to the finite evolution time $t_\textrm{max} = 10$.   
 
 
 

docs/src/tutorial/Observable.md

@@ -34,6 +34,6 @@ Here `Tree` is an `ObservableTree` object that contains all observables to be ca
 ```@example Observable
 Obs["SzSz"][(1, 2)]
 ```
-is the correlation `\langle S_1^z S_2^z\rangle`. One can perform a simple quantum mechanics calculation to check this result.
+is the correlation $\langle S_1^z S_2^z\rangle$. One can perform a simple quantum mechanics calculation to check this result.
 
 Here is just a simple example to show the basic usage, more complex examples that contain fermion correlations and multi-site correlations (e.g. pairing correlations) can be found in the concrete example for [Hubbard model](@ref Hubbard).