Compatibility for MPSKit 0.13 (#183)

* CompatHelper: bump compat for MPSKit to 0.13, (keep existing compat)

* Actually import `MPSKit.physicalspace`

* Refactor derivative contractions

* Add twosite derivative contraction

* Include PEPO contractions

* Formatter

* Drop compat for older MPSKit versions

* Fix typo

* Update VUMPS environments calling syntax (#167)

* `bose_hubbard_model` charge convention

* Update `add_physical_charge` implementation for `LocalOperator`

* Remove TODO

* Move Bose-Hubbard example

* And actually point to moved example

* Add test for twosite derivative operator

---------

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Lukas Devos <ldevos98@gmail.com>
Co-authored-by: Lander Burgelman <39218680+leburgel@users.noreply.github.com>
Co-authored-by: leburgel <lander.burgelman@ugent.be>

Author: 4 people

Project.toml

@@ -31,7 +31,7 @@ FiniteDifferences = "0.12"
 KrylovKit = "0.9.5"
 LinearAlgebra = "1"
 LoggingExtras = "1"
-MPSKit = "0.12.4"
+MPSKit = "0.13"
 MPSKitModels = "0.4"
 OhMyThreads = "0.7"
 OptimKit = "0.3, 0.4"

src/PEPSKit.jl

@@ -7,7 +7,7 @@ using VectorInterface
 using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations
 using ChainRulesCore, Zygote
 using LoggingExtras
-import MPSKit: leading_boundary, loginit!, logiter!, logfinish!, logcancel!
+import MPSKit: leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 using MPSKitModels
 using FiniteDifferences
 using OhMyThreads: tmap

src/algorithms/contractions/vumps_contractions.jl

@@ -134,52 +134,84 @@ function MPSKit.contract_mpo_expval(
     return environment_overlap(GL´, GR)
 end
 
+# PEPS Derivative contractions
+# ----------------------------
+# This is appropriating the MPSKit MPO derivative structures, which might not be the best
+# idea in the long run.
+
+const PEPS_C_Hamiltonian{S,N} = MPSKit.MPO_C_Hamiltonian{
+    <:GenericMPSTensor{S,N},<:GenericMPSTensor{S,N}
+} # this one is technically type-piracy
+PEPS_C_Hamiltonian(GL, GR) = MPSKit.MPODerivativeOperator(GL, (), GR)
+
+const PEPS_AC_Hamiltonian{S,N} = MPSKit.MPO_AC_Hamiltonian{
+    <:GenericMPSTensor{S,N},<:PEPSSandwich,<:GenericMPSTensor{S,N}
+}
+PEPS_AC_Hamiltonian(GL, O, GR) = MPSKit.MPODerivativeOperator(GL, (O,), GR)
+
+const PEPS_AC2_Hamiltonian{S,N} = MPSKit.MPO_AC2_Hamiltonian{
+    <:GenericMPSTensor{S,N},<:PEPSSandwich,<:PEPSSandwich,<:GenericMPSTensor{S,N}
+}
+PEPS_AC2_Hamiltonian(GL, O1, O2, GR) = MPSKit.MPODerivativeOperator(GL, (O1, O2), GR)
+
+# Constructors
 #
-# Derivative contractions
-#
-
-function MPSKit.∂C(
-    C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N}
-) where {S,N}
+function MPSKit.C_hamiltonian(site::Int, below, ::InfiniteTransferMatrix, above, envs)
+    GL = leftenv(envs, site + 1, below)
     GL = twistdual(GL, 1)
+    GR = rightenv(envs, site, below)
     GR = twistdual(GR, numind(GR))
-    return _∂C(C, GL, GR)
+    return PEPS_C_Hamiltonian(GL, GR)
 end
-@generated function _∂C(
-    C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N}
-) where {S,N}
-    C´_e = tensorexpr(:C´, -1, -2)
-    C_e = tensorexpr(:C, 1, 2)
-    GL_e = tensorexpr(:GL, (-1, (3:(N + 1))...), 1)
-    GR_e = tensorexpr(:GR, 2:(N + 1), -2)
-    return macroexpand(@__MODULE__, :(return @tensor $C´_e := $GL_e * $C_e * $GR_e))
+
+function MPSKit.AC_hamiltonian(
+    site::Int, below, operator::InfiniteTransferPEPS, above, envs
+)
+    GL = leftenv(envs, site, below)
+    GL = twistdual(GL, 1)
+    GR = rightenv(envs, site, below)
+    GR = twistdual(GR, numind(GR))
+    return PEPS_AC_Hamiltonian(GL, operator[site], GR)
 end
 
-function MPSKit.∂AC(
-    AC::GenericMPSTensor{S,N},
-    O::Union{PEPSSandwich,PEPOSandwich},
-    GL::GenericMPSTensor{S,N},
-    GR::GenericMPSTensor{S,N},
-) where {S,N}
+function MPSKit.AC2_hamiltonian(
+    site::Int, below, operator::InfiniteTransferPEPS, above, envs
+)
+    GL = leftenv(envs, site, below)
     GL = twistdual(GL, 1)
+    GR = rightenv(envs, site + 1, below)
     GR = twistdual(GR, numind(GR))
-    return _∂AC(AC, O, GL, GR)
+    return PEPS_AC2_Hamiltonian(GL, operator[site], operator[site + 1], GR)
 end
 
-## PEPS
+# Actions
+#
+@generated function (h::PEPS_C_Hamiltonian{S,N})(C::MPSBondTensor{S}) where {S,N}
+    C´_e = tensorexpr(:C´, -1, -2)
+    C_e = tensorexpr(:C, 1, 2)
+    GL_e = tensorexpr(:(h.leftenv), (-1, (3:(N + 1))...), 1)
+    GR_e = tensorexpr(:(h.rightenv), (2:(N + 1)...,), -2)
+    return macroexpand(@__MODULE__, :(return @tensor $C´_e := $GL_e * $C_e * $GR_e))
+end
 
-function _∂AC(
-    AC::GenericMPSTensor{S,3},
-    O::PEPSSandwich,
-    GL::GenericMPSTensor{S,3},
-    GR::GenericMPSTensor{S,3},
-) where {S}
+function (h::PEPS_AC_Hamiltonian{S,N})(AC::GenericMPSTensor{S,N}) where {S,N}
     return @autoopt @tensor AC′[χ_SW D_S_above D_S_below; χ_SE] :=
-        GL[χ_SW D_W_above D_W_below; χ_NW] *
+        h.leftenv[χ_SW D_W_above D_W_below; χ_NW] *
         AC[χ_NW D_N_above D_N_below; χ_NE] *
-        GR[χ_NE D_E_above D_E_below; χ_SE] *
-        ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
-        conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below])
+        h.rightenv[χ_NE D_E_above D_E_below; χ_SE] *
+        ket(h.operators[1])[d; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra(h.operators[1])[d; D_N_below D_E_below D_S_below D_W_below])
+end
+
+function (h::PEPS_AC2_Hamiltonian{S,3})(AC2::AbstractTensorMap{<:Any,S,3,3}) where {S}
+    return @autoopt @tensor AC2′[χ_SW D_S_above1 D_S_below1; χ_SE D_S_below2 D_S_above2] :=
+        h.leftenv[χ_SW D_W_above1 D_W_below1; χ_NW] *
+        AC2[χ_NW D_N_above1 D_N_below1; χ_NE D_N_below2 D_N_above2] *
+        h.rightenv[χ_NE D_E_above2 D_E_below2; χ_SE] *
+        ket(h.operators[1])[d1; D_N_above1 D_E_above1 D_S_above1 D_W_above1] *
+        conj(bra(h.operators[1])[d1; D_N_below1 D_E_below1 D_S_below1 D_W_below1]) *
+        ket(h.operators[2])[d2; D_N_above2 D_E_above2 D_S_above2 D_E_above1] *
+        conj(bra(h.operators[2])[d2; D_N_below2 D_E_below2 D_S_below2 D_E_below1])
 end
 
 # PEPS derivative
@@ -198,22 +230,32 @@ function ∂peps(
         conj(ĀC[χ_SW D_S_above D_S_below; χ_SE])
 end
 
-## PEPO
+# PEPO Derivative contractions
+# ----------------------------
+const PEPO_AC_Hamiltonian{S,N,H} = MPSKit.MPO_AC_Hamiltonian{
+    <:GenericMPSTensor{S,N},<:PEPOSandwich{H},<:GenericMPSTensor{S,N}
+}
+PEPO_AC_Hamiltonian(GL, O, GR) = MPSKit.MPODerivativeOperator(GL, (O,), GR)
+
+function MPSKit.AC_hamiltonian(
+    site::Int, below, operator::InfiniteTransferPEPO, above, envs
+)
+    GL = leftenv(envs, site, below)
+    GL = twistdual(GL, 1)
+    GR = rightenv(envs, site, below)
+    GR = twistdual(GR, numind(GR))
+    return PEPO_AC_Hamiltonian(GL, operator[site], GR)
+end
 
-@generated function _∂AC(
-    AC::GenericMPSTensor{S,N},
-    O::PEPOSandwich{H},
-    GL::GenericMPSTensor{S,N},
-    GR::GenericMPSTensor{S,N},
-) where {S,N,H}
+@generated function (h::PEPO_AC_Hamiltonian{S,N,H})(AC::GenericMPSTensor{S,N}) where {S,N,H}
     # sanity check
-    @assert H == N - 3
+    @assert H == N - 3 "Incompatible number of legs and layers"
 
     AC´_e = _pepo_edge_expr(:AC´, :SW, :SE, :S, H)
     AC_e = _pepo_edge_expr(:AC, :NW, :NE, :N, H)
-    GL_e = _pepo_edge_expr(:GL, :SW, :NW, :W, H)
-    GR_e = _pepo_edge_expr(:GR, :NE, :SE, :E, H)
-    ket_e, bra_e, pepo_es = _pepo_sandwich_expr(:O, H)
+    GL_e = _pepo_edge_expr(:(h.leftenv), :SW, :NW, :W, H)
+    GR_e = _pepo_edge_expr(:(h.rightenv), :NE, :SE, :E, H)
+    ket_e, bra_e, pepo_es = _pepo_sandwich_expr(:(h.operators[1]), H)
 
     rhs = Expr(:call, :*, AC_e, GL_e, GR_e, ket_e, Expr(:call, :conj, bra_e), pepo_es...)
 

src/environments/vumps_environments.jl

@@ -5,9 +5,9 @@ using MPSKit: InfiniteEnvironments
 
 function MPSKit.issamespace(
     env::InfiniteEnvironments,
-    above::InfiniteMPS,
-    operator::InfiniteTransferMatrix,
     below::InfiniteMPS,
+    operator::InfiniteTransferMatrix,
+    above::InfiniteMPS,
 )
     L = MPSKit.check_length(above, operator, below)
     for i in 1:L

src/operators/localoperator.jl

@@ -264,17 +264,19 @@ end
 """
     MPSKit.add_physical_charge(H::LocalOperator, charges::AbstractMatrix{<:Sector}) where {S}
 
-Change the spaces of a `LocalOperator` by fusing in an auxiliary charge on every site,
-according to a given matrix of 'auxiliary' physical charges.
+Change the spaces of a `LocalOperator` by fusing in an auxiliary charge into the domain of
+the operator on every site, according to a given matrix of 'auxiliary' physical charges.
 """
 function MPSKit.add_physical_charge(H::LocalOperator, charges::AbstractMatrix{<:Sector})
     size(H.lattice) == size(charges) ||
         throw(ArgumentError("Incompatible lattice and auxiliary charge sizes"))
     sectortype(H) === eltype(charges) ||
         throw(SectorMismatch("Incompatible lattice and auxiliary charge sizes"))
 
-    # make indexing periodic, for convenience
-    Paux = PeriodicArray(map(c -> Vect[typeof(c)](c => 1), charges))
+    # auxiliary spaces will be fused into codomain, so need to dualize the space to fuse
+    # the charge into the domain as desired
+    # also, make indexing periodic for convenience
+    Paux = PeriodicArray(map(c -> Vect[typeof(c)](c => 1)', charges))
 
     # new physical spaces
     Pspaces = map(fuse, H.lattice, Paux)

src/operators/models.jl

@@ -198,7 +198,7 @@ function MPSKitModels.bose_hubbard_model(
     elseif symmetry === U1Irrep
         isinteger(2n) ||
             throw(ArgumentError("`U₁` symmetry requires halfinteger particle number"))
-        H = MPSKit.add_physical_charge(H, fill(U1Irrep(-n), size(spaces)...)) # TODO: settle on convention here?
+        H = MPSKit.add_physical_charge(H, fill(U1Irrep(n), size(spaces)...))
     else
         throw(ArgumentError("symmetry not implemented"))
     end

test/boundarymps/vumps.jl

@@ -10,6 +10,7 @@ Random.seed!(29384293742893)
 const vumps_alg = VUMPS(;
     tol=1e-6, alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitian=false), verbosity=2
 )
+
 @testset "(1, 1) PEPS" begin
     psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2))
 
@@ -19,6 +20,11 @@ const vumps_alg = VUMPS(;
     mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     N = abs(sum(expectation_value(mps, T)))
 
+    mps2, = changebonds(mps, T, OptimalExpand(; trscheme=truncdim(30)))
+    mps2, env2, ϵ = leading_boundary(mps2, T, vumps_alg)
+    N2 = abs(sum(expectation_value(mps2, T)))
+    @test N ≈ N2 rtol = 1e-2
+
     ctm, = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi)
     N´ = abs(norm(psi, ctm))
 

test/bose_hubbard.jl test/examples/bose_hubbard.jltest/bose_hubbard.jl renamed to test/examples/bose_hubbard.jl

@@ -76,7 +76,7 @@ end
             include("examples/pwave.jl")
         end
         @time @safetestset "U1-symmetric Bose-Hubbard model" begin
-            include("bose_hubbard.jl")
+            include("examples/bose_hubbard.jl")
         end
     end
 end