diff --git a/docs/make.jl b/docs/make.jl index f8a741973..0b3a17ba3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -54,6 +54,9 @@ makedocs(; "tutorials/installation.md", "tutorials/first_groundstate.md", "tutorials/thermodynamic_limit.md", + "tutorials/time_evolution.md", + "tutorials/excitations.md", + "tutorials/using_symmetries.md", ], "Manual" => [ "man/intro.md", @@ -64,11 +67,18 @@ makedocs(; "man/parallelism.md", ], "How-to" => [ + "howto/index.md", "howto/states.md", "howto/hamiltonians.md", + "howto/groundstate_algorithms.md", "howto/bond_dimension.md", + "howto/time_evolution.md", "howto/observables.md", "howto/entanglement.md", + "howto/excitations.md", + ], + "Concepts" => [ + "concepts/algorithm_landscape.md", ], "Examples" => "examples/index.md", "Library" => [ diff --git a/docs/src/concepts/algorithm_landscape.md b/docs/src/concepts/algorithm_landscape.md new file mode 100644 index 000000000..161e3a31d --- /dev/null +++ b/docs/src/concepts/algorithm_landscape.md @@ -0,0 +1,141 @@ +# [The algorithm landscape](@id concept_algorithm_landscape) + +MPSKit deliberately separates *what* you want to compute from *how* it gets computed. +Entry points such as [`find_groundstate`](@ref), [`timestep`](@ref), [`excitations`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref) each accept several interchangeable algorithm structs, and the package ships more than a dozen of them. +That flexibility exists because no single algorithm wins everywhere: some only apply to finite or only to infinite systems, some can grow the bond dimension while others cannot, and their relative performance depends on the model at hand. + + +This page is the decision guide. +It starts from a table that maps each task onto the algorithm(s) of choice, and then walks through the reasoning behind each row. +It explains *why* you would pick one algorithm over another; for the *how* — the actual calls, keywords, and worked recipes — follow the links into the how-to pages. + +## The decision table + +| Task | Finite system | Infinite system | +|:-----|:--------------|:----------------| +| **Ground state** ([`find_groundstate`](@ref)) | [`DMRG`](@ref) (workhorse, fixed bond dimension); [`DMRG2`](@ref) (grows bond dimension, requires `trscheme`); [`GradientGrassmann`](@ref) (final polish) | [`VUMPS`](@ref) (workhorse, needs a unique ground state); [`IDMRG`](@ref) / [`IDMRG2`](@ref) (two-site requires `trscheme` and a unit cell of at least two sites); [`GradientGrassmann`](@ref) (final polish) | +| **Time evolution** ([`timestep`](@ref) / [`time_evolve`](@ref)) | [`TDVP`](@ref) (fixed bond dimension); [`TDVP2`](@ref) (grows bond dimension, requires `trscheme`); or [`make_time_mpo`](@ref) ([`WI`](@ref) / [`WII`](@ref) / [`TaylorCluster`](@ref)) applied with [`approximate`](@ref) | [`TDVP`](@ref) (no two-site variant exists); or [`make_time_mpo`](@ref) applied with [`approximate`](@ref) | +| **Excitations** ([`excitations`](@ref)) | [`QuasiparticleAnsatz`](@ref) (the only one supporting charged `sector`s); [`FiniteExcited`](@ref) (penalty method); [`ChepigaAnsatz`](@ref) / [`ChepigaAnsatz2`](@ref) (cheap, from ground-state environments) | [`QuasiparticleAnsatz`](@ref) (momentum-resolved, the only choice) | +| **Boundary / statistical mechanics** ([`leading_boundary`](@ref)) | apply the transfer MPO row by row with [`approximate`](@ref) | [`VUMPS`](@ref); [`VOMPS`](@ref) (power method); [`IDMRG`](@ref) / [`IDMRG2`](@ref); [`GradientGrassmann`](@ref) (hermitian, positive transfer matrices) | +| **Compression / approximation** ([`approximate`](@ref), [`changebonds`](@ref)) | [`approximate`](@ref) with [`DMRG`](@ref) / [`DMRG2`](@ref); [`SvdCut`](@ref) via [`changebonds`](@ref) for local truncation | [`approximate`](@ref) with [`IDMRG`](@ref) / [`IDMRG2`](@ref) / [`VOMPS`](@ref); [`SvdCut`](@ref) via [`changebonds`](@ref) for local truncation | + + + +A few structural facts hold across the whole table and are worth internalizing early. +Every two-site algorithm (`DMRG2`, `IDMRG2`, `TDVP2`) requires an explicit `trscheme` keyword and can change the bond dimension as it runs; the single-site variants with their default settings cannot. +`IDMRG2` additionally needs a unit cell of at least two sites, and `TDVP2` exists only for finite MPS. +Finally, algorithms compose: the `&` operator chains two algorithms into one, running the first to completion and handing its result to the second, which is how two-site warm-up passes and gradient-descent polishing stages are combined with a workhorse algorithm in a single call. + +## Ground states + +The classic approach is alternating local optimization: [`DMRG`](@ref) sweeps through a finite chain, optimizing one site while all others are held fixed. + +The catch is the fixed bond dimension: a single-site update can never enlarge the virtual spaces, so the precision of the calculation is locked in by the initial state. +This bites hardest when symmetries are involved, because then not just the total bond dimension but its distribution over charge sectors is frozen, and a poor initial distribution cannot be repaired. + +[`DMRG2`](@ref) fixes this by optimizing two neighbouring sites jointly and truncating back down, which lets the bond dimension (and its sector distribution) adapt, at a higher cost per sweep. + + +For infinite systems, two philosophies compete. +[`IDMRG`](@ref) grows the system from the middle outwards, repeatedly inserting and optimizing new sites until the boundary is no longer felt; [`IDMRG2`](@ref) is its two-site, bond-growing variant. + +Because convergence requires the effective system to outgrow the correlation length, IDMRG can be slow to converge for critical systems, where that length diverges. + +[`VUMPS`](@ref) instead works with a genuinely uniform state: each local update is followed by a re-gauging step that replaces *every* tensor in the infinite chain with the updated one, so the effect of an update is felt throughout the system immediately. + +This often gives VUMPS a higher convergence rate than IDMRG, which is why it is the default infinite-system workhorse. + +The price is an injectivity requirement: VUMPS assumes a unique ground state, and it is not the right tool when the state it should converge to is non-injective. + +Like DMRG, VUMPS is single-site and cannot alter the bond dimension. + +[`GradientGrassmann`](@ref) approaches the problem from a third direction: the MPS tensors form a Riemannian manifold (a Grassmann manifold), and one can run gradient descent directly on it, for finite and infinite states alike. +Its niche is the tail of the optimization: close to convergence its rate is often the best of the lot, while far from convergence the sweeping algorithms tend to make faster progress. + +The practical consequence is the chaining pattern: run a cheap workhorse first, then hand over to gradient descent, e.g. `VUMPS(...) & GradientGrassmann(...)`. + +This pattern is baked into `find_groundstate` itself: called with only keywords, it picks `DMRG` for a finite state and `VUMPS` for an infinite one, appends a `GradientGrassmann` stage on infinite states when the requested tolerance is tighter than `1e-4`, and prepends a two-site pass (`DMRG2` or `IDMRG2`) whenever you supply a `trscheme`. +Since gradient descent is also a single-site method, growing the bond dimension remains the job of that two-site pre-pass or of [`changebonds`](@ref). + +For call syntax, keyword tables, and worked chaining examples, see [Ground-state algorithms](@ref howto_groundstate_algorithms). + +## Time evolution + +MPSKit solves the time-dependent Schrödinger equation along two distinct routes, and the choice between them is a genuine trade-off rather than a finite/infinite split. + +The first route, [`TDVP`](@ref), never builds the evolution operator at all. +It projects the Schrödinger equation onto the tangent space of the current MPS, solves the projected equation for a small time step, and repeats. + +Its two-site variant [`TDVP2`](@ref) plays the same role as `DMRG2` does for `DMRG`: it lets the bond dimension grow to absorb the entanglement generated by the evolution, at extra cost, and it exists only for finite systems. + + +The second route splits the problem in two: first approximate the evolution operator ``\exp(-iH\,dt)`` itself as an MPO using [`make_time_mpo`](@ref) — with [`WI`](@ref), [`WII`](@ref), or [`TaylorCluster`](@ref) as the approximation scheme — and then apply that MPO to the state with [`approximate`](@ref). +The appeal is amortization: for a time-independent Hamiltonian and a fixed step size the MPO is built once and reused for every step, and the accuracy of the operator approximation is controlled independently of the accuracy of its application. + + + +Both routes accept an `imaginary_evolution` keyword for evolution in imaginary time. +For step-by-step recipes along either route, see [Time evolution](@ref howto_time_evolution). + +## Excitations + +Resolving states deep in the spectrum is generally out of reach, but three families of algorithms target the low-lying part, each with a distinct character. + + +The [`QuasiparticleAnsatz`](@ref) is the most broadly applicable: it works for finite and infinite systems, and it is the only algorithm that can target excitations carrying a nontrivial symmetry charge, via the `sector` keyword. +It builds an excited state by replacing a single tensor of the ground-state MPS — summed over all positions on a finite chain, or in a momentum-carrying plane-wave superposition on an infinite one — and solves the resulting eigenvalue problem. + +Because the variational class consists of local perturbations on top of the ground state, it is the natural choice for quasiparticle-like excitations, and on infinite systems it is the only option, giving direct access to dispersion relations. + + +[`FiniteExcited`](@ref) takes a brute-force approach available only on finite chains: it reruns a full ground-state optimization on a modified Hamiltonian that carries an energy penalty for overlapping with all previously found states. +Each new excited state therefore costs another complete ground-state search, and the orthogonality to earlier states is only approximate (enforced by the penalty `weight`, not exactly). + +Its advantage is that it makes no assumption about the *form* of the excited state: since each state is a fully variational `FiniteMPS`, it can capture excitations that are not well described as a local perturbation of the ground state. + + +The [`ChepigaAnsatz`](@ref) (and its two-site refinement [`ChepigaAnsatz2`](@ref)) is the cheapest of the three, also finite-only. +It observes that the gauged ground-state MPS tensors act as isometries projecting the Hamiltonian into a low-energy subspace, so the low-lying spectrum can be read off by diagonalizing the effective Hamiltonian already available from the ground-state environments, with no additional sweeping. + +This works best precisely where excitations are hard for the other methods: in critical systems with long-range correlations, where the excitation weight is spread across the whole chain. + + +For the call signatures, momentum scans, and sector-targeting recipes, see [Excited states](@ref howto_excitations). + +## Boundaries and statistical mechanics + +MPS algorithms are not limited to Hamiltonian problems. +A two-dimensional classical partition function can be written as an infinite power of a row-to-row transfer MPO, and contracting the network amounts to finding that operator's dominant eigenvector — a boundary MPS. + +This is the job of [`leading_boundary`](@ref), which accepts a familiar cast: [`VUMPS`](@ref) and [`IDMRG`](@ref)/[`IDMRG2`](@ref) carry over directly from the ground-state problem, [`GradientGrassmann`](@ref) applies when the transfer MPO is hermitian and positive, and [`VOMPS`](@ref) is a power method specific to this setting, which iteratively approximates the operator-times-state product by a new state of the same bond dimension. + + + +## Compression and changing bond dimension + +Two mechanisms round out the landscape by manipulating states rather than solving for new ones. +[`approximate`](@ref) variationally fits a new MPS, typically of different bond dimension, to the result of applying an MPO to a state; the sweeping ground-state algorithms (`DMRG`/`DMRG2` for finite, `IDMRG`/`IDMRG2`/`VOMPS` for infinite) double as its optimization engines. +This is the same machinery that applies time-evolution MPOs, and combined with [`SvdCut`](@ref) it yields a globally optimal truncation of a state. + +[`changebonds`](@ref), by contrast, performs direct local surgery on a state: truncating with [`SvdCut`](@ref), or expanding with [`OptimalExpand`](@ref), [`RandExpand`](@ref), or [`VUMPSSvdCut`](@ref) so that the single-site algorithms above have room to work with. +The trade-offs between those expansion schemes, and recipes for when to grow, are covered in [Controlling bond dimension](@ref howto_bond_dimension). + +## Where to go next + +The how-to pages turn each row of the table into runnable recipes: [Ground-state algorithms](@ref howto_groundstate_algorithms), [Time evolution](@ref howto_time_evolution), and [Excited states](@ref howto_excitations), with [Controlling bond dimension](@ref howto_bond_dimension) supporting all three. +For the complete signatures, keyword lists, and docstrings of every algorithm named here, see the library reference: [Ground-state algorithms](@ref lib_groundstate), [Time evolution](@ref lib_time_evolution), and [Excitations](@ref lib_excitations). + + diff --git a/docs/src/howto/excitations.md b/docs/src/howto/excitations.md new file mode 100644 index 000000000..39b3b03ec --- /dev/null +++ b/docs/src/howto/excitations.md @@ -0,0 +1,168 @@ +# [Excited states](@id howto_excitations) + +The examples on this page use MPSKit.jl, MPSKitModels.jl, and TensorKit.jl. +See [Installation](@ref tutorial_installation) for how to add these packages to your environment. + +[`excitations`](@ref) is the single entry point for computing energy eigenstates beyond the ground state. +This page shows how to call it for a gap, a full dispersion relation, a charged excitation, and a handful of excited states on a finite chain. +For what each algorithm actually does and why you would choose one over another, see [Excitations](@ref lib_excitations). +All examples share a single namespace: + +```@example excitations_howto +using MPSKit, MPSKitModels, TensorKit +``` + +--- + +## 1. Get a single excitation gap on an infinite chain + +On an `InfiniteMPS`, [`QuasiparticleAnsatz`](@ref) perturbs every site of the unit cell in a plane-wave superposition with a fixed `momentum`, given as a `Real` in radians per unit cell. +Pass the ground state and (optionally) its environments straight through from [`find_groundstate`](@ref): + +```@example excitations_howto +g = 2.0 +H_inf = transverse_field_ising(; g) +ψ₀_inf = InfiniteMPS(ℂ^2, ℂ^12) +ψ_inf, envs_inf, = find_groundstate(ψ₀_inf, H_inf; verbosity = 0) + +Es_inf, ϕs_inf = excitations(H_inf, QuasiparticleAnsatz(), 0.0, ψ_inf, envs_inf; num = 1) +Es_inf[1] +``` + +The values in `Es_inf` are excitation *gaps* above the ground-state energy density, not total energies: internally the ground-state energy per site is subtracted before diagonalizing. +`ϕs_inf` holds the corresponding quasiparticle states (`num` of them), which behave like normal vectors for `eigsolve`-style post-processing but are not `FiniteMPS`/`InfiniteMPS` objects themselves. +Raise `num` to get more than one state at the same momentum, e.g. `num = 3` for the three lowest excitations at that momentum. + +--- + +## 2. Scan the dispersion relation + +Pass a range (or any vector) of momenta instead of a single number to sweep the whole Brillouin zone in one call: + +```@example excitations_howto +momenta = range(0, π, 5) +Es_disp, ϕs_disp = excitations( + H_inf, QuasiparticleAnsatz(), momenta, ψ_inf, envs_inf; + num = 1, verbosity = 0 +) +size(Es_disp) +``` + +With a vector of `length(momenta)` momenta, `Es_disp` and `ϕs_disp` come back as `(length(momenta), num)` matrices rather than plain vectors — index `Es_disp[:, n]` for the dispersion of the `n`-th branch, or use `vec(Es_disp)` when `num = 1`. +`verbosity = 0` silences the per-momentum `@info` line that this method otherwise prints; raise it to see progress on a longer scan. +Momenta are independent of each other, so this method also accepts `parallel = true` (the default) to distribute them over available threads/workers; pass `parallel = false` to force sequential evaluation. + +--- + +## 3. Target a symmetry sector + +By default the optimization looks for the lowest excitation with trivial (vacuum) total charge, `sector = leftunit(ψ)`. +Passing a different `TensorKit` sector targets a quasiparticle with that charge instead — only `QuasiparticleAnsatz` supports this keyword. +Build the ground state with symmetric tensors first, then request the sector on the excitation call: + +```@example excitations_howto +g = 10.0 +L = 12 +H_Z2 = transverse_field_ising(Z2Irrep, FiniteChain(L); g) +ψ₀_Z2 = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => 8, 1 => 8)) +ψ_Z2, envs_Z2, = find_groundstate(ψ₀_Z2, H_Z2; verbosity = 0) + +Es_triv, = excitations(H_Z2, QuasiparticleAnsatz(), ψ_Z2, envs_Z2; num = 1) +Es_charged, ϕs_charged = excitations( + H_Z2, QuasiparticleAnsatz(), ψ_Z2, envs_Z2; + num = 1, sector = Z2Irrep(1) +) +Es_triv[1], Es_charged[1] +``` + +Here the `Z2Irrep(1)` excitation corresponds to a single flipped spin, the lowest physical excitation of the transverse-field Ising model. +[`ChepigaAnsatz`](@ref)/[`ChepigaAnsatz2`](@ref) do not support charged excitations at all: passing a nontrivial `sector` to either raises an error, and [`FiniteExcited`](@ref) has no `sector` keyword in the first place. + +--- + +## 4. Excited states on a finite chain + +Momentum is not a conserved quantity on a finite chain, so the finite method of `excitations` has no momentum argument; drop it entirely and call `QuasiparticleAnsatz` on the ground state directly: + +```@example excitations_howto +L = 12 +H_fin = transverse_field_ising(FiniteChain(L); g) +ψ₀_fin = FiniteMPS(L, ℂ^2, ℂ^16) +ψ_fin, envs_fin, = find_groundstate(ψ₀_fin, H_fin; verbosity = 0) + +Es_qp, ϕs_qp = excitations(H_fin, QuasiparticleAnsatz(), ψ_fin, envs_fin; num = 1) +Es_qp[1] +``` + +[`FiniteExcited`](@ref) takes a different approach: it repeatedly finds the ground state of `H + weight * Σᵢ |ψᵢ⟩⟨ψᵢ|`, penalizing overlap with the ground state and any excited states already found, and returns full `FiniteMPS` objects instead of quasiparticle states: + +```@example excitations_howto +fe_alg = FiniteExcited(; gsalg = DMRG(; verbosity = 0), weight = 10.0) +Es_fe, ψs_fe = excitations(H_fin, fe_alg, ψ_fin; num = 2) +Es_fe +``` + +Unlike `QuasiparticleAnsatz`, the values in `Es_fe` are total energies of the excited states, directly comparable to `expectation_value(ψ_fin, H_fin)` on the ground state. + +Because each call to `FiniteExcited` reruns a full ground-state optimization under the hood, it scales worse with `num` than the other methods here; reach for it when you need excited states of a genuinely different character than the ground state (so the projector penalty, rather than a local perturbation, is what finds them). + +[`ChepigaAnsatz`](@ref) is a cheaper alternative for excitations that are qualitatively similar to the ground state: it diagonalizes the effective Hamiltonian at a single site `pos` (default the middle of the chain) using the ground-state environments, with no extra sweeping: + +```@example excitations_howto +Es_ch, ψs_ch = excitations(H_fin, ChepigaAnsatz(), ψ_fin, envs_fin; num = 1, pos = L ÷ 2) +Es_ch[1] +``` + +[`ChepigaAnsatz2`](@ref) does the same with a two-site block at `pos, pos + 1`, which costs more but is typically more accurate; it truncates the optimized two-site tensor back down with a `trscheme` keyword (`notrunc()` by default): + +```@example excitations_howto +Es_ch2, ψs_ch2 = excitations(H_fin, ChepigaAnsatz2(; trscheme = truncrank(16)), ψ_fin, envs_fin; num = 1) +Es_ch2[1] +``` + +Like `FiniteExcited`, the energies returned by both Chepiga variants are total energies, not gaps. + + +--- + +## 5. Check excitation quality + +[`variance`](@ref) accepts a quasiparticle state directly and reports the variance of the energy, with smaller values indicating a better-converged excitation: + +```@example excitations_howto +variance(ϕs_qp[1], H_fin) +``` + +It also works on the infinite quasiparticle states from §1–§3: + +```@example excitations_howto +variance(ϕs_inf[1], H_inf) +``` + +!!! warning "Variance of infinite quasiparticle states" + `variance` on an infinite quasiparticle state carries an unresolved implementation note in `src/algorithms/toolbox.jl` and may be unreliable; verify its output before relying on it as a convergence diagnostic. + It also throws an `ArgumentError` for domain-wall (topological) excitations, where it is not implemented at all. + +To measure other observables on a finite quasiparticle state, convert it to a plain `FiniteMPS` first: + +```@example excitations_howto +excited_state = convert(FiniteMPS, ϕs_qp[1]) +real(expectation_value(excited_state, H_fin)) +``` + +--- + +## 6. Domain-wall excitations + +`excitations` also accepts two *different* ground states, `excitations(H, QuasiparticleAnsatz(), momentum, ψ_left, envs_left, ψ_right, envs_right; ...)`, which builds a quasiparticle that interpolates between them — a domain-wall (topological) excitation rather than a local perturbation on top of a single ground state. +This is real, exported functionality, but it has no dedicated test or example in the repository at the time of writing, and constructing two genuinely distinct, well-converged ground states to feed it (for instance the two symmetry-broken ground states of an ordered phase) is itself nontrivial to set up reliably in a short recipe. + + +--- + +## Where to go next + +For growing the bond dimension of the ground states these recipes start from, see [Controlling bond dimension](@ref howto_bond_dimension). +For background on the quasiparticle ansatz and the other algorithms used here, see [Excitations](@ref lib_excitations). +For general expectation values and correlators, including on the converted excited states from §5, see [Computing observables](@ref howto_observables). +Spectral functions built from these excitations (`propagator`, `DynamicalDMRG`, and related solvers) are a separate topic covered in the "Linear problems and spectral functions" section of the [Public API](@ref public_api) reference. diff --git a/docs/src/howto/groundstate_algorithms.md b/docs/src/howto/groundstate_algorithms.md new file mode 100644 index 000000000..4cb753a73 --- /dev/null +++ b/docs/src/howto/groundstate_algorithms.md @@ -0,0 +1,207 @@ +# [Ground-state algorithms](@id howto_groundstate_algorithms) + +The examples on this page use MPSKit.jl, MPSKitModels.jl, TensorKit.jl, and TensorKitTensors.jl. +See [Installation](@ref tutorial_installation) for how to add these packages to your environment. + +[`find_groundstate`](@ref) is the single entry point for optimizing an MPS towards the ground state of a Hamiltonian. +This page shows how to pick and configure the algorithm it runs, for both finite and infinite systems. +For what each algorithm actually does and why you would choose one over another, see [Ground-state algorithms](@ref lib_groundstate). +All examples share a single namespace: + +```@example groundstate_algs +using MPSKit, MPSKitModels, TensorKit +using TensorKitTensors.SpinOperators: σˣ, σᶻ +``` + +--- + +## 1. Get a ground state with defaults + +Called with just a state and a Hamiltonian, `find_groundstate` inspects the type of the initial state and picks a matching algorithm for you. +For a `FiniteMPS` it runs [`DMRG`](@ref) with the keywords you pass through (`tol`, `maxiter`, `verbosity`): + +```@example groundstate_algs +L = 8 +ψ₀ = FiniteMPS(L, ℂ^2, ℂ^8) +H = transverse_field_ising(FiniteChain(L); g = 0.5) + +ψ, envs, ϵ = find_groundstate(ψ₀, H; tol = 1.0e-8, maxiter = 50, verbosity = 0) +ϵ +``` + +For an `InfiniteMPS` it instead runs [`VUMPS`](@ref), and if the requested `tol` is tighter than `1e-4` it chains a [`GradientGrassmann`](@ref) pass afterwards to polish the last few digits (see [§4](#4-refine-convergence-with-gradientgrassmann)): + +```@example groundstate_algs +ψ₀_inf = InfiniteMPS(ℂ^2, ℂ^6) +H_inf = transverse_field_ising(; g = 0.5) + +ψ_inf, envs_inf, ϵ_inf = find_groundstate(ψ₀_inf, H_inf; verbosity = 0) +ϵ_inf +``` + +Passing a `trscheme` keyword switches on a two-site pre-pass that can grow the bond dimension before the single-site algorithm takes over. +On a `FiniteMPS` this prepends [`DMRG2`](@ref); on an `InfiniteMPS` it prepends [`IDMRG2`](@ref) (which needs a unit cell of at least two sites, see [§3](#3-configure-infinite-system-algorithms)). + +```@example groundstate_algs +ψ_auto, envs_auto, ϵ_auto = find_groundstate( + ψ₀, H; + trscheme = truncrank(16), verbosity = 0 +) +ϵ_auto +``` + +!!! note "When to reach for an explicit algorithm" + The keyword form above covers the common cases. + Reach for an explicit algorithm struct (`DMRG`, `DMRG2`, `VUMPS`, `IDMRG`, `IDMRG2`, `GradientGrassmann`), or a chain of them with `&`, whenever you need finer control than the heuristic provides — the rest of this page shows how. + +--- + +## 2. Configure finite-system DMRG + +Pass a [`DMRG`](@ref) struct explicitly to set `tol`, `maxiter`, and `verbosity` directly: + +```@example groundstate_algs +ψ_dmrg, envs_dmrg, ϵ_dmrg = find_groundstate( + ψ₀, H, + DMRG(; tol = 1.0e-8, maxiter = 50, verbosity = 0) +) +ϵ_dmrg +``` + +`DMRG` updates one site at a time, so with its defaults it cannot change the bond dimension: whatever bond dimension `ψ₀` starts with is what it keeps. + +[`DMRG2`](@ref) optimizes two sites at once and truncates back down, which lets it grow (or shrink) the bond dimension as it sweeps, at extra cost per step. +Unlike `DMRG`, `DMRG2` has no default truncation scheme, so `trscheme` is required: + +```@example groundstate_algs +ψ_dmrg2, envs_dmrg2, ϵ_dmrg2 = find_groundstate( + ψ₀, H, + DMRG2(; trscheme = truncrank(16), maxiter = 5, verbosity = 0) +) +ϵ_dmrg2 +``` + +A common pattern is to warm up with `DMRG2` to grow the bond dimension, then refine with the cheaper single-site `DMRG`. +The `&` chaining operator runs the first algorithm to completion, then feeds its result into the second: + +```@example groundstate_algs +warmup_then_refine = DMRG2(; trscheme = truncrank(16), maxiter = 3, verbosity = 0) & + DMRG(; tol = 1.0e-8, maxiter = 30, verbosity = 0) + +ψ_c, envs_c, ϵ_c = find_groundstate(ψ₀, H, warmup_then_refine) +ϵ_c +``` + +For more on choosing `trscheme` and growing bond dimension in general, see [Controlling bond dimension](@ref howto_bond_dimension). + +--- + +## 3. Configure infinite-system algorithms + +[`VUMPS`](@ref) is the default single-site algorithm for an `InfiniteMPS`, and takes the same `tol`/`maxiter`/`verbosity` keywords: + +```@example groundstate_algs +ψ_v, envs_v, ϵ_v = find_groundstate( + ψ₀_inf, H_inf, + VUMPS(; tol = 1.0e-8, maxiter = 50, verbosity = 0) +) +ϵ_v +``` + + +[`IDMRG`](@ref) is the infinite analogue of `DMRG`: it grows the system by repeatedly inserting sites in the middle and re-optimizing, until boundary effects wash out. + +```@example groundstate_algs +ψ_i, envs_i, ϵ_i = find_groundstate( + ψ₀_inf, H_inf, + IDMRG(; tol = 1.0e-8, maxiter = 50, verbosity = 0) +) +ϵ_i +``` + + + +In practice, prefer `VUMPS` unless you specifically need `IDMRG`'s ability to change the bond dimension one site at a time. + +[`IDMRG2`](@ref) is the two-site, bond-dimension-changing variant, and mirrors `DMRG2`: `trscheme` is required, and it needs a unit cell of at least two sites. + +!!! warning "Unit cell size" + `IDMRG2` throws an `ArgumentError` on a single-site `InfiniteMPS`. + Build the initial state and Hamiltonian with a unit cell of two (or more) sites instead. + +```@example groundstate_algs +J = 1.0 +g = 0.5 +X = σˣ() +Z = σᶻ() + +lattice_2 = PeriodicVector([ℂ^2, ℂ^2]) +H_inf_2 = InfiniteMPOHamiltonian( + lattice_2, + (1, 2) => -J * X ⊗ X, + (2, 3) => -J * X ⊗ X, + (1,) => -g * Z, + (2,) => -g * Z, +) +ψ₀_2 = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^2, ℂ^2]) + +ψ_i2, envs_i2, ϵ_i2 = find_groundstate( + ψ₀_2, H_inf_2, + IDMRG2(; trscheme = truncrank(16), maxiter = 5, verbosity = 0) +) +ϵ_i2 +``` + +--- + +## 4. Refine convergence with GradientGrassmann + +[`GradientGrassmann`](@ref) performs Riemannian gradient descent directly on the manifold of (finite or infinite) MPS, using an optimizer from OptimKit (`ConjugateGradient` by default via the `method` keyword). + +Chain it after `VUMPS` (or `DMRG`) with `&` to combine both regimes in one call — this is exactly what `find_groundstate`'s heuristic does once `tol` is tighter than `1e-4`: + +```@example groundstate_algs +refine = VUMPS(; tol = 1.0e-6, maxiter = 20, verbosity = 0) & + GradientGrassmann(; tol = 1.0e-10, maxiter = 50, verbosity = 0) + +ψ_g, envs_g, ϵ_g = find_groundstate(ψ₀_inf, H_inf, refine) +ϵ_g +``` + +Since `GradientGrassmann` is also a single-site algorithm, it cannot change the bond dimension either: grow it beforehand with `DMRG2`/`IDMRG2` or the `changebonds` recipes in [Controlling bond dimension](@ref howto_bond_dimension). + +--- + +## 5. Control output and tolerances + +Every algorithm accepts a `verbosity` keyword as a plain integer: + +| `verbosity` | Output | +|:-----------:|:-------| +| `0` | nothing | +| `1` | warnings only | +| `2` | convergence information | +| `3` | per-iteration information (the default) | +| `4` | everything | + +`find_groundstate` returns `(ψ, envs, ϵ)`. +`ϵ` is the final convergence-error measure (a Galerkin residual) of whichever algorithm ran last — it quantifies how well the sweeps converged, not the error in the energy itself. + +The optional third positional argument to `find_groundstate` lets you reuse `envs` from a previous call instead of recomputing it, which is useful when tightening the tolerance on a state you already optimized: + +```@example groundstate_algs +ψ_v2, envs_v2, ϵ_v2 = find_groundstate( + ψ_v, H_inf, + VUMPS(; tol = 1.0e-10, maxiter = 50, verbosity = 0), + envs_v +) +ϵ_v2 +``` + +--- + +## Where to go next + +For growing, shrinking, and inspecting bond dimension during or between these calculations, see [Controlling bond dimension](@ref howto_bond_dimension). +For background on when each algorithm applies and how it relates to the others, see [Ground-state algorithms](@ref lib_groundstate). +To see `find_groundstate` used end to end on a finite chain, start from [Your first ground state](@ref tutorial_first_groundstate); for the infinite-system counterpart, see [The thermodynamic limit](@ref tutorial_thermodynamic_limit). diff --git a/docs/src/howto/index.md b/docs/src/howto/index.md new file mode 100644 index 000000000..bc821b6ec --- /dev/null +++ b/docs/src/howto/index.md @@ -0,0 +1,84 @@ +# [How-to guides](@id howto_index) + +These pages are task recipes: short, runnable answers to "how do I do X?". +They assume you already know the basics — if you are new to MPSKit, start with [Your first ground state](@ref tutorial_first_groundstate) instead. +Each recipe below stands on its own, so feel free to jump straight to the one you need. + +## States and operators + +**[Constructing states](@ref howto_states)** — building `FiniteMPS`, `InfiniteMPS`, `WindowMPS`, and `MultilineMPS` objects. +- A finite MPS from length, physical space, and maximum bond dimension — the default random-state constructor. +- Choosing the initializer and element type — `rand`/`randn` and real vs. complex tensors. +- Per-site physical and virtual spaces — heterogeneous lattices via vectors of spaces. +- A product state (trivial virtual space) — an unentangled bond-dimension-1 state. +- From your own site tensors — wrap existing `TensorMap`s into a canonical-form MPS. +- An infinite MPS — single- and multi-site unit cells, from spaces or from tensors. +- A window MPS — a mutable finite region embedded in infinite environments. +- A multiline MPS — stacking `InfiniteMPS` rows for boundary-MPS methods. +- States with symmetries — building MPS with `Rep[G]` graded spaces. + +**[Building Hamiltonians](@ref howto_hamiltonians)** — assembling MPO Hamiltonians from local operators. +- Finite Hamiltonian from local terms — `FiniteMPOHamiltonian` from `inds => operator` pairs. +- Infinite (translation-invariant) Hamiltonian — `InfiniteMPOHamiltonian` on a unit cell. +- Converting between boundary conditions — open vs. periodic finite chains from an infinite Hamiltonian. +- A window Hamiltonian — carving a finite interval out of an infinite Hamiltonian with `WindowMPOHamiltonian`. + +## Finding ground states + +**[Ground-state algorithms](@ref howto_groundstate_algorithms)** — configuring `find_groundstate`. +- Get a ground state with defaults — letting `find_groundstate` pick an algorithm for you. +- Configure finite-system DMRG — explicit `DMRG`/`DMRG2`, and chaining them with `&`. +- Configure infinite-system algorithms — `VUMPS`, `IDMRG`, and `IDMRG2`. +- Refine convergence with GradientGrassmann — Riemannian gradient descent after a cheaper warm-up. +- Control output and tolerances — `verbosity` levels and reusing `envs` between calls. + +**[Controlling bond dimension](@ref howto_bond_dimension)** — inspecting, growing, and shrinking bond dimension. +- Inspecting the current bond dimension — `left_virtualspace`/`right_virtualspace` and `dim`. +- Growing bond dimension — `RandExpand` (no Hamiltonian needed) and `OptimalExpand`. +- Reducing bond dimension — `SvdCut` and the in-place `changebonds!`. +- Truncation schemes — `truncrank`, `trunctol`, `notrunc`, `truncspace`, and combining them with `&`. +- Growing during finite MPS optimization — `DMRG2` and the `trscheme` keyword of `find_groundstate`. +- Growing during infinite MPS optimization — `IDMRG2` and `VUMPSSvdCut`. +- Chaining algorithms — composing bond-change and ground-state algorithms with `&`. + +## Dynamics + +**[Time evolution](@ref howto_time_evolution)** — real- and imaginary-time evolution of an MPS. +- Evolve a state through one time step — `timestep` with `TDVP`. +- Evolve over a time span — `time_evolve` across a vector of time points. +- Grow the bond dimension while evolving — `TDVP2` with a mandatory `trscheme`. +- Evolve an infinite state — single-site `TDVP` on an `InfiniteMPS`. +- Imaginary-time evolution — `imaginary_evolution = true` to cool towards the ground state. +- Build a time-evolution MPO — `make_time_mpo` (`WII`, `TaylorCluster`, `WI`) plus `approximate`. + +## Measurements + +**[Computing observables](@ref howto_observables)** — extracting physical quantities from an MPS. +- Local (one-site) expectation value — `expectation_value(ψ, i => O)`. +- Multi-site (contiguous) expectation value — tensor-product operators on an index tuple. +- Energy (full-MPO expectation value) — `expectation_value(ψ, H)` for a Hamiltonian MPO. +- Two-point correlators — `correlator`, including a full correlation profile over a range. +- Energy variance as a convergence check — `variance` as a diagnostic after a ground-state search. + +**[Entanglement entropy and spectrum](@ref howto_entanglement)** — reading off entanglement from the gauge tensors. +- Entanglement entropy at a single cut — `entropy(ψ, site)`. +- Entropy profile across every cut — collecting `entropy` over all sites. +- The entanglement spectrum — `entanglement_spectrum` as a sector-resolved vector. +- Sector-resolved spectrum — indexing the spectrum by symmetry sector with `keys`/`pairs`. +- Entanglement of an infinite MPS — `entropy`/`entanglement_spectrum` per unit-cell site. +- Plotting the spectrum — the `entanglementplot` recipe (requires Plots.jl). + +## Excitations + +**[Excited states](@ref howto_excitations)** — computing energy eigenstates beyond the ground state. +- Get a single excitation gap on an infinite chain — `QuasiparticleAnsatz` at a fixed momentum. +- Scan the dispersion relation — passing a range of momenta in one call. +- Target a symmetry sector — a charged quasiparticle via the `sector` keyword. +- Excited states on a finite chain — `QuasiparticleAnsatz`, `FiniteExcited`, and the Chepiga ansätze. +- Check excitation quality — `variance` on a quasiparticle state. +- Domain-wall excitations — quasiparticles interpolating between two distinct ground states. + +## Missing a recipe? + +If the task you're after isn't listed here, please open an issue at [QuantumKitHub/MPSKit.jl](https://github.com/QuantumKitHub/MPSKit.jl/issues) describing what you're trying to do. +Concrete task descriptions make the best new recipes. diff --git a/docs/src/howto/time_evolution.md b/docs/src/howto/time_evolution.md new file mode 100644 index 000000000..36e5061d3 --- /dev/null +++ b/docs/src/howto/time_evolution.md @@ -0,0 +1,158 @@ +# [Time evolution](@id howto_time_evolution) + +The examples on this page use MPSKit.jl, MPSKitModels.jl, TensorKit.jl, and TensorKitTensors.jl. +See [Installation](@ref tutorial_installation) for how to add these packages to your environment. + +MPSKit solves the (real- or imaginary-time) Schrödinger equation `i ∂ψ/∂t = H ψ` in two ways: by projecting the equation onto the MPS tangent space at every step ([`TDVP`](@ref)/[`TDVP2`](@ref)), or by first building an approximate evolution operator as an MPO ([`make_time_mpo`](@ref)) and repeatedly applying it. +This page gives task recipes for both routes. +All examples share a single namespace: + +```@example time_evo +using MPSKit, MPSKitModels, TensorKit +using TensorKitTensors.SpinOperators: σˣ, σᶻ +``` + +--- + +## 1. Evolve a state through one time step + +[`timestep`](@ref) advances a state by a single `dt` under a Hamiltonian, using whichever [`TDVP`](@ref)-family algorithm you pass. +Its argument order is state, Hamiltonian, current time, time step, algorithm: + +```@example time_evo +L = 8 +g₀ = 0.5 +H₀ = transverse_field_ising(FiniteChain(L); g = g₀) + +ψ₀ = FiniteMPS(L, ℂ^2, ℂ^8) +ψ₀, = find_groundstate(ψ₀, H₀, DMRG(; verbosity = 0)) +expectation_value(ψ₀, 4 => σᶻ()) +``` + +Now quench to a different transverse field `g₁` and take a single step: + +```@example time_evo +g₁ = 2.0 +H₁ = transverse_field_ising(FiniteChain(L); g = g₁) + +dt = 0.05 +ψ₁, envs₁ = timestep(ψ₀, H₁, 0.0, dt, TDVP()) +expectation_value(ψ₁, 4 => σᶻ()) +``` + +`timestep` returns the updated state together with an `envs` cache; reuse `envs₁` in the next call to avoid recomputation. +An in-place `timestep!` also exists, but only for finite MPS. + +--- + +## 2. Evolve over a time span + +[`time_evolve`](@ref) steps through an explicit vector of time points instead of a single `dt`, carrying the algorithm and environments through the whole span: + +```@example time_evo +t_span = 0:dt:(4dt) +ψ_span, envs_span = time_evolve(ψ₀, H₁, t_span, TDVP(); verbosity = 0) +expectation_value(ψ_span, 4 => σᶻ()) +``` + +`t_span` need not be uniformly spaced — `time_evolve` steps pairwise between consecutive entries, so any `AbstractVector` of increasing times works. +There is no exported `time_evolve!`; use `time_evolve` and rebind the result. + +--- + +## 3. Grow the bond dimension while evolving + +Single-site `TDVP` cannot change the bond dimension: whatever `ψ₀` starts with is what it keeps. +After a quench, entanglement typically grows and a fixed bond dimension eventually becomes insufficient. + +[`TDVP2`](@ref) updates two sites at a time and truncates back down, so it can grow (or shrink) the bond dimension as it evolves. +Unlike `TDVP`, `TDVP2` requires `trscheme` — there is no default: + +```@example time_evo +ψ_tdvp2, envs_tdvp2 = timestep(ψ₀, H₁, 0.0, dt, TDVP2(; trscheme = truncrank(16))) +dim(left_virtualspace(ψ_tdvp2, 4)) +``` + +`TDVP2` only has a finite-MPS method. +For finite systems, an alternative to switching algorithms entirely is single-site `TDVP` with `alg_expand` set to a bond-expansion algorithm such as [`OptimalExpand`](@ref) (Controlled Bond Expansion, "CBE-TDVP"); see [Controlling bond dimension](@ref howto_bond_dimension) for `OptimalExpand` and other `changebonds` recipes. + +--- + +## 4. Evolve an infinite state + +Single-site `TDVP` also works directly on an `InfiniteMPS`: + +```@example time_evo +ψ₀_inf = InfiniteMPS(ℂ^2, ℂ^8) +H₀_inf = transverse_field_ising(; g = g₀) +ψ₀_inf, = find_groundstate(ψ₀_inf, H₀_inf, VUMPS(; verbosity = 0)) + +H₁_inf = transverse_field_ising(; g = g₁) +ψ₁_inf, envs₁_inf = timestep(ψ₀_inf, H₁_inf, 0.0, dt, TDVP()) +expectation_value(ψ₁_inf, 1 => σᶻ()) +``` + +`TDVP2` has no `InfiniteMPS` method, and single-site `TDVP` cannot change the bond dimension on an infinite state either. +Grow the bond dimension beforehand with `changebonds` (e.g. `OptimalExpand` or `VUMPSSvdCut`) — see [Controlling bond dimension](@ref howto_bond_dimension) — then evolve at the fixed, larger bond dimension. + +--- + +## 5. Imaginary-time evolution + +Passing `imaginary_evolution = true` evolves under `exp(-H dt)` instead of `exp(-iH dt)`, using the same real `dt`: + +```@example time_evo +ψ_im, envs_im = timestep(ψ₀, H₁, 0.0, dt, TDVP(); imaginary_evolution = true) +norm(ψ_im) +``` + +Imaginary-time evolution renormalizes the state at every step, so `norm(ψ_im)` stays `1` regardless of how the un-normalized weight would otherwise change. +Repeated imaginary-time steps damp excited-state components faster than the ground state, so this is often used as a (slower) alternative to `find_groundstate` for driving a state towards the ground state of `H₁`. + +`time_evolve` accepts the same `imaginary_evolution` keyword for a span of imaginary-time steps. + +--- + +## 6. Build a time-evolution MPO + +When `H` is time-independent and `dt` is fixed, an alternative to repeated `timestep` calls is to build the evolution operator once as an MPO with [`make_time_mpo`](@ref), then apply it repeatedly with [`approximate`](@ref). + + +[`WII`](@ref) builds a low-order MPO approximation and works for both finite and infinite Hamiltonians: + +```@example time_evo +O = make_time_mpo(H₁, dt, WII()) +``` + +[`TaylorCluster`](@ref) gives higher-order control via its `N` keyword (and accepts an additional `tol` keyword in `make_time_mpo`): + +```@example time_evo +O_taylor = make_time_mpo(H₁, dt, TaylorCluster(; N = 2); tol = 1.0e-10) +``` + +[`WI`](@ref) is a ready-made first-order `TaylorCluster` constant, so it is used as a value, not called as a constructor: + +```@example time_evo +O_wi = make_time_mpo(H₁, dt, WI) +``` + +Applying the MPO to a finite state uses [`approximate`](@ref) with a finite ground-state-style algorithm such as [`DMRG2`](@ref), which returns a 3-tuple including the final convergence error: + +```@example time_evo +ψ_mpo, envs_mpo, ϵ_mpo = approximate( + ψ₀, (O, ψ₀), DMRG2(; trscheme = truncrank(16), verbosity = 0) +) +expectation_value(ψ_mpo, 4 => σᶻ()) +``` + +Repeat the `approximate` call with the same `O` for successive time steps to build up a longer evolution. +Imaginary-time MPOs are built the same way, by passing `imaginary_evolution = true` to `make_time_mpo`; a real `dt` is promoted internally, so no manual complex conversion is needed. + +--- + +## Where to go next + +For choosing and configuring ground-state algorithms to prepare the pre-quench state, see [Ground-state algorithms](@ref howto_groundstate_algorithms). +For growing or shrinking bond dimension between or during evolution steps, see [Controlling bond dimension](@ref howto_bond_dimension). +For extracting expectation values and correlators from the evolved state, see [Computing observables](@ref howto_observables). +For background on the TDVP and time-evolution-MPO approaches and how they relate, see [Time evolution](@ref lib_time_evolution). diff --git a/docs/src/tutorials/excitations.md b/docs/src/tutorials/excitations.md new file mode 100644 index 000000000..fa77470da --- /dev/null +++ b/docs/src/tutorials/excitations.md @@ -0,0 +1,116 @@ +# [Quasiparticle excitations](@id tutorial_excitations) + +The previous tutorials ended with a ground state: the lowest-energy state of the transverse-field Ising model, first on a finite chain and then directly in [the thermodynamic limit](@ref tutorial_thermodynamic_limit). +The natural next question is what lies *above* it: how much energy does it cost to excite the system? +For a translation-invariant chain the answer is organized by momentum — for each momentum ``k`` there is a lowest excitation energy ``\Delta E(k)``, and the resulting curve is the **dispersion relation** of the model. +Its minimum over all momenta is the **energy gap**, one of the most basic characterizations of a quantum phase. + + +In this tutorial we compute the dispersion relation of the infinite transverse-field Ising chain with MPSKit's quasiparticle ansatz, and finish with a plot of ``\Delta E(k)`` across the Brillouin zone — compared against the exact solution. + +## Loading the packages + +As in the previous tutorials, every code block on this page shares one Julia session, so we load the packages once. + +```@example excitations +using MPSKit, MPSKitModels, TensorKit +using Plots +``` + +## 1. Find the ground state + +Excitations are computed *on top of* a ground state, so the first step is the calculation you already know from [The thermodynamic limit](@ref tutorial_thermodynamic_limit): build the infinite Hamiltonian, make a random `InfiniteMPS`, and converge it with `VUMPS`. + +This time we set the field to `g = 2.0`, deep in the paramagnetic phase, where the model is **gapped**: the lowest excitation costs a finite amount of energy, which is exactly what we want to measure. + + +```@example excitations +g = 2.0 +H = transverse_field_ising(; g) +ψ₀ = InfiniteMPS(ℂ^2, ℂ^12) +ψ, envs, ϵ = find_groundstate(ψ₀, H, VUMPS(; verbosity = 0)) +``` + +We keep all three return values this time: the optimized state `ψ` and the environments `envs` both feed directly into the excitation calculation below, so nothing has to be recomputed. + +## 2. One excitation at one momentum + +The **quasiparticle ansatz** builds an excited state directly on top of the uniform ground state. +The idea is simple to picture: take the converged ground state and perturb it locally, replacing the tensor at one site with a new one that we get to optimize. +Because the chain is infinite and translation invariant, we do not place this perturbation at any particular site; instead we superpose it across *all* sites with a plane-wave phase, which gives the excitation a definite momentum ``k``. +Optimizing the perturbation then yields the lowest excited state at that momentum. + + +The call is [`excitations`](@ref) with the [`QuasiparticleAnsatz`](@ref) algorithm, a momentum (a real number, in radians per site), and the ground state with its environments. +Let us ask for the excitation at the edge of the Brillouin zone, ``k = \pi``: + +```@example excitations +E, ϕ = excitations(H, QuasiparticleAnsatz(), π, ψ, envs) +E +``` + +Two things to note about the return values: + +- `E` is a *vector* of excitation energies, of length `num` — the keyword controlling how many excitations to compute at this momentum, which defaults to `num = 1`, so here it has a single entry. +- The entries of `E` are energies **above the ground state** — gaps at this momentum, not total energies. The ground-state energy is subtracted internally, so you can read them off directly. + +The second return value `ϕ` holds the corresponding quasiparticle states, which can be used for further post-processing; we will not need them in this tutorial. + +## 3. The full dispersion + +To trace out the whole dispersion relation we simply pass a *range* of momenta instead of a single number. +By symmetry it is enough to scan from ``0`` to ``\pi``, and we use 16 points to keep the runtime modest. + + +```@example excitations +momenta = range(0, π, 16) +Es, ϕs = excitations(H, QuasiparticleAnsatz(), momenta, ψ, envs; verbosity = 0) +size(Es) +``` + +With a range of momenta the energies come back as a matrix of size `(length(momenta), num)` — here `(16, 1)`, one row per momentum and one column because we kept the default `num = 1`. +We pass `verbosity = 0` to silence the progress line this method otherwise prints for every momentum. +The momenta are independent of one another, so MPSKit works on them in parallel by default. + +## 4. Plot the dispersion + +Now for the payoff. +This particular model is exactly solvable, so we can plot our numerical dispersion right on top of the known answer: + +```math +\Delta E(k) = 2\sqrt{1 + g^2 - 2 g \cos k}. +``` + + + +For this Hermitian problem the computed energies come back as real numbers; the `real.(...)` below is a harmless safeguard for the general case, where the eigenvalue solver may return a complex number type with numerically vanishing imaginary parts. + +```@example excitations +k_exact = range(0, π, 200) +ΔE_exact = @. 2 * sqrt(1 + g^2 - 2g * cos(k_exact)) +plot(k_exact, ΔE_exact; label = "exact", xlabel = "momentum k", ylabel = "ΔE(k)", title = "TFIM dispersion (g = $g)") +scatter!(momenta, real.(Es); label = "quasiparticle ansatz (D = 12)") +``` + +The 16 computed points fall right on the exact curve. +The dispersion rises monotonically from ``k = 0`` to ``k = \pi``, so its minimum — the gap — sits at zero momentum, where the exact value is ``\Delta E(0) = 2(g - 1)``. +Our first matrix entry is precisely that point, so we can close with a numerical check: + +```@example excitations +real(Es[1, 1]), 2 * (g - 1) +``` + +A ground state at bond dimension 12 plus a variational quasiparticle on top reproduces the exact gap of the model — that is the quasiparticle ansatz working as intended. + + +## Where to go next + +You have computed a full dispersion relation on top of an infinite ground state and read off the energy gap. + +The [`excitations`](@ref) entry point can do considerably more than what we used here: it can target excitations carrying a nontrivial symmetry charge, build topological (domain-wall) excitations that interpolate between two different ground states, and compute excited states of *finite* chains, where momentum is no longer a good quantum number and different algorithms take over. +All of these are recipes in [Excited states](@ref howto_excitations). +For what each excitation algorithm actually does and when to choose it, see the library reference [Excitations](@ref lib_excitations). + + diff --git a/docs/src/tutorials/thermodynamic_limit.md b/docs/src/tutorials/thermodynamic_limit.md index 3f51bb745..0be68f69e 100644 --- a/docs/src/tutorials/thermodynamic_limit.md +++ b/docs/src/tutorials/thermodynamic_limit.md @@ -129,9 +129,6 @@ In the thermodynamic limit the symmetry breaking is genuine — the two opposite You have now run the same TFIM calculation twice — once at finite size, once directly at `L = ∞` — and seen how little the code had to change. From here you can go beyond ground states. -A natural next step is to compute the excitations above this infinite ground state (the model's quasiparticle spectrum), or to exploit the symmetries of the model to make the calculation cheaper and more accurate. - - - +A natural next step is to [compute the excitations above this infinite ground state](@ref tutorial_excitations) (the model's quasiparticle spectrum), or to [exploit the symmetries of the model](@ref tutorial_using_symmetries) to make the calculation cheaper and more accurate. To go deeper on the individual steps used here, see [Constructing states](@ref howto_states), [Controlling bond dimension](@ref howto_bond_dimension), and [Entanglement entropy and spectrum](@ref howto_entanglement); the algorithm reference is [Ground-state algorithms](@ref lib_groundstate). diff --git a/docs/src/tutorials/time_evolution.md b/docs/src/tutorials/time_evolution.md new file mode 100644 index 000000000..cff01897a --- /dev/null +++ b/docs/src/tutorials/time_evolution.md @@ -0,0 +1,138 @@ +# [A quantum quench](@id tutorial_time_evolution) + +The previous tutorials computed ground states — static snapshots of a model at its lowest energy. +This tutorial adds the time axis: we take a state that is *not* an eigenstate of its Hamiltonian and watch it evolve under the Schrödinger equation, + +```math +|\psi(t)\rangle = e^{-iHt}\,|\psi(0)\rangle , +``` + +tracking one local observable as a function of time. +The protocol we use is the simplest and most common one in the field, a *global quench*, and the workhorse algorithm is [`TDVP`](@ref), the time-dependent variational principle, driven one step at a time through [`timestep`](@ref). + +We stay with the transverse-field Ising model from [Your first ground state](@ref tutorial_first_groundstate), so the model-building and ground-state steps below should look familiar. + +## Loading the packages + +Every code block on this page shares one Julia session, so we load the packages once. +As before, the model comes from MPSKitModels, the local spin operators from TensorKitTensors, and `Plots` draws the final figure. + +```@example time-evolution +using MPSKit, MPSKitModels, TensorKit +using TensorKitTensors.SpinOperators: σˣ, σᶻ +using Plots +``` + +## 1. Prepare the initial state + +Time evolution needs a definite starting state, and the standard choice is the ground state of some Hamiltonian. +We take a chain of `L = 12` sites with a transverse field `g₀ = 0.5` — the ordered side of the model — and find its ground state exactly as in [the first tutorial](@ref tutorial_first_groundstate), silencing the convergence log with `verbosity = 0`. +The bond dimension `D = 16` is comfortably large for a ground state of this size; we will see below why time evolution wants more headroom than a ground-state calculation. + +```@example time-evolution +L = 12 +D = 16 +g₀ = 0.5 +H₀ = transverse_field_ising(FiniteChain(L); g = g₀) +ψ₀ = FiniteMPS(L, ℂ^2, ℂ^D) +ψ, = find_groundstate(ψ₀, H₀, DMRG(; verbosity = 0)) +nothing # hide +``` + +The observable we will track through the evolution is the transverse magnetization ``\langle\sigma^x\rangle`` at the middle of the chain, away from the open ends. + +We measure its baseline value in the pre-quench ground state: + +```@example time-evolution +i_mid = L ÷ 2 +real(expectation_value(ψ, i_mid => σˣ())) +``` + +## 2. The quench + +A *global quench* is the sudden change of a parameter of the Hamiltonian, everywhere at once: we prepare the ground state of `H₀`, then at `t = 0` switch the Hamiltonian to a different `H₁` and let the state evolve under it. +Because `ψ` is an eigenstate of `H₀` but not of `H₁`, it is no longer stationary — the quench injects energy into the system, and nontrivial dynamics follows. + +Our quench takes the transverse field from `g₀ = 0.5` all the way across the phase transition to `g₁ = 2.0`: + +```@example time-evolution +g₁ = 2.0 +H₁ = transverse_field_ising(FiniteChain(L); g = g₁) +``` + + + +## 3. A single time step + +The elementary move of real-time evolution in MPSKit is [`timestep`](@ref), which advances a state by one small increment `dt`. +Its arguments are, in order: the state, the Hamiltonian to evolve under, the current time, the step size, and the algorithm. +Here we take the very first step, from `t = 0.0` to `t = dt`, with single-site [`TDVP`](@ref): + +```@example time-evolution +dt = 0.05 +ψ_t, envs = timestep(ψ, H₁, 0.0, dt, TDVP()) +real(expectation_value(ψ_t, i_mid => σˣ())) +``` + +TDVP integrates the Schrödinger equation projected onto the space of MPS with the current bond dimension, which is why it slots so naturally into an MPS workflow. + +`timestep` returns two things: + +- `ψ_t` — the evolved state at time `dt` (a new state; `ψ` is left untouched). +- `envs` — the environments, cached partial contractions belonging to the new state and `H₁`. + +The `envs` are the reason evolution loops are cheap to keep running: passing them back into the next `timestep` call lets it start from the cached contractions instead of recomputing them from scratch. + +The transverse magnetization has already moved slightly away from its `t = 0` value — the state is on its way. + +## 4. Evolving in a loop + +Real-time evolution is nothing more than this single step, repeated. +We already took step 1 above, so the loop below performs the remaining steps, up to `n_steps = 40` in total (a final time of `t = 2.0`), recording the transverse magnetization at the middle of the chain after every step. +Note how each iteration feeds the previous `envs` back in as the optional last argument, and how the current time `(n - 1) * dt` advances with the loop. + +```@example time-evolution +n_steps = 40 +times = (0:n_steps) .* dt +m = zeros(n_steps + 1) +m[1] = real(expectation_value(ψ, i_mid => σˣ())) # t = 0, before the quench dynamics +m[2] = real(expectation_value(ψ_t, i_mid => σˣ())) # t = dt, from the single step above +for n in 2:n_steps + global ψ_t, envs + ψ_t, envs = timestep(ψ_t, H₁, (n - 1) * dt, dt, TDVP(), envs) + m[n + 1] = real(expectation_value(ψ_t, i_mid => σˣ())) +end +m[end] +``` + +(The `global` keyword is needed because the loop rebinds `ψ_t` and `envs`, which live outside it; inside a function you would not need it.) + +Writing the loop by hand like this keeps every moving part visible, which is the point of a tutorial. +For production use, [`time_evolve`](@ref) wraps exactly this loop and steps through a whole vector of time points in one call — see [Time evolution](@ref howto_time_evolution). + +## 5. Magnetization over time + +The payoff: the transverse magnetization at the middle of the chain, as a function of time after the quench. + +```@example time-evolution +plot(times, m; + xlabel = "t", ylabel = "⟨σˣ⟩ at site $i_mid", + label = "TDVP, D = $D", title = "TFIM transverse magnetization after a quench") +``` + +The curve shows how the observable responds to the sudden change in the field: starting from its pre-quench value, ``\langle\sigma^x\rangle`` relaxes towards a new value set by `H₁`, with oscillations along the way. + + +!!! warning "Fixed bond dimension means finite reach in time" + Single-site TDVP keeps the bond dimension fixed at whatever the initial state has. + After a quench, however, the entanglement of the evolving state grows with time, so a fixed bond dimension can only follow the true dynamics faithfully up to some finite time — beyond it, the simulation quietly loses accuracy rather than failing loudly. + + The practical checks and remedies — two-site [`TDVP2`](@ref), which grows the bond dimension as it truncates, and bond-expansion options for single-site TDVP — are collected in [Time evolution](@ref howto_time_evolution). + +## Where to go next + +You have run your first dynamics simulation: prepare a ground state, quench the Hamiltonian, step the state forward in time, and read off an observable at every step. + +The natural reference for everything this page glossed over is the [Time evolution](@ref howto_time_evolution) how-to: evolving over a time span in one call, growing the bond dimension during evolution, imaginary time, and evolving infinite states. +Speaking of which — everything here was done on a finite chain, but `timestep` works just as well on the `InfiniteMPS` states introduced in [The thermodynamic limit](@ref tutorial_thermodynamic_limit). +And for measuring more than a single local magnetization on the evolved states, see [Computing observables](@ref howto_observables). diff --git a/docs/src/tutorials/using_symmetries.md b/docs/src/tutorials/using_symmetries.md new file mode 100644 index 000000000..6a394f676 --- /dev/null +++ b/docs/src/tutorials/using_symmetries.md @@ -0,0 +1,197 @@ + + +# [Using symmetries](@id tutorial_using_symmetries) + +In [Your first ground state](@ref tutorial_first_groundstate) and [The thermodynamic limit](@ref tutorial_thermodynamic_limit) we treated the transverse-field Ising model (TFIM) as a generic spin chain. +But the TFIM is not generic: it has a symmetry, and in this tutorial we teach MPSKit about it. + +Recall the Hamiltonian, + +```math +H = -J\left(\sum_{\langle i,j\rangle} \sigma^z_i\,\sigma^z_j + g\sum_i \sigma^x_i\right), +``` + +and consider the *global spin flip* ``P = \prod_i \sigma^x_i``, which flips every spin at once. +Conjugating by ``P`` sends ``\sigma^z_i \to -\sigma^z_i``, so the interaction term ``\sigma^z_i\sigma^z_j`` picks up two minus signs and is unchanged, while the field term ``\sigma^x_i`` commutes with ``P`` trivially. +Hence ``H`` commutes with ``P``. +Since ``P^2 = 1``, this is a ``\mathbb{Z}_2`` symmetry, and every eigenstate of ``H`` can be labelled by a parity quantum number: *even* (``P = +1``) or *odd* (``P = -1``). + +MPSKit, through the TensorKit tensor backend, can bake this symmetry directly into the tensors of the MPS. +Doing so buys you two things. +First, the tensors become **block-sparse**: at the same total bond dimension the computer multiplies smaller dense blocks, which is faster. +Second, every state you compute carries an explicit **sector label**, so "the lowest odd-parity excitation" becomes something you can ask for directly. +This tutorial demonstrates both, by redoing the finite-chain TFIM calculation once without and once with the symmetry. + +## Loading the packages + +Every code block on this page shares one Julia session, so we load the packages once. +`Z2Irrep` and `Z2Space`, the symmetry-aware building blocks used below, come from TensorKit. + +```@example using-symmetries +using MPSKit, MPSKitModels, TensorKit +``` + +## 1. Recap: the ground state without symmetry + +We start from the workflow of the first tutorial: a chain of `L = 16` sites, a random `FiniteMPS`, and a DMRG ground-state search. +Two small changes from before: we set the field to `g = 2.0`, and we use a total bond dimension of 16. + + +Why `g = 2.0`? +This puts us deep in the paramagnetic phase, where the ground state respects the spin-flip symmetry. +That matters for what comes next: an MPS built from symmetric tensors lives in exactly one parity sector and *cannot* spontaneously break the symmetry, so a fair comparison needs a point where the true ground state is symmetric to begin with. + +```@example using-symmetries +L = 16 +H = transverse_field_ising(FiniteChain(L); g = 2.0) +ψ₀ = FiniteMPS(L, ℂ^2, ℂ^16) +ψ, envs, ϵ = find_groundstate(ψ₀, H, DMRG(; verbosity = 0)) +E = expectation_value(ψ, H) +``` + +This is our reference number: the ground-state energy computed with plain, symmetry-oblivious tensors. + +## 2. The same model, with the symmetry made explicit + +To exploit the symmetry we change two lines: the Hamiltonian and the initial state. + +For the Hamiltonian, we pass the symmetry as an extra first argument. +`transverse_field_ising(Z2Irrep, ...)` builds the *same* Hamiltonian as before, but out of tensors that manifestly commute with the spin flip: + +```@example using-symmetries +H_Z2 = transverse_field_ising(Z2Irrep, FiniteChain(L); g = 2.0) +``` + +For the state, the plain spaces `ℂ^2` and `ℂ^16` are replaced by *graded* spaces that keep track of parity: + +```@example using-symmetries +ψ₀_Z2 = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => 8, 1 => 8)) +``` + +The syntax reads as a list of `sector => dimension` pairs, where sector `0` is the even (``P = +1``) irrep of ``\mathbb{Z}_2`` and sector `1` is the odd (``P = -1``) one: + +- The physical space `Z2Space(0 => 1, 1 => 1)` is the familiar two-dimensional spin-1/2 site, now split into its symmetry content: one even state and one odd state. + +- The virtual space `Z2Space(0 => 8, 1 => 8)` says the bond carries 8 states of even parity and 8 of odd parity — 16 in total, matching the `ℂ^16` of the plain run, so the two calculations have exactly the same variational power. + +From here the workflow is unchanged: + +```@example using-symmetries +ψ_Z2, envs_Z2, ϵ_Z2 = find_groundstate(ψ₀_Z2, H_Z2, DMRG(; verbosity = 0)) +E_Z2 = expectation_value(ψ_Z2, H_Z2) +``` + +Both runs found the same ground state, and the two energies agree to numerical precision: + +```@example using-symmetries +E, E_Z2 +``` + +!!! note "Same physics, different bookkeeping" + Nothing about the model changed — only the way its tensors are stored. + The symmetric calculation restricts the search to states of definite (here: even) parity, and stores only the tensor blocks the symmetry allows to be nonzero. + +## 3. The payoff, part 1: block-sparse tensors + +Where did the symmetry go? +Into the *structure* of the state. +Ask for the virtual space at the central bond and you no longer get an anonymous `ℂ^16`, but a space that knows its sector decomposition: + +```@example using-symmetries +V = left_virtualspace(ψ_Z2, L ÷ 2) +``` + +Its total dimension is still 16: + +```@example using-symmetries +dim(V) +``` + +The same sector labels show up in every quantity derived from the state. +The entanglement spectrum at the central cut, for instance, now comes back resolved by sector — compare [Entanglement entropy and spectrum](@ref howto_entanglement), where the same call on an unsymmetric state produced a single `Trivial()` block: + +```@example using-symmetries +spectrum = entanglement_spectrum(ψ_Z2, L ÷ 2) +collect(keys(spectrum)) +``` + +Iterating `pairs` gives each sector together with its singular values: + +```@example using-symmetries +collect(pairs(spectrum)) +``` + + +This block structure is where the speedup comes from: instead of multiplying one dense 16-dimensional bond index, the computer multiplies two independent blocks of roughly half that size, and the forbidden matrix elements between the sectors are never stored or touched at all. +At bond dimension 16 the difference is negligible, but the saving grows with the bond dimension and with the size of the symmetry group. + +## 4. The payoff, part 2: sectors label the physics + +The sector labels are not just an implementation detail — they classify the eigenstates of ``H``, and MPSKit lets you target a sector directly. + + +In the paramagnetic phase the lowest excitation of the TFIM is, roughly speaking, a single flipped spin. +Flipping one spin changes the parity of the state, so this excitation lives in the *odd* sector — a different sector than the (even) ground state. + +The [`excitations`](@ref) function computes excited states on top of a converged ground state; on a finite chain it takes the Hamiltonian, an algorithm, the ground state, and its environments, and returns energies measured *above* the ground state. +By default it searches the trivial (even) sector: + +```@example using-symmetries +Es_even, ϕs_even = excitations(H_Z2, QuasiparticleAnsatz(), ψ_Z2, envs_Z2; num = 1) +Es_even[1] +``` + +The `sector` keyword redirects the search to the odd sector: + +```@example using-symmetries +Es_odd, ϕs_odd = excitations( + H_Z2, QuasiparticleAnsatz(), ψ_Z2, envs_Z2; + num = 1, sector = Z2Irrep(1) +) +Es_odd[1] +``` + +The odd-sector excitation is indeed the lower one: + +```@example using-symmetries +Es_odd[1] < Es_even[1] +``` + +Without the symmetry built into the tensors, this question could not even be posed: the plain calculation of Section 1 has no notion of parity to select on. +More ways to use `excitations` — dispersion relations, other algorithms, infinite chains — are collected in [Excited states](@ref howto_excitations). + +## Where to go next + +You have run the flagship TFIM calculation with its ``\mathbb{Z}_2`` symmetry made explicit: the same physics at the same total bond dimension, but with block-sparse tensors and sector labels on everything the calculation produces. + +The same syntax scales up to larger symmetry groups, where the payoff grows. +For a U(1) symmetry (particle number, magnetization) the graded spaces list integer or half-integer charges instead of parities — worked constructions are in [Constructing states](@ref howto_states), Section 9. + +For non-abelian symmetries such as SU(2) the gains are more dramatic still, because each symmetric block then represents an entire multiplet of states; the spin-1 Haldane chain and XXZ Heisenberg pages in the examples gallery show this in action. + +To continue the tutorial track, [Quasiparticle excitations](@ref tutorial_excitations) develops the excitation calculation of Section 4 into a full dispersion relation; the recipe collection for excited states is [Excited states](@ref howto_excitations), and the one for building symmetric states is [Constructing states](@ref howto_states). + +