diff --git a/docs/make.jl b/docs/make.jl index 3e7b253eb..63ea39962 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -46,7 +46,7 @@ mathengine = MathJax3( # examples pages examples_optimization = joinpath.( - ["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard", "c4v_ctmrg"], Ref("index.md") + ["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard", "c4v_ctmrg", "triangular_ctmrg"], Ref("index.md") ) examples_time_evolution = joinpath.( ["heisenberg_su", "hubbard_su", "j1j2_su"], Ref("index.md") diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 7196cfa2a..5433cdec7 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -22,6 +22,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- Triangular CTMRG for PEPS and classical partition functions, with onsite and nearest-neighbor PEPS observables and native optimization using `FixedPointGradient` - `TensorKit.storagetype` for `InfinitePEPS`, `InfinitePEPO`, `InfinitePartitionFunction`, `InfiniteSquareNetwork` and `CTMRGEnv` (#429) ### Changed diff --git a/docs/src/examples/triangular_ctmrg/index.md b/docs/src/examples/triangular_ctmrg/index.md new file mode 100644 index 000000000..e05f36ff9 --- /dev/null +++ b/docs/src/examples/triangular_ctmrg/index.md @@ -0,0 +1,165 @@ +```@meta +EditURL = "../../../../examples/triangular_ctmrg/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/PEPSKit.jl/gh-pages?filepath=dev/examples/triangular_ctmrg/main.ipynb) +[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](https://nbviewer.jupyter.org/github/QuantumKitHub/PEPSKit.jl/blob/gh-pages/dev/examples/triangular_ctmrg/main.ipynb) +[![](https://img.shields.io/badge/download-project-orange)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumKitHub/PEPSKit.jl/examples/tree/gh-pages/dev/examples/triangular_ctmrg) + +````julia +using Markdown #hide +```` + +# [Triangular CTMRG and optimization](@id examples_triangular_ctmrg) + +This example contracts and optimizes an antiferromagnetic spin-1/2 Heisenberg PEPS +on the triangular lattice using [`InfinitePEPSTriangular`](@ref), +[`CTMRGEnvTriangular`](@ref), and the native [`fixedpoint`](@ref) interface. +A single-site tensor and fixed 120° sublattice rotations represent a three-sublattice +state. We start from a seeded random real tensor with `D = 2` and a complex +environment with `χ = 32`. Check convergence in `χ` before drawing quantitative +conclusions from the energy. + +From the repository root, prepare the example environment once with +`julia --project=examples -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'`, +then run `julia --project=examples --threads=2 examples/triangular_ctmrg/main.jl`. + +````julia +using Random +using TensorKit, PEPSKit +using TensorKitTensors.SpinOperators +Random.seed!(123456789); +nothing #hide +```` + +## Unit cell and nearest neighbors + +PEPSKit indexes the unit cell as a Julia matrix: `[r, c]` means row `r`, column `c`, +starting at 1. For `[A B C; D E F; G H I]`, `[2, 1]` is `D` and `[1, 2]` is `B`. +The six virtual legs are ordered `NW, NE, E, SE, SW, W`, with offsets +`(-1, 0), (-1, 1), (0, 1), (1, 0), (1, -1), (0, -1)` respectively. + +To count every nearest-neighbor bond once, we use the following three offsets: + +| Bond family | Offset `(Δr, Δc)` | Virtual leg | Relative rotation | +|:------------|:------------------|:------------|:------------------| +| 0° | `(0, 1)` | E | `2π/3` | +| 60° | `(1, -1)` | SW | `-4π/3` | +| 120° | `(1, 0)` | SE | `-2π/3` | + +The bond-family angles are unoriented; the relative rotations refer to the directed +offsets in the table. With `R(θ) = exp(-im * θ * S_y)`, the physical tensor at +`[r, c]` is `R((r - 1)q[1] + (c - 1)q[2]) * A`, where `q = (-2π/3, 2π/3)`. +This gives `[A B C; C A B; B C A]`, with `B = R(2π/3) * A` and +`C = R(4π/3) * A`, up to irrelevant spin-1/2 tensor phases. It is the PEPSKit +version of the variPEPS pattern `[A B C; B C A; C A B]` with `q = (2π/3, 2π/3)`. + +We optimize `A` in a single-site unit cell. All three sublattices share its +variational parameters, and the wavevector stays fixed during optimization. + +````julia +D = 2 +χ = 32 +P, V = ℂ^2, ℂ^D +peps₀ = InfinitePEPSTriangular(randn, Float64, P, V; unitcell = (1, 1)); +nothing #hide +```` + +## Hamiltonian and CTMRG + +The Heisenberg interaction is invariant under a common rotation of both spins, so +only their relative rotation enters each bond. We absorb it into the gate on the +second site: `hδ = (id(P) ⊗ R(q⋅δ))' * h * (id(P) ⊗ R(q⋅δ))`. The diagonal +angle `-4π/3` gives the same gate as `2π/3`. This lets the existing +[`LocalOperator`](@ref) and `expectation_value` interfaces evaluate the physical +three-sublattice energy using one tensor and three bond gates. + +We use antiferromagnetic coupling `J = 1`. Keep the neighbor coordinates unwrapped, +including those outside the unit cell: the contractions apply periodicity while +preserving which bond is intended. Triangular expectation values support onsite and nearest-neighbor +terms in all three bond families. + +````julia +h = S_x_S_x() + S_y_S_y() + S_z_S_z() +offsets = (CartesianIndex(0, 1), CartesianIndex(1, -1), CartesianIndex(1, 0)) +q = (-2π / 3, 2π / 3) +rotation(θ) = TensorMap([cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)], P, P) +gates = map(offsets) do δ + U = id(P) ⊗ rotation(q[1] * δ[1] + q[2] * δ[2]) + return U' * h * U +end +H = LocalOperator( + physicalspace(peps₀), + [(I, I + δ) => gate for I in CartesianIndices(size(peps₀)) for (δ, gate) in zip(offsets, gates)], +); +nothing #hide +```` + +The environment constructor takes the doubled virtual space `V ⊗ V'` of the PEPS +norm network and an environment space `ℂ^χ`. We initialize it with `ones`. +The PEPS tensor and the rotated Heisenberg gates are real, but the current gauge +fixing for fixed-point differentiation requires a complex environment. + +[`leading_boundary`](@ref) selects triangular CTMRG from the environment; here we +also name it explicitly, so the same settings can be passed to the optimizer. +We use the default `projector_alg = :full` for the corner construction. +The alternative `:twothirds` uses two thirds of the environment. `conditioning = true` (the default) conditions the direct +Ta/Tb edge split. It does not toggle corner conditioning. + +````julia +env_initial = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', ℂ^χ; unitcell = size(peps₀)) +boundary_alg = (; + alg = :SimultaneousCTMRGTriangular, tol = 1.0e-7, maxiter = 100, + projector_alg = :full, trunctype = :FixedSpaceTruncation, +) +env₀, info_ctmrg = leading_boundary(env_initial, peps₀; boundary_alg...); +@show info_ctmrg.converged info_ctmrg.convergence_error info_ctmrg.contraction_metrics; +nothing #hide +```` + +Convergence compares the corner and edge singular-value spectra. The truncation +metrics describe the local SVD errors, not the error in the energy. Check CTMRG +convergence and increase `χ` when assessing the accuracy of an entangled calculation. + +`expectation_value` sums all terms in `H`, so this Hamiltonian gives the energy per +unit cell. Divide by the number of sites to obtain the energy per site: + +````julia +E₀ = real(expectation_value(peps₀, H, env₀)) +@show E₀ / length(peps₀); +nothing #hide +```` + +## Native optimization + +We use L-BFGS with [`FixedPointGradient`](@ref) and a GMRES linear solver to +differentiate the CTMRG fixed point. Triangular optimization currently requires +`precondition_alg = nothing` +and does not support symmetrization. Alternatively, `gradient_alg = nothing` +differentiates through the CTMRG iterations. + +PEPSKit adapts the CTMRG and gradient tolerances during optimization, using looser +solves when the gradient is large. Check the final gradient norm and increase the +optimization budget if necessary. +Convergence of CTMRG at a fixed `χ` does not establish convergence in `χ`; repeat +the calculation at larger environment dimensions to assess the energy. + +````julia +peps, env, E, info_opt = fixedpoint( + H, peps₀, env₀; + boundary_alg, + gradient_alg = (; alg = :FixedPointGradient, solver_alg = (; alg = :GMRES)), + precondition_alg = nothing, + optimizer_alg = (; alg = :LBFGS, tol = 1.0e-6, maxiter = 100), verbosity = 3, +); +@show E / length(peps) info_opt.gradnorms[end]; +nothing #hide +```` + +The `D = 2` ansatz allows entanglement beyond the 120° product state, whose +energy is `-3/8` per site. The returned `E` is the energy per unit cell, which here +contains a single site. Increase both `D` and `χ` to study the quantum ground state. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* diff --git a/docs/src/examples/triangular_ctmrg/main.ipynb b/docs/src/examples/triangular_ctmrg/main.ipynb new file mode 100644 index 000000000..6d2c3df97 --- /dev/null +++ b/docs/src/examples/triangular_ctmrg/main.ipynb @@ -0,0 +1,246 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Markdown #hide" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Triangular CTMRG and optimization\n", + "\n", + "This example contracts and optimizes an antiferromagnetic spin-1/2 Heisenberg PEPS\n", + "on the triangular lattice using `InfinitePEPSTriangular`,\n", + "`CTMRGEnvTriangular`, and the native `fixedpoint` interface.\n", + "A single-site tensor and fixed 120° sublattice rotations represent a three-sublattice\n", + "state. We start from a seeded random real tensor with `D = 2` and a complex\n", + "environment with `χ = 32`. Check convergence in `χ` before drawing quantitative\n", + "conclusions from the energy.\n", + "\n", + "From the repository root, prepare the example environment once with\n", + "`julia --project=examples -e 'using Pkg; Pkg.develop(path=\".\"); Pkg.instantiate()'`,\n", + "then run `julia --project=examples --threads=2 examples/triangular_ctmrg/main.jl`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Random\n", + "using TensorKit, PEPSKit\n", + "using TensorKitTensors.SpinOperators\n", + "Random.seed!(123456789);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Unit cell and nearest neighbors\n", + "\n", + "PEPSKit indexes the unit cell as a Julia matrix: `[r, c]` means row `r`, column `c`,\n", + "starting at 1. For `[A B C; D E F; G H I]`, `[2, 1]` is `D` and `[1, 2]` is `B`.\n", + "The six virtual legs are ordered `NW, NE, E, SE, SW, W`, with offsets\n", + "`(-1, 0), (-1, 1), (0, 1), (1, 0), (1, -1), (0, -1)` respectively.\n", + "\n", + "To count every nearest-neighbor bond once, we use the following three offsets:\n", + "\n", + "| Bond family | Offset `(Δr, Δc)` | Virtual leg | Relative rotation |\n", + "|:------------|:------------------|:------------|:------------------|\n", + "| 0° | `(0, 1)` | E | `2π/3` |\n", + "| 60° | `(1, -1)` | SW | `-4π/3` |\n", + "| 120° | `(1, 0)` | SE | `-2π/3` |\n", + "\n", + "The bond-family angles are unoriented; the relative rotations refer to the directed\n", + "offsets in the table. With `R(θ) = exp(-im * θ * S_y)`, the physical tensor at\n", + "`[r, c]` is `R((r - 1)q[1] + (c - 1)q[2]) * A`, where `q = (-2π/3, 2π/3)`.\n", + "This gives `[A B C; C A B; B C A]`, with `B = R(2π/3) * A` and\n", + "`C = R(4π/3) * A`, up to irrelevant spin-1/2 tensor phases. It is the PEPSKit\n", + "version of the variPEPS pattern `[A B C; B C A; C A B]` with `q = (2π/3, 2π/3)`.\n", + "\n", + "We optimize `A` in a single-site unit cell. All three sublattices share its\n", + "variational parameters, and the wavevector stays fixed during optimization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "D = 2\n", + "χ = 32\n", + "P, V = ℂ^2, ℂ^D\n", + "peps₀ = InfinitePEPSTriangular(randn, Float64, P, V; unitcell = (1, 1));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Hamiltonian and CTMRG\n", + "\n", + "The Heisenberg interaction is invariant under a common rotation of both spins, so\n", + "only their relative rotation enters each bond. We absorb it into the gate on the\n", + "second site: `hδ = (id(P) ⊗ R(q⋅δ))' * h * (id(P) ⊗ R(q⋅δ))`. The diagonal\n", + "angle `-4π/3` gives the same gate as `2π/3`. This lets the existing\n", + "`LocalOperator` and `expectation_value` interfaces evaluate the physical\n", + "three-sublattice energy using one tensor and three bond gates.\n", + "\n", + "We use antiferromagnetic coupling `J = 1`. Keep the neighbor coordinates unwrapped,\n", + "including those outside the unit cell: the contractions apply periodicity while\n", + "preserving which bond is intended. Triangular expectation values support onsite and nearest-neighbor\n", + "terms in all three bond families." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "h = S_x_S_x() + S_y_S_y() + S_z_S_z()\n", + "offsets = (CartesianIndex(0, 1), CartesianIndex(1, -1), CartesianIndex(1, 0))\n", + "q = (-2π / 3, 2π / 3)\n", + "rotation(θ) = TensorMap([cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)], P, P)\n", + "gates = map(offsets) do δ\n", + " U = id(P) ⊗ rotation(q[1] * δ[1] + q[2] * δ[2])\n", + " return U' * h * U\n", + "end\n", + "H = LocalOperator(\n", + " physicalspace(peps₀),\n", + " [(I, I + δ) => gate for I in CartesianIndices(size(peps₀)) for (δ, gate) in zip(offsets, gates)],\n", + ");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The environment constructor takes the doubled virtual space `V ⊗ V'` of the PEPS\n", + "norm network and an environment space `ℂ^χ`. We initialize it with `ones`.\n", + "The PEPS tensor and the rotated Heisenberg gates are real, but the current gauge\n", + "fixing for fixed-point differentiation requires a complex environment.\n", + "\n", + "`leading_boundary` selects triangular CTMRG from the environment; here we\n", + "also name it explicitly, so the same settings can be passed to the optimizer.\n", + "We use the default `projector_alg = :full` for the corner construction.\n", + "The alternative `:twothirds` uses two thirds of the environment. `conditioning = true` (the default) conditions the direct\n", + "Ta/Tb edge split. It does not toggle corner conditioning." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "env_initial = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', ℂ^χ; unitcell = size(peps₀))\n", + "boundary_alg = (;\n", + " alg = :SimultaneousCTMRGTriangular, tol = 1.0e-7, maxiter = 100,\n", + " projector_alg = :full, trunctype = :FixedSpaceTruncation,\n", + ")\n", + "env₀, info_ctmrg = leading_boundary(env_initial, peps₀; boundary_alg...);\n", + "@show info_ctmrg.converged info_ctmrg.convergence_error info_ctmrg.contraction_metrics;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convergence compares the corner and edge singular-value spectra. The truncation\n", + "metrics describe the local SVD errors, not the error in the energy. Check CTMRG\n", + "convergence and increase `χ` when assessing the accuracy of an entangled calculation.\n", + "\n", + "`expectation_value` sums all terms in `H`, so this Hamiltonian gives the energy per\n", + "unit cell. Divide by the number of sites to obtain the energy per site:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "E₀ = real(expectation_value(peps₀, H, env₀))\n", + "@show E₀ / length(peps₀);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Native optimization\n", + "\n", + "We use L-BFGS with `FixedPointGradient` and a GMRES linear solver to\n", + "differentiate the CTMRG fixed point. Triangular optimization currently requires\n", + "`precondition_alg = nothing`\n", + "and does not support symmetrization. Alternatively, `gradient_alg = nothing`\n", + "differentiates through the CTMRG iterations.\n", + "\n", + "PEPSKit adapts the CTMRG and gradient tolerances during optimization, using looser\n", + "solves when the gradient is large. Check the final gradient norm and increase the\n", + "optimization budget if necessary.\n", + "Convergence of CTMRG at a fixed `χ` does not establish convergence in `χ`; repeat\n", + "the calculation at larger environment dimensions to assess the energy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "peps, env, E, info_opt = fixedpoint(\n", + " H, peps₀, env₀;\n", + " boundary_alg,\n", + " gradient_alg = (; alg = :FixedPointGradient, solver_alg = (; alg = :GMRES)),\n", + " precondition_alg = nothing,\n", + " optimizer_alg = (; alg = :LBFGS, tol = 1.0e-6, maxiter = 100), verbosity = 3,\n", + ");\n", + "@show E / length(peps) info_opt.gradnorms[end];" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `D = 2` ansatz allows entanglement beyond the 120° product state, whose\n", + "energy is `-3/8` per site. The returned `E` is the energy per unit cell, which here\n", + "contains a single site. Increase both `D` and `χ` to study the quantum ground state." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.11.6", + "language": "julia", + "name": "julia-1.11" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} diff --git a/examples/triangular_ctmrg/main.jl b/examples/triangular_ctmrg/main.jl new file mode 100644 index 000000000..6f416b570 --- /dev/null +++ b/examples/triangular_ctmrg/main.jl @@ -0,0 +1,146 @@ +using Markdown #hide +md""" +# [Triangular CTMRG and optimization](@id examples_triangular_ctmrg) + +This example contracts and optimizes an antiferromagnetic spin-1/2 Heisenberg PEPS +on the triangular lattice using [`InfinitePEPSTriangular`](@ref), +[`CTMRGEnvTriangular`](@ref), and the native [`fixedpoint`](@ref) interface. +A single-site tensor and fixed 120° sublattice rotations represent a three-sublattice +state. We start from a seeded random real tensor with `D = 2` and a complex +environment with `χ = 32`. Check convergence in `χ` before drawing quantitative +conclusions from the energy. + +From the repository root, prepare the example environment once with +`julia --project=examples -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'`, +then run `julia --project=examples --threads=2 examples/triangular_ctmrg/main.jl`. +""" + +using Random +using TensorKit, PEPSKit +using TensorKitTensors.SpinOperators +Random.seed!(123456789); + +md""" +## Unit cell and nearest neighbors + +PEPSKit indexes the unit cell as a Julia matrix: `[r, c]` means row `r`, column `c`, +starting at 1. For `[A B C; D E F; G H I]`, `[2, 1]` is `D` and `[1, 2]` is `B`. +The six virtual legs are ordered `NW, NE, E, SE, SW, W`, with offsets +`(-1, 0), (-1, 1), (0, 1), (1, 0), (1, -1), (0, -1)` respectively. + +To count every nearest-neighbor bond once, we use the following three offsets: + +| Bond family | Offset `(Δr, Δc)` | Virtual leg | Relative rotation | +|:------------|:------------------|:------------|:------------------| +| 0° | `(0, 1)` | E | `2π/3` | +| 60° | `(1, -1)` | SW | `-4π/3` | +| 120° | `(1, 0)` | SE | `-2π/3` | + +The bond-family angles are unoriented; the relative rotations refer to the directed +offsets in the table. With `R(θ) = exp(-im * θ * S_y)`, the physical tensor at +`[r, c]` is `R((r - 1)q[1] + (c - 1)q[2]) * A`, where `q = (-2π/3, 2π/3)`. +This gives `[A B C; C A B; B C A]`, with `B = R(2π/3) * A` and +`C = R(4π/3) * A`, up to irrelevant spin-1/2 tensor phases. It is the PEPSKit +version of the variPEPS pattern `[A B C; B C A; C A B]` with `q = (2π/3, 2π/3)`. + +We optimize `A` in a single-site unit cell. All three sublattices share its +variational parameters, and the wavevector stays fixed during optimization. +""" + +D = 2 +χ = 32 +P, V = ℂ^2, ℂ^D +peps₀ = InfinitePEPSTriangular(randn, Float64, P, V; unitcell = (1, 1)); + +md""" +## Hamiltonian and CTMRG + +The Heisenberg interaction is invariant under a common rotation of both spins, so +only their relative rotation enters each bond. We absorb it into the gate on the +second site: `hδ = (id(P) ⊗ R(q⋅δ))' * h * (id(P) ⊗ R(q⋅δ))`. The diagonal +angle `-4π/3` gives the same gate as `2π/3`. This lets the existing +[`LocalOperator`](@ref) and `expectation_value` interfaces evaluate the physical +three-sublattice energy using one tensor and three bond gates. + +We use antiferromagnetic coupling `J = 1`. Keep the neighbor coordinates unwrapped, +including those outside the unit cell: the contractions apply periodicity while +preserving which bond is intended. Triangular expectation values support onsite and nearest-neighbor +terms in all three bond families. +""" + +h = S_x_S_x() + S_y_S_y() + S_z_S_z() +offsets = (CartesianIndex(0, 1), CartesianIndex(1, -1), CartesianIndex(1, 0)) +q = (-2π / 3, 2π / 3) +rotation(θ) = TensorMap([cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)], P, P) +gates = map(offsets) do δ + U = id(P) ⊗ rotation(q[1] * δ[1] + q[2] * δ[2]) + return U' * h * U +end +H = LocalOperator( + physicalspace(peps₀), + [(I, I + δ) => gate for I in CartesianIndices(size(peps₀)) for (δ, gate) in zip(offsets, gates)], +); + +md""" +The environment constructor takes the doubled virtual space `V ⊗ V'` of the PEPS +norm network and an environment space `ℂ^χ`. We initialize it with `ones`. +The PEPS tensor and the rotated Heisenberg gates are real, but the current gauge +fixing for fixed-point differentiation requires a complex environment. + +[`leading_boundary`](@ref) selects triangular CTMRG from the environment; here we +also name it explicitly, so the same settings can be passed to the optimizer. +We use the default `projector_alg = :full` for the corner construction. +The alternative `:twothirds` uses two thirds of the environment. `conditioning = true` (the default) conditions the direct +Ta/Tb edge split. It does not toggle corner conditioning. +""" + +env_initial = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', ℂ^χ; unitcell = size(peps₀)) +boundary_alg = (; + alg = :SimultaneousCTMRGTriangular, tol = 1.0e-7, maxiter = 100, + projector_alg = :full, trunctype = :FixedSpaceTruncation, +) +env₀, info_ctmrg = leading_boundary(env_initial, peps₀; boundary_alg...); +@show info_ctmrg.converged info_ctmrg.convergence_error info_ctmrg.contraction_metrics; + +md""" +Convergence compares the corner and edge singular-value spectra. The truncation +metrics describe the local SVD errors, not the error in the energy. Check CTMRG +convergence and increase `χ` when assessing the accuracy of an entangled calculation. + +`expectation_value` sums all terms in `H`, so this Hamiltonian gives the energy per +unit cell. Divide by the number of sites to obtain the energy per site: +""" + +E₀ = real(expectation_value(peps₀, H, env₀)) +@show E₀ / length(peps₀); + +md""" +## Native optimization + +We use L-BFGS with [`FixedPointGradient`](@ref) and a GMRES linear solver to +differentiate the CTMRG fixed point. Triangular optimization currently requires +`precondition_alg = nothing` +and does not support symmetrization. Alternatively, `gradient_alg = nothing` +differentiates through the CTMRG iterations. + +PEPSKit adapts the CTMRG and gradient tolerances during optimization, using looser +solves when the gradient is large. Check the final gradient norm and increase the +optimization budget if necessary. +Convergence of CTMRG at a fixed `χ` does not establish convergence in `χ`; repeat +the calculation at larger environment dimensions to assess the energy. +""" + +peps, env, E, info_opt = fixedpoint( + H, peps₀, env₀; + boundary_alg, + gradient_alg = (; alg = :FixedPointGradient, solver_alg = (; alg = :GMRES)), + precondition_alg = nothing, + optimizer_alg = (; alg = :LBFGS, tol = 1.0e-6, maxiter = 100), verbosity = 3, +); +@show E / length(peps) info_opt.gradnorms[end]; + +md""" +The `D = 2` ansatz allows entanglement beyond the 120° product state, whose +energy is `-3/8` per site. The returned `E` is the energy per unit cell, which here +contains a single site. Increase both `D` and `χ` to study the quantum ground state. +""" diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b715dc2ec..466a66add 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -75,8 +75,12 @@ include("utility/initialization.jl") include("networks/tensors.jl") include("networks/local_sandwich.jl") include("networks/infinitesquarenetwork.jl") +include("networks/tensors_triangular.jl") +include("networks/local_triangular_sandwich.jl") +include("networks/infinitetriangularnetwork.jl") include("states/infinitepeps.jl") +include("states/infinitepepstriangular.jl") include("states/infinitepartitionfunction.jl") include("utility/symmetrization.jl") @@ -92,6 +96,7 @@ include("operators/lattices/squarelattice.jl") include("operators/models.jl") include("environments/ctmrg_environments.jl") +include("environments/ctmrg_environments_triangular.jl") include("environments/c4v_ctmrg_environments.jl") include("environments/vumps_environments.jl") include("environments/suweight.jl") @@ -148,6 +153,17 @@ include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_tensor.jl") include("algorithms/truncation/bond_truncation.jl") +include("algorithms/contractions/triangular_ctmrg/double_corner.jl") +include("algorithms/contractions/triangular_ctmrg/renormalize_edge.jl") +include("algorithms/contractions/triangular_ctmrg/halfinf_env.jl") +include("algorithms/contractions/triangular_ctmrg/renormalize_corner.jl") +include("algorithms/contractions/triangular_ctmrg/network_value.jl") +include("algorithms/contractions/triangular_ctmrg/expectation_value.jl") +include("algorithms/ctmrg_triangular/ctmrg.jl") +include("algorithms/ctmrg_triangular/utility.jl") +include("algorithms/ctmrg_triangular/simultaneous.jl") +include("algorithms/ctmrg_triangular/gaugefix.jl") + include("algorithms/compress/local.jl") include("algorithms/time_evolution/apply_gate.jl") @@ -170,6 +186,7 @@ include("algorithms/expectation_value/network_value.jl") include("algorithms/expectation_value/correlator_adapters.jl") include("algorithms/expectation_value/correlators.jl") include("algorithms/toolbox.jl") +include("algorithms/toolbox_triangular.jl") include("algorithms/optimization/implicit_differentiation.jl") include("algorithms/optimization/preconditioning.jl") @@ -180,6 +197,7 @@ include("algorithms/select_algorithm.jl") using .Defaults: set_scheduler! export set_scheduler! export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, OrthAdjoint, QRAdjoint +export CTMRGEnvTriangular, SimultaneousCTMRGTriangular export C4vCTMRGEnv, CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG export initialize_ctmrg_environment, RandomInitialization, ProductStateInitialization, ApplicationInitialization, IdentityInitialization @@ -204,9 +222,10 @@ export ALSTruncation, FullEnvTruncation export SimpleUpdate export TimeEvolver, timestep, time_evolve -export InfiniteSquareNetwork +export InfiniteSquareNetwork, InfiniteTriangularNetwork export InfinitePartitionFunction export InfinitePEPS, InfiniteTransferPEPS +export InfinitePEPSTriangular export SUWeight export InfinitePEPO, InfiniteTransferPEPO diff --git a/src/algorithms/contractions/triangular_ctmrg/double_corner.jl b/src/algorithms/contractions/triangular_ctmrg/double_corner.jl new file mode 100644 index 000000000..649d4b5ab --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/double_corner.jl @@ -0,0 +1,23 @@ +# Directions refer to the local frame after rotl60; t/b denote ket/bra layers. + +function build_double_corner_matrix_triangular(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir::Int, r::Int, c::Int) where {P <: PFTensorTriangular} + @tensor opt = true mat[χW D240; χE D300] := + env.C[_coordinates(dir, 0, r, c, size(network))...][χNW D120; χN] * + env.C[_coordinates(dir + 1, 0, r, c, size(network))...][χN D60; χNE] * + env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][χW D180; χNW] * + env.Eb[_coordinates(dir + 1, 1, r, c, size(network))...][χNE D0; χE] * + rotl60(network[r, c], mod(dir - 1, 6))[D180 D240 D300; D120 D60 D0] + return mat +end + +function build_double_corner_matrix_triangular(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir::Int, r::Int, c::Int) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + @tensor opt = true mat[χW Dt240 Db240; χE Dt300 Db300] := + env.C[_coordinates(dir, 0, r, c, size(network))...][χNW Dt120 Db120; χN] * + env.C[_coordinates(dir + 1, 0, r, c, size(network))...][χN Dt60 Db60; χNE] * + env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][χW Dt180 Db180; χNW] * + env.Eb[_coordinates(dir + 1, 1, r, c, size(network))...][χNE Dt0 Db0; χE] * + rotl60(ket(O_r_c), mod(dir - 1, 6))[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * + conj(rotl60(bra(O_r_c), mod(dir - 1, 6))[d; Db120 Db60 Db0 Db300 Db240 Db180]) + return mat +end diff --git a/src/algorithms/contractions/triangular_ctmrg/expectation_value.jl b/src/algorithms/contractions/triangular_ctmrg/expectation_value.jl new file mode 100644 index 000000000..7a5ff04a2 --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/expectation_value.jl @@ -0,0 +1,42 @@ +function _contract_corners((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, op::AbstractTensorMap{E, S, 1, 1}) where {P <: PEPSSandwichTriangular, E, S} + O_r_c = network[r, c] + return @tensor opt = true env.C[1, _prev(r, end), c][χNW Dt120 Db120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN Dt60 Db60; χNE] * env.C[3, r, _next(c, end)][χNE Dt0 Db0; χSE] * + env.C[4, _next(r, end), c][χSE Dt300 Db300; χS] * env.C[5, _next(r, end), _prev(c, end)][χS Dt240 Db240; χSW] * env.C[6, r, _prev(c, end)][χSW Dt180 Db180; χNW] * + ket(O_r_c)[dt; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * conj(bra(O_r_c)[db; Db120 Db60 Db0 Db300 Db240 Db180]) * op[db; dt] +end + +function _contract_edges_0((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, op::AbstractTensorMap{E, S, 2, 2}) where {P <: PEPSSandwichTriangular, E, S} + O_r_c = network[r, c] + O_r_cn = network[r, _next(c, end)] + return @tensor opt = true ket(O_r_c)[dtL; DLt120 DLt60 DLt0 DLt300 DLt240 DLt180] * ket(O_r_cn)[dtR; DRt120 DRt60 DRt0 DRt300 DRt240 DLt0] * + conj(bra(O_r_c)[dbL; DLb120 DLb60 DLb0 DLb300 DLb240 DLb180]) * conj(bra(O_r_cn)[dbR; DRb120 DRb60 DRb0 DRb300 DRb240 DLb0]) * + env.C[1, _prev(r, end), c][χNW DLt120 DLb120; χaN] * env.C[2, _prev(r, end), _next(c + 1, end)][χbN DRt60 DRb60; χNE] * env.C[3, r, _next(c + 1, end)][χNE DRt0 DRb0; χSE] * + env.C[4, _next(r, end), _next(c, end)][χSE DRt300 DRb300; χaS] * env.C[5, _next(r, end), _prev(c, end)][χbS DLt240 DLb240; χSW] * env.C[6, r, _prev(c, end)][χSW DLt180 DLb180; χNW] * + env.Eb[1, _prev(r, end), _next(c, end)][χaN DLt60 DLb60; χCN] * env.Ea[1, _prev(r, end), _next(c, end)][χCN DRt120 DRb120; χbN] * + env.Eb[4, _next(r, end), c][χaS DRt240 DRb240; χCS] * env.Ea[4, _next(r, end), c][χCS DLt300 DLb300; χbS] * + op[dbL dbR; dtL dtR] +end + +function _contract_edges_60((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, op::AbstractTensorMap{E, S, 2, 2}) where {P <: PEPSSandwichTriangular, E, S} + O_r_c = network[r, c] + O_rn_cp = network[_next(r, end), _prev(c, end)] + return @tensor opt = true ket(O_r_c)[dtTR; DTRt120 DTRt60 DTRt0 DTRt300 DBLt60 DTRt180] * ket(O_rn_cp)[dtBL; DBLt120 DBLt60 DBLt0 DBLt300 DBLt240 DBLt180] * + conj(bra(O_r_c)[dbTR; DTRb120 DTRb60 DTRb0 DTRb300 DBLb60 DTRb180]) * conj(bra(O_rn_cp)[dbBL; DBLb120 DBLb60 DBLb0 DBLb300 DBLb240 DBLb180]) * + env.C[1, _prev(r, end), c][χbNW DTRt120 DTRb120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTRt60 DTRb60; χNE] * env.C[3, r, _next(c, end)][χNE DTRt0 DTRb0; χaSE] * + env.C[4, _next(r + 1, end), _prev(c, end)][χbSE DBLt300 DBLb300; χS] * env.C[5, _next(r + 1, end), _prev(c - 1, end)][χS DBLt240 DBLb240; χSW] * env.C[6, _next(r, end), _prev(c - 1, end)][χSW DBLt180 DBLb180; χaNW] * + env.Eb[3, _next(r, end), c][χaSE DTRt300 DTRb300; χCSE] * env.Ea[3, _next(r, end), c][χCSE DBLt0 DBLb0; χbSE] * + env.Eb[6, r, _prev(c, end)][χaNW DBLt120 DBLb120; χCNW] * env.Ea[6, r, _prev(c, end)][χCNW DTRt180 DTRb180; χbNW] * + op[dbTR dbBL; dtTR dtBL] +end + +function _contract_edges_120((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, op::AbstractTensorMap{E, S, 2, 2}) where {P <: PEPSSandwichTriangular, E, S} + O_r_c = network[r, c] + O_rn_c = network[_next(r, end), c] + return @tensor opt = true ket(O_r_c)[dtTL; DTLt120 DTLt60 DTLt0 DTLt300 DTLt240 DTLt180] * ket(O_rn_c)[dtBR; DTLt300 DBRt60 DBRt0 DBRt300 DBRt240 DBRt180] * + conj(bra(O_r_c)[dbTL; DTLb120 DTLb60 DTLb0 DTLb300 DTLb240 DTLb180]) * conj(bra(O_rn_c)[dbBR; DTLb300 DBRb60 DBRb0 DBRb300 DBRb240 DBRb180]) * + env.C[1, _prev(r, end), c][χNW DTLt120 DTLb120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTLt60 DTLb60; χaNE] * env.C[3, _next(r, end), _next(c, end)][χbNE DBRt0 DBRb0; χSE] * + env.C[4, _next(r + 1, end), c][χSE DBRt300 DBRb300; χS] * env.C[5, _next(r + 1, end), _prev(c, end)][χS DBRt240 DBRb240; χaSW] * env.C[6, r, _prev(c, end)][χbSW DTLt180 DTLb180; χNW] * + env.Eb[2, r, _next(c, end)][χaNE DTLt0 DTLb0; χCNE] * env.Ea[2, r, _next(c, end)][χCNE DBRt60 DBRb60; χbNE] * + env.Eb[5, _next(r, end), _prev(c, end)][χaSW DBRt180 DBRb180; χCSW] * env.Ea[5, _next(r, end), _prev(c, end)][χCSW DTLt240 DTLb240; χbSW] * + op[dbTL dbBR; dtTL dtBR] +end diff --git a/src/algorithms/contractions/triangular_ctmrg/halfinf_env.jl b/src/algorithms/contractions/triangular_ctmrg/halfinf_env.jl new file mode 100644 index 000000000..f93e2baaa --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/halfinf_env.jl @@ -0,0 +1,41 @@ +# Directions refer to the local frame after rotl60; t/b denote ket/bra layers. + +function build_halfinfinite_environments(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir, r, c) where {P <: PEPSSandwichTriangular} + rL, cL = _coordinates(dir + 4, 0, r, c, size(network))[2:3] + rR, cR = _coordinates(dir + 3, 0, r, c, size(network))[2:3] + O_L = network[rL, cL] + O_R = network[rR, cR] + @tensor opt = true ρL[χSE Dt0 Db0; χN Dt60 Db60] := + env.C[_coordinates(dir, 0, rL, cL, size(network))...][χNW Dt120 Db120; χN] * + env.C[_coordinates(dir - 1, 0, rL, cL, size(network))...][χSW Dt180 Db180; χNW] * + env.C[_coordinates(dir - 2, 0, rL, cL, size(network))...][χS Dt240 Db240; χSW] * + env.Ea[_coordinates(dir + 3, 0, rL, cL, size(network))...][χSE Dt300 Db300; χS] * + rotl60(ket(O_L), mod(dir - 1, 6))[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * + conj(rotl60(bra(O_L), mod(dir - 1, 6))[d; Db120 Db60 Db0 Db300 Db240 Db180]) + @tensor opt = true ρR[χN Dt120 Db120; χSW Dt180 Db180] := + env.C[_coordinates(dir + 1, 0, rR, cR, size(network))...][χN Dt60 Db60; χNE] * + env.C[_coordinates(dir + 2, 0, rR, cR, size(network))...][χNE Dt0 Db0; χSE] * + env.C[_coordinates(dir + 3, 0, rR, cR, size(network))...][χSE Dt300 Db300; χS] * + env.Eb[_coordinates(dir + 3, 1, rR, cR, size(network))...][χS Dt240 Db240; χSW] * + rotl60(ket(O_R), mod(dir - 1, 6))[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * + conj(rotl60(bra(O_R), mod(dir - 1, 6))[d; Db120 Db60 Db0 Db300 Db240 Db180]) + return ρL, ρR +end + +function build_halfinfinite_environments(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir, r, c) where {P <: PFTensorTriangular} + rL, cL = _coordinates(dir + 4, 0, r, c, size(network))[2:3] + rR, cR = _coordinates(dir + 3, 0, r, c, size(network))[2:3] + @tensor opt = true ρL[χSE D0; χN D60] := + env.C[_coordinates(dir, 0, rL, cL, size(network))...][χNW D120; χN] * + env.C[_coordinates(dir - 1, 0, rL, cL, size(network))...][χSW D180; χNW] * + env.C[_coordinates(dir - 2, 0, rL, cL, size(network))...][χS D240; χSW] * + env.Ea[_coordinates(dir + 3, 0, rL, cL, size(network))...][χSE D300; χS] * + rotl60(network[rL, cL], mod(dir - 1, 6))[D180 D240 D300; D120 D60 D0] + @tensor opt = true ρR[χN D120; χSW D180] := + env.C[_coordinates(dir + 1, 0, rR, cR, size(network))...][χN D60; χNE] * + env.C[_coordinates(dir + 2, 0, rR, cR, size(network))...][χNE D0; χSE] * + env.C[_coordinates(dir + 3, 0, rR, cR, size(network))...][χSE D300; χS] * + env.Eb[_coordinates(dir + 3, 1, rR, cR, size(network))...][χS D240; χSW] * + rotl60(network[rR, cR], mod(dir - 1, 6))[D180 D240 D300; D120 D60 D0] + return ρL, ρR +end diff --git a/src/algorithms/contractions/triangular_ctmrg/network_value.jl b/src/algorithms/contractions/triangular_ctmrg/network_value.jl new file mode 100644 index 000000000..b0c537145 --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/network_value.jl @@ -0,0 +1,103 @@ +function _contract_edges_0((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PFTensorTriangular} + return @tensor opt = true network[r, c][DL180 DL240 DL300; DL120 DL60 DL0] * network[r, _next(c, end)][DL0 DR240 DR300; DR120 DR60 DR0] * + env.C[1, _prev(r, end), c][χNW DL120; χaN] * env.C[2, _prev(r, end), _next(c + 1, end)][χbN DR60; χNE] * env.C[3, r, _next(c + 1, end)][χNE DR0; χSE] * + env.C[4, _next(r, end), _next(c, end)][χSE DR300; χaS] * env.C[5, _next(r, end), _prev(c, end)][χbS DL240; χSW] * env.C[6, r, _prev(c, end)][χSW DL180; χNW] * + env.Eb[1, _prev(r, end), _next(c, end)][χaN DL60; χCN] * env.Ea[1, _prev(r, end), _next(c, end)][χCN DR120; χbN] * + env.Eb[4, _next(r, end), c][χaS DR240; χCS] * env.Ea[4, _next(r, end), c][χCS DL300; χbS] +end + +function _contract_edges_60((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PFTensorTriangular} + return @tensor opt = true network[r, c][DTR180 DBL60 DTR300; DTR120 DTR60 DTR0] * network[_next(r, end), _prev(c, end)][DBL180 DBL240 DBL300; DBL120 DBL60 DBL0] * + env.C[1, _prev(r, end), c][χbNW DTR120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTR60; χNE] * env.C[3, r, _next(c, end)][χNE DTR0; χaSE] * + env.C[4, _next(r + 1, end), _prev(c, end)][χbSE DBL300; χS] * env.C[5, _next(r + 1, end), _prev(c - 1, end)][χS DBL240; χSW] * env.C[6, _next(r, end), _prev(c - 1, end)][χSW DBL180; χaNW] * + env.Eb[3, _next(r, end), c][χaSE DTR300; χCSE] * env.Ea[3, _next(r, end), c][χCSE DBL0; χbSE] * + env.Eb[6, r, _prev(c, end)][χaNW DBL120; χCNW] * env.Ea[6, r, _prev(c, end)][χCNW DTR180; χbNW] +end + +function _contract_edges_120((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PFTensorTriangular} + return @tensor opt = true network[r, c][DTL180 DTL240 DTL300; DTL120 DTL60 DTL0] * network[_next(r, end), c][DBR180 DBR240 DBR300; DTL300 DBR60 DBR0] * + env.C[1, _prev(r, end), c][χNW DTL120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTL60; χaNE] * env.C[3, _next(r, end), _next(c, end)][χbNE DBR0; χSE] * + env.C[4, _next(r + 1, end), c][χSE DBR300; χS] * env.C[5, _next(r + 1, end), _prev(c, end)][χS DBR240; χaSW] * env.C[6, r, _prev(c, end)][χbSW DTL180; χNW] * + env.Eb[2, r, _next(c, end)][χaNE DTL0; χCNE] * env.Ea[2, r, _next(c, end)][χCNE DBR60; χbNE] * + env.Eb[5, _next(r, end), _prev(c, end)][χaSW DBR180; χCSW] * env.Ea[5, _next(r, end), _prev(c, end)][χCSW DTL240; χbSW] +end + +function _contract_site_large((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PFTensorTriangular} + return @tensor opt = true network[_prev(r, end), c][DNW180 DW60 DNW300; DNW120 DNW60 DNW0] * network[_prev(r, end), _next(c, end)][DNW0 DNE240 DNE300; DNE120 DNE60 DNE0] * + network[r, _next(c, end)][DE180 DE240 DE300; DNE300 DE60 DE0] * network[_next(r, end), c][DSE180 DSE240 DSE300; DSE120 DE240 DSE0] * network[_next(r, end), _prev(c, end)][DSW180 DSW240 DSW300; DSW120 DSW60 DSE180] * + network[r, _prev(c, end)][DW180 DW240 DSW120; DW120 DW60 DW0] * network[r, c][DW0 DSW60 DSE120; DNW300 DNE240 DE180] * + env.C[1, _prev(r - 1, end), c][χaNW DNW120; χbN] * env.Eb[1, _prev(r - 1, end), _next(c, end)][χbN DNW60; χCN] * env.Ea[1, _prev(r - 1, end), _next(c, end)][χCN DNE120; χaN] * + env.C[2, _prev(r - 1, end), _next(c + 1, end)][χaN DNE60; χbNE] * env.Eb[2, _prev(r, end), _next(c + 1, end)][χbNE DNE0; χCNE] * env.Ea[2, _prev(r, end), _next(c + 1, end)][χCNE DE60; χaNE] * + env.C[3, r, _next(c + 1, end)][χaNE DE0; χbSE] * env.Eb[3, _next(r, end), _next(c, end)][χbSE DE300; χCSE] * env.Ea[3, _next(r, end), _next(c, end)][χCSE DSE0; χaSE] * + env.C[4, _next(r + 1, end), c][χaSE DSE300; χbS] * env.Eb[4, _next(r + 1, end), _prev(c, end)][χbS DSE240; χCS] * env.Ea[4, _next(r + 1, end), _prev(c, end)][χCS DSW300; χaS] * + env.C[5, _next(r + 1, end), _prev(c - 1, end)][χaS DSW240; χbSW] * env.Eb[5, _next(r, end), _prev(c - 1, end)][χbSW DSW180; χCSW] * env.Ea[5, _next(r, end), _prev(c - 1, end)][χCSW DW240; χaSW] * + env.C[6, r, _prev(c - 1, end)][χaSW DW180; χbNW] * env.Eb[6, _prev(r, end), _prev(c, end)][χbNW DW120; χCNW] * env.Ea[6, _prev(r, end), _prev(c, end)][χCNW DNW180; χaNW] +end + +function _contract_corners((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PFTensorTriangular} + return @tensor opt = true env.C[1, _prev(r, end), c][χNW D120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN D60; χNE] * env.C[3, r, _next(c, end)][χNE D0; χSE] * + env.C[4, _next(r, end), c][χSE D300; χS] * env.C[5, _next(r, end), _prev(c, end)][χS D240; χSW] * env.C[6, r, _prev(c, end)][χSW D180; χNW] * + network[r, c][D180 D240 D300; D120 D60 D0] +end + +function _contract_edges_0((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + O_r_cn = network[r, _next(c, end)] + return @tensor opt = true ket(O_r_c)[dL; DLt120 DLt60 DLt0 DLt300 DLt240 DLt180] * ket(O_r_cn)[dR; DRt120 DRt60 DRt0 DRt300 DRt240 DLt0] * + conj(bra(O_r_c)[dL; DLb120 DLb60 DLb0 DLb300 DLb240 DLb180]) * conj(bra(O_r_cn)[dR; DRb120 DRb60 DRb0 DRb300 DRb240 DLb0]) * + env.C[1, _prev(r, end), c][χNW DLt120 DLb120; χaN] * env.C[2, _prev(r, end), _next(c + 1, end)][χbN DRt60 DRb60; χNE] * env.C[3, r, _next(c + 1, end)][χNE DRt0 DRb0; χSE] * + env.C[4, _next(r, end), _next(c, end)][χSE DRt300 DRb300; χaS] * env.C[5, _next(r, end), _prev(c, end)][χbS DLt240 DLb240; χSW] * env.C[6, r, _prev(c, end)][χSW DLt180 DLb180; χNW] * + env.Eb[1, _prev(r, end), _next(c, end)][χaN DLt60 DLb60; χCN] * env.Ea[1, _prev(r, end), _next(c, end)][χCN DRt120 DRb120; χbN] * + env.Eb[4, _next(r, end), c][χaS DRt240 DRb240; χCS] * env.Ea[4, _next(r, end), c][χCS DLt300 DLb300; χbS] +end + +function _contract_edges_60((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + O_rn_cp = network[_next(r, end), _prev(c, end)] + return @tensor opt = true ket(O_r_c)[dTR; DTRt120 DTRt60 DTRt0 DTRt300 DBLt60 DTRt180] * ket(O_rn_cp)[dBL; DBLt120 DBLt60 DBLt0 DBLt300 DBLt240 DBLt180] * + conj(bra(O_r_c)[dTR; DTRb120 DTRb60 DTRb0 DTRb300 DBLb60 DTRb180]) * conj(bra(O_rn_cp)[dBL; DBLb120 DBLb60 DBLb0 DBLb300 DBLb240 DBLb180]) * + env.C[1, _prev(r, end), c][χbNW DTRt120 DTRb120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTRt60 DTRb60; χNE] * env.C[3, r, _next(c, end)][χNE DTRt0 DTRb0; χaSE] * + env.C[4, _next(r + 1, end), _prev(c, end)][χbSE DBLt300 DBLb300; χS] * env.C[5, _next(r + 1, end), _prev(c - 1, end)][χS DBLt240 DBLb240; χSW] * env.C[6, _next(r, end), _prev(c - 1, end)][χSW DBLt180 DBLb180; χaNW] * + env.Eb[3, _next(r, end), c][χaSE DTRt300 DTRb300; χCSE] * env.Ea[3, _next(r, end), c][χCSE DBLt0 DBLb0; χbSE] * + env.Eb[6, r, _prev(c, end)][χaNW DBLt120 DBLb120; χCNW] * env.Ea[6, r, _prev(c, end)][χCNW DTRt180 DTRb180; χbNW] +end + +function _contract_edges_120((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + O_rn_c = network[_next(r, end), c] + return @tensor opt = true ket(O_r_c)[dTL; DTLt120 DTLt60 DTLt0 DTLt300 DTLt240 DTLt180] * ket(O_rn_c)[dBR; DTLt300 DBRt60 DBRt0 DBRt300 DBRt240 DBRt180] * + conj(bra(O_r_c)[dTL; DTLb120 DTLb60 DTLb0 DTLb300 DTLb240 DTLb180]) * conj(bra(O_rn_c)[dBR; DTLb300 DBRb60 DBRb0 DBRb300 DBRb240 DBRb180]) * + env.C[1, _prev(r, end), c][χNW DTLt120 DTLb120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN DTLt60 DTLb60; χaNE] * env.C[3, _next(r, end), _next(c, end)][χbNE DBRt0 DBRb0; χSE] * + env.C[4, _next(r + 1, end), c][χSE DBRt300 DBRb300; χS] * env.C[5, _next(r + 1, end), _prev(c, end)][χS DBRt240 DBRb240; χaSW] * env.C[6, r, _prev(c, end)][χbSW DTLt180 DTLb180; χNW] * + env.Eb[2, r, _next(c, end)][χaNE DTLt0 DTLb0; χCNE] * env.Ea[2, r, _next(c, end)][χCNE DBRt60 DBRb60; χbNE] * + env.Eb[5, _next(r, end), _prev(c, end)][χaSW DBRt180 DBRb180; χCSW] * env.Ea[5, _next(r, end), _prev(c, end)][χCSW DTLt240 DTLb240; χbSW] +end + +function _contract_site_large((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + O_rp_c = network[_prev(r, end), c] + O_rp_cn = network[_prev(r, end), _next(c, end)] + O_rn_cp = network[_next(r, end), _prev(c, end)] + O_rn_c = network[_next(r, end), c] + O_r_cn = network[r, _next(c, end)] + O_r_cp = network[r, _prev(c, end)] + return @tensor opt = true ket(O_rp_c)[dNW; DNWt120 DNWt60 DNWt0 DNWt300 DWt60 DNWt180] * ket(O_rp_cn)[dNE; DNEt120 DNEt60 DNEt0 DNEt300 DNEt240 DNWt0] * + ket(O_r_cn)[dE; DNEt300 DEt60 DEt0 DEt300 DEt240 DEt180] * ket(O_rn_c)[dSE; DSEt120 DEt240 DSEt0 DSEt300 DSEt240 DSEt180] * ket(O_rn_cp)[dSW; DSWt120 DSWt60 DSEt180 DSWt300 DSWt240 DSWt180] * + ket(O_r_cp)[dW; DWt120 DWt60 DWt0 DSWt120 DWt240 DWt180] * ket(O_r_c)[dCenter; DNWt300 DNEt240 DEt180 DSEt120 DSWt60 DWt0] * + conj(bra(O_rp_c)[dNW; DNWb120 DNWb60 DNWb0 DNWb300 DWb60 DNWb180]) * conj(bra(O_rp_cn)[dNE; DNEb120 DNEb60 DNEb0 DNEb300 DNEb240 DNWb0]) * + conj(bra(O_r_cn)[dE; DNEb300 DEb60 DEb0 DEb300 DEb240 DEb180]) * conj(bra(O_rn_c)[dSE; DSEb120 DEb240 DSEb0 DSEb300 DSEb240 DSEb180]) * conj(bra(O_rn_cp)[dSW; DSWb120 DSWb60 DSEb180 DSWb300 DSWb240 DSWb180]) * + conj(bra(O_r_cp)[dW; DWb120 DWb60 DWb0 DSWb120 DWb240 DWb180]) * conj(bra(O_r_c)[dCenter; DNWb300 DNEb240 DEb180 DSEb120 DSWb60 DWb0]) * + env.C[1, _prev(r - 1, end), c][χaNW DNWt120 DNWb120; χbN] * env.Eb[1, _prev(r - 1, end), _next(c, end)][χbN DNWt60 DNWb60; χCN] * env.Ea[1, _prev(r - 1, end), _next(c, end)][χCN DNEt120 DNEb120; χaN] * + env.C[2, _prev(r - 1, end), _next(c + 1, end)][χaN DNEt60 DNEb60; χbNE] * env.Eb[2, _prev(r, end), _next(c + 1, end)][χbNE DNEt0 DNEb0; χCNE] * env.Ea[2, _prev(r, end), _next(c + 1, end)][χCNE DEt60 DEb60; χaNE] * + env.C[3, r, _next(c + 1, end)][χaNE DEt0 DEb0; χbSE] * env.Eb[3, _next(r, end), _next(c, end)][χbSE DEt300 DEb300; χCSE] * env.Ea[3, _next(r, end), _next(c, end)][χCSE DSEt0 DSEb0; χaSE] * + env.C[4, _next(r + 1, end), c][χaSE DSEt300 DSEb300; χbS] * env.Eb[4, _next(r + 1, end), _prev(c, end)][χbS DSEt240 DSEb240; χCS] * env.Ea[4, _next(r + 1, end), _prev(c, end)][χCS DSWt300 DSWb300; χaS] * + env.C[5, _next(r + 1, end), _prev(c - 1, end)][χaS DSWt240 DSWb240; χbSW] * env.Eb[5, _next(r, end), _prev(c - 1, end)][χbSW DSWt180 DSWb180; χCSW] * env.Ea[5, _next(r, end), _prev(c - 1, end)][χCSW DWt240 DWb240; χaSW] * + env.C[6, r, _prev(c - 1, end)][χaSW DWt180 DWb180; χbNW] * env.Eb[6, _prev(r, end), _prev(c, end)][χbNW DWt120 DWb120; χCNW] * env.Ea[6, _prev(r, end), _prev(c, end)][χCNW DNWt180 DNWb180; χaNW] +end + +function _contract_corners((r, c)::Tuple{Int, Int}, network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + return @tensor opt = true env.C[1, _prev(r, end), c][χNW Dt120 Db120; χN] * env.C[2, _prev(r, end), _next(c, end)][χN Dt60 Db60; χNE] * env.C[3, r, _next(c, end)][χNE Dt0 Db0; χSE] * + env.C[4, _next(r, end), c][χSE Dt300 Db300; χS] * env.C[5, _next(r, end), _prev(c, end)][χS Dt240 Db240; χSW] * env.C[6, r, _prev(c, end)][χSW Dt180 Db180; χNW] * + ket(O_r_c)[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * conj(bra(O_r_c)[d; Db120 Db60 Db0 Db300 Db240 Db180]) +end diff --git a/src/algorithms/contractions/triangular_ctmrg/renormalize_corner.jl b/src/algorithms/contractions/triangular_ctmrg/renormalize_corner.jl new file mode 100644 index 000000000..e662e4091 --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/renormalize_corner.jl @@ -0,0 +1,23 @@ +# Directions refer to the local frame after rotl60; t/b denote ket/bra layers. + +function renormalize_corner_triangular((dir, r, c), network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PFTensorTriangular} + @tensor opt = true env_C_new[χSW D300; χSE] := + env.C[_coordinates(dir, 0, r, c, size(network))...][χNW D120; χN] * + env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][χW D180; χNW] * + env.Eb[_coordinates(dir, 1, r, c, size(network))...][χN D60; χNE] * + rotl60(network[r, c], mod(dir - 1, 6))[D180 D240 D300; D120 D60 D0] * + Pas[_coordinates(dir - 1, 5, r, c, size(network))...][χSW; χW D240] * Pbs[dir, r, c][χNE D0; χSE] + return env_C_new +end + +function renormalize_corner_triangular((dir, r, c), network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + @tensor opt = true env_C_new[χSW Dt300 Db300; χSE] := + env.C[_coordinates(dir, 0, r, c, size(network))...][χNW Dt120 Db120; χN] * + env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][χW Dt180 Db180; χNW] * + env.Eb[_coordinates(dir, 1, r, c, size(network))...][χN Dt60 Db60; χNE] * + rotl60(ket(O_r_c), mod(dir - 1, 6))[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * + conj(rotl60(bra(O_r_c), mod(dir - 1, 6))[d; Db120 Db60 Db0 Db300 Db240 Db180]) * + Pas[_coordinates(dir - 1, 5, r, c, size(network))...][χSW; χW Dt240 Db240] * Pbs[dir, r, c][χNE Dt0 Db0; χSE] + return env_C_new +end diff --git a/src/algorithms/contractions/triangular_ctmrg/renormalize_edge.jl b/src/algorithms/contractions/triangular_ctmrg/renormalize_edge.jl new file mode 100644 index 000000000..6b579f0e2 --- /dev/null +++ b/src/algorithms/contractions/triangular_ctmrg/renormalize_edge.jl @@ -0,0 +1,22 @@ +# Directions refer to the local frame after rotl60; t/b denote ket/bra layers. + +function semi_renormalize_edge(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Pas, Pbs, dir, r, c) where {P <: PEPSSandwichTriangular} + O_r_c = network[r, c] + @tensor opt = true mat[χW Dt240 Db240; χE Dt300 Db300] := + Pas[_coordinates(dir, 5, r, c, size(network))...][χW; χNW Dt180 Db180] * Pbs[dir, r, c][χNE Dt0 Db0; χE] * + env.Ea[_coordinates(dir, 0, r, c, size(network))...][χNW Dt120 Db120; χN] * + env.Eb[_coordinates(dir, 1, r, c, size(network))...][χN Dt60 Db60; χNE] * + rotl60(ket(O_r_c), mod(dir - 1, 6))[d; Dt120 Dt60 Dt0 Dt300 Dt240 Dt180] * + conj(rotl60(bra(O_r_c), mod(dir - 1, 6))[d; Db120 Db60 Db0 Db300 Db240 Db180]) + return mat / norm(mat) +end + +function semi_renormalize_edge(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Pas, Pbs, dir, r, c) where {P <: PFTensorTriangular} + @tensor opt = true mat[χW D240; χE D300] := + Pas[_coordinates(dir, 5, r, c, size(network))...][χW; χNW D180] * + Pbs[dir, r, c][χNE D0; χE] * + env.Ea[_coordinates(dir, 0, r, c, size(network))...][χNW D120; χN] * + env.Eb[_coordinates(dir, 1, r, c, size(network))...][χN D60; χNE] * + rotl60(network[r, c], mod(dir - 1, 6))[D180 D240 D300; D120 D60 D0] + return mat / norm(mat) +end diff --git a/src/algorithms/ctmrg_triangular/ctmrg.jl b/src/algorithms/ctmrg_triangular/ctmrg.jl new file mode 100644 index 000000000..235ef8834 --- /dev/null +++ b/src/algorithms/ctmrg_triangular/ctmrg.jl @@ -0,0 +1,77 @@ +abstract type CTMRGAlgorithmTriangular end + +function CTMRGAlgorithmTriangular(; alg = :SimultaneousCTMRGTriangular, kwargs...) + alg == :SimultaneousCTMRGTriangular || throw(ArgumentError("unknown triangular CTMRG algorithm: $alg")) + return SimultaneousCTMRGTriangular(; kwargs...) +end + +""" + leading_boundary(env₀::CTMRGEnvTriangular, network; alg=:SimultaneousCTMRGTriangular, kwargs...) + leading_boundary(env₀::CTMRGEnvTriangular, network, alg::CTMRGAlgorithmTriangular) + +Keyword arguments are forwarded to `SimultaneousCTMRGTriangular`: `tol`, `maxiter`, +`miniter`, `verbosity`, `decomposition_alg`, `conditioning`, `projector_alg`, and `trunctype`. + +Return the contracted environment and convergence information. `info.converged` and +`info.convergence_error` describe convergence; `info.convergence_metric` is retained as +an alias for the latter. Convergence is the maximum singular-value distance of all +six corners and twelve edge halves at every site, following PEPSKit's native tensor +bipartitions and sector-aware distance. Spectral convergence does not imply an +elementwise fixed point. `info.contraction_metrics` contains `corner_truncation_error` +and `edge_truncation_error`, maximized over all sites and directions of the final +iteration, and their maximum `truncation_error`. + +Errors are discarded norms divided by the retained SVD spectrum norm, following +PEPSKit's corner-projector convention. They include singular values suppressed by +the inverse-square-root cutoff, but exclude preliminary seed and conditioning SVDs. +These are local truncation diagnostics, not errors on the final observables. +""" +function leading_boundary(env₀::CTMRGEnvTriangular, network::InfiniteTriangularNetwork; kwargs...) + alg = select_algorithm(leading_boundary, env₀; kwargs...) + return leading_boundary(env₀, network, alg) +end +function leading_boundary( + env₀::CTMRGEnvTriangular, network::InfiniteTriangularNetwork, alg::CTMRGAlgorithmTriangular + ) + log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) + return with_verbosity(; alg.verbosity) do + env = deepcopy(env₀) + CS, TS = ignore_derivatives() do + return convergence_spectra(env₀, alg) + end + η = one(real(scalartype(network))) + ctmrg_loginit!(log, η, network, env₀) + local info_iter + converged = false + for iter in 1:(alg.maxiter) + env, info_iter = ctmrg_iteration(network, env, alg) + η, CS, TS = calc_convergence(env, CS, TS, alg) + + if η ≤ alg.tol && iter ≥ alg.miniter + ctmrg_logfinish!(log, iter, η, network, env) + converged = true + break + end + if iter == alg.maxiter + ctmrg_logcancel!(log, iter, η, network, env) + else + ctmrg_logiter!(log, iter, η, network, env) + end + end + info = (; + converged, convergence_error = η, convergence_metric = η, + contraction_metrics = info_iter.contraction_metrics, + ) + return env, info + end +end +function leading_boundary(env₀::CTMRGEnvTriangular, state, args...; kwargs...) + return leading_boundary(env₀, InfiniteTriangularNetwork(state), args...; kwargs...) +end + +convergence_tensors(env::CTMRGEnvTriangular, alg) = (env.C, [env.Ea; env.Eb]) + +function convergence_spectra(env::CTMRGEnvTriangular, alg) + corners, edges = convergence_tensors(env, alg) + return corner_spectra(corners, alg), edge_spectra(edges, alg) +end diff --git a/src/algorithms/ctmrg_triangular/gaugefix.jl b/src/algorithms/ctmrg_triangular/gaugefix.jl new file mode 100644 index 000000000..e7e250b74 --- /dev/null +++ b/src/algorithms/ctmrg_triangular/gaugefix.jl @@ -0,0 +1,77 @@ +function gauge_fix(envfinal::CTMRGEnvTriangular, envprev::CTMRGEnvTriangular, alg::ScramblingEnvGauge) + signs, corner_phases, edge_phases = compute_gauge_fix_gauge(envfinal, envprev, alg) + return fix_phases(envfinal, signs, corner_phases, edge_phases) +end + +function compute_gauge_fix_gauge(envfinal::CTMRGEnvTriangular, envprev::CTMRGEnvTriangular, alg::ScramblingEnvGauge) + same_spaces = all(zip((envfinal.C, envfinal.Ea, envfinal.Eb), (envprev.C, envprev.Ea, envprev.Eb))) do (final, prev) + return size(final) == size(prev) && all(space.(final) .== space.(prev)) + end + same_spaces || throw(ArgumentError("Spaces of envfinal and envprev are not the same")) + signs = compute_relative_phases(envfinal, envprev, alg) + corner_phases, edge_phases = compute_global_phases(fix_relative_phases(envfinal, signs), envprev) + return signs, corner_phases, edge_phases +end + +# Follow the boundary tangent, alternating the two edge halves at each site. +function _triangular_edge_chain(env::CTMRGEnvTriangular, dir, r, c) + tensors = [env.Eb[dir, r, c], env.Ea[dir, r, c]] + unitcell = size(env)[2:3] + _, r′, c′ = _coordinates(dir, 2, r, c, unitcell) + while (r′, c′) != (r, c) + push!(tensors, env.Eb[dir, r′, c′], env.Ea[dir, r′, c′]) + _, r′, c′ = _coordinates(dir, 2, r′, c′, unitcell) + end + return tensors +end + +function compute_relative_phases(envfinal::CTMRGEnvTriangular, envprev::CTMRGEnvTriangular, ::ScramblingEnvGauge) + signs = map((1, 0)) do shift + return map(eachcoordinate(envfinal, 1:6)) do (dir, r, c) + Tsprev = circshift(_triangular_edge_chain(envprev, dir, r, c), -shift) + Tsfinal = circshift(_triangular_edge_chain(envfinal, dir, r, c), -shift) + M = map(t -> randn(storagetype(t), space(t)), Tsfinal) + ρprev = right_transfermatrix_fixedpoint(Tsprev, M, Arnoldi()) + ρfinal = right_transfermatrix_fixedpoint(Tsfinal, M, Arnoldi()) + Qprev, = left_orth!(ρprev; positive = true) + Qfinal, = left_orth!(ρfinal; positive = true) + return Qprev * Qfinal' + end + end + return (; Ea = signs[1], Eb = signs[2]) +end + +function fix_relative_phases(env::CTMRGEnvTriangular, signs) + unitcell = size(env)[2:3] + # Triangular corners have the same tensor layout as the edge halves. + C = map(eachcoordinate(env, 1:6)) do (dir, r, c) + next = _coordinates(dir, 2, r, c, unitcell) + return fix_gauge_edge(env.C[dir, r, c], signs.Eb[mod1(dir - 1, 6), r, c], signs.Eb[next...]) + end + Ea = map(eachcoordinate(env, 1:6)) do (dir, r, c) + next = _coordinates(dir, 2, r, c, unitcell) + return fix_gauge_edge(env.Ea[dir, r, c], signs.Ea[dir, r, c], signs.Eb[next...]) + end + Eb = map(eachcoordinate(env, 1:6)) do (dir, r, c) + return fix_gauge_edge(env.Eb[dir, r, c], signs.Eb[dir, r, c], signs.Ea[dir, r, c]) + end + return CTMRGEnvTriangular(C, Ea, Eb) +end + +function compute_global_phases(envfix::CTMRGEnvTriangular, envprev::CTMRGEnvTriangular) + corner_phases = dot.(envprev.C, envfix.C) + edge_phases = (; Ea = dot.(envprev.Ea, envfix.Ea), Eb = dot.(envprev.Eb, envfix.Eb)) + return corner_phases, edge_phases +end + +function fix_global_phases(env::CTMRGEnvTriangular, corner_phases, edge_phases) + return CTMRGEnvTriangular(env.C .* inv.(corner_phases), env.Ea .* inv.(edge_phases.Ea), env.Eb .* inv.(edge_phases.Eb)) +end + +function fix_phases(env::CTMRGEnvTriangular, signs, corner_phases, edge_phases) + return fix_global_phases(fix_relative_phases(env, signs), corner_phases, edge_phases) +end + +function calc_elementwise_convergence(envfinal::CTMRGEnvTriangular, envfix::CTMRGEnvTriangular; atol::Real = 1.0e-6) + return max(norm(envfinal.C .- envfix.C, Inf), norm(envfinal.Ea .- envfix.Ea, Inf), norm(envfinal.Eb .- envfix.Eb, Inf)) +end diff --git a/src/algorithms/ctmrg_triangular/simultaneous.jl b/src/algorithms/ctmrg_triangular/simultaneous.jl new file mode 100644 index 000000000..816e44b5c --- /dev/null +++ b/src/algorithms/ctmrg_triangular/simultaneous.jl @@ -0,0 +1,258 @@ +""" +$(TYPEDEF) + +Simultaneous CTMRG for triangular PEPS and classical networks, with direct Ta/Tb +updates of all six edge directions. See [`leading_boundary`](@ref) for convergence and truncation diagnostics. + +## Constructors + + SimultaneousCTMRGTriangular(; kwargs...) + +* `tol::Real=$(Defaults.ctmrg_tol)` : Convergence tolerance for corner and edge spectra. +* `maxiter::Int=$(Defaults.ctmrg_maxiter)` : Maximum number of CTMRG iterations. +* `miniter::Int=$(Defaults.ctmrg_miniter)` : Minimum number of CTMRG iterations. +* `verbosity::Int=$(Defaults.ctmrg_verbosity)` : Output verbosity, as in [`leading_boundary`](@ref). +* `decomposition_alg::Union{SVDAdjoint,NamedTuple}=(;)` : SVD algorithm and reverse rule, + as in square CTMRG. A `NamedTuple` specifies the forward algorithm; use [`SVDAdjoint`](@ref) + to configure its pullback. Full corner and edge conditioning require a dense SVD. +* `projector_alg::Symbol=:full` : Corner-projector construction, either `:twothirds` or `:full`. +* `conditioning::Bool=true` : Square-root SVD conditioning of the Ta/Tb edge split. + This does not change the corner projectors; `:full` always splits its opposite + corner factor with a square-root SVD, while `:twothirds` uses its factors directly. +* `trunctype::Symbol=:truncrank` : Keep the initial environment dimension. + Use `:FixedSpaceTruncation` to preserve each environment bond space, including its sectors. +""" +struct SimultaneousCTMRGTriangular{S <: SVDAdjoint} <: CTMRGAlgorithmTriangular + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + conditioning::Bool + projector_alg::Symbol + trunctype::Symbol + decomposition_alg::S +end +function SimultaneousCTMRGTriangular(; + tol = Defaults.ctmrg_tol, + maxiter = Defaults.ctmrg_maxiter, + miniter = Defaults.ctmrg_miniter, + verbosity = Defaults.ctmrg_verbosity, + conditioning = true, + projector_alg = :full, + trunctype = :truncrank, + decomposition_alg = (;) + ) + if decomposition_alg isa NamedTuple + decomposition_alg = SVDAdjoint(; fwd_alg = decomposition_alg) + end + return SimultaneousCTMRGTriangular(tol, maxiter, miniter, verbosity, conditioning, projector_alg, trunctype, decomposition_alg) +end + +# Based on +# https://arxiv.org/pdf/2510.04907 + +function ctmrg_iteration(network::InfiniteTriangularNetwork, env::CTMRGEnvTriangular, alg::SimultaneousCTMRGTriangular) + trunc_type = getfield(@__MODULE__, alg.trunctype) + if trunc_type == FixedSpaceTruncation + trunc = FixedSpaceTruncation() + else + trunc = trunc_type(dim(domain(env.C[1, 1, 1])[1])) + end + Pas, Pbs, S, corner_truncation_error = calculate_projectors(network, env, trunc, alg.projector_alg, alg.decomposition_alg) + + env = renormalize_corners!(network, env, Pas, Pbs) + env = normalize_corners(env) + + mats, env = semi_renormalize(network, env, Pas, Pbs, trunc, alg.decomposition_alg) + env, edge_truncation_error = renormalize_edges(network, env, mats, trunc, alg.decomposition_alg; alg.conditioning) + truncation_error = max(corner_truncation_error, edge_truncation_error) + info = (; + contraction_metrics = (; truncation_error, corner_truncation_error, edge_truncation_error), S, + ) + return env, info +end + +function calculate_projectors(network, env, trunc, projector_alg, decomposition_alg = SVDAdjoint()) + if projector_alg == :full + return calculate_full_projectors(network, env, trunc, decomposition_alg) + elseif projector_alg == :twothirds + return calculate_twothirds_projectors(network, env, trunc, decomposition_alg) + else + @error "projector_alg = $projector_alg not defined" + end +end + +# TODO: remove once this once everything is plumbed through a proper projector algorithm +_truncation_strategy(trunc::TruncationStrategy, _) = trunc +function _truncation_strategy(::FixedSpaceTruncation, edge) + tspace = space(edge, 1) + return isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace) +end + +# Include retained values suppressed by the inverse square root, without subtracting +# nearly equal norms. Normalize by norm(S), as in the square CTMRG projectors. +function _relative_truncation_error(S, isqrtS, truncation_error) + return ignore_derivatives() do + discarded_S = DiagonalTensorMap(S.data .* iszero.(isqrtS.data), space(S, 1)) + discarded_norm = hypot(truncation_error, norm(discarded_S)) + retained_norm = norm(S) + iszero(retained_norm) && return iszero(discarded_norm) ? zero(discarded_norm) : oftype(discarded_norm, Inf) + return discarded_norm / retained_norm + end +end + +_normalize_corner_environment(ρ) = ρ / norm(ρ) + +# The double corners around the elementary triangle of projector `[dir, r, c]` sit at +# `(r, c)`, its neighbor in direction `dir + 2` and its neighbor in direction `dir + 3`. +# `Pb` acts on the leg of `(r, c)`, `Pa` on the leg of its neighbor in direction `dir + 2`. +function calculate_twothirds_projectors(network::InfiniteTriangularNetwork{P}, env, trunc, decomposition_alg = SVDAdjoint()) where {P} + projectors = dtmap(eachcoordinate(network, 1:6)) do (dir, r, c) + ρL = build_double_corner_matrix_triangular(network, env, _prev(dir, 6), r, c) + ρR = build_double_corner_matrix_triangular(network, env, _next(dir, 6), _coordinates(dir, 2, r, c, size(network))[2:3]...) + ρL = _normalize_corner_environment(ρL) + ρR = _normalize_corner_environment(ρR) + ρρ = _normalize_corner_environment(ρL * ρR) + + trunc′ = _truncation_strategy(trunc, env.Ea[dir, r, c]) + svd_alg = _set_decomposition_truncation(decomposition_alg, trunc′) + U, S, V, truncation_error = svd_trunc(ρρ, svd_alg) + iqsrtS = sdiag_pow(S, -0.5) + truncation_error = _relative_truncation_error(S, iqsrtS, truncation_error) + + Pb = ρR * V' * iqsrtS + Pa = iqsrtS * U' * ρL + return Pa, Pb, S, truncation_error + end + return getindex.(projectors, 1), getindex.(projectors, 2), getindex.(projectors, 3), maximum(x -> x[4], projectors) +end + +function calculate_full_projectors(network::InfiniteTriangularNetwork{P}, env, trunc, decomposition_alg = SVDAdjoint()) where {P} + projectors = dtmap(eachcoordinate(network, 1:6)) do (dir, r, c) + ρL = build_double_corner_matrix_triangular(network, env, mod1(dir - 1, 6), r, c) + ρR = build_double_corner_matrix_triangular(network, env, mod1(dir + 1, 6), _coordinates(dir, 2, r, c, size(network))[2:3]...) + ρ̄ = build_double_corner_matrix_triangular(network, env, mod1(dir + 3, 6), _coordinates(dir, 3, r, c, size(network))[2:3]...) + ρ̄ = _normalize_corner_environment(ρ̄) + Ū, S̄, V̄ᴴ = svd_full(ρ̄, decomposition_alg) + sqrtS̄ = sdiag_pow(S̄, 0.5) + ρ̄ᴿ = Ū * sqrtS̄ + ρ̄ᴸ = sqrtS̄ * V̄ᴴ + ρL = _normalize_corner_environment(ρ̄ᴸ * ρL) + ρR = _normalize_corner_environment(ρR * ρ̄ᴿ) + ρρ = _normalize_corner_environment(ρL * ρR) + + trunc′ = _truncation_strategy(trunc, env.Ea[dir, r, c]) + svd_alg = _set_decomposition_truncation(decomposition_alg, trunc′) + U, S, Vᴴ, truncation_error = svd_trunc(ρρ, svd_alg) + isqrtS = sdiag_pow(S, -0.5) + truncation_error = _relative_truncation_error(S, isqrtS, truncation_error) + Pb = ρR * Vᴴ' * isqrtS + Pa = isqrtS * U' * ρL + + return Pa, Pb, S, truncation_error + end + return getindex.(projectors, 1), getindex.(projectors, 2), getindex.(projectors, 3), maximum(x -> x[4], projectors) +end + +function renormalize_corners!(network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PEPSSandwichTriangular} + coordinates = eachcoordinate(network, 1:6) + new_corners′ = similar(env.C) + new_corners = dtmap!!(new_corners′, coordinates) do (dir, r, c) + return renormalize_corner_triangular((dir, r, c), network, env, Pas, Pbs) + end + return CTMRGEnvTriangular(new_corners, env.Ea, env.Eb) +end + +function renormalize_corners!(network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PFTensorTriangular} + coordinates = eachcoordinate(network, 1:6) + new_corners′ = similar(env.C) + new_corners = dtmap!!(new_corners′, coordinates) do (dir, r, c) + return renormalize_corner_triangular((dir, r, c), network, env, Pas, Pbs) + end + return CTMRGEnvTriangular(new_corners, env.Ea, env.Eb) +end + +function _permute_edge(t::Union{T1, T2}) where {E, S, T1 <: AbstractTensorMap{E, S, 2, 1}, T2 <: AbstractTensorMap{E, S, 1, 2}} + return permute(t, ((1, 3), (2,))) +end + +function _permute_edge(t::Union{T1, T2}) where {E, S, T1 <: AbstractTensorMap{E, S, 3, 1}, T2 <: AbstractTensorMap{E, S, 1, 3}} + return permute(t, ((1, 3, 4), (2,))) +end + +function semi_renormalize(network::InfiniteTriangularNetwork, env::CTMRGEnvTriangular, Pas, Pbs, trunc, decomposition_alg = SVDAdjoint()) + edges = dtmap(eachcoordinate(network, 1:6)) do (dir, r, c) + mat = semi_renormalize_edge(network, env, Pas, Pbs, dir, r, c) + trunc′ = _truncation_strategy(trunc, env.Ea[dir, r, c]) + svd_alg = _set_decomposition_truncation(decomposition_alg, trunc′) + U, S, V = svd_trunc(mat, svd_alg) + sqrtS = sdiag_pow(S, 0.5) + Ẽb = U * sqrtS + Ẽa = _permute_edge(sqrtS * V) + return mat, Ẽa, Ẽb + end + return getindex.(edges, 1), CTMRGEnvTriangular(env.C, getindex.(edges, 2), getindex.(edges, 3)) +end + +function _truncate_edge(mat, ρL, ρR, trunc, decomposition_alg = SVDAdjoint(); conditioning = true) + ρL /= norm(ρL) + ρR /= norm(ρR) + if conditioning + _, SL, VLᴴ = svd_full(ρL, decomposition_alg) + UR, SR, _ = svd_full(ρR, decomposition_alg) + sqrtSL = sdiag_pow(SL, 0.5) + sqrtSR = sdiag_pow(SR, 0.5) + ρL = sqrtSL * VLᴴ + ρR = UR * sqrtSR + end + + ρ = ρL * mat * ρR + + ρ /= norm(ρ) + svd_alg = _set_decomposition_truncation(decomposition_alg, trunc) + U, S, Vᴴ, truncation_error = svd_trunc(ρ, svd_alg) + isqrtS = sdiag_pow(S, -0.5) + truncation_error = _relative_truncation_error(S, isqrtS, truncation_error) + Ea = _permute_edge(isqrtS * U' * ρL * mat) + Eb = mat * ρR * Vᴴ' * isqrtS + return Ea, Eb, truncation_error +end + +function renormalize_edges(network::InfiniteTriangularNetwork, env::CTMRGEnvTriangular, mats, trunc, decomposition_alg = SVDAdjoint(); conditioning = true) + edges = dtmap(eachcoordinate(network, 1:6)) do (dir, r, c) + ρL, ρR = build_halfinfinite_environments(network, env, dir, r, c) + trunc′ = _truncation_strategy(trunc, env.Ea[dir, r, c]) + return _truncate_edge(mats[dir, r, c], ρL, ρR, trunc′, decomposition_alg; conditioning) + end + env′ = normalize_edges(CTMRGEnvTriangular(env.C, getindex.(edges, 1), getindex.(edges, 2))) + return env′, maximum(x -> x[3], edges) +end + +function normalize_corners(env::CTMRGEnvTriangular) + C_normalized = map(env.C) do C + return C / norm(C) + end + return CTMRGEnvTriangular(C_normalized, env.Ea, env.Eb) +end + +function normalize_edges(env::CTMRGEnvTriangular) + Ea_normalized = map(env.Ea) do e + return e / norm(e) + end + Eb_normalized = map(env.Eb) do e + return e / norm(e) + end + return CTMRGEnvTriangular(env.C, Ea_normalized, Eb_normalized) +end + +function calculate_error(Ss, Ss_prev) + ε = Inf + for (S, S_prev) in zip(Ss, Ss_prev) + if space(S) == space(S_prev) + ε = norm(S^4 - S_prev^4) + else + return Inf + end + end + return ε +end diff --git a/src/algorithms/ctmrg_triangular/utility.jl b/src/algorithms/ctmrg_triangular/utility.jl new file mode 100644 index 000000000..fc294266c --- /dev/null +++ b/src/algorithms/ctmrg_triangular/utility.jl @@ -0,0 +1,60 @@ +function _coordinates(dir::Int, rot::Int, r::Int, c::Int, unitcell::Tuple{Int, Int}) + rows, cols = unitcell + if mod1(dir + rot, 6) == 1 + return (mod1(dir, 6), mod1(r - 1, rows), mod1(c, cols)) + elseif mod1(dir + rot, 6) == 2 + return (mod1(dir, 6), mod1(r - 1, rows), mod1(c + 1, cols)) + elseif mod1(dir + rot, 6) == 3 + return (mod1(dir, 6), mod1(r, rows), mod1(c + 1, cols)) + elseif mod1(dir + rot, 6) == 4 + return (mod1(dir, 6), mod1(r + 1, rows), mod1(c, cols)) + elseif mod1(dir + rot, 6) == 5 + return (mod1(dir, 6), mod1(r + 1, rows), mod1(c - 1, cols)) + elseif mod1(dir + rot, 6) == 6 + return (mod1(dir, 6), mod1(r, rows), mod1(c - 1, cols)) + end +end + +function _coordinates(dir::Int, r::Int, c::Int, unitcell::Tuple{Int, Int}) + return _coordinates(dir, 0, r, c, unitcell)[2:3] +end + +function _truncway(trunctype::Symbol, D::E) where {E <: ElementarySpace} + if trunctype == :truncdim + return truncdim(dim(D)) + else + return notrunc() + end +end + +function rotr60(A::T, dir::Int) where {E, S, T <: AbstractTensorMap{E, S, 1, 6}} + if dir == 0 + return A + else + return rotr60(permute(A, ((1,), (7, 2, 3, 4, 5, 6))), dir - 1) + end +end + +function rotr60(A::T, dir::Int) where {E, S, T <: AbstractTensorMap{E, S, 3, 3}} + if dir == 0 + return A + else + return rotr60(permute(A, ((2, 3, 6), (1, 4, 5))), dir - 1) + end +end + +function rotl60(A::T, dir::Int) where {E, S, T <: AbstractTensorMap{E, S, 1, 6}} + if dir == 0 + return A + else + return rotl60(permute(A, ((1,), (3, 4, 5, 6, 7, 2))), dir - 1) + end +end + +function rotl60(A::T, dir::Int) where {E, S, T <: AbstractTensorMap{E, S, 3, 3}} + if dir == 0 + return A + else + return rotl60(permute(A, ((4, 1, 2), (5, 6, 3))), dir - 1) + end +end diff --git a/src/algorithms/expectation_value/expectation_value.jl b/src/algorithms/expectation_value/expectation_value.jl index e8d25e3eb..322aec7bd 100644 --- a/src/algorithms/expectation_value/expectation_value.jl +++ b/src/algorithms/expectation_value/expectation_value.jl @@ -9,17 +9,21 @@ Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / This can be done either for a PEPS, or alternatively for a density matrix PEPO. In the latter case the first signature corresponds to a single layer PEPO contraction, while the second signature yields a bilayer contraction instead. + +For an `InfinitePEPSTriangular`, dense one-site and nearest-neighbor terms are supported. +The three bond families have offsets `(0, 1)`, `(1, -1)`, and `(1, 0)`. """ function MPSKit.expectation_value( bra::S, O::LocalOperator, ket::S, env - ) where {S <: InfiniteState} + ) where {S <: Union{InfiniteState, InfinitePEPSTriangular}} checklattice(bra, O, ket) + isempty(O.terms) && return zero(scalartype(ket)) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly return local_expectation_value(inds, bra, operator, ket, env) end return sum(term_vals) end -MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) +MPSKit.expectation_value(peps::Union{InfinitePEPS, InfinitePEPSTriangular}, O::LocalOperator, env) = expectation_value(peps, O, peps, env) function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) checklattice(state, O) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl index 37ae322b4..ccf1f3ebf 100644 --- a/src/algorithms/optimization/implicit_differentiation.jl +++ b/src/algorithms/optimization/implicit_differentiation.jl @@ -202,21 +202,29 @@ Evaluating the gradient of the cost function for CTMRG: - With explicit evaluation of the geometric sum, the gradient is computed by differentiating the cost function with the environment kept fixed, and then manually adding the gradient contributions from the environments. =# -_scrambling_env_gauge(::CTMRGAlgorithm) = ScramblingEnvGauge() +_scrambling_env_gauge(::Union{CTMRGAlgorithm, CTMRGAlgorithmTriangular}) = ScramblingEnvGauge() +_scrambling_env_gauge(::C4vCTMRG) = ScramblingEnvGaugeC4v() function _set_fixed_truncation(alg::CTMRGAlgorithm) alg_fixed = @set alg.projector_alg = _set_truncation(alg.projector_alg, FixedSpaceTruncation()) return alg_fixed end +function _set_fixed_truncation(alg::SimultaneousCTMRGTriangular) + return @set alg.trunctype = :FixedSpaceTruncation +end + +_ctmrg_network(state, ::CTMRGAlgorithm) = InfiniteSquareNetwork(state) +_ctmrg_network(state, ::CTMRGAlgorithmTriangular) = InfiniteTriangularNetwork(state) + # compute the CTMRG gradient through fixed-point differentiation function _rrule( gradmode::FixedPointGradient, config::RuleConfig, ::typeof(MPSKit.leading_boundary), - envinit::AbstractCTMRGEnv, + envinit::Union{CTMRGEnv, CTMRGEnvTriangular}, state, - alg::CTMRGAlgorithm, + alg::Union{CTMRGAlgorithm, CTMRGAlgorithmTriangular}, ) _check_algorithm_combination(alg, gradmode) @@ -225,13 +233,12 @@ function _rrule( # prepare iterating function corresponding to a single gauge-fixed CTMRG iteration alg_fixed = _set_fixed_truncation(alg) # fix spaces during differentiation alg_gauge = _scrambling_env_gauge(alg) # select appropriate gauge-fixing algorithm - env_conv, info = ctmrg_iteration(InfiniteSquareNetwork(state), env, alg_fixed) - signs, = compute_gauge_fix_gauge(env_conv, env, alg_gauge) - # fix global phases relative to the fixed point inside the iteration + env_conv, info = ctmrg_iteration(_ctmrg_network(state, alg), env, alg_fixed) + signs, corner_phases, edge_phases = compute_gauge_fix_gauge(env_conv, env, alg_gauge) function gauge_fixed_iteration(A, x) - return fix_global_phases( - fix_relative_phases(ctmrg_iteration(InfiniteSquareNetwork(A), x, alg_fixed)[1], signs), - env, + return fix_phases( + ctmrg_iteration(_ctmrg_network(A, alg), x, alg_fixed)[1], + signs, corner_phases, edge_phases, ) end # prepare its pullback @@ -413,6 +420,17 @@ function _gauge_fix_c4v_projector(::C4vCTMRG{<:C4vQRProjector}, signs, info) return info.Q * signs[1]' end +function _check_algorithm_combination(::CTMRGAlgorithmTriangular, ::ImplicitGradient) + throw(ArgumentError("triangular CTMRG supports FixedPointGradient; characteristic equations for ImplicitGradient are not implemented")) +end + +function _rrule( + gradmode::ImplicitGradient, config::RuleConfig, + ::typeof(MPSKit.leading_boundary), envinit, state, alg::CTMRGAlgorithmTriangular, + ) + return _check_algorithm_combination(alg, gradmode) +end + function _check_algorithm_combination(::SequentialCTMRG, ::ImplicitGradient) msg = "The `:ImplicitGradient` algorithm is not yet compatible with the `SequentialCTMRG` algorithm; select \ SimultaneousCTMRG instead to use :ImplicitGradient" diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 808ed04f5..04a9eadb9 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -1,10 +1,10 @@ """ - cost_function(peps::InfinitePEPS, env, O::LocalOperator) + cost_function(peps::Union{InfinitePEPS, InfinitePEPSTriangular}, env, O::LocalOperator) Real part of expectation value of `O`, used as the cost function in variational PEPS optimization. Prints a warning if the expectation value yields a finite imaginary part (up to a tolerance). """ -function cost_function(peps::InfinitePEPS, env, O::LocalOperator) +function cost_function(peps::Union{InfinitePEPS, InfinitePEPSTriangular}, env, O::LocalOperator) E = MPSKit.expectation_value(peps, O, env) ignore_derivatives() do isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || @@ -62,7 +62,13 @@ function PEPSOptimize(; boundary_alg = (;), gradient_alg = (;), optimizer_alg = (;), precondition_alg = (;), reuse_env = Defaults.reuse_env, symmetrization = nothing, ) - boundary_algorithm = _alg_or_nt(CTMRGAlgorithm, boundary_alg) + boundary_type = if parent_alg(boundary_alg) isa CTMRGAlgorithmTriangular || + (boundary_alg isa NamedTuple && get(boundary_alg, :alg, nothing) == :SimultaneousCTMRGTriangular) + CTMRGAlgorithmTriangular + else + CTMRGAlgorithm + end + boundary_algorithm = _alg_or_nt(boundary_type, boundary_alg) gradient_algorithm = _alg_or_nt(GradientAlgorithm, gradient_alg) optimizer_algorithm = _alg_or_nt(OptimKit.OptimizationAlgorithm, optimizer_alg) precondition_algorithm = _alg_or_nt(PreconditionAlgorithm, precondition_alg) @@ -108,9 +114,9 @@ function _OptimizationAlgorithm(; end """ - fixedpoint(operator, peps₀::InfinitePEPS, env₀; kwargs...) -> peps_final, env_final, cost_final, info + fixedpoint(operator, peps₀::Union{InfinitePEPS, InfinitePEPSTriangular}, env₀; kwargs...) -> peps_final, env_final, cost_final, info # expert version: - fixedpoint(operator, peps₀::InfinitePEPS, env₀, alg::PEPSOptimize; finalize!=OptimKit._finalize!) + fixedpoint(operator, peps₀::Union{InfinitePEPS, InfinitePEPSTriangular}, env₀, alg::PEPSOptimize; finalize!=OptimKit._finalize!) Find the fixed point of `operator` (i.e. the ground state) starting from `peps₀` according to the supplied optimization parameters. The initial environment `env₀` serves as an @@ -118,7 +124,11 @@ initial guess for the first boundary contraction run. By default, a random initi Each local tensor is normalized initially, and optimization preserves its Euclidean norm by projecting gradients and preconditioned directions onto its real tangent space. The optimization parameters can be supplied via the keyword arguments or directly as a -`PEPSOptimize` struct. The following keyword arguments are supported: +`PEPSOptimize` struct. Triangular PEPS support `FixedPointGradient` with a complex +environment, or `gradient_alg = nothing` to differentiate the CTMRG iterations. +They require `precondition_alg = nothing`; symmetrization is not yet supported. + +The following keyword arguments are supported: ## Keyword arguments @@ -223,13 +233,13 @@ information `NamedTuple` which contains the following entries: * `gradnorms_unitcell` : History of gradient norms for each respective unit cell entry. * `times` : History of optimization step execution times. """ -function fixedpoint(operator, peps₀::InfinitePEPS, env₀; kwargs...) +function fixedpoint(operator, peps₀::Union{InfinitePEPS, InfinitePEPSTriangular}, env₀; kwargs...) extra_kwarg_keys = (:finalize!, :hasconverged, :shouldstop) # these will not be passed to `select_algorithm`, only to the 2nd `fixedpoint` call alg = select_algorithm(fixedpoint, env₀; filter(kw -> !(first(kw) in extra_kwarg_keys), kwargs)...) return fixedpoint(operator, peps₀, env₀, alg; filter(kw -> first(kw) in extra_kwarg_keys, kwargs)...) end function fixedpoint( - operator, peps₀::InfinitePEPS, env₀, alg::PEPSOptimize; + operator, peps₀::Union{InfinitePEPS, InfinitePEPSTriangular}, env₀, alg::PEPSOptimize; (finalize!) = OptimKit._finalize!, hasconverged = OptimKit.DefaultHasConverged(alg.optimizer_alg.gradtol), shouldstop = OptimKit.DefaultShouldStop(alg.optimizer_alg.maxiter), @@ -347,15 +357,30 @@ function check_input(::typeof(fixedpoint), peps₀, env₀, alg::PEPSOptimize) return nothing end +function check_input(::typeof(fixedpoint), peps₀::InfinitePEPSTriangular, env₀, alg::PEPSOptimize) + env₀ isa CTMRGEnvTriangular || + throw(ArgumentError("triangular PEPS optimization requires a CTMRGEnvTriangular environment")) + parent_alg(alg.boundary_alg) isa CTMRGAlgorithmTriangular || + throw(ArgumentError("triangular PEPS optimization requires a triangular CTMRG algorithm")) + if parent_alg(alg.gradient_alg) isa FixedPointGradient && scalartype(env₀) <: Real + throw(ArgumentError("triangular fixed-point differentiation requires a complex environment because gauge fixing generally produces complex gauges")) + end + isnothing(alg.precondition_alg) || + throw(ArgumentError("triangular preconditioning is not supported; set precondition_alg = nothing")) + isnothing(alg.symmetrization) || + throw(ArgumentError("triangular symmetrization is not supported; set symmetrization = nothing")) + return nothing +end + """ - peps_normalize(A::InfinitePEPS) + peps_normalize(A::Union{InfinitePEPS, InfinitePEPSTriangular}) -Normalize the individual tensors in the unit cell of an `InfinitePEPS` such that they each +Normalize the individual tensors in the PEPS unit cell such that they each have unit Euclidean norm. """ -function peps_normalize(A::InfinitePEPS) +function peps_normalize(A::Union{InfinitePEPS, InfinitePEPSTriangular}) normalized_tensors = normalize.(unitcell(A)) - return InfinitePEPS(normalized_tensors) + return typeof(A)(normalized_tensors) end """Real part of the Hermitian inner product.""" @@ -368,11 +393,11 @@ $(SIGNATURES) Project each direction onto the tangent space of its fixed-norm PEPS tensor using [`project_tangent`](@ref) with the real part of the Hermitian inner product. Complex tensors are treated as a real vector space, retaining the tangent phase direction ``i A``. """ -function peps_project_tangent(η::InfinitePEPS, peps::InfinitePEPS) +function peps_project_tangent(η::Union{InfinitePEPS, InfinitePEPSTriangular}, peps::Union{InfinitePEPS, InfinitePEPSTriangular}) tangents = map(unitcell(η), unitcell(peps)) do direction, A return project_tangent(direction, A; inner = real_inner) end - return InfinitePEPS(tangents) + return typeof(peps)(tangents) end """ @@ -396,8 +421,8 @@ function peps_retract(x, η, α) env = deepcopy(x[2]) retractions = norm_preserving_retract.(unitcell(peps), unitcell(η), α; inner = real_inner) - peps′ = InfinitePEPS(map(first, retractions)) - ξ = InfinitePEPS(map(last, retractions)) + peps′ = typeof(peps)(map(first, retractions)) + ξ = typeof(peps)(map(last, retractions)) return (peps′, env), ξ end diff --git a/src/algorithms/select_algorithm.jl b/src/algorithms/select_algorithm.jl index 2671a8053..34d128ba2 100644 --- a/src/algorithms/select_algorithm.jl +++ b/src/algorithms/select_algorithm.jl @@ -158,3 +158,10 @@ function select_algorithm( return CTMRGAlgorithm(; alg, tol, verbosity, decomposition_alg, kwargs...) end + +function select_algorithm( + ::typeof(leading_boundary), env₀::CTMRGEnvTriangular; + alg = :SimultaneousCTMRGTriangular, kwargs..., + ) + return CTMRGAlgorithmTriangular(; alg, kwargs...) +end diff --git a/src/algorithms/toolbox_triangular.jl b/src/algorithms/toolbox_triangular.jl new file mode 100644 index 000000000..f767fb0c5 --- /dev/null +++ b/src/algorithms/toolbox_triangular.jl @@ -0,0 +1,86 @@ +""" + network_value(network::InfiniteTriangularNetwork, env::CTMRGEnvTriangular) + +Return the value per unit cell using the principal cube root of the contraction-ratio +product. Boundary phases cancel before taking the cube root. For general complex +networks, this estimator determines the phase only modulo `2π/3`. +""" +function network_value(network::InfiniteTriangularNetwork, env::CTMRGEnvTriangular) + value = prod(Iterators.product(axes(network)...)) do (r, c) + nw_corners = complex(_contract_corners((r, c), network, env)) + nw_full = complex(_contract_site_large((r, c), network, env)) + nw_0 = complex(_contract_edges_0((r, c), network, env)) + nw_60 = complex(_contract_edges_60((r, c), network, env)) + nw_120 = complex(_contract_edges_120((r, c), network, env)) + ratio = nw_full * nw_corners^2 / (nw_0 * nw_60 * nw_120) + # Keep the full phase until the unit-cell product, without cubing its magnitude. + return cbrt(abs(ratio)) * sign(ratio) + end + return abs(value) * cis(angle(value) / 3) +end + +function local_expectation_value( + inds, bra::InfinitePEPSTriangular, operator::AbstractTensorMap, + ket::InfinitePEPSTriangular, env::CTMRGEnvTriangular, + ) + site = first(inds) + if length(inds) == 1 + contract = _contract_corners + elseif length(inds) == 2 + offsets = (CartesianIndex(0, 1), CartesianIndex(1, -1), CartesianIndex(1, 0)) + offset = inds[2] - inds[1] + dir = findfirst(==(offset), offsets) + if isnothing(dir) + dir = findfirst(==(-offset), offsets) + isnothing(dir) && throw(ArgumentError("only nearest-neighbor triangular terms are supported")) + site = inds[2] + operator = permute(operator, ((2, 1), (4, 3))) + end + contract = (_contract_edges_0, _contract_edges_60, _contract_edges_120)[dir] + else + throw(ArgumentError("only one-site and two-site triangular terms are supported")) + end + network = InfiniteTriangularNetwork(ket, bra) + site = mod1.(Tuple(site), size(network)) + return contract(site, network, env, operator) / contract(site, network, env) +end + +function energy(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, onesite_op, twosite_op) where {P <: PEPSSandwichTriangular} + E = sum(Iterators.product(axes(network)...)) do (r, c) + expval_onesite = _contract_corners((r, c), network, env, onesite_op) / _contract_corners((r, c), network, env) + expval_twosite = _contract_edges_0((r, c), network, env, twosite_op) / _contract_edges_0((r, c), network, env) + return expval_onesite + expval_twosite + end + ignore_derivatives() do + isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || + @warn "Expectation value is not real: $E." + end + return real(E) +end +function energy(state::InfinitePEPSTriangular, env::CTMRGEnvTriangular, args...) + return energy(InfiniteTriangularNetwork(state), env, args...) +end + +""" + energy(network::InfiniteTriangularNetwork, env, onesite_op, twosite_ops::NTuple{3}) + +Compute the energy per unit cell with a one-site term and nearest-neighbor terms +in the 0°, 60°, and 120° bond families. The oriented bonds connect `(r, c)` to +`(r, c + 1)`, `(r + 1, c - 1)`, and `(r + 1, c)`, respectively. +Divide by `length(network)` for the energy per site. +""" +function energy(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, onesite_op, twosite_ops::NTuple{3, AbstractTensorMap}) where {P <: PEPSSandwichTriangular} + contractions = (_contract_edges_0, _contract_edges_60, _contract_edges_120) + E = sum(Iterators.product(axes(network)...)) do site + expval_onesite = _contract_corners(site, network, env, onesite_op) / _contract_corners(site, network, env) + expval_twosite = sum(zip(contractions, twosite_ops)) do (contract, op) + return contract(site, network, env, op) / contract(site, network, env) + end + return expval_onesite + expval_twosite + end + ignore_derivatives() do + isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || + @warn "Expectation value is not real: $E." + end + return real(E) +end diff --git a/src/environments/ctmrg_environments_triangular.jl b/src/environments/ctmrg_environments_triangular.jl new file mode 100644 index 000000000..454bdd9d9 --- /dev/null +++ b/src/environments/ctmrg_environments_triangular.jl @@ -0,0 +1,197 @@ +""" +$(TYPEDEF) + +Triangular corner transfer-matrix environment, with one corner and two edge halves +(`Ea` and `Eb`, also called Ta and Tb) for each of six directions at every site. The arrays +are indexed as `[direction, row, column]`, with directions ordered `NW, NE, E, SE, SW, W`. + +## Fields + +$(TYPEDFIELDS) + +## Constructors + + CTMRGEnvTriangular(f, T, D, χspace; unitcell=(1, 1)) + +Initialize normalized tensors using `f` and scalar type `T`, with uniform virtual +space `D` and environment space `χspace`. For a PEPS with virtual space `V`, use +`D = V ⊗ V'`; for a classical network, use its single-layer virtual space. +The unit-cell size must match the state or network. Use a complex scalar type for +optimization with [`FixedPointGradient`](@ref). +""" +struct CTMRGEnvTriangular{TC, TA, TB} + "6 x rows x cols array of corner C tensors, where the first dimension specifies the spatial direction" + C::Array{TC, 3} + "6 x rows x cols array of edge Ta tensors, where the first dimension specifies the spatial direction" + Ea::Array{TA, 3} + "6 x rows x cols array of edge Tb tensors, where the first dimension specifies the spatial direction" + Eb::Array{TB, 3} +end + +function get_Ds(D::A) where {A <: ProductSpace} + return [dir > 3 ? _elementwise_dual(D) : D for dir in 1:6] +end + +function get_Ds(D::A) where {A <: ElementarySpace} + return [dir > 3 ? ProductSpace(D') : ProductSpace(D) for dir in 1:6] +end + +function CTMRGEnvTriangular( + f, T, D::Union{A, B}, chis::B; unitcell::Tuple{Int, Int} = (1, 1) + ) where {A <: ProductSpace, B <: ElementarySpace} + Ds = get_Ds(D) + N = length(Ds[1]) + st = spacetype(Ds[1]) + + T_type = tensormaptype(st, N + 1, 1, T) + + Cs = Array{T_type}(undef, 6, unitcell...) + Eas = Array{T_type}(undef, 6, unitcell...) + Ebs = Array{T_type}(undef, 6, unitcell...) + + for dir in 1:6, r in 1:unitcell[1], c in 1:unitcell[2] + C = _edge_tensor(f, T, chis, Ds[dir], chis) + Ea = _edge_tensor(f, T, chis, Ds[dir], chis) + Eb = _edge_tensor(f, T, chis, Ds[mod1(dir + 1, 6)], chis) + + C /= norm(C) + Ea /= norm(Ea) + Eb /= norm(Eb) + + Cs[dir, r, c] = C + Eas[dir, r, c] = Ea + Ebs[dir, r, c] = Eb + end + return CTMRGEnvTriangular(Cs, Eas, Ebs) +end + +Base.size(env::CTMRGEnvTriangular, args...) = size(env.C, args...) +Base.axes(x::CTMRGEnvTriangular, args...) = axes(x.C, args...) +function eachcoordinate(x::CTMRGEnvTriangular) + return collect(Iterators.product(axes(x, 2), axes(x, 3))) +end +function eachcoordinate(x::CTMRGEnvTriangular, dirs) + return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) +end +Base.real(env::CTMRGEnvTriangular) = CTMRGEnvTriangular(real.(env.C), real.(env.Ea), real.(env.Eb)) +Base.complex(env::CTMRGEnvTriangular) = CTMRGEnvTriangular(complex.(env.C), complex.(env.Ea), complex.(env.Eb)) + +cornertype(env::CTMRGEnvTriangular) = cornertype(typeof(env)) +cornertype(::Type{CTMRGEnvTriangular{C, A, B}}) where {C, A, B} = C +edgetype(env::CTMRGEnvTriangular) = edgetype(typeof(env)) +edgetype(::Type{CTMRGEnvTriangular{C, A, B}}) where {C, A, B} = promote_type(A, B) + +TensorKit.spacetype(::Type{E}) where {E <: CTMRGEnvTriangular} = spacetype(cornertype(E)) + +function VI.scalartype(::Type{CTMRGEnvTriangular{C, A, B}}) where {C, A, B} + return promote_type(scalartype(C), scalartype(A), scalartype(B)) +end + +# In-place update of environment +function update!(env::CTMRGEnvTriangular{T}, env´::CTMRGEnvTriangular{T}) where {T} + env.C .= env´.C + env.Ea .= env´.Ea + env.Eb .= env´.Eb + return env +end + +# Custom adjoint for CTMRGEnvTriangular constructor, needed for fixed-point differentiation +function ChainRulesCore.rrule( + ::Type{CTMRGEnvTriangular}, C::Array{T, 3}, Ea::Array{T, 3}, Eb::Array{T, 3} + ) where {T} + function triangularenv_pullback(Δenv_) + Δenv = unthunk(Δenv_) + return NoTangent(), Δenv.C, Δenv.Ea, Δenv.Eb + end + return CTMRGEnvTriangular(C, Ea, Eb), triangularenv_pullback +end + +# math basics used for tangent accumulation in backwards pass +function Base.:+(e₁::CTMRGEnvTriangular, e₂::CTMRGEnvTriangular) + return CTMRGEnvTriangular(e₁.C + e₂.C, e₁.Ea + e₂.Ea, e₁.Eb + e₂.Eb) +end +function Base.:-(e₁::CTMRGEnvTriangular, e₂::CTMRGEnvTriangular) + return CTMRGEnvTriangular(e₁.C - e₂.C, e₁.Ea - e₂.Ea, e₁.Eb - e₂.Eb) +end +Base.:*(α::Number, e::CTMRGEnvTriangular) = CTMRGEnvTriangular(α * e.C, α * e.Ea, α * e.Eb) + +# Custom adjoint for CTMRGEnvTriangular getproperty, to avoid creating named tuples in backward pass +function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnvTriangular, name::Symbol) + result = getproperty(e, name) + if name === :C + function C_pullback(ΔC_) + ΔC = unthunk(ΔC_) + return NoTangent(), CTMRGEnvTriangular(ΔC, zerovector(e.Ea), zerovector(e.Eb)), NoTangent() + end + return result, C_pullback + elseif name === :Ea + function Ea_pullback(ΔEa_) + ΔEa = unthunk(ΔEa_) + return NoTangent(), CTMRGEnvTriangular(zerovector(e.C), ΔEa, zerovector(e.Eb)), NoTangent() + end + return result, Ea_pullback + elseif name === :Eb + function Eb_pullback(ΔEb_) + ΔEb = unthunk(ΔEb_) + return NoTangent(), CTMRGEnvTriangular(zerovector(e.C), zerovector(e.Ea), ΔEb), NoTangent() + end + return result, Eb_pullback + else + # this should never happen because already errored in forwards pass + throw(ArgumentError("No rrule for getproperty of $name")) + end +end + +function VI.zerovector(env::CTMRGEnvTriangular, ::Type{S}) where {S <: Number} + _zerovector = Base.Fix2(zerovector, S) + return CTMRGEnvTriangular(map(_zerovector, env.C), map(_zerovector, env.Ea), map(_zerovector, env.Eb)) +end +function VI.zerovector!(env::CTMRGEnvTriangular) + foreach(zerovector!, env.C) + foreach(zerovector!, env.Ea) + foreach(zerovector!, env.Eb) + return env +end +VI.zerovector!!(env::CTMRGEnvTriangular) = zerovector!(env) + +function VI.scale(env::CTMRGEnvTriangular, α::Number) + _scale = Base.Fix2(scale, α) + return CTMRGEnvTriangular(map(_scale, env.C), map(_scale, env.Ea), map(_scale, env.Eb)) +end +function VI.scale!(env::CTMRGEnvTriangular, α::Number) + _scale! = Base.Fix2(scale!, α) + foreach(_scale!, env.C) + foreach(_scale!, env.Ea) + foreach(_scale!, env.Eb) + return env +end +function VI.scale!(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular, α::Number) + _scale!(x, y) = scale!(x, y, α) + foreach(_scale!, env₁.C, env₂.C) + foreach(_scale!, env₁.Ea, env₂.Ea) + foreach(_scale!, env₁.Eb, env₂.Eb) + return env₁ +end +VI.scale!!(env::CTMRGEnvTriangular, α::Number) = scale!(env, α) +VI.scale!!(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular, α::Number) = scale!(env₁, env₂, α) + +function VI.add(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular, α::Number, β::Number) + _add(x, y) = add(x, y, α, β) + return CTMRGEnvTriangular( + map(_add, env₁.C, env₂.C), map(_add, env₁.Ea, env₂.Ea), map(_add, env₁.Eb, env₂.Eb) + ) +end +function VI.add!(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular, α::Number, β::Number) + _add!(x, y) = add!(x, y, α, β) + foreach(_add!, env₁.C, env₂.C) + foreach(_add!, env₁.Ea, env₂.Ea) + foreach(_add!, env₁.Eb, env₂.Eb) + return env₁ +end +VI.add!!(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular, α::Number, β::Number) = add!(env₁, env₂, α, β) + +# Exploiting the fact that VectorInterface works for tuples: +function VI.inner(env₁::CTMRGEnvTriangular, env₂::CTMRGEnvTriangular) + return inner((env₁.C, env₁.Ea, env₁.Eb), (env₂.C, env₂.Ea, env₂.Eb)) +end +VI.norm(env::CTMRGEnvTriangular) = norm((env.C, env.Ea, env.Eb)) diff --git a/src/networks/infinitetriangularnetwork.jl b/src/networks/infinitetriangularnetwork.jl new file mode 100644 index 000000000..ac3071c9f --- /dev/null +++ b/src/networks/infinitetriangularnetwork.jl @@ -0,0 +1,155 @@ +const PEPSTensorTriangular{S <: ElementarySpace} = AbstractTensorMap{<:Any, S, 1, 6} + +""" +$(TYPEDEF) + +Contractible triangular network. Wraps a matrix of 'rank-6-tensor-like' objects. + +## Fields + +$(TYPEDFIELDS) +""" +struct InfiniteTriangularNetwork{O} + A::Matrix{O} + InfiniteTriangularNetwork{O}(A::Matrix{O}) where {O} = new{O}(A) + function InfiniteTriangularNetwork(A::Matrix) + for I in eachindex(IndexCartesian(), A) + r, c = Tuple(I) + virtualspace(A[r, c], 1) == + _elementwise_dual(virtualspace(A[_prev(r, end), c], 4)) || + throw(SpaceMismatch("Northwest virtual space at site $((r, c)) does not match.")) + virtualspace(A[r, c], 2) == + _elementwise_dual(virtualspace(A[_prev(r, end), _next(c, end)], 5)) || + throw(SpaceMismatch("Northeast virtual space at site $((r, c)) does not match.")) + virtualspace(A[r, c], 3) == + _elementwise_dual(virtualspace(A[r, _next(c, end)], 6)) || + throw(SpaceMismatch("East virtual space at site $((r, c)) does not match.")) + end + return InfiniteTriangularNetwork{eltype(A)}(A) + end +end +InfiniteTriangularNetwork(n::InfiniteTriangularNetwork) = n + +## Unit cell interface + +unitcell(n::InfiniteTriangularNetwork) = n.A +Base.size(n::InfiniteTriangularNetwork, args...) = size(unitcell(n), args...) +Base.length(n::InfiniteTriangularNetwork) = length(unitcell(n)) +Base.eltype(n::InfiniteTriangularNetwork) = eltype(typeof(n)) +Base.eltype(::Type{InfiniteTriangularNetwork{O}}) where {O} = O + +Base.copy(n::InfiniteTriangularNetwork) = InfiniteTriangularNetwork(copy(unitcell(n))) +function Base.similar(n::InfiniteTriangularNetwork, T::Type{TorA} = scalartype(n)) where {TorA} + return InfiniteTriangularNetwork(map(t -> similar(t, T), unitcell(n))) +end +function Base.repeat(n::InfiniteTriangularNetwork, counts...) + return InfiniteTriangularNetwork(repeat(unitcell(n), counts...)) +end + +## Indexing +Base.getindex(n::InfiniteTriangularNetwork, args...) = Base.getindex(unitcell(n), args...) +function Base.setindex!(n::InfiniteTriangularNetwork, args...) + return (Base.setindex!(unitcell(n), args...); n) +end +Base.axes(n::InfiniteTriangularNetwork, args...) = axes(unitcell(n), args...) +eachcoordinate(n::InfiniteTriangularNetwork) = collect(Iterators.product(axes(n)...)) +function eachcoordinate(n::InfiniteTriangularNetwork, dirs) + return collect(Iterators.product(dirs, axes(n, 1), axes(n, 2))) +end + +## Spaces + +TensorKit.spacetype(::Type{T}) where {T <: InfiniteTriangularNetwork} = spacetype(eltype(T)) +function virtualspace(n::InfiniteTriangularNetwork, r::Int, c::Int, dir) + Nr, Nc = size(n) + return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir) +end + +## Vector interface + +function VectorInterface.scalartype(::Type{T}) where {T <: InfiniteTriangularNetwork} + return scalartype(eltype(T)) +end +function VectorInterface.zerovector(A::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(zerovector(unitcell(A))) +end + +## Math (for Zygote accumulation) + +function Base.:+(A₁::InfiniteTriangularNetwork, A₂::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(_add_localsandwich.(unitcell(A₁), unitcell(A₂))) +end +function Base.:-(A₁::InfiniteTriangularNetwork, A₂::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(_subtract_localsandwich.(unitcell(A₁), unitcell(A₂))) +end +function Base.:*(α::Number, A::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(_mul_localsandwich.(α, unitcell(A))) +end +Base.:*(A::InfiniteTriangularNetwork, α::Number) = α * A +function Base.:/(A::InfiniteTriangularNetwork, α::Number) + return A * inv(α) +end +function LinearAlgebra.dot(A₁::InfiniteTriangularNetwork, A₂::InfiniteTriangularNetwork) + return dot(unitcell(A₁), unitcell(A₂)) +end +LinearAlgebra.norm(A::InfiniteTriangularNetwork) = norm(unitcell(A)) + +## (Approximate) equality + +function Base.:(==)(A₁::InfiniteTriangularNetwork, A₂::InfiniteTriangularNetwork) + return all(zip(unitcell(A₁), unitcell(A₂))) do (p₁, p₂) + return p₁ == p₂ + end +end +function Base.isapprox(A₁::InfiniteTriangularNetwork, A₂::InfiniteTriangularNetwork; kwargs...) + return all(zip(unitcell(A₁), unitcell(A₂))) do (p₁, p₂) + return _isapprox_localsandwich(p₁, p₂; kwargs...) + end +end + +## Rotations + +function rotl60(n::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(rotl60(_rotl60_localsandwich.(unitcell(n)))) +end +function rotr60(n::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(rotr60(_rotr60_localsandwich.(unitcell(n)))) +end +function Base.rot180(n::InfiniteTriangularNetwork) + return InfiniteTriangularNetwork(rot180(_rot180_localsandwich.(unitcell(n)))) +end + +## Chainrules + +# generic implementation +function ChainRulesCore.rrule( + ::typeof(Base.getindex), network::InfiniteTriangularNetwork, r::Int, c::Int + ) + O = network[r, c] + + function getindex_pullback(ΔO_) + ΔO = map(unthunk, ΔO_) + if ΔO isa Tangent + ΔO = ChainRulesCore.construct(typeof(O), ChainRulesCore.backing(ΔO)) + end + Δnetwork = zerovector(network) + Δnetwork[r, c] = ΔO + return NoTangent(), Δnetwork, NoTangent(), NoTangent() + end + return O, getindex_pullback +end + +# specialized PFTensorTriangular implementation +function ChainRulesCore.rrule( + ::typeof(Base.getindex), network::InfiniteTriangularNetwork{<:PFTensorTriangular}, r::Int, c::Int + ) + O = network[r, c] + + function getindex_pullback(ΔO_) + ΔO = unthunk(ΔO_) + Δnetwork = zerovector(network) + Δnetwork[r, c] = ΔO + return NoTangent(), Δnetwork, NoTangent(), NoTangent() + end + return O, getindex_pullback +end diff --git a/src/networks/local_triangular_sandwich.jl b/src/networks/local_triangular_sandwich.jl new file mode 100644 index 000000000..a5fe0ab99 --- /dev/null +++ b/src/networks/local_triangular_sandwich.jl @@ -0,0 +1,43 @@ +# +# Interface for local effective-rank-6 tensor sandwiches +# + +# route all virtualspace getters through a single method for convenience +northwest_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[1] +northeast_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[2] +east_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[3] +southeast_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[4] +southwest_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[5] +west_virtualspace(O::AbstractTensorMap{E, S, 1, 6}, args...) where {E, S} = domain(O)[6] + +## Rotations + +# generic local interface +_rotl60_localsandwich(O) = rotl60.(O) +_rotr60_localsandwich(O) = rotr60.(O) + +## PartitionFunction + +# specialized local rotation interface +_rotl60_localsandwich(O::PFTensorTriangular) = rotl60(O) +_rotr60_localsandwich(O::PFTensorTriangular) = rotr60(O) +_rot180_localsandwich(O::PFTensorTriangular) = rot180(O) + +# specialized local math interface +_add_localsandwich(O1::PFTensorTriangular, O2::PFTensorTriangular) = O1 + O2 +_subtract_localsandwich(O1::PFTensorTriangular, O2::PFTensorTriangular) = O1 - O2 +_mul_localsandwich(α::Number, O::PFTensorTriangular) = α * O +_isapprox_localsandwich(O1::PFTensorTriangular, O2::PFTensorTriangular; kwargs...) = isapprox(O1, O2; kwargs...) + +## PEPS + +const PEPSSandwichTriangular{T <: PEPSTensorTriangular} = Tuple{T, T} + +ket(O::PEPSSandwichTriangular) = O[1] +bra(O::PEPSSandwichTriangular) = O[2] + +function virtualspace(O::PEPSSandwichTriangular, dir) + return virtualspace(ket(O), dir) ⊗ virtualspace(bra(O), dir)' +end + +TensorKit.spacetype(::Type{P}) where {P <: PEPSSandwichTriangular} = spacetype(eltype(P)) diff --git a/src/networks/tensors_triangular.jl b/src/networks/tensors_triangular.jl new file mode 100644 index 000000000..9366feb7c --- /dev/null +++ b/src/networks/tensors_triangular.jl @@ -0,0 +1,101 @@ +# +# Partition function +# + +""" + const PartitionFunctionTensorTriangular{S} + +Default type for partition function tensors with 6 virtual indices, conventionally ordered +as: ``T : W ⊗ SW ⊗ SE ← NW ⊗ NE ⊗ E``. Here ``NW``, ``NE``, ``E``, ``SE``, ``SW`` and ``W`` denote the +north-west (120°), north-east (60°), east (0°), south-east (300°), south-west (240°) and west (180°) spaces, respectively, +where the angles denote the directions of the legs with respect to the positive horizontal axis. + +``` + NW NE + ╲ ╱ + ╲ ╱ + ╲ ╱ + W ----- T ----- E + ╱ ╲ + ╱ ╲ + ╱ ╲ + SW SE +``` +""" +const PartitionFunctionTensorTriangular{S <: ElementarySpace} = AbstractTensorMap{<:Any, S, 3, 3} +const PFTensorTriangular = PartitionFunctionTensorTriangular + +""" + PartitionFunctionTensorTriangular(f, ::Type{T}, NWspace::S, + [NEspace::S], [Espace::S], [SEspace::S], + [SWspace::S], [Wspace::S]) where {T,S<:ElementarySpace} + +Construct a `PartitionFunctionTensorTriangular` tensor based on the +north-west, north-east, east, south-east, south-west and west spaces. +The tensor elements are generated based on `f` and the element type is specified in `T`. +""" +function PartitionFunctionTensorTriangular( + f, ::Type{T}, + NWspace::S, NEspace::S = NWspace, Espace::S = NWspace, SEspace::S = NWspace, + SWspace::S = NWspace, Wspace::S = NWspace, + ) where {T, S <: ElementarySpace} + return f(T, Wspace ⊗ SWspace ⊗ SEspace ← NWspace ⊗ NEspace ⊗ Espace) +end + +Base.rot180(t::PFTensorTriangular) = permute(t, ((6, 5, 4), (3, 2, 1))) + +function virtualspace(t::PFTensorTriangular, dir) + invp = (4, 5, 6, 3, 2, 1) # internally, virtual directions are ordered as NW, NE, E, SE, SW... + return space(t, invp[dir]) +end + +# +# PEPS +# + +""" + const PEPSTensorTriangular{S} + +Default type for PEPS tensors with a single physical index, and 6 virtual indices, +conventionally ordered as: ``T : P ← NW ⊗ NE ⊗ E ⊗ SE ⊗ SW ⊗ W``. Here, ``P`` denotes the physical space +and ``NW``, ``NE``, ``E``, ``SE``, ``SW`` and ``W`` denote the north-west (120°), north-east (60°), +east (0°), south-east (300°), south-west (240°) and west (180°) spaces, respectively, +where the angles denote the directions of the legs with respect to the positive horizontal axis. + +``` + NW NE + ╲ ╱ + ╲ ╱ + ╲ ╱ + W ----- T ----- E + ╱|╲ + ╱ | ╲ + ╱ | ╲ + SW P SE +``` +""" +const PEPSTensorTriangular{S <: ElementarySpace} = AbstractTensorMap{<:Any, S, 1, 6} + +""" + PEPSTensorTriangular(f, ::Type{T}, Pspace::S, NWspace::S, + [NEspace::S], [Espace::S], [SEspace::S], + [SWspace::S], [Wspace::S]) where {T,S<:ElementarySpace} + +Construct a PEPS tensor based on the physical, north-west, north-east, east, south-east, south-west and west spaces. +The tensor elements are generated based on `f` and the element type is specified in `T`. +""" +function PEPSTensorTriangular( + f, ::Type{T}, + Pspace::S, + NWspace::S, NEspace::S = NWspace, Espace::S = NWspace, SEspace::S = NWspace, + SWspace::S = NWspace, Wspace::S = NWspace, + ) where {T, S <: ElementarySpace} + return f(T, Pspace ← NWspace ⊗ NEspace ⊗ Espace ⊗ SEspace ⊗ SWspace ⊗ Wspace) +end + +rotl60(t::PEPSTensorTriangular) = permute(t, ((1,), (3, 4, 5, 6, 7, 2))) +rotr60(t::PEPSTensorTriangular) = permute(t, ((1,), (7, 2, 3, 4, 5, 6))) +Base.rot180(t::PEPSTensorTriangular) = permute(t, ((1,), (5, 6, 7, 2, 3, 4))) + +physicalspace(t::PEPSTensorTriangular) = space(t, 1) +virtualspace(t::PEPSTensorTriangular, dir) = space(t, dir + 1) diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 137cc0441..67d1ae254 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -89,10 +89,10 @@ end function checklattice(::Type{Bool}, H1::LocalOperator, H2::LocalOperator) return physicalspace(H1) == physicalspace(H2) end -function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalOperator) +function checklattice(::Type{Bool}, peps::Union{InfinitePEPS, InfinitePEPSTriangular}, O::LocalOperator) return physicalspace(peps) == physicalspace(O) end -function checklattice(::Type{Bool}, H::LocalOperator, peps::InfinitePEPS) +function checklattice(::Type{Bool}, H::LocalOperator, peps::Union{InfinitePEPS, InfinitePEPSTriangular}) return checklattice(Bool, peps, H) end function checklattice(::Type{Bool}, pepo::InfinitePEPO, O::LocalOperator) diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl index fd94825e1..b1162fb6a 100644 --- a/src/operators/transfermatrix.jl +++ b/src/operators/transfermatrix.jl @@ -153,10 +153,10 @@ function initialize_mps( return InfiniteMPS( [ f( - T, - virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)), - virtualspaces[mod1(i, end)], - ) for i in 1:length(O) + T, + virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)), + virtualspaces[mod1(i, end)], + ) for i in 1:length(O) ]; kwargs... ) end diff --git a/src/states/infinitepepstriangular.jl b/src/states/infinitepepstriangular.jl new file mode 100644 index 000000000..1ace2ae81 --- /dev/null +++ b/src/states/infinitepepstriangular.jl @@ -0,0 +1,241 @@ +""" + struct InfinitePEPSTriangular{T<:PEPSTensorTriangular} + +Represents an infinite projected entangled-pair state on a 2D triangular lattice. + +## Fields + +$(TYPEDFIELDS) +""" +struct InfinitePEPSTriangular{T <: PEPSTensorTriangular} + A::Matrix{T} + InfinitePEPSTriangular{T}(A::Matrix{T}) where {T <: PEPSTensorTriangular} = new{T}(A) + function InfinitePEPSTriangular(A::Array{T, 2}) where {T <: PEPSTensorTriangular} + bosonic_braiding = BraidingStyle(sectortype(T)) === Bosonic() + for (r, c) in Tuple.(CartesianIndices(A)) + (bosonic_braiding || !isdual(physicalspace(A[r, c]))) || + throw(ArgumentError("Dual physical spaces for symmetry sectors with non-trivial twists are not allowed (for now).")) + northwest_virtualspace(A[r, c]) == southeast_virtualspace(A[_prev(r, end), c])' || + throw(SpaceMismatch("Northwest virtual space at site $((r, c)) does not match.")) + northeast_virtualspace(A[r, c]) == southwest_virtualspace(A[_prev(r, end), _next(c, end)])' || + throw(SpaceMismatch("Northeast virtual space at site $((r, c)) does not match.")) + east_virtualspace(A[r, c]) == west_virtualspace(A[r, _next(c, end)])' || + throw(SpaceMismatch("East virtual space at site $((r, c)) does not match.")) + dim(space(A[r, c])) > 0 || @warn "no fusion channels at site ($r, $c)" + end + return new{T}(A) + end +end + +## Constructors + +""" + InfinitePEPSTriangular(A::AbstractMatrix{T}) + +Create an `InfinitePEPSTriangular` by specifying a matrix containing the PEPS tensors at each site in +the unit cell. +""" +function InfinitePEPSTriangular(A::AbstractMatrix{<:PEPSTensorTriangular}) + return InfinitePEPSTriangular(Array(deepcopy(A))) # TODO: find better way to copy +end + +""" + InfinitePEPSTriangular([f=randn, T=ComplexF64,] Pspaces::A, NWspaces::A, [NEspaces::A], [Espaces::A]) where {A<:AbstractMatrix{ElementarySpace}} + +Create an `InfinitePEPSTriangular` by specifying matrices of physical, northwest, northeast +and east virtual spaces. Omitted virtual spaces default to `NWspaces`; the opposite +spaces are the duals of the corresponding neighboring spaces, with periodic boundaries. +""" +function InfinitePEPSTriangular( + f, T::Type{<:Number}, Pspaces::M, NWspaces::M, NEspaces::M = NWspaces, Espaces::M = NWspaces + ) where {M <: AbstractMatrix{<:ElementarySpace}} + size(Pspaces) == size(NWspaces) == size(NEspaces) == size(Espaces) || + throw(ArgumentError("Input spaces should have equal sizes.")) + + SEspaces = adjoint.(circshift(NWspaces, (-1, 0))) + SWspaces = adjoint.(circshift(NEspaces, (-1, 1))) + Wspaces = adjoint.(circshift(Espaces, (0, 1))) + + A = map(Pspaces, NWspaces, NEspaces, Espaces, SEspaces, SWspaces, Wspaces) do P, NW, NE, E, SE, SW, W + return PEPSTensorTriangular(f, T, P, NW, NE, E, SE, SW, W) + end + + return InfinitePEPSTriangular(A) +end +function InfinitePEPSTriangular( + Pspaces::A, virtual_spaces...; kwargs... + ) where {A <: Union{AbstractMatrix{<:ElementarySpace}, ElementarySpace}} + return InfinitePEPSTriangular(randn, ComplexF64, Pspaces, virtual_spaces...; kwargs...) +end + +""" + InfinitePEPSTriangular(A::PEPSTensorTriangular; unitcell=(1, 1)) + +Create an `InfinitePEPSTriangular` by specifying a tensor and unit cell. + +The unit cell uses one-based matrix coordinates `[r, c]`. The virtual legs are ordered +`NW, NE, E, SE, SW, W`, connecting to neighbors at offsets +`(-1, 0), (-1, 1), (0, 1), (1, 0), (1, -1), (0, -1)`, respectively. +Contractions wrap these coordinates periodically; direct `peps[r, c]` indexing follows +ordinary matrix bounds. A three-sublattice cell is `[A B C; C A B; B C A]`. +""" +function InfinitePEPSTriangular(A::T; unitcell::Tuple{Int, Int} = (1, 1)) where {T <: PEPSTensorTriangular} + return InfinitePEPSTriangular(fill(A, unitcell)) +end + +# expand PEPS spaces to unit cell size +function _fill_state_virtual_spaces_triangular( + NWspace::S, NEspace::S = NWspace, Espace::S = NWspace; unitcell::Tuple{Int, Int} = (1, 1) + ) where {S <: ElementarySpace} + return fill(NWspace, unitcell), fill(NEspace, unitcell), fill(Espace, unitcell) +end + +""" + InfinitePEPSTriangular([f=randn, T=ComplexF64,] Pspace, NWspace, [NEspace], [Espace]; unitcell=(1,1)) + +Create an `InfinitePEPSTriangular` with uniform physical, northwest, northeast and east +virtual spaces and the given unit-cell size. Omitted virtual spaces default to `NWspace`. +""" +function InfinitePEPSTriangular( + f, T::Type{<:Number}, Pspace::S, vspaces...; unitcell::Tuple{Int, Int} = (1, 1) + ) where {S <: ElementarySpace} + return InfinitePEPSTriangular( + f, T, + _fill_state_physical_spaces(Pspace; unitcell), + _fill_state_virtual_spaces_triangular(vspaces...; unitcell)..., + ) +end + +## Unit cell interface + +unitcell(t::InfinitePEPSTriangular) = t.A +Base.size(A::InfinitePEPSTriangular, args...) = size(unitcell(A), args...) +Base.length(A::InfinitePEPSTriangular) = length(unitcell(A)) +Base.eltype(::Type{InfinitePEPSTriangular{T}}) where {T} = T +Base.eltype(A::InfinitePEPSTriangular) = eltype(typeof(A)) + +Base.copy(A::InfinitePEPSTriangular) = InfinitePEPSTriangular(copy(unitcell(A))) +function Base.similar(A::InfinitePEPSTriangular, T::Type{TorA} = scalartype(A)) where {TorA} + return InfinitePEPSTriangular(map(t -> similar(t, T), unitcell(A))) +end +Base.repeat(A::InfinitePEPSTriangular, counts...) = InfinitePEPSTriangular(repeat(unitcell(A), counts...)) + +Base.getindex(A::InfinitePEPSTriangular, args...) = Base.getindex(unitcell(A), args...) +Base.setindex!(A::InfinitePEPSTriangular, args...) = (Base.setindex!(unitcell(A), args...); A) +Base.axes(A::InfinitePEPSTriangular, args...) = axes(unitcell(A), args...) +eachcoordinate(A::InfinitePEPSTriangular) = collect(Iterators.product(axes(A)...)) +function eachcoordinate(A::InfinitePEPSTriangular, dirs) + return collect(Iterators.product(dirs, axes(A, 1), axes(A, 2))) +end + +## Spaces + +TensorKit.spacetype(::Type{T}) where {T <: InfinitePEPSTriangular} = spacetype(eltype(T)) +virtualspace(n::InfinitePEPSTriangular, dir) = virtualspace.(unitcell(n), dir) +function virtualspace(n::InfinitePEPSTriangular, r::Int, c::Int, dir) + Nr, Nc = size(n) + return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir) +end +physicalspace(n::InfinitePEPSTriangular) = physicalspace.(unitcell(n)) +function physicalspace(n::InfinitePEPSTriangular, r::Int, c::Int) + Nr, Nc = size(n) + return physicalspace(n[mod1(r, Nr), mod1(c, Nc)]) +end + +## InfiniteTriangularNetwork interface + +function InfiniteTriangularNetwork(top::InfinitePEPSTriangular, bot::InfinitePEPSTriangular = top) + size(top) == size(bot) || throw( + ArgumentError("Top PEPS, bottom PEPS and PEPO rows should have the same length") + ) + return InfiniteTriangularNetwork(map(tuple, unitcell(top), unitcell(bot))) +end + +## Vector interface + +VI.scalartype(::Type{NT}) where {NT <: InfinitePEPSTriangular} = scalartype(eltype(NT)) +VI.zerovector(A::InfinitePEPSTriangular) = InfinitePEPSTriangular(zerovector(unitcell(A))) + +function VI.scale(ψ::InfinitePEPSTriangular, α::Number) + _scale = Base.Fix2(scale, α) + return InfinitePEPSTriangular(map(_scale, unitcell(ψ))) +end +function VI.scale!(ψ::InfinitePEPSTriangular, α::Number) + _scale! = Base.Fix2(scale!, α) + foreach(_scale!, unitcell(ψ)) + return ψ +end +function VI.scale!(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular, α::Number) + _scale!(x, y) = scale!(x, y, α) + foreach(_scale!, unitcell(ψ₁), unitcell(ψ₂)) + return ψ₁ +end +VI.scale!!(ψ::InfinitePEPSTriangular, α::Number) = scale!(ψ, α) +VI.scale!!(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular, α::Number) = scale!(ψ₁, ψ₂, α) + +function VI.add(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular, α::Number, β::Number) + _add(x, y) = add(x, y, α, β) + return InfinitePEPSTriangular(map(_add, unitcell(ψ₁), unitcell(ψ₂))) +end +function VI.add!(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular, α::Number, β::Number) + _add!(x, y) = add!(x, y, α, β) + foreach(_add!, unitcell(ψ₁), unitcell(ψ₂)) + return ψ₁ +end +VI.add!!(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular, α::Number, β::Number) = add!(ψ₁, ψ₂, α, β) +VI.inner(ψ₁::InfinitePEPSTriangular, ψ₂::InfinitePEPSTriangular) = inner(unitcell(ψ₁), unitcell(ψ₂)) + +## Math + +function Base.:+(A₁::InfinitePEPSTriangular, A₂::InfinitePEPSTriangular) + return InfinitePEPSTriangular(unitcell(A₁) + unitcell(A₂)) +end +function Base.:-(A₁::InfinitePEPSTriangular, A₂::InfinitePEPSTriangular) + return InfinitePEPSTriangular(unitcell(A₁) - unitcell(A₂)) +end +Base.:*(α::Number, A::InfinitePEPSTriangular) = InfinitePEPSTriangular(α * unitcell(A)) +Base.:*(A::InfinitePEPSTriangular, α::Number) = α * A +Base.:/(A::InfinitePEPSTriangular, α::Number) = InfinitePEPSTriangular(unitcell(A) / α) +LinearAlgebra.dot(A₁::InfinitePEPSTriangular, A₂::InfinitePEPSTriangular) = inner(A₁, A₂) +LinearAlgebra.norm(A::InfinitePEPSTriangular) = norm(unitcell(A)) + +## (Approximate) equality +function Base.:(==)(A₁::InfinitePEPSTriangular, A₂::InfinitePEPSTriangular) + return all(zip(unitcell(A₁), unitcell(A₂))) do (p₁, p₂) + return p₁ == p₂ + end +end +function Base.isapprox(A₁::InfinitePEPSTriangular, A₂::InfinitePEPSTriangular; kwargs...) + return all(zip(unitcell(A₁), unitcell(A₂))) do (p₁, p₂) + return isapprox(p₁, p₂; kwargs...) + end +end + + +## Chainrules + +function ChainRulesCore.rrule(::typeof(Base.getindex), network::InfinitePEPSTriangular, args...) + tensor = network[args...] + + function getindex_pullback(Δtensor_) + Δtensor = unthunk(Δtensor_) + Δnetwork = zerovector(network) + Δnetwork[args...] = Δtensor + return NoTangent(), Δnetwork, NoTangent(), NoTangent() + end + return tensor, getindex_pullback +end + +function ChainRulesCore.rrule( + ::Type{InfiniteTriangularNetwork}, top::InfinitePEPSTriangular, bot::InfinitePEPSTriangular + ) + network = InfiniteTriangularNetwork(top, bot) + + function InfiniteTriangularNetwork_pullback(Δnetwork_) + Δnetwork = unthunk(Δnetwork_) + Δtop = InfinitePEPSTriangular(map(ket, unitcell(Δnetwork))) + Δbot = InfinitePEPSTriangular(map(bra, unitcell(Δnetwork))) + return NoTangent(), Δtop, Δbot + end + return network, InfiniteTriangularNetwork_pullback +end diff --git a/src/utility/svd.jl b/src/utility/svd.jl index 8a39ebf1d..8d42d9d34 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -51,6 +51,10 @@ struct SVDAdjoint{F, R} rrule_alg::R end +function _set_decomposition_truncation(alg::SVDAdjoint, trunc::TruncationStrategy) + return @set alg.fwd_alg = TruncatedAlgorithm(alg.fwd_alg, trunc) +end + const SVD_FWD_SYMBOLS = IdDict{Symbol, Any}( :DefaultAlgorithm => DefaultAlgorithm, :DivideAndConquer => DivideAndConquer, @@ -116,6 +120,25 @@ function SVDAdjoint(; fwd_alg = (;), rrule_alg = (;)) return SVDAdjoint(fwd_algorithm, rrule_algorithm) end +""" + svd_full(t, alg::SVDAdjoint) + svd_full!(t, alg::SVDAdjoint) + +Full SVD with a dense forward algorithm and the gauge verbosity and degeneracy tolerance +of `alg.rrule_alg`. Both `FullPullback` and `TruncPullback` use the full pullback here, +since no singular values are discarded. +""" +MatrixAlgebraKit.svd_full(t, alg::SVDAdjoint) = svd_full!(copy(t), alg) +function MatrixAlgebraKit.svd_full!(t, alg::SVDAdjoint) + alg.fwd_alg isa MatrixAlgebraKit.Algorithm && alg.rrule_alg isa Union{FullPullback, TruncPullback} || + throw(ArgumentError("full SVD requires a dense forward algorithm and FullPullback or TruncPullback")) + return svd_full!(t, alg.fwd_alg) +end +function MatrixAlgebraKit.svd_full!(t::AdjointTensorMap, alg::SVDAdjoint) + u, s, vt = svd_full!(adjoint(t), alg) + return adjoint(vt), adjoint(s), adjoint(u) +end + """ svd_trunc(t, alg::SVDAdjoint) svd_trunc!(t, alg::SVDAdjoint) @@ -285,6 +308,26 @@ end ## Reverse-rule algorithms # +# svd_full! rrule wrapping MatrixAlgebraKit's svd_pullback! +function ChainRulesCore.rrule( + ::typeof(svd_full!), t::AbstractTensorMap, alg::SVDAdjoint{F, R} + ) where {F <: MatrixAlgebraKit.Algorithm, R <: Union{FullPullback, TruncPullback}} + USV = svd_full!(t, alg.fwd_alg) + gtol = _get_pullback_gauge_tol(alg.rrule_alg.verbosity) + function svd_full!_pullback(ΔUSV′) + ΔUSV = unthunk.(ΔUSV′) + Δt = svd_pullback!( + zeros(storagetype(t), space(t)), t, USV, ΔUSV; + gauge_atol = gtol(ΔUSV), degeneracy_atol = alg.rrule_alg.degeneracy_atol, + ) + return NoTangent(), Δt, NoTangent() + end + function svd_full!_pullback(::Tuple{ZeroTangent, ZeroTangent, ZeroTangent}) + return NoTangent(), ZeroTangent(), NoTangent() + end + return USV, svd_full!_pullback +end + # svd_trunc! rrule wrapping MatrixAlgebraKit's svd_pullback! # https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/blob/b76c7bb60014ecfead6925d0df6cb4b8d7c2668a/src/pullbacks/svd.jl#L33 function ChainRulesCore.rrule( diff --git a/src/utility/util.jl b/src/utility/util.jl index df25d2ae5..0da7bd4c5 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -23,6 +23,7 @@ function _maxabs(data::AbstractArray) return norm(data, Inf) end _maxabs(s::AbstractTensorMap) = _maxabs(s.data) +_maxabs(s::AdjointTensorMap) = _maxabs(adjoint(s)) """ sdiag_pow(s, pow::Real; tol::Real=eps(real(scalartype(s)))^(3 / 4)) diff --git a/test/ctmrg/triangular.jl b/test/ctmrg/triangular.jl new file mode 100644 index 000000000..3fa413a6b --- /dev/null +++ b/test/ctmrg/triangular.jl @@ -0,0 +1,742 @@ +using Test +using Random +using StableRNGs +using TensorKit +using MPSKit +using PEPSKit +using OptimKit +using Zygote +using LinearAlgebra +using KrylovKit + +using PEPSKit: unitcell, norm_preserving_retract + +sd = 123456 +Random.seed!(sd) + +D = 2 +χ = 7 +T = ComplexF64 + +ctmrg_tol = 1.0e-8 +ctmrg_maxiter = 300 +ctmrg_verbosity = 2 + +const ising_βc_triangular = BigFloat(BigFloat(asinh(BigFloat(sqrt(BigFloat(1.0) / BigFloat(3.0))))) / BigFloat(2.0)) +const f_onsager_triangular::BigFloat = -3.20253248660790791834355252025862951439 + +function classical_ising_triangular(β) + t = Float64[exp(β) exp(-β); exp(-β) exp(β)] + + r = eigen(t) + nt = r.vectors * sqrt(LinearAlgebra.Diagonal(r.values)) * r.vectors + + O = zeros(2, 2, 2, 2, 2, 2) + O[1, 1, 1, 1, 1, 1] = 1 + O[2, 2, 2, 2, 2, 2] = 1 + + H = [1 1; 1 -1] / sqrt(2) + + @tensor o[-1 -2 -3; -4 -5 -6] := O[1 2 3; 4 5 6] * nt[-1; 1] * nt[-2; 2] * nt[-3; 3] * nt[-4; 4] * nt[-5; 5] * nt[-6; 6] + @tensor o2[-1 -2 -3; -4 -5 -6] := o[1 2 3; 4 5 6] * H[-1; 1] * H[-2; 2] * H[-3; 3] * H[-4; 4] * H[-5; 5] * H[-6; 6] + return TensorMap(o2, ℂ^2 * ℂ^2 * ℂ^2, ℂ^2 * ℂ^2 * ℂ^2) +end + +@testset "CTM_triangular - Coordinates" begin + neighbors = ((3, 1), (3, 2), (1, 2), (2, 1), (2, 4), (1, 4)) + for dir in 1:6 + @test PEPSKit._coordinates(dir, 1, 1, (3, 4)) == neighbors[dir] + end +end + +@testset "CTM_triangular - Network value phases" begin + V = ComplexSpace(1) + for sandwich in (false, true) + network = if sandwich + A = normalize(ones(Float64, ComplexSpace(2), V^3 ⊗ (V')^3)) + InfiniteTriangularNetwork(InfinitePEPSTriangular(A; unitcell = (2, 2))) + else + InfiniteTriangularNetwork(fill(ones(Float64, V^3, V^3), 2, 2)) + end + env = CTMRGEnvTriangular(ones, ComplexF64, sandwich ? V ⊗ V' : V, V; unitcell = size(network)) + @test network_value(network, env) ≈ 1 + # Arbitrary boundary phases must cancel across the unit cell. + for field in (:C, :Ea, :Eb), phase in (-1, cis(2π / 3)) + env′ = deepcopy(env) + getproperty(env′, field)[3, 1, 1] *= phase + @test network_value(network, env′) ≈ 1 + end + end + + # Two negative scalar tensors give a positive value per unit cell. + network = InfiniteTriangularNetwork(fill(-ones(Float64, V^3, V^3), 1, 2)) + env = CTMRGEnvTriangular(ones, ComplexF64, V, V; unitcell = size(network)) + @test network_value(network, env) ≈ 1 + + for α in (1.1, 1.1 + 0.1im) + g = only(Zygote.gradient(n -> real(network_value(n, env)), α * network)) + @test all(t -> t ≈ -conj(α) * ones(Float64, V^3, V^3), unitcell(g)) + end + + # Taking the cube root must not require cubing the full unit-cell magnitude. + for scale in (1.0e-30, 1.0e30) + network = InfiniteTriangularNetwork(fill(scale * ones(Float64, V^3, V^3), 2, 2)) + env = CTMRGEnvTriangular(ones, ComplexF64, V, V; unitcell = size(network)) + @test network_value(network, env) ≈ scale^length(network) + end +end + +@testset "CTM_triangular - Direct edge split" begin + Random.seed!(19461) + vspace = ComplexSpace(2) + envspace = ComplexSpace(3) + for (V, inds) in ((envspace ⊗ vspace, (3, 2)), (envspace ⊗ vspace ⊗ vspace', (4, 2, 3))), conditioning in (false, true), ill_conditioned in (false, true) + mat, ρL, ρR = (randn(ComplexF64, V, V) + 5 * id(V) for _ in 1:3) + if ill_conditioned + block(ρL, Trivial()) .= Diagonal([ones(dim(V) - 1); 1.0e-14]) + block(ρR, Trivial()) .= Diagonal([ones(dim(V) - 1); 1.0e-14]) + end + L, R = ρL / norm(ρL), ρR / norm(ρR) + if conditioning + _, SL, VL = svd_full(L) + UR, SR, _ = svd_full(R) + L = PEPSKit.sdiag_pow(SL, 0.5) * VL + R = UR * PEPSKit.sdiag_pow(SR, 0.5) + end + ρ = L * mat * R + for trunc in (notrunc(), truncrank(3), truncspace(envspace)) + Ea, Eb = PEPSKit._truncate_edge(mat, ρL, ρR, trunc; conditioning) + reconstructed = Eb * permute(Ea, ((1,), inds)) + U, S, Vᴴ = svd_trunc(ρ / norm(ρ); trunc) + @test space(Ea, 1) == space(S, 1) + @test L * reconstructed * R ≈ norm(ρ)^2 * U * S * Vᴴ + if trunc == notrunc() && !ill_conditioned + @test reconstructed / norm(reconstructed) ≈ mat / norm(mat) + end + end + end +end + +@testset "CTM_triangular - SVD pullback settings" begin + Random.seed!(21735) + V = ℂ^2 ⊗ ℂ^2 + for T in (Float64, ComplexF64), pullback in (:FullPullback, :TruncPullback) + tensors = ntuple(_ -> randn(T, V, V) + 5 * id(V), 3) + directions = map(t -> normalize(randn(T, space(t))), tensors) + target = randn(T, V, V) + decomposition_alg = SVDAdjoint(; fwd_alg = (; alg = :QRIteration), rrule_alg = (; alg = pullback, verbosity = 0)) + loss = function (mat, L, R) + Ea, Eb, = PEPSKit._truncate_edge(mat, L, R, truncrank(2), decomposition_alg) + return real(dot(target, Eb * permute(Ea, ((1,), (3, 2))))) + end + result = @test_logs min_level = Base.CoreLogging.Warn withgradient(loss, tensors...) + ε = 1.0e-5 + for i in 1:3 + plus = Base.setindex(tensors, tensors[i] + ε * directions[i], i) + minus = Base.setindex(tensors, tensors[i] - ε * directions[i], i) + @test real(dot(result.grad[i], directions[i])) ≈ (loss(plus...) - loss(minus...)) / (2ε) rtol = 1.0e-5 atol = 1.0e-7 + end + @test decomposition_alg.fwd_alg isa PEPSKit.MatrixAlgebraKit.QRIteration + end +end + +@testset "CTM_triangular - Direct edge coordinates" begin + Random.seed!(68191) + sz = (3, 4) + vspace = ComplexSpace(1) + offsets = ((-1, 0), (-1, 1), (0, 1), (1, 0), (1, -1), (0, -1)) + for sandwich in (false, true) + tensors = map(Iterators.product(1:sz[1], 1:sz[2])) do _ + if sandwich + V = vspace^3 ⊗ (vspace')^3 + return randn(ComplexF64, vspace, V), randn(ComplexF64, vspace, V) + end + return randn(ComplexF64, vspace^3, vspace^3) + end + network = InfiniteTriangularNetwork(tensors) + env = CTMRGEnvTriangular(randn, ComplexF64, sandwich ? vspace ⊗ vspace' : vspace, vspace; unitcell = sz) + at(A, d, r, c) = only(A[mod1(d, 6), mod1(r, sz[1]), mod1(c, sz[2])].data) + function site(r, c) + t = tensors[mod1(r, sz[1]), mod1(c, sz[2])] + return sandwich ? only(t[1].data) * conj(only(t[2].data)) : only(t.data) + end + for dir in 1:6, r in 1:sz[1], c in 1:sz[2] + δL, δR = offsets[mod1(dir + 4, 6)], offsets[mod1(dir + 3, 6)] + rL, cL = r + δL[1], c + δL[2] + rR, cR = r + δR[1], c + δR[2] + left = site(rL, cL) * at(env.Ea, dir + 3, rL + δR[1], cL + δR[2]) + right = site(rR, cR) * at(env.Eb, dir + 3, rR + δL[1], cR + δL[2]) + for i in 0:2 + δ = offsets[mod1(dir - i, 6)] + left *= at(env.C, dir - i, rL + δ[1], cL + δ[2]) + δ = offsets[mod1(dir + i + 1, 6)] + right *= at(env.C, dir + i + 1, rR + δ[1], cR + δ[2]) + end + ρL, ρR = PEPSKit.build_halfinfinite_environments(network, env, dir, r, c) + @test only(ρL.data) ≈ left + @test only(ρR.data) ≈ right + end + end +end + +@testset "CTM_triangular - Simultaneous edge update" begin + Random.seed!(38471) + vspace = ComplexSpace(2) + ket = randn(ComplexF64, ComplexSpace(2), vspace^3 ⊗ (vspace')^3) + pf = randn(ComplexF64, vspace^3, vspace^3) + for (sandwich, V) in ((pf, vspace), ((ket, ket), vspace ⊗ vspace')) + network = InfiniteTriangularNetwork(fill(sandwich, (1, 1))) + env = CTMRGEnvTriangular(randn, ComplexF64, V, ComplexSpace(2)) + trunc = PEPSKit.FixedSpaceTruncation() + Pas, Pbs, = PEPSKit.calculate_projectors(network, env, trunc, :twothirds) + mats, preliminary = PEPSKit.semi_renormalize(network, env, Pas, Pbs, trunc) + for conditioning in (false, true) + previous = deepcopy(preliminary) + edges = map(PEPSKit.eachcoordinate(network, 1:6)) do (dir, r, c) + ρL, ρR = PEPSKit.build_halfinfinite_environments(network, preliminary, dir, r, c) + trunc′ = PEPSKit._truncation_strategy(trunc, preliminary.Ea[dir, r, c]) + return PEPSKit._truncate_edge(mats[dir, r, c], ρL, ρR, trunc′; conditioning) + end + expected = PEPSKit.normalize_edges(CTMRGEnvTriangular(preliminary.C, getindex.(edges, 1), getindex.(edges, 2))) + updated, = PEPSKit.renormalize_edges(network, preliminary, mats, trunc; conditioning) + @test all(isapprox.(updated.Ea, expected.Ea)) + @test all(isapprox.(updated.Eb, expected.Eb)) + @test preliminary.Ea == previous.Ea + @test preliminary.Eb == previous.Eb + end + end +end + +@testset "CTM_triangular - Random tensor" begin + for conditioning in [true false] + for projector_alg in [:twothirds :full] + Random.seed!(79413165445) + alg = SimultaneousCTMRGTriangular(; + tol = ctmrg_tol, maxiter = ctmrg_maxiter, verbosity = ctmrg_verbosity, + conditioning, projector_alg + ) + pspace = ComplexSpace(2) + vspace = ComplexSpace(D) + envspace = ComplexSpace(χ) + + ket = randn(T, pspace, vspace ⊗ vspace ⊗ vspace ⊗ vspace' ⊗ vspace' ⊗ vspace') + bra = copy(ket) + pf = randn(T, vspace ⊗ vspace ⊗ vspace, vspace ⊗ vspace ⊗ vspace) + sandwiches = [pf, (ket, bra)] + # sandwiches = [(ket, bra), pf] + unitcell = (2, 2) + + for (sandwich, V) in zip(sandwiches, [vspace, vspace ⊗ vspace']) + # for (sandwich, vspace) in zip(sandwiches, [vspace ⊗ vspace', vspace]) + network = InfiniteTriangularNetwork(fill(sandwich, unitcell)) + env₀ = CTMRGEnvTriangular(randn, T, V, envspace; unitcell) + env, info = leading_boundary(env₀, network, alg) + @test info.convergence_metric < 1.0 # this is not much of a test + end + end + end +end + +@testset "CTM_triangular - Classical Ising" begin + χ_local = 20 + T_local = Float64 + + for conditioning in [true false] + for projector_alg in [:twothirds :full] + Random.seed!(156484561351) + + alg = SimultaneousCTMRGTriangular(; + tol = ctmrg_tol, maxiter = ctmrg_maxiter, verbosity = ctmrg_verbosity, + conditioning, projector_alg + ) + sz = (1, 1) + T = classical_ising_triangular(ising_βc_triangular) + pf = InfiniteTriangularNetwork(fill(T, sz)) + + vspace = codomain(T)[1] + envspace = ComplexSpace(χ_local) + env₀ = CTMRGEnvTriangular(randn, T_local, vspace, envspace; unitcell = sz) + env, info = leading_boundary(env₀, pf, alg) + + nw_value = network_value(pf, env) + lz = real(log(nw_value)) + fs = lz * -1 / ising_βc_triangular + @test fs ≈ f_onsager_triangular rtol = 1.0e-4 + end + end +end + + +## Gradient tests for triangular CTMRG +# ------------------------------------ + +Pspaces = [ComplexSpace(2)] +Vspaces = [ComplexSpace(D)] +Espaces = [ComplexSpace(χ)] +# TODO: actually add triangular lattice models... +models = [transverse_field_ising(InfiniteSquare(); J = 1.0, g = 1.0)] +names = ["Heisenberg"] + +gradtol = 1.0e-4 +ctmrg_verbosity = 0 +ctmrg_algs = [[:triangular]] +projector_algs = [[:twothirds, :full]] +conditionings = [[true, false]] +# svd_rrule_algs = nothing # TODO: allow this +gradient_algs = [[nothing]] +# gradient_iterschemes = [[:fixed, :diffgauge]] # TODO: support this +steps = -0.01:0.005:0.01 + +# have to use a custom retraction, since the state type of hardcoded... +function peps_triangular_retract(x, η, α) + peps = x[1] + env = deepcopy(x[2]) + + retractions = norm_preserving_retract.(unitcell(peps), unitcell(η), α) + peps´ = typeof(peps)(map(first, retractions)) + ξ = typeof(peps)(map(last, retractions)) + + return (peps´, env), ξ +end + + +@testset "Triangular CTMRG energy gradient test for $(names[i]) model" verbose = true for i in + eachindex( + models + ) + Pspace = Pspaces[i] + Vspace = Vspaces[i] + Espace = Espaces[i] + # calgs = ctmrg_algs[i] + palgs = projector_algs[i] + conds = conditionings[i] + # salgs = svd_rrule_algs[i] + galgs = gradient_algs[i] + # gischemes = gradient_iterschemes[i] + + # TODO: allow passing actual LocalOperator + H = models[i] + O1 = last(first(filter(t -> length(first(t)) == 1, H.terms))) + O2 = last(first(filter(t -> length(first(t)) == 2, H.terms))) + + @testset "ctmrg_alg=:simultaneous_triangular, projector_alg=:$projector_alg, conditioning=$conditioning, gradient_alg=:$gradient_alg" for ( + projector_alg, conditioning, gradient_alg, + ) in Iterators.product(palgs, conds, galgs) + + @info "optimtest of ctmrg_alg=:simultaneous_triangular, projector_alg=:$projector_alg, conditioning=$conditioning, gradient_alg=:$gradient_alg on $(names[i])" + Random.seed!(42039482030) + dir = InfinitePEPSTriangular(Pspace, Vspace) + psi = InfinitePEPSTriangular(Pspace, Vspace) + + # TODO: actually go through CTMRGAlgorithm selection + contrete_ctmrg_alg = SimultaneousCTMRGTriangular(; + projector_alg, conditioning, verbosity = ctmrg_verbosity, + trunctype = :FixedSpaceTruncation + ) + # instantiate because hook_pullback doesn't go through the keyword selector... + concrete_gradient_alg = if isnothing(gradient_alg) + nothing + else + PEPSKit.GradientAlgorithm(; alg = gradient_alg, tol = gradtol) + end + env, = leading_boundary( + CTMRGEnvTriangular(randn, ComplexF64, Vspace ⊗ Vspace', Espace), + psi, contrete_ctmrg_alg + ) + alphas, fs, dfs1, dfs2 = OptimKit.optimtest( + (psi, env), + dir; + alpha = steps, + retract = peps_triangular_retract, + inner = PEPSKit.real_inner, + ) do (peps, env) + E, g = Zygote.withgradient(peps) do psi + env2, = PEPSKit.hook_pullback( + leading_boundary, + env, + psi, + contrete_ctmrg_alg; + alg_rrule = concrete_gradient_alg, + ) + return PEPSKit.energy(psi, env2, O1, O2) + end + + return E, only(g) + end + @test dfs1 ≈ dfs2 atol = 1.0e-2 + end +end + +@testset "CTM_triangular - Truncation diagnostics" begin + V = ComplexSpace(2) ⊗ ComplexSpace(2) + for T in (Float64, ComplexF64), conditioning in (false, true), scale in (1.0e-30, 1.0, 1.0e30) + for values in ([4.0, 3.0, 2.0, 1.0], [4.0, 3.0, 1.0e-9, 2.0e-10], [4.0, 3.0, 1.0e-14, 2.0e-15], [4.0, 3.0, 0.0, 0.0]) + mat = TensorMap(Matrix(Diagonal(T.(scale .* values))), V, V) + for rank in 1:5 + _, _, error = PEPSKit._truncate_edge(mat, id(V), id(V), truncrank(rank); conditioning) + kept = values[1:min(rank, end)] + discarded = [values[(min(rank, end) + 1):end]; filter(<(eps(Float64)^(3 / 4) * maximum(kept)), kept)] + @test error ≈ norm(discarded) / norm(kept) rtol = 1.0e-12 + end + _, _, error = PEPSKit._truncate_edge(mat, id(V), id(V), notrunc(); conditioning) + discarded = filter(<(eps(Float64)^(3 / 4) * maximum(values)), values) + @test error ≈ norm(discarded) / norm(values) rtol = 1.0e-12 + end + end + + # Non-Abelian sector dimensions must weight both the cutoff loss and norm(S). + S = DiagonalTensorMap([4.0, 3.0, 1.0e-14], SU2Space(0 => 1, 1 // 2 => 2)) + isqrtS = PEPSKit.sdiag_pow(S, -0.5) + error = PEPSKit._relative_truncation_error(S, isqrtS, 2.0e-15) + @test error ≈ sqrt((2.0e-15)^2 + 2 * (1.0e-14)^2) / sqrt(4.0^2 + 2 * 3.0^2) + for S in (DiagonalTensorMap(zeros(2), ComplexSpace(2)), DiagonalTensorMap(Float64[], ComplexSpace(0))) + @test iszero(PEPSKit._relative_truncation_error(S, S, 0.0)) + @test isinf(PEPSKit._relative_truncation_error(S, S, 1.0)) + end +end + +@testset "CTM_triangular - Convergence spectra" begin + Random.seed!(71834) + vspace = ComplexSpace(2) + alg = SimultaneousCTMRGTriangular(; verbosity = 0) + for V in (vspace, vspace ⊗ vspace') + env = CTMRGEnvTriangular(randn, ComplexF64, V, vspace; unitcell = (1, 2)) + CS, TS = PEPSKit.convergence_spectra(env, alg) + @test size(CS) == (6, 1, 2) + @test size(TS) == (12, 1, 2) + @test first(PEPSKit.calc_convergence(env, CS, TS, alg)) == 0 + for field in (:C, :Ea, :Eb) + changed = deepcopy(env) + tensors = getproperty(changed, field) + expected = norm(tensors[6, 1, 2]) + tensors[6, 1, 2] *= 2 + @test first(PEPSKit.calc_convergence(changed, CS, TS, alg)) ≈ expected + end + end +end + +@testset "CTM_triangular - Iteration diagnostics" begin + Random.seed!(88193) + vspace = ComplexSpace(2) + sz = (1, 2) + for sandwich in (false, true) + tensors = map(Iterators.product(1:sz[1], 1:sz[2])) do _ + if sandwich + ket = randn(ComplexF64, vspace, vspace^3 ⊗ (vspace')^3) + return ket, copy(ket) + end + return randn(ComplexF64, vspace^3, vspace^3) + end + network = InfiniteTriangularNetwork(tensors) + env₀ = CTMRGEnvTriangular(randn, ComplexF64, sandwich ? vspace ⊗ vspace' : vspace, vspace; unitcell = sz) + for projector_alg in (:twothirds, :full), conditioning in (false, true), trunctype in (:truncrank, :FixedSpaceTruncation) + alg = SimultaneousCTMRGTriangular(; projector_alg, conditioning, trunctype, verbosity = 0) + env₁, info₁ = PEPSKit.ctmrg_iteration(network, env₀, alg) + _, info₂ = PEPSKit.ctmrg_iteration(network, env₁, alg) + for (maxiter, miniter, tol, expected) in ((1, 1, Inf, info₁), (2, 2, Inf, info₂), (1, 1, 0.0, info₁)) + alg = SimultaneousCTMRGTriangular(; projector_alg, conditioning, trunctype, verbosity = 0, maxiter, miniter, tol) + env, info = leading_boundary(env₀, network, alg) + CS, TS = PEPSKit.convergence_spectra(maxiter == 1 ? env₀ : env₁, alg) + @test info.convergence_error ≈ first(PEPSKit.calc_convergence(env, CS, TS, alg)) + metrics = info.contraction_metrics + @test metrics == expected.contraction_metrics + @test all(isfinite, metrics) + @test metrics.corner_truncation_error > 0 + @test metrics.edge_truncation_error > 0 + @test metrics.truncation_error == max(metrics.corner_truncation_error, metrics.edge_truncation_error) + @test info.converged == isinf(tol) + @test info.convergence_error == info.convergence_metric + end + end + end +end + +@testset "Triangular algorithm selection" begin + V = ComplexSpace(1) + state = InfinitePEPSTriangular(ones(ComplexF64, ComplexSpace(2), V^3 ⊗ (V')^3)) + env₀ = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', V) + kwargs = (; tol = 1.0e-9, maxiter = 2, miniter = 2, verbosity = 0, conditioning = false, trunctype = :FixedSpaceTruncation) + alg = SimultaneousCTMRGTriangular(; kwargs...) + @test PEPSKit.select_algorithm(leading_boundary, env₀; kwargs...) == alg + env, info = leading_boundary(env₀, state, alg) + for network in (state, InfiniteTriangularNetwork(state)) + env′, info′ = leading_boundary(env₀, network; alg = :SimultaneousCTMRGTriangular, kwargs...) + @test (env′.C, env′.Ea, env′.Eb) == (env.C, env.Ea, env.Eb) + @test info′ == info + end + @test_throws ArgumentError leading_boundary(env₀, state; alg = :unknown) + + decomposition_alg = SVDAdjoint(; fwd_alg = (; alg = :QRIteration), rrule_alg = (; verbosity = 1)) + configured = PEPSKit.select_algorithm(leading_boundary, env₀; kwargs..., decomposition_alg) + @test configured.decomposition_alg === decomposition_alg + @test SimultaneousCTMRGTriangular(; decomposition_alg = (; alg = :QRIteration)).decomposition_alg.fwd_alg isa PEPSKit.MatrixAlgebraKit.QRIteration + @test PEPSKit.PEPSOptimize(; boundary_alg = (; alg = :SimultaneousCTMRGTriangular, decomposition_alg), precondition_alg = nothing).boundary_alg.decomposition_alg === decomposition_alg + env′, info′ = leading_boundary(env₀, state, configured) + @test all(isapprox.(env′.C, env.C)) + @test all(isapprox.(env′.Ea, env.Ea)) + @test all(isapprox.(env′.Eb, env.Eb)) + @test info′.converged == info.converged + @test SimultaneousCTMRGTriangular().projector_alg === :full + @test SimultaneousCTMRGTriangular(; projector_alg = :twothirds).projector_alg === :twothirds + @test PEPSKit._set_fixed_truncation(configured).decomposition_alg === decomposition_alg +end + +@testset "Triangular fixedpoint interface" begin + P, V = ComplexSpace(2), ComplexSpace(1) + state = InfinitePEPSTriangular(ones(Float64, P, V^3 ⊗ (V')^3)) + env = CTMRGEnvTriangular(ones, Float64, V ⊗ V', V) + alg = PEPSOptimize(SimultaneousCTMRGTriangular(), nothing, OptimKit.LBFGS(), nothing, true, nothing) + @test isnothing(PEPSKit.check_input(fixedpoint, state, env, alg)) + @test_throws ArgumentError fixedpoint(LocalOperator(physicalspace(state)), state, nothing, alg) +end + +@testset "Triangular native optimization" begin + P, V = ComplexSpace(2), ComplexSpace(1) + Sz = TensorMap([1.0 0; 0 -1] / 2, P, P) + for T in (Float64, ComplexF64) + tensors = [TensorMap(reshape(T[1, a], 2, 1, 1, 1, 1, 1, 1), P, V^3 ⊗ (V')^3) for a in [0.6 1.3]] + state = InfinitePEPSTriangular(tensors) + original = deepcopy(state) + env = CTMRGEnvTriangular(ones, T, V ⊗ V', V; unitcell = size(state)) + H = LocalOperator(physicalspace(state), [((1, c),) => Sz for c in 1:2]) + boundary_alg = (; miniter = 1, maxiter = 1) + kwargs = (; boundary_alg, gradient_alg = nothing, precondition_alg = nothing, optimizer_alg = (; maxiter = 2, tol = 0), verbosity = 0) + alg = PEPSKit.select_algorithm(fixedpoint, env; kwargs...) + @test PEPSKit.parent_alg(alg.boundary_alg) isa SimultaneousCTMRGTriangular + direct_alg = PEPSOptimize(; boundary_alg = (; alg = :SimultaneousCTMRGTriangular, miniter = 1, maxiter = 1), gradient_alg = nothing, precondition_alg = nothing) + @test direct_alg.boundary_alg isa SimultaneousCTMRGTriangular + @test direct_alg.boundary_alg.maxiter == 1 + @test PEPSKit.adapt_solver(alg.boundary_alg; iter = 2, g_global = 0.1) isa SimultaneousCTMRGTriangular + iterations = Int[] + finalize! = (x, f, g, iter) -> (push!(iterations, iter); (x, f, g)) + initial_cost = cost_function(state, env, H) + ψ, env′, E, info = if T === Float64 + fixedpoint(H, state, env, alg; finalize!) + else + fixedpoint(H, state, env; kwargs..., finalize!) + end + @test ψ isa typeof(state) + @test env′ isa typeof(env) + @test E < initial_cost + @test E ≈ cost_function(ψ, env′, H) + @test all(isapprox.(norm.(unitcell(ψ)), 1)) + @test state == original + @test iterations == [1, 2] + @test info.costs[end] ≈ E + @test info.last_gradient isa typeof(state) + @test all(isfinite, info.gradnorms) + @test length(info.contraction_metrics) == length(info.gradnorms_unitcell) == length(info.times) == 2 + + if T === ComplexF64 + ψfixed, envfixed, Efixed, infofixed = fixedpoint( + H, state, env; kwargs..., gradient_alg = (; tol = 1.0e-8), + boundary_alg = (; miniter = 2, maxiter = 10, tol = 1.0e-10), + ) + @test Efixed ≈ E + @test infofixed.costs ≈ info.costs + @test infofixed.gradnorms ≈ info.gradnorms + @test Efixed ≈ cost_function(ψfixed, envfixed, H) + else + @test_throws ArgumentError fixedpoint(H, state, env; kwargs..., gradient_alg = (;)) + end + @test_throws ArgumentError fixedpoint(H, state, env; kwargs..., gradient_alg = (; alg = :ImplicitGradient)) + @test_throws ArgumentError fixedpoint(H, state, nothing, alg) + @test_throws ArgumentError fixedpoint(H, state, env; boundary_alg) + @test_throws ArgumentError fixedpoint(H, state, env; boundary_alg, gradient_alg = nothing) + @test_throws ArgumentError fixedpoint(H, state, env; kwargs..., symmetrization = Rotate()) + @test_throws ArgumentError fixedpoint(H, state, env; kwargs..., boundary_alg = SimultaneousCTMRG()) + end +end + +@testset "Triangular AF Heisenberg optimization" begin + rng = StableRNG(48631) + V = ComplexSpace(2) + state = InfinitePEPSTriangular((args...) -> randn(rng, args...), ComplexF64, V, V) + Sx = TensorMap(ComplexF64[0 1; 1 0] / 2, V, V) + Sy = TensorMap(ComplexF64[0 -im; im 0] / 2, V, V) + Sz = TensorMap(ComplexF64[1 0; 0 -1] / 2, V, V) + # Ferromagnetic nearest-neighbor Heisenberg energy per site. + h = -(Sx ⊗ Sx + Sy ⊗ Sy + Sz ⊗ Sz) + offsets = (CartesianIndex(0, 1), CartesianIndex(1, -1), CartesianIndex(1, 0)) + H = LocalOperator(physicalspace(state), [(CartesianIndex(1, 1), CartesianIndex(1, 1) + δ) => h for δ in offsets]) + env₀ = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', ComplexSpace(4)) + boundary_alg = SimultaneousCTMRGTriangular(; projector_alg = :full, tol = 1.0e-10, maxiter = 100, miniter = 4, verbosity = 0) + gradient_alg = PEPSKit.FixedPointGradient(; solver_alg = (; alg = :GMRES, krylovdim = 30), tol = 1.0e-8, maxiter = 10, verbosity = 0) + env₀, info₀ = leading_boundary(env₀, state, boundary_alg) + @test info₀.converged + # A mixed one-site density matrix rules out a padded product state. + m = [real(expectation_value(state, LocalOperator(physicalspace(state), ((1, 1),) => S), env₀)) for S in (Sx, Sy, Sz)] + @test 0.5 ≤ 0.5 + 2sum(abs2, m) < 1 - 1.0e-4 + E₀ = cost_function(state, env₀, H) + ψ, env, E, info = fixedpoint( + H, state, env₀; boundary_alg, gradient_alg, precondition_alg = nothing, + optimizer_alg = (; maxiter = 1, tol = 0), verbosity = 0, + ) + @test E < E₀ + @test info.costs[end] ≈ E + @test length(info.contraction_metrics) == 1 + @test all(isfinite, info.gradnorms) + @test info.last_gradient isa typeof(ψ) + @test isfinite(norm(info.last_gradient)) && norm(info.last_gradient) > 0 + @test all(isapprox.(norm.(unitcell(ψ)), 1)) + env′, info′ = leading_boundary(env, ψ, boundary_alg) + @test info′.converged + @test E ≈ cost_function(ψ, env′, H) atol = 1.0e-8 + # Check the returned nearest-neighbor gradient along one normalized direction. + direction = info.last_gradient / norm(info.last_gradient) + ε = 1.0e-4 + energies = map((-ε, ε)) do step + ψ′ = ψ + step * direction + env′, info′ = leading_boundary(env, ψ′, boundary_alg) + @test info′.converged + cost_function(ψ′, env′, H) + end + @test real(dot(info.last_gradient, direction)) ≈ (energies[2] - energies[1]) / (2ε) rtol = 1.0e-4 atol = 1.0e-6 +end + +@testset "Triangular gauge fixing" begin + Random.seed!(63412) + for (V, χ, sandwich) in ( + (ComplexSpace(2), ComplexSpace(2), false), + (ComplexSpace(2), ComplexSpace(2), true), + (Z2Space(0 => 1, 1 => 1), Z2Space(0 => 1, 1 => 1), true), + ) + sz = (2, 3) # diagonal boundaries traverse both unit-cell axes + env = CTMRGEnvTriangular(randn, ComplexF64, sandwich ? V ⊗ V' : V, χ; unitcell = sz) + signs = map((:Ea, :Eb)) do _ + map(env.C) do t + first(left_orth(randn(ComplexF64, space(t, 1), space(t, 1)); positive = true)) + end + end + rotated = PEPSKit.fix_relative_phases(env, (; Ea = signs[1], Eb = signs[2])) + original = deepcopy(env) + # Check that the gauges cancel in the physical boundary contractions. + network = if sandwich + InfiniteTriangularNetwork(InfinitePEPSTriangular([randn(ComplexF64, V, V^3 ⊗ (V')^3) for _ in 1:sz[1], _ in 1:sz[2]])) + else + InfiniteTriangularNetwork([randn(ComplexF64, V^3, V^3) for _ in 1:sz[1], _ in 1:sz[2]]) + end + for contract in (PEPSKit._contract_corners, PEPSKit._contract_edges_0, PEPSKit._contract_edges_60, PEPSKit._contract_edges_120) + @test contract((1, 2), network, rotated) ≈ contract((1, 2), network, env) + end + rotated = CTMRGEnvTriangular(im .* rotated.C, -rotated.Ea, -im .* rotated.Eb) + fixed = PEPSKit.gauge_fix(rotated, env, PEPSKit.ScramblingEnvGauge()) + @test PEPSKit.calc_elementwise_convergence(env, fixed) < 1.0e-10 + @test (env.C, env.Ea, env.Eb) == (original.C, original.Ea, original.Eb) + different = CTMRGEnvTriangular(ones, ComplexF64, sandwich ? V ⊗ V' : V, χ; unitcell = (1, 1)) + @test_throws ArgumentError PEPSKit.gauge_fix(different, env, PEPSKit.ScramblingEnvGauge()) + end +end + +@testset "Triangular fixed-point differentiation" begin + Random.seed!(48531) + V = ComplexSpace(2) + fields = [0.02 -0.03] + β = 0.1 + W = sqrt(Symmetric([exp(β / 2) exp(-β / 2); exp(-β / 2) exp(β / 2)])) + state = InfinitePEPSTriangular( + map(fields) do h + data = [exp(h * (3 - 2s) / 2 + im * h * s) * prod(W[s, i] for i in Tuple(I)) for s in 1:2, I in CartesianIndices((2, 2, 2, 2, 2, 2))] + TensorMap(data, V, V^3 ⊗ (V')^3) + end + ) + Sz = TensorMap(ComplexF64[1 0; 0 -1] / 2, V, V) + H = LocalOperator(physicalspace(state), [((1, c),) => Sz for c in 1:2]) + direction = InfinitePEPSTriangular(map(t -> Sz * t, unitcell(state))) + env₀ = CTMRGEnvTriangular(ones, ComplexF64, V ⊗ V', V; unitcell = size(state)) + gradalg = PEPSKit.FixedPointGradient(; solver_alg = (; alg = :GMRES, krylovdim = 30), tol = 1.0e-8, maxiter = 10, verbosity = 0) + for (projector_alg, conditioning) in ((:twothirds, false), (:full, true)) + alg = SimultaneousCTMRGTriangular(; projector_alg, conditioning, maxiter = 200, miniter = 4, tol = 1.0e-10, verbosity = 0) + env, info = leading_boundary(env₀, state, alg) + @test info.converged + env′, = PEPSKit.ctmrg_iteration(InfiniteTriangularNetwork(state), env, PEPSKit._set_fixed_truncation(alg)) + fixed = PEPSKit.gauge_fix(env′, env, PEPSKit.ScramblingEnvGauge()) + @test PEPSKit.calc_elementwise_convergence(env, fixed) < 1.0e-7 + E, gs = withgradient(state) do A + env′, = PEPSKit.hook_pullback(leading_boundary, env, A, alg; alg_rrule = gradalg) + cost_function(A, env′, H) + end + g = only(gs) + @test g isa typeof(state) + # The custom rule reconverges from env before evaluating the energy. + env_ref, = leading_boundary(env, state, alg) + @test E ≈ cost_function(state, env_ref, H) + ε = 1.0e-4 + energies = map((-ε, ε)) do step + A = state + step * direction + env′, info′ = leading_boundary(env, A, alg) + @test info′.converged + cost_function(A, env′, H) + end + @test real(dot(g, direction)) ≈ (energies[2] - energies[1]) / (2ε) rtol = 1.0e-5 + end +end + +@testset "Triangular nearest-neighbor energy" begin + Random.seed!(28514) + pspace, vspace = ComplexSpace(2), ComplexSpace(1) + vspaces = vspace^3 ⊗ (vspace')^3 + Sx = TensorMap(ComplexF64[0 1; 1 0] / 2, pspace, pspace) + Sy = TensorMap(ComplexF64[0 -im; im 0] / 2, pspace, pspace) + Sz = TensorMap(ComplexF64[1 0; 0 -1] / 2, pspace, pspace) + spins = (Sx, Sy, Sz) + # Unequal site norms expose neighbor-index errors even in identity contractions. + vectors = [randn(ComplexF64, 2) for r in 1:3, c in 1:4] + tensors = map(v -> TensorMap(reshape(v, 2, 1, 1, 1, 1, 1, 1), pspace, vspaces), vectors) + state = InfinitePEPSTriangular(tensors) + network = InfiniteTriangularNetwork(state) + env = CTMRGEnvTriangular(ones, ComplexF64, vspace ⊗ vspace', vspace; unitcell = size(state)) + contractions = (PEPSKit._contract_edges_0, PEPSKit._contract_edges_60, PEPSKit._contract_edges_120) + offsets = ((0, 1), (1, -1), (1, 0)) + pairs = ((Sx, Sz), (Sy, Sx), (Sz, Sy)) + ops = map(p -> p[1] ⊗ p[2], pairs) + expectation(v, op) = dot(v, reshape(convert(Array, op), 2, 2) * v) / dot(v, v) + expected = 0.0 + expected_horizontal = 0.0 + terms = Pair[] + for r in 1:3, c in 1:4 + v = vectors[r, c] + onsite = real(expectation(v, Sz)) + expected += onsite + expected_horizontal += onsite + push!(terms, ((r, c),) => Sz) + for (i, (contract, (δr, δc))) in enumerate(zip(contractions, offsets)) + w = vectors[mod1(r + δr, 3), mod1(c + δc, 4)] + n = dot(v, v) * dot(w, w) + bond = expectation(v, pairs[i][1]) * expectation(w, pairs[i][2]) + @test contract((r, c), network, env) ≈ n + @test contract((r, c), network, env, id(pspace ⊗ pspace)) ≈ n + @test contract((r, c), network, env, ops[i]) / n ≈ bond + expected += real(bond) + i == 1 && (expected_horizontal += real(bond)) + push!(terms, ((r, c), (r + δr, c + δc)) => ops[i]) + end + end + @test PEPSKit.energy(state, env, Sz, ops) ≈ expected + @test PEPSKit.energy(state, env, Sz, ops[1]) ≈ expected_horizontal + + H = LocalOperator(physicalspace(state), terms) + @test expectation_value(state, H, env) ≈ expected + @test expectation_value(state, im * H, env) ≈ im * expected + @test cost_function(state, env, H) ≈ expected + @test iszero(expectation_value(state, LocalOperator(physicalspace(state)), env)) + @test iszero(expectation_value(InfinitePEPS(pspace, vspace), LocalOperator(fill(pspace, 1, 1)), nothing)) + + # Distinct bra/ket and an asymmetric diagonal bond check operator ordering. + bra = InfinitePEPSTriangular(map(t -> 2Sx * t, tensors)) + bond = LocalOperator(physicalspace(state), ((1, 1), (2, 0)) => ops[2]) + transition(v, op) = dot(reverse(v), reshape(convert(Array, op), 2, 2) * v) / dot(reverse(v), v) + @test expectation_value(bra, bond, state, env) ≈ transition(vectors[1, 1], Sy) * transition(vectors[2, 4], Sx) + @test_throws ArgumentError expectation_value(state, LocalOperator(fill(pspace, 1, 1), ((1, 1),) => Sz), env) + @test_throws ArgumentError expectation_value(state, LocalOperator(physicalspace(state), ((1, 1), (1, 3)) => ops[1]), env) + @test_throws ArgumentError expectation_value(state, LocalOperator(physicalspace(state), ((1, 1), (1, 2), (2, 1)) => Sz ⊗ Sz ⊗ Sz), env) + + # Spin-1/2 120° product state: each of the three bonds contributes -1/8 per site. + rotation(θ) = TensorMap([cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)], pspace, pspace) + A = TensorMap(reshape([1.0, 0.0], 2, 1, 1, 1, 1, 1, 1), pspace, vspaces) + B, C = rotation(2π / 3) * A, rotation(4π / 3) * A + state = InfinitePEPSTriangular([A B C; C A B; B C A]) + env = CTMRGEnvTriangular(ones, Float64, vspace ⊗ vspace', vspace; unitcell = (3, 3)) + h = sum(S ⊗ S for S in spins) + @test PEPSKit.energy(state, env, zero(Sz), (h, h, h)) / length(state) ≈ -3 / 8 +end diff --git a/test/testsuite/ctmrg/unitcell.jl b/test/testsuite/ctmrg/unitcell.jl index fb7b24e70..abcd74ebe 100644 --- a/test/testsuite/ctmrg/unitcell.jl +++ b/test/testsuite/ctmrg/unitcell.jl @@ -31,11 +31,11 @@ function test_unitcell( Pspaces, [ (c,) => adapt( - AT, randn( - scalartype(peps), - Pspaces[c], Pspaces[c], - ) - ) for c in CartesianIndices(unitcell) + AT, randn( + scalartype(peps), + Pspaces[c], Pspaces[c], + ) + ) for c in CartesianIndices(unitcell) ]..., ) @test expectation_value(peps, random_op, env) isa Number diff --git a/test/testsuite/utility/diff_maps.jl b/test/testsuite/utility/diff_maps.jl index a1a904f87..d986327c3 100644 --- a/test/testsuite/utility/diff_maps.jl +++ b/test/testsuite/utility/diff_maps.jl @@ -17,6 +17,8 @@ function utility_diff_maps(AT) gradient = only(Zygote.gradient(X -> real(sum(sdiag_pow(X, power).data)), s)) @test storagetype(gradient) === storagetype(s) @test gradient.data ≈ power .* s.data .^ (power - 1) + dense_gradient = only(Zygote.gradient(X -> real(sum(sdiag_pow(X, power).data)), TensorMap(s))) + @test dense_gradient ≈ TensorMap(gradient) end end end diff --git a/test/testsuite/utility/svd_wrapper.jl b/test/testsuite/utility/svd_wrapper.jl index 7259b347a..6c98e5a0a 100644 --- a/test/testsuite/utility/svd_wrapper.jl +++ b/test/testsuite/utility/svd_wrapper.jl @@ -10,6 +10,52 @@ using MatrixAlgebraKit: TruncatedAlgorithm, diagview, svd_trunc_no_error function utility_svd_wrapper(AT) return @testset "svd_wrapper ($AT)" begin + @testset "Full SVD" begin + Random.seed!(41279) + for T in (Float64, ComplexF64), (V, W) in ( + (ℂ^2, ℂ^3), (ℂ^3, ℂ^2), + (Z2Space(0 => 2, 1 => 1), Z2Space(0 => 1, 1 => 2)), + ), pullback in (:FullPullback, :TruncPullback) + A = adapt(AT, randn(T, V, W)) + R = adapt(AT, randn(T, V, W)) + direction = normalize(adapt(AT, randn(T, V, W))) + alg = SVDAdjoint(; rrule_alg = (; alg = pullback)) + loss = function (A) + U, S, V = svd_full(A, alg) + return real(dot(R, U * (S * S' * S) * V)) + end + original = copy(A) + U, S, V = svd_full(A, alg) + @test U * S * V ≈ A + @test space.((U, S, V)) == space.(svd_full(A)) + @test A == original + @test prod(svd_full(A', alg)) ≈ A' + g = only(gradient(loss, A)) + g′ = only( + gradient(A) do A + U, S, V = svd_full(A', alg) + real(dot(R', U * (S * S' * S) * V)) + end + ) + @test g′ ≈ g + ε = 1.0e-5 + @test real(dot(g, direction)) ≈ (loss(A + ε * direction) - loss(A - ε * direction)) / (2ε) rtol = 1.0e-6 atol = 1.0e-8 + end + + A = adapt(AT, randn(ComplexF64, ℂ^2, ℂ^3)) + @test_throws ArgumentError svd_full(A, SVDAdjoint(; fwd_alg = (; alg = :GKL))) + for verbosity in (0, 1) + alg = SVDAdjoint(; rrule_alg = (; verbosity)) + (U, S, V), back = rrule(svd_full!, copy(A), alg) + if verbosity == 0 + @test_logs min_level = Base.CoreLogging.Warn back((im * U, zero(S), zero(V))) + else + @test_logs (:warn, r"`svd` cotangents sensitive to gauge choice") back((im * U, zero(S), zero(V))) + end + @test back((ZeroTangent(), ZeroTangent(), ZeroTangent()))[2] isa ZeroTangent + end + end + # Gauge-invariant loss function function lossfun(svd_trunc_f, A, alg, R = randn(space(A)), trunc = notrunc()) alg = @set alg.fwd_alg = TruncatedAlgorithm(alg.fwd_alg, trunc) diff --git a/test/types/infinitetriangularnetwork.jl b/test/types/infinitetriangularnetwork.jl new file mode 100644 index 000000000..bf2797673 --- /dev/null +++ b/test/types/infinitetriangularnetwork.jl @@ -0,0 +1,85 @@ +using Test +using TensorKit +using PEPSKit +using Zygote +using LinearAlgebra +using PEPSKit: unitcell, virtualspace + +@testset "Triangular constructors" begin + for T in (Float64, ComplexF64), (V, W) in ( + (ComplexSpace(1), ComplexSpace(2)), + (U1Space(0 => 1), U1Space(0 => 1, 1 => 1)), + ) + NW = [V W V; W V W] + NE = [W V V; V W W] + E = [V V W; W V V] + P = fill(W, 2, 3) + state = InfinitePEPSTriangular(ones, T, P, NW, NE, E) + A = unitcell(state) + for (dir, δr, δc) in ((1, -1, 0), (2, -1, 1), (3, 0, 1)), r in 1:2, c in 1:3 + @test virtualspace(state, r, c, dir) == virtualspace(state, r + δr, c + δc, dir + 3)' + end + @test InfinitePEPSTriangular(copy(A)) == state + @test unitcell(InfiniteTriangularNetwork(state)) == map(tuple, A, A) + @test_throws ArgumentError InfinitePEPSTriangular(ones, T, P, NW[:, 1:2], NE, E) + + # Same six virtual legs in the classical tensor layout W SW SE ← NW NE E. + pf(t) = ones(T, domain(t)[6]' ⊗ domain(t)[5]' ⊗ domain(t)[4]', domain(t)[1] ⊗ domain(t)[2] ⊗ domain(t)[3]) + tensors = map(pf, A) + @test unitcell(InfiniteTriangularNetwork(tensors)) == tensors + + # Break each leg at a periodic boundary without changing its dimension. + for dir in 1:6 + invalid = copy(A) + spaces = collect(domain(A[1, 3])) + spaces[dir] = spaces[dir]' + invalid[1, 3] = ones(T, codomain(A[1, 3]), ProductSpace(spaces...)) + @test_throws SpaceMismatch InfinitePEPSTriangular(invalid) + @test_throws SpaceMismatch InfiniteTriangularNetwork(map(tuple, invalid, A)) + @test_throws SpaceMismatch InfiniteTriangularNetwork(map(tuple, A, invalid)) + @test_throws SpaceMismatch InfiniteTriangularNetwork(map(pf, invalid)) + end + + uniform = InfinitePEPSTriangular(ones, T, W, V; unitcell = (1, 1)) + @test size(InfinitePEPSTriangular(uniform[1, 1]; unitcell = (2, 3))) == (2, 3) + end + + P = Vect[FermionParity](0 => 1, 1 => 1) + V = Vect[FermionParity](0 => 1) + @test size(InfinitePEPSTriangular(ones, Float64, P, V)) == (1, 1) + @test_throws ArgumentError InfinitePEPSTriangular(ones(Float64, P', V^3 ⊗ (V')^3)) +end + +@testset "Triangular constructor differentiation" begin + V = ComplexSpace(1) + A = unitcell(InfinitePEPSTriangular(ones, Float64, V, V; unitcell = (1, 2))) + for construct in (InfinitePEPSTriangular, a -> InfiniteTriangularNetwork(map(tuple, a, a))) + f(a) = sum(unitcell(construct(a))) do t + x = t isa Tuple ? first(t) : t + return real(dot(x, x)) + end + g = only(Zygote.gradient(f, A)) + @test all(isapprox.(g, 2 .* A)) + end +end + +@testset "Triangular environment types" begin + V = ComplexSpace(1) + for sandwich in (false, true) + env = CTMRGEnvTriangular(ones, Float64, sandwich ? V ⊗ V' : V, V) + complex_edges = complex.(env.Ea) + for (Ea, Eb, E) in ( + (env.Ea, env.Eb, eltype(env.Ea)), + (complex_edges, env.Eb, eltype(complex_edges)), + (env.Ea, complex_edges, eltype(complex_edges)), + ) + env′ = CTMRGEnvTriangular(env.C, Ea, Eb) + @test PEPSKit.cornertype(env′) === eltype(env.C) + @test PEPSKit.cornertype(typeof(env′)) === eltype(env.C) + @test PEPSKit.edgetype(env′) === E + @test PEPSKit.edgetype(typeof(env′)) === E + @test spacetype(env′) === typeof(V) + @test spacetype(typeof(env′)) === typeof(V) + end + end +end