[Breaking] Redefine functions to take symmetries as argument

Project.toml

@@ -1,6 +1,6 @@
 name = "TNRKit"
 uuid = "7dfc3ef0-df9b-475f-b8e2-b91f34f5d84d"
-version = "0.4.0"
+version = "0.5.0"
 authors = ["Victor Vanthilt, Adwait Naravane, Atsushi Ueda and contributors"]
 
 [deps]

README.md

@@ -74,7 +74,7 @@ For example:
 ```julia
 using TNRKit, TensorKit
 
-T = classical_ising_symmetric(ising_βc) # partition function of classical Ising model at the critical point
+T = classical_ising(ising_βc) # partition function of classical Ising model at the critical point
 scheme = BTRG(T) # Bond-weighted TRG (excellent choice)
 data = run!(scheme, truncrank(16), maxiter(25)) # max bond-dimension of 16, for 25 iterations
 ```
@@ -106,14 +106,23 @@ To choose the verbosity level, simply use `run!(...; verbosity=n)`. The default
 
 ## Included Models on the square lattice
 TNRKit includes several common models out of the box.
-- Ising model: `classical_ising(β; h=0)` and `classical_ising_symmetric(β)`, which has a $\mathbb{Z}_2$ grading on each leg.
-- Potts model: `classical_potts(q, β)` and `classical_potts_symetric(q, β)`, which has a $\mathbb{Z}_q$ grading on each leg.
-- Six Vertex model: `sixvertex(scalartype, spacetype; a=1.0, b=1.0, c=1.0)`
-- Clock model: `classical_clock` and `classical_clock_symmetric`, which has a $\mathbb{Z}_q$ grading on each leg.
-- XY model: `classical_XY_U1_symmetric` and `classical_XY_O2_symmetric`
-- Real $\phi^4$ model: `phi4_real` and  `phi4_real_Z2`, which has a $\mathbb{Z}_2$ grading on each leg.
-- Complex $\phi^4$ model: `phi4_complex`,  `phi4_complex_U1`, which has a $U(1)$ grading on each leg and `phi4_complex_Z2Z2`, which has a $\mathbb{Z}_2 \times \mathbb{Z}_2$ grading on each leg.
+- Ising model in 2D: `classical_ising(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Ising model in 2D with impurities: `classical_ising_impurity(β; h=0)`.
+- Ising model in 3D: `classical_ising_3D(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Potts model in 2D: `classical_potts(S, q, β)`, where `S` can be `Trivial` or `ZNIrrep{q}` to specify the symmetry.
+- Potts model in 2D with impurities: `classical_potts_impurity(q, β)`.
+- Six Vertex model: `sixvertex(S, elt; a=1.0, b=1.0, c=1.0)` where `S` can be `Trivial`, `U1Irrep` or `CU1Irrep` to specify the symmetry and `elt` can be any number type (default is `Float64`).
+- Clock model: `classical_clock(S, q, β)` where `S` can be `Trivial` or `ZNIrrep{q}` to specify the symmetry.
+- XY model in 2D: `classical_XY(S, β, charge_trunc)` where `S` can be `U1Irrep` or `CU1Irrep` to specify the symmetry.
+- Real $\phi^4$ model: `phi4_real(S, K, μ0, λ, h)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Real $\phi^4$ model with impurities: `phi4_real_imp1(S, K, μ0, λ, h)` and `phi4_real_imp2(S, K, μ0, λ, h)` where `S` can be `Trivial`.
+- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep` or `U1Irrep` to specify the symmetry.
+- Gross-Neveu model: `gross_neveu_start(S, μ, m, g)` where `S` can be `FermionParity` to specify the symmetry.
 
 ## Included Models on the triangular lattice
 TNRKit includes several common models out of the box.
-- Ising model: `classical_ising_triangular` and `classical_ising_triangular_symmetric`, which has a $ℤ_2$ grading on each leg.
+- Ising model: `classical_ising_triangular(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+
+## Included Models on the honeycomb lattice
+TNRKit includes several common models out of the box.
+- Ising model: `classical_ising_honeycomb(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.

docs/src/finalizers.md

@@ -60,17 +60,17 @@ TNRKit exports the following pre-built `Finalizer` instances:
 using TNRKit
 
 # Default finalization (simple norm)
-T = classical_ising_symmetric(ising_βc)
+T = classical_ising(ising_βc)
 scheme = TRG(T)
 data = run!(scheme, truncrank(16), maxiter(25))
 
 # Use the two-by-two normalizer (more stable)
-T = classical_ising_symmetric(ising_βc)
+T = classical_ising(ising_βc)
 scheme = TRG(T)
 data = run!(scheme, truncrank(16), maxiter(25); finalizer=two_by_two_Finalizer)
 
 # Track ground state degeneracy throughout the simulation
-T = classical_ising_symmetric(ising_βc)
+T = classical_ising(ising_βc)
 scheme = TRG(T)
 gsd_data = run!(scheme, truncrank(16), maxiter(25); finalizer=GSDegeneracy_Finalizer)
 ```

docs/src/index.md

@@ -52,7 +52,7 @@ For example:
 ```julia
 using TNRKit, TensorKit
 
-T = classical_ising_symmetric(ising_βc) # partition function of classical Ising model at the critical point
+T = classical_ising(ising_βc) # partition function of classical Ising model at the critical point
 scheme = BTRG(T) # Bond-weighted TRG (excellent choice)
 data = run!(scheme, truncrank(16), maxiter(25)) # max bond-dimension of 16, for 25 iterations
 ```
@@ -80,15 +80,26 @@ to choose the verbosity level, simply use `run!(...; verbosity=n)`. The default
 
 ## Included Models on the square lattice
 TNRKit includes several common models out of the box.
-- Ising model: [`classical_ising`](@ref) and [`classical_ising_symmetric`](@ref), which has a Z2 grading on each leg.
-- Potts model: [`classical_potts`](@ref) and [`classical_potts_symmetric`](@ref), which has a Zq grading on each leg.
-- Six Vertex model: [`sixvertex`](@ref)
-- Clock model: [`classical_clock`](@ref) and [`classical_clock_symmetric`](@ref) which has a Zq grading on each leg.
-- XY model: [`classical_XY_U1_symmetric`](@ref) and [`classical_XY_O2_symmetric`](@ref).
+- Ising model in 2D: `classical_ising(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Ising model in 2D with impurities: `classical_ising_impurity(β; h=0)`.
+- Ising model in 3D: `classical_ising_3D(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Potts model in 2D: `classical_potts(S, q, β)`, where `S` can be `Trivial` or `ZNIrrep{q}` to specify the symmetry.
+- Potts model in 2D with impurities: `classical_potts_impurity(q, β)`.
+- Six Vertex model: `sixvertex(S, elt; a=1.0, b=1.0, c=1.0)` where `S` can be `Trivial`, `U1Irrep` or `CU1Irrep` to specify the symmetry and `elt` can be any number type (default is `Float64`).
+- Clock model: `classical_clock(S, q, β)` where `S` can be `Trivial` or `ZNIrrep{q}` to specify the symmetry.
+- XY model in 2D: `classical_XY(S, β, charge_trunc)` where `S` can be `U1Irrep` or `CU1Irrep` to specify the symmetry.
+- Real $\phi^4$ model: `phi4_real(S, K, μ0, λ, h)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+- Real $\phi^4$ model with impurities: `phi4_real_imp1(S, K, μ0, λ, h)` and `phi4_real_imp2(S, K, μ0, λ, h)` where `S` can be `Trivial`.
+- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep` or `U1Irrep` to specify the symmetry.
+- Gross-Neveu model: `gross_neveu_start(S, μ, m, g)` where `S` can be `FermionParity` to specify the symmetry.
 
 ## Included Models on the triangular lattice
 TNRKit includes several common models out of the box.
-- Ising model: [`classical_ising_triangular`](@ref) and [`classical_ising_triangular_symmetric`](@ref), which has a Z2 grading on each leg.
+- Ising model: `classical_ising_triangular(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
+
+## Included Models on the honeycomb lattice
+TNRKit includes several common models out of the box.
+- Ising model: `classical_ising_honeycomb(S, β; h=0)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry.
 
 If you want to implement your own model you must respect the leg-convention assumed by all TNRKit schemes.
 

src/TNRKit.jl

@@ -84,12 +84,10 @@ export run!
 include("models/ising.jl")
 include("models/ising_triangular.jl")
 include("models/ising_honeycomb.jl")
-export classical_ising, classical_ising_symmetric, ising_βc, f_onsager, ising_cft_exact,
-    ising_βc_3D, classical_ising_symmetric_3D, classical_ising_3D, classical_ising_impurity,
-    classical_ising_triangular, classical_ising_triangular_symmetric,
-    ising_βc_triangular, f_onsager_triangular,
-    classical_ising_honeycomb, classical_ising_honeycomb_symmetric,
-    ising_βc_honeycomb, f_onsager_honeycomb, XY_βc
+export classical_ising, ising_βc, f_onsager, ising_cft_exact,
+    ising_βc_3D, classical_ising_3D, classical_ising_impurity,
+    classical_ising_triangular, ising_βc_triangular, f_onsager_triangular,
+    classical_ising_honeycomb, ising_βc_honeycomb, f_onsager_honeycomb
 
 include("models/gross-neveu.jl")
 export gross_neveu_start
@@ -98,22 +96,19 @@ include("models/sixvertex.jl")
 export sixvertex
 
 include("models/potts.jl")
-export classical_potts, classical_potts_symmetric, potts_βc, classical_potts_impurity
+export classical_potts, potts_βc, classical_potts_impurity
 
 include("models/clock.jl")
-export classical_clock, classical_clock_symmetric
+export classical_clock
 
 include("models/XY.jl")
-export classical_XY_U1_symmetric
-export classical_XY_O2_symmetric
+export classical_XY, XY_βc
 
 include("models/phi4_real.jl")
 export phi4_real, phi4_real_imp1, phi4_real_imp2
-export phi4_real_Z2
 
 include("models/phi4_complex.jl")
 export phi4_complex, phi4_complex_impϕ, phi4_complex_impϕdag, phi4_complex_impϕabs, phi4_complex_impϕ2, phi4_complex_all
-export phi4_complex_U1, phi4_complex_Z2Z2
 
 # utility functions
 include("utility/free_energy.jl")

src/models/XY.jl

@@ -14,75 +14,64 @@ end
 const XY_βc = 1.1199 # This is an approximation!
 
 """
-$(SIGNATURES)
+    classical_XY(charge_trunc::Int; kwargs...)
+    classical_XY(beta::Float64, charge_trunc::Int; kwargs...)
+    classical_XY(::Type{U1Irrep}, beta::Float64, charge_trunc::Int; T::Type{<:Number} = Float64)
+    classical_XY(::Type{CU1Irrep}, beta::Float64, charge_trunc::Int; T::Type{<:Number} = Float64)
 
 Constructs the partition function tensor for a symmetric 2D square lattice
-for the classical XY model with U(1) symmetry, using inverse temperature `beta`
+for the classical XY model using inverse temperature `beta`
 and charge truncation `charge_trunc`.
+Compatible with U(1) symmetry or CU(1) = O(2) symmetry on each of its spaces.
+Defaults to CU(1) symmetry if the symmetry type is not provided.
 
 ### Examples
 ```julia
-    classical_XY_U1_symmetric(0.9, 6)
+    classical_XY(U1Irrep, 0.9, 6)
+    classical_XY(CU1Irrep, 0.9, 4)
 ```
 
 ### References
 * [Yu et. al. 10.1103/PhysRevE.89.013308 (2014)](@cite Yu_2014)
-
-See also: [`classical_XY_O2_symmetric`](@ref).
 """
-function classical_XY_U1_symmetric(beta::Float64, charge_trunc::Int)
+function classical_XY(beta::Float64, charge_trunc::Int; kwargs...)
+    return classical_XY(CU1Irrep, beta, charge_trunc; kwargs...)
+end
+classical_XY(charge_trunc::Int; kwargs...) = classical_XY(XY_βc, charge_trunc; kwargs...)
+classical_XY(::Type{U1Irrep}, charge_trunc::Int; kwargs...) = classical_XY(U1Irrep, XY_βc, charge_trunc; kwargs...)
+function classical_XY(::Type{U1Irrep}, beta::Float64, charge_trunc::Int; T::Type{<:Number} = Float64)
     FunU1 = U1Space(map(x -> (x => 1), (-charge_trunc):charge_trunc))
 
-    m = ones(Float64, FunU1 ← FunU1 ⊗ FunU1)
-
-    bond = zeros(Float64, FunU1 ← FunU1)
+    m = ones(T, FunU1 ← FunU1 ⊗ FunU1)
 
-    for sector in fusiontrees(bond)
-        charge = sector[1].uncoupled[1].charge
-        bond[sector...] .= besseli(charge, beta)
+    bond = zeros(T, FunU1 ← FunU1)
+    for (s, f) in fusiontrees(bond)
+        charge = s.uncoupled[1].charge
+        bond[s, f] .= besseli(charge, beta)
     end
 
     return algebraic_initialization(m, bond)
 end
-
-"""
-$(SIGNATURES)
-
-Constructs the partition function tensor for a symmetric 2D square lattice
-for the classical XY model with O(2) symmetry, using inverse temperature `beta`
-and charge truncation `charge_trunc`.
-
-### Examples
-```julia
-    classical_XY_O2_symmetric(0.9, 6)
-```
-
-### References
-* [Yu et. al. 10.1103/PhysRevE.89.013308 (2014)](@cite Yu_2014)
-
-See also: [`classical_XY_U1_symmetric`](@ref).
-"""
-function classical_XY_O2_symmetric(beta::Float64, charge_trunc::Int)
+function classical_XY(::Type{CU1Irrep}, beta::Float64, charge_trunc::Int; T::Type{<:Number} = Float64)
     FunU1_0 = CU1Space((0, 0) => 1)
     FunU1_1 = CU1Space(((i, 2) => 1 for i in 1:charge_trunc))
     FunU1 = FunU1_0 ⊕ FunU1_1
 
-    m = zeros(Float64, FunU1 ← FunU1 ⊗ FunU1)
-
+    m = zeros(T, FunU1 ← FunU1 ⊗ FunU1)
     for (to, from) in fusiontrees(m)
         left, right = from.uncoupled
-        if (left == right) && left != CU1Irrep(0, 0) && from.coupled == CU1Irrep(0, 0)
+        if (left == right) && !isunit(left) && isunit(from.coupled)
             m[to, from] .= sqrt(2)
         else
             m[to, from] .= 1
         end
     end
 
-    bond = zeros(Float64, FunU1 ← FunU1)
+    bond = zeros(T, FunU1 ← FunU1)
 
-    for sector in fusiontrees(bond)
-        charge = sector[1].uncoupled[1].j
-        bond[sector...] .= besseli(charge, beta)
+    for (s, f) in fusiontrees(bond)
+        charge = s.uncoupled[1].j
+        bond[s, f] .= besseli(charge, beta)
     end
 
     return algebraic_initialization(m, bond)

src/models/clock.jl

@@ -1,48 +1,48 @@
-"""
-$(SIGNATURES)
-
-Constructs the partition function tensor for the classical clock model with `q` states
-and a given inverse temperature `β`.
-"""
-function classical_clock(q::Int, β::Float64)
+function clock_tensor(q::Int, β::Real; T::Type{<:Number} = Float64)
     V = ℂ^q
-    A_clock = zeros(Float64, V ⊗ V ← V ⊗ V)
+    A_clock = zeros(T, V ⊗ V ← V ⊗ V)
     clock(i, j) = -cos(2π / q * (i - j))
 
-    for i in 1:q
-        for j in 1:q
-            for k in 1:q
-                for l in 1:q
-                    E = clock(i, j) + clock(j, l) + clock(l, k) + clock(k, i)
-                    A_clock[i, j, k, l] = exp(-β * E)
-                end
-            end
-        end
+    for i in 1:q, j in 1:q, k in 1:q, l in 1:q
+        E = clock(i, j) + clock(j, l) + clock(l, k) + clock(k, i)
+        A_clock[i, j, k, l] = exp(-β * E)
     end
+
     return A_clock
 end
 
 """
-$(SIGNATURES)
+    classical_clock(q::Int, β::Real; kwargs...)
+    classical_clock(::Type{Trivial}, q::Int, β::Real; T::Type{<:Number} = Float64)
+    classical_clock(::Type{ZNIrrep{N}}, q::Int, β::Real; T::Type{<:Number} = Float64) where {N}
 
 Constructs the partition function tensor for the classical clock model with `q` states
 and a given inverse temperature `β`.
 
-This tensor has explicit ℤq symmetry on each of it spaces.
+Compatible with no symmetry or with explicit ℤq symmetry on each of its spaces.
+Defaults to ℤq symmetry if the symmetry type is not provided.
 """
-function classical_clock_symmetric(q::Int, β::Float64)
-    A = classical_clock(q, β)
+function classical_clock(q::Int, β::Real; kwargs...)
+    return classical_clock(ZNIrrep{q}, q, β; kwargs...)
+end
+function classical_clock(::Type{Trivial}, q::Int, β::Real; kwargs...)
+    return clock_tensor(q, β; kwargs...)
+end
+function classical_clock(::Type{ZNIrrep{N}}, q::Int, β::Real; T::Type{<:Number} = Float64) where {N}
+    @assert N == q "number of irreps must match the number of states"
+    A = classical_clock(Trivial, q, β; T = T)
 
     # Construct the Fourier matrix for the clock model
-    U = zeros(ComplexF64, q, q)
+    Udat = zeros(ComplexF64, q, q)
     for i in 0:(q - 1)
         for j in 0:(q - 1)
-            U[i + 1, j + 1] = exp(2im * π / q * i * j) / sqrt(q)
+            Udat[i + 1, j + 1] = cispi(2 / q * i * j) / sqrt(q)
         end
     end
-    U = TensorMap(U, ℂ^q ← ℂ^q)
+    U = TensorMap(Udat, ℂ^q ← ℂ^q)
 
     @tensor Anew[-1 -2;-3 -4] := A[1 2; 3 4] * U[4; -4] * conj(U[1; -1]) * U[3; -3] * conj(U[2; -2])
     V = ZNSpace{q}(i => 1 for i in 0:(q - 1))
-    return real(TensorMap(convert(Array, Anew), V ⊗ V ← V ⊗ V))
+    t = TensorMap(convert(Array, Anew), V ⊗ V ← V ⊗ V)
+    return T <: Real ? real(t) : t
 end

src/models/gross-neveu.jl

@@ -7,17 +7,17 @@ function P_tensor()
     return P
 end
 
-function gross_neveu_8_leg_tensor(μ::Number, m::Number, g::Number)
+function gross_neveu_8_leg_tensor(μ::Number, m::Number, g::Number; T::Type{<:Complex} = ComplexF64)
     # Manually defining the A and A_bar tensors
-    A = zeros(ComplexF64, 2, 2, 2, 2)
+    A = zeros(T, 2, 2, 2, 2)
     A[2, 2, 1, 1] = 1 + 1im
     A[1, 1, 2, 2] = -1 - 1im
     A[2, 1, 1, 2] = 1 - 1im
     A[2, 1, 2, 1] = 2
     A[1, 2, 1, 2] = -2im
     A[1, 2, 2, 1] = 1 - 1im
 
-    A_bar = zeros(ComplexF64, 2, 2, 2, 2)
+    A_bar = zeros(T, 2, 2, 2, 2)
     A_bar[2, 2, 1, 1] = -1 + 1im
     A_bar[1, 1, 2, 2] = 1 - 1im
     A_bar[2, 1, 1, 2] = -1 - 1im
@@ -28,12 +28,12 @@ function gross_neveu_8_leg_tensor(μ::Number, m::Number, g::Number)
     # Utility Kronecker delta function
     δ(x, y) = ==(x, y)
 
-    T = zeros(ComplexF64, 2, 2, 2, 2, 2, 2, 2, 2)
+    t = zeros(T, 2, 2, 2, 2, 2, 2, 2, 2)
     P = P_tensor()
     V = Vect[FermionParity](0 => 1, 1 => 1)
     for (pi1, pj1, pi2, pj2, i1, j1, i2, j2) in Iterators.product([0:1 for _ in 1:8]...)
         p = P[pi1 + 1, pj1 + 1, pi2 + 1, pj2 + 1, i1 + 1, j1 + 1, i2 + 1, j2 + 1]
-        T[pi1 + 1, pj1 + 1, pi2 + 1, pj2 + 1, i2 + 1, j2 + 1, i1 + 1, j1 + 1] =
+        t[pi1 + 1, pj1 + 1, pi2 + 1, pj2 + 1, i2 + 1, j2 + 1, i1 + 1, j1 + 1] =
             ((-1)^p) * exp(0.5 * μ * (i2 - j2 + pi2 - pj2)) * ((1 / sqrt(2))^(i1 + i2 + j1 + j2 + pi1 + pi2 + pj1 + pj2)) *
             (
             ((m + 2)^2 + 2 * g^2) * δ(i1 + i2 + pj1 + pj2, 0) * δ(j1 + j2 + pi1 + pi2, 0) -
@@ -43,20 +43,24 @@ function gross_neveu_8_leg_tensor(μ::Number, m::Number, g::Number)
         )
 
     end
-    return TensorMap(T, V ⊗ V ⊗ V ⊗ V ← V ⊗ V ⊗ V ⊗ V)
+    return TensorMap(t, V ⊗ V ⊗ V ⊗ V ← V ⊗ V ⊗ V ⊗ V)
 end
 
 
 """
-$(SIGNATURES)
+    gross_neveu_start([::Type{FermionParity}], μ::Number, m::Number, g::Number; T::Type{<:Complex} = ComplexF64)
 
 Constructs the partition function tensor for the Gross-Neveu model with given parameters `μ`, `m`, and `g`.
+Compatible with explicit fermion parity symmetry on each of its spaces.
 
 ### References
 * [Akiyama et. al. J. Phys.: Condens. Matter 36 (2024) 343002](@cite akiyama2024)
 """
-function gross_neveu_start(μ::Number, m::Number, g::Number)
-    T_unfused = gross_neveu_8_leg_tensor(μ, m, g)
+function gross_neveu_start(μ::Number, m::Number, g::Number; kwargs...)
+    return gross_neveu_start(FermionParity, μ, m, g; kwargs...)
+end
+function gross_neveu_start(::Type{FermionParity}, μ::Number, m::Number, g::Number; T::Type{<:Complex} = ComplexF64)
+    T_unfused = gross_neveu_8_leg_tensor(μ, m, g; T = T)
     V = Vect[FermionParity](0 => 1, 1 => 1)
     U = isometry(fuse(V, V), V ⊗ V)
     Udg = adjoint(U)