Add potts model

src/MPSKitModels.jl

@@ -37,6 +37,7 @@ export e⁺, e⁻, e⁺⁺, e⁻⁻, e⁺⁻, e⁻⁺
 
 export transverse_field_ising
 export kitaev_model
+export quantum_potts
 export heisenberg_XXX, heisenberg_XXZ, heisenberg_XYZ
 export bilinear_biquadratic_model
 export hubbard_model, bose_hubbard_model

src/models/hamiltonians.jl

@@ -1,7 +1,6 @@
 #===========================================================================================
     Ising model
 ===========================================================================================#
-ising_kwargs = (; J=1.0, g=1.0)
 
 """
     transverse_field_ising([elt::Type{<:Number}], [symmetry::Type{<:Sector}],
@@ -234,6 +233,44 @@ function bilinear_biquadratic_model(elt::Type{<:Number}=ComplexF64,
     end
 end
 
+#===========================================================================================
+    Potts models
+===========================================================================================#
+"""
+    quantum_potts([elt::Type{<:Number}], [symmetry::Type{<:Sector}],
+                  [lattice::AbstractLattice]; q=3, J=1.0, g=1.0)
+
+MPO for the hamiltonian of the quantum Potts model, as defined by
+```math
+H = - J \\sum_{\\langle i,j \\rangle} Z_i^\\dagger Z_j + Z_i Z_j^\\dagger
+- g \\sum_i (X_i + X_i^\\dagger)
+```
+where the operators ``Z`` and ``X`` are the ``q``-rotation operators satisfying
+``Z^q = X^q = 1`` and ``ZX = \\omega XZ`` where ``\\omega = e^{2πi/q}``.
+"""
+function quantum_potts end
+function quantum_potts(lattice::AbstractLattice; kwargs...)
+    return quantum_potts(ComplexF64, Trivial, lattice; kwargs...)
+end
+function quantum_potts(symmetry::Type{<:Sector},
+                       lattice::AbstractLattice=InfiniteChain(1); kwargs...)
+    return quantum_potts(ComplexF64, symmetry, lattice; kwargs...)
+end
+function quantum_potts(elt::Type{<:Number}, lattice::AbstractLattice;
+                       kwargs...)
+    return quantum_potts(elt, Trivial, lattice; kwargs...)
+end
+function quantum_potts(elt::Type{<:Number}=ComplexF64,
+                       symmetry::Type{<:Sector}=Trivial,
+                       lattice::AbstractLattice=InfiniteChain(1);
+                       q=3, J=1.0, g=1.0)
+    return @mpoham sum(nearest_neighbours(lattice)) do (i, j)
+        return -J * potts_exchange(elt, symmetry; q){i,j}
+    end - sum(vertices(lattice)) do i
+          return g * potts_field(elt, symmetry; q){i}
+          end
+end
+
 #===========================================================================================
     Hubbard models
 ===========================================================================================#

src/operators/spinoperators.jl

@@ -444,3 +444,64 @@ const SS = S_exchange
 
 """Pauli exchange operator."""
 σσ(args...; kwargs...) = 4 * S_exchange(args...; kwargs...)
+
+"""
+    potts_exchange([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; q=3)
+
+The Potts exchange operator ``Z ⊗ Z' + Z' ⊗ Z``, where ``Z^q = 1``.
+"""
+function potts_exchange end
+potts_exchange(; kwargs...) = potts_exchange(ComplexF64, Trivial; kwargs...)
+potts_exchange(elt::Type{<:Number}; kwargs...) = potts_exchange(elt, Trivial; kwargs...)
+function potts_exchange(symmetry::Type{<:Sector}; kwargs...)
+    return potts_exchange(ComplexF64, symmetry; kwargs...)
+end
+
+function potts_exchange(elt::Type{<:Number}, ::Type{Trivial}; q=3)
+    pspace = ComplexSpace(q)
+    Z = TensorMap(zeros, elt, pspace ← pspace)
+    for i in 1:q
+        Z[i, i] = cis(2π * (i - 1) / q)
+    end
+    return Z ⊗ Z' + Z' ⊗ Z
+end
+function potts_exchange(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
+    @assert q == Q "q must match the irrep charge"
+    pspace = Vect[ZNIrrep{q}](i => 1 for i in 0:(q - 1))
+    aspace = Vect[ZNIrrep{q}](1 => 1, -1 => 1)
+    Z_left = TensorMap(ones, elt, pspace ← pspace ⊗ aspace)
+    Z_right = TensorMap(ones, elt, aspace ⊗ pspace ← pspace)
+    return contract_twosite(Z_left, Z_right)
+end
+
+"""
+    potts_field([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; q=3) 
+
+The Potts field operator ``X + X'``, where ``X^q = 1``.
+"""
+function potts_field end
+potts_field(; kwargs...) = potts_field(ComplexF64, Trivial; kwargs...)
+potts_field(elt::Type{<:Number}; kwargs...) = potts_field(elt, Trivial; kwargs...)
+function potts_field(symmetry::Type{<:Sector}; kwargs...)
+    return potts_field(ComplexF64, symmetry; kwargs...)
+end
+
+function potts_field(elt::Type{<:Number}, ::Type{Trivial}; q=3)
+    pspace = ComplexSpace(q)
+    X = TensorMap(zeros, elt, pspace ← pspace)
+    for i in 1:q
+        X[mod1(i - 1, q), i] = one(elt)
+    end
+    return X + X'
+end
+# TODO: generalize to arbitrary q
+function potts_field(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
+    @assert q == Q "q must match the irrep charge"
+    @assert q == 3 "only q = 3 is implemented"
+    pspace = Vect[ZNIrrep{q}](i => 1 for i in 0:(q - 1))
+    X = TensorMap(zeros, elt, pspace ← pspace)
+    for (c, b) in blocks(X)
+        b .= isone(c) ? 2 : -1
+    end
+    return X
+end

test/potts.jl

@@ -0,0 +1,21 @@
+using MPSKit
+using TensorKit
+
+alg = VUMPS(; maxiter=25, verbosity=0)
+E₀ = -(4 / 3 + 2sqrt(3) / π)
+
+@testset "Trivial" begin
+    H = quantum_potts(; q=3)
+    ψ = InfiniteMPS(3, 32)
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end
+
+@testset "Z3Irrep" begin
+    H = quantum_potts(Z3Irrep; q=3)
+    ψ = InfiniteMPS(Rep[ℤ₃](i => 1 for i in 0:2), Rep[ℤ₃](i => 12 for i in 0:2))
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end