Add quantum Ashkin-Teller model

docs/src/man/models.md

@@ -15,6 +15,7 @@ heisenberg_XXZ
 heisenberg_XYZ
 bilinear_biquadratic_model
 tj_model
+ashkin_teller
 hubbard_model
 bose_hubbard_model
 quantum_chemistry_hamiltonian

src/MPSKitModels.jl

@@ -45,6 +45,7 @@ export heisenberg_XXX, heisenberg_XXZ, heisenberg_XYZ
 export bilinear_biquadratic_model
 export hubbard_model, bose_hubbard_model
 export tj_model
+export ashkin_teller
 export quantum_chemistry_hamiltonian
 
 export classical_ising

src/models/hamiltonians.jl

@@ -461,4 +461,79 @@ function tj_model(
     end
 end
 
+#===========================================================================================
+    Quantum Ashkin-Teller model
+===========================================================================================#
+"""
+    ashkin_teller([T::Type{<:Number} = ComplexF64],
+                  [S = ProductSector{Tuple{Z2Irrep, Z2Irrep}}],
+                  [lattice::AbstractLattice = InfiniteChain(1)];
+                  h = 1.0, J = 1.0, λ = 1.0)
+
+MPO for the hamiltonian of the quantum Ashkin-Teller model. The model is
+defined on a chain of two qubits per site. Writing Pauli operators on each of these qubits
+as ``\\sigma `` and ``\\tau ``, the Hamiltonian reads:
+```math
+H = -h \\sum_i\\bigg(\\sigma_i^x + \\tau_i^x + \\lambda \\sigma_i^x \\tau_i^x \\bigg)
+-J \\sum_{\\langle i,j \\rangle} \\bigg( \\sigma_i^z \\sigma_j^z + \\tau_j^z \\tau_j^z
++\\lambda \\sigma_i^z \\sigma_j^z \\tau_i^z \\tau_j^z \\bigg).
+
+```
+Currently the Hamiltonian is only supported with ```Z2Irrep ⊠  Z2Irrep``` symmetry.
+"""
+function ashkin_teller end
+function ashkin_teller(lattice::AbstractLattice; kwargs...)
+    return ashkin_teller(ComplexF64, ProductSector{Tuple{Z2Irrep, Z2Irrep}}, lattice; kwargs...)
+end
+function ashkin_teller(T::Type{<:Number}, lattice::AbstractLattice; kwargs...)
+    return ashkin_teller(T, ProductSector{Tuple{Z2Irrep, Z2Irrep}}, lattice; kwargs...)
+end
+function ashkin_teller(
+        T::Type{<:Number} = ComplexF64,
+        S = ProductSector{Tuple{Z2Irrep, Z2Irrep}},
+        lattice::AbstractLattice = InfiniteChain(1);
+        h = 1.0, J = 1.0, λ = 1.0
+    )
+
+    S == ProductSector{Tuple{Z2Irrep, Z2Irrep}} || error("Only implemented with ℤ₂×ℤ₂ symmetry")
+
+    V = Vect[S](c => 1 for c in values(S))
+
+    # Single site operators
+    XI = ones(T, V ← V)
+    block(XI, S(1, 0)) .*= -1
+    block(XI, S(1, 1)) .*= -1
+    IX = ones(T, V ← V)
+    block(IX, S(0, 1)) .*= -1
+    block(IX, S(1, 1)) .*= -1
+    XX = ones(T, V ← V)
+    block(XX, S(0, 1)) .*= -1
+    block(XX, S(1, 0)) .*= -1
+
+    # Nearest-neighbour terms
+    ZIZI = zeros(T, V ⊗ V ← V ⊗ V)
+    IZIZ = zeros(T, V ⊗ V ← V ⊗ V)
+    ZZZZ = zeros(T, V ⊗ V ← V ⊗ V)
+    for (s, f) in fusiontrees(ZIZI)
+        if s.uncoupled == map(x -> flip_charge(x, (1, 0)), f.uncoupled)
+            ZIZI[s, f] .= 1
+        end
+        if s.uncoupled == map(x -> flip_charge(x, (0, 1)), f.uncoupled)
+            IZIZ[s, f] .= 1
+        end
+        if s.uncoupled == map(x -> flip_charge(x, (1, 1)), f.uncoupled)
+            ZZZZ[s, f] .= 1
+        end
+    end
+
+    return @mpoham begin
+        sum(vertices(lattice)) do i
+            return -h * (XI{i} + IX{i} + λ * XX{i})
+        end +
+            sum(nearest_neighbours(lattice)) do (i, j)
+            -J * (ZIZI{i, j} + IZIZ{i, j} + λ * ZZZZ{i, j})
+        end
+    end
+end
+
 # TODO: add (hardcore) bosonic t-J model (https://arxiv.org/abs/2409.15424)

src/utility.jl

@@ -57,3 +57,15 @@ function split_twosite(O::AbstractTensorMap{<:Any, <:Any, 2, 2})
     @plansor R[a p'; p] := sqrtS[a; 1] * V[1; p p']
     return L, R
 end
+
+#===========================================================================================
+    Other utility
+===========================================================================================#
+
+"""
+    flip_charge(charge, flip)
+
+    Take the product of charge with flip.
+    Works only for abelian symmetries.
+"""
+flip_charge(charge, flip) = only(charge ⊗ sectortype(charge)(flip))

test/ashkin_teller.jl

@@ -0,0 +1,34 @@
+using Test
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+# Setup
+
+gammas = [
+    pi / 2, pi / 3, pi / 4, pi / 6, 0,
+]
+# Exact GS energy density via mapping to spin 1/2 XXZ and Bethe ansatz
+# https://www.sciencedirect.com/science/article/pii/0003491688900152
+E0s = [
+    -8 / pi, -12 / 4, -4 * (sqrt(2) / pi + 1 / (2 * sqrt(2))),
+    -4 * (1 / pi + 11 / (12 * sqrt(3))), 2 - 8 * log(2),
+]
+alg = VUMPS(; maxiter = 100, verbosity = 0)
+
+S = Z2Irrep ⊠ Z2Irrep
+V = Vect[S](c => 1 for c in values(S))
+W = Vect[S](c => 6 for c in values(S))
+
+# Test
+
+@testset "Ashkin-Teller" for (gamma, E0) in zip(gammas, E0s)
+
+    H = ashkin_teller(h = 1, J = 1, λ = cos(gamma))
+    Ψ = InfiniteMPS(V, W)
+
+    @test imag(expectation_value(Ψ, H)) ≈ 0 atol = 1.0e-12
+
+    Ψ0, _ = find_groundstate(Ψ, H, alg)
+    @test real(expectation_value(Ψ0, H)) ≈ E0 atol = 1.0e-3
+end