Merge branch 'pb-lorentzian-broadening' of github.com:quantumghent/PEPSKit.jl into pb-lorentzian-broadening

Author: pbrehmer

examples/boundary_mps.jl

@@ -23,7 +23,7 @@ T = InfiniteTransferPEPS(peps, 1, 1)
 # encodes the norm of the state
 
 # Fist we build an initial guess for the boundary MPS, choosing a bond dimension of 20
-mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)])
+mps = initialize_mps(T, [ComplexSpace(20)])
 
 # We then find the leading boundary MPS fixed point using the VUMPS algorithm
 mps, env, ϵ = leading_boundary(mps, T, VUMPS())
@@ -51,7 +51,7 @@ N´ = abs(norm(peps, ctm))
 peps2 = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2))
 T2 = PEPSKit.MultilineTransferPEPS(peps2, 1)
 
-mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2))
+mps2 = initialize_mps(T2, fill(ComplexSpace(20), 2, 2))
 mps2, env2, ϵ = leading_boundary(mps2, T2, VUMPS())
 N2 = abs(prod(expectation_value(mps2, T2)))
 
@@ -76,7 +76,7 @@ N2´ = abs(norm(peps2, ctm2))
 pepo = ising_pepo(1)
 T3 = InfiniteTransferPEPO(peps, pepo, 1, 1)
 
-mps3 = PEPSKit.initializeMPS(T3, [ComplexSpace(20)])
+mps3 = initialize_mps(T3, [ComplexSpace(20)])
 mps3, env3, ϵ = leading_boundary(mps3, T3, VUMPS())
 @show N3 = abs(prod(expectation_value(mps3, T3)))
 

src/PEPSKit.jl

@@ -89,7 +89,7 @@ export InfinitePartitionFunction
 export InfinitePEPS, InfiniteTransferPEPS
 export SUWeight, InfiniteWeightPEPS
 export InfinitePEPO, InfiniteTransferPEPO
-export initializeMPS, initializePEPS
+export initialize_mps, initializePEPS
 export ReflectDepth, ReflectWidth, Rotate, RotateReflect
 export symmetrize!, symmetrize_retract_and_finalize!
 export showtypeofgrad

src/algorithms/contractions/vumps_contractions.jl

@@ -4,9 +4,29 @@ using MPSKit: GenericMPSTensor, MPSBondTensor
 # Environment transfer functions
 #
 
+function MPSKit.transfer_left(
+    GL::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    A::GenericMPSTensor{S,N},
+    Ā::GenericMPSTensor{S,N},
+) where {S,N}
+    Ā = twistdual(Ā, 2:N)
+    return _transfer_left(GL, O, A, Ā)
+end
+
+function MPSKit.transfer_right(
+    GR::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    A::GenericMPSTensor{S,N},
+    Ā::GenericMPSTensor{S,N},
+) where {S,N}
+    Ā = twistdual(Ā, 2:N)
+    return _transfer_right(GR, O, A, Ā)
+end
+
 ## PEPS
 
-function MPSKit.transfer_left(
+function _transfer_left(
     GL::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     A::GenericMPSTensor{S,3},
@@ -20,7 +40,7 @@ function MPSKit.transfer_left(
         A[χ_NW D_N_above D_N_below; χ_NE]
 end
 
-function MPSKit.transfer_right(
+function _transfer_right(
     GR::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     A::GenericMPSTensor{S,3},
@@ -36,7 +56,7 @@ end
 
 ## PEPO
 
-@generated function MPSKit.transfer_left(
+@generated function _transfer_left(
     GL::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     A::GenericMPSTensor{S,N},
@@ -65,7 +85,7 @@ end
     return macroexpand(@__MODULE__, :(return @autoopt @tensor $GL´_e := $rhs))
 end
 
-@generated function MPSKit.transfer_right(
+@generated function _transfer_right(
     GR::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     A::GenericMPSTensor{S,N},
@@ -118,7 +138,14 @@ end
 # Derivative contractions
 #
 
-@generated function MPSKit.∂C(
+function MPSKit.∂C(
+    C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N}
+) where {S,N}
+    GL = twistdual(GL, 1)
+    GR = twistdual(GR, numind(GR))
+    return _∂C(C, 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)
@@ -128,9 +155,20 @@ end
     return macroexpand(@__MODULE__, :(return @tensor $C´_e := $GL_e * $C_e * $GR_e))
 end
 
+function MPSKit.∂AC(
+    AC::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    GL::GenericMPSTensor{S,N},
+    GR::GenericMPSTensor{S,N},
+) where {S,N}
+    GL = twistdual(GL, 1)
+    GR = twistdual(GR, numind(GR))
+    return _∂AC(AC, O, GL, GR)
+end
+
 ## PEPS
 
-function MPSKit.∂AC(
+function _∂AC(
     AC::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     GL::GenericMPSTensor{S,3},
@@ -162,7 +200,7 @@ end
 
 ## PEPO
 
-@generated function MPSKit.∂AC(
+@generated function _∂AC(
     AC::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     GL::GenericMPSTensor{S,N},

src/algorithms/toolbox.jl

@@ -142,12 +142,11 @@ environment `env`.
 """
 function _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv)
     r, c = ind
-    return tr(
-        env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] *
-        env.corners[NORTHEAST, _prev(r, end), c] *
-        env.corners[SOUTHEAST, r, c] *
-        env.corners[SOUTHWEST, r, _prev(c, end)],
-    )
+    C_NW = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)]
+    C_NE = env.corners[NORTHEAST, _prev(r, end), c]
+    C_SE = env.corners[SOUTHEAST, r, c]
+    C_SW = env.corners[SOUTHWEST, r, _prev(c, end)]
+    return @tensor C_NW[1; 2] * C_NE[2; 3] * C_SE[3; 4] * C_SW[4; 1]
 end
 
 """

src/networks/local_sandwich.jl

@@ -55,6 +55,39 @@ function virtualspace(O::PEPSSandwich, dir)
     return virtualspace(ket(O), dir) ⊗ virtualspace(bra(O), dir)'
 end
 
+# not overloading MPOTensor because that defines AbstractTensorMap{<:Any,S,2,2}(::PEPSTensor, ::PEPSTensor)
+# ie type piracy
+mpotensor(top::PEPSTensor) = mpotensor((top, top))
+function mpotensor((top, bot)::PEPSSandwich)
+    @assert virtualspace(top, NORTH) == dual(virtualspace(top, SOUTH)) &&
+        virtualspace(bot, NORTH) == dual(virtualspace(bot, SOUTH)) &&
+        virtualspace(top, EAST) == dual(virtualspace(top, WEST)) &&
+        virtualspace(bot, EAST) == dual(virtualspace(bot, WEST)) &&
+        isdual(virtualspace(top, NORTH)) &&
+        isdual(virtualspace(bot, NORTH)) &&
+        isdual(virtualspace(top, EAST)) &&
+        isdual(virtualspace(bot, EAST)) "Method not yet implemented for given virtual spaces"
+
+    F_west = isomorphism(
+        storagetype(top),
+        fuse(virtualspace(top, WEST), virtualspace(bot, WEST)'),
+        virtualspace(top, WEST) ⊗ virtualspace(bot, WEST)',
+    )
+    F_south = isomorphism(
+        storagetype(top),
+        fuse(virtualspace(top, SOUTH), virtualspace(bot, SOUTH)'),
+        virtualspace(top, SOUTH) ⊗ virtualspace(bot, SOUTH)',
+    )
+    @tensor O[west south; north east] :=
+        top[phys; top_north top_east top_south top_west] *
+        conj(bot[phys; bot_north bot_east bot_south bot_west]) *
+        twist(F_west, 3)[west; top_west bot_west] *
+        twist(F_south, 3)[south; top_south bot_south] *
+        conj(F_west[east; top_east bot_east]) *
+        conj(F_south[north; top_north bot_north])
+    return O
+end
+
 ## PEPO
 
 const PEPOSandwich{N,T<:PEPSTensor,P<:PEPOTensor} = Tuple{T,T,Vararg{P,N}}

src/operators/transfermatrix.jl

@@ -125,53 +125,59 @@ const MultilineTransferMatrix = Union{MultilineTransferPEPS,MultilineTransferPEP
 virtualspace(O::InfiniteTransferMatrix, i, dir) = virtualspace(O[i], dir)
 
 """
-    initializeMPS(
+    initialize_mps(
+        f=randn,
+        T=scalartype(O),
         O::Union{InfiniteTransferPEPS,InfiniteTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,1}
     )
-    initializeMPS(
+    initialize_mps(
+        f=randn,
+        T=scalartype(O),
         O::Union{MultilineTransferPEPS,MultilineTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,2}
     )
 
 Inialize a boundary MPS for the transfer operator `O` by specifying an array of virtual
 spaces consistent with the unit cell.
 """
-function initializeMPS(
-    O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S,1}
+function initialize_mps(O::Union{InfiniteTransferMatrix,MultilineTransferMatrix}, arg) # initialize(f=randn, T=scalartype(O), O, ...)
+    return initialize_mps(randn, scalartype(O), O, arg)
+end
+function initialize_mps(
+    f, T, O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S,1}
 ) where {S}
     return InfiniteMPS([
-        randn(
-            scalartype(O),
+        f(
+            T,
             virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)),
             virtualspaces[mod1(i, end)],
         ) for i in 1:length(O)
     ])
 end
-function initializeMPS(O::InfiniteTransferMatrix, χ::Int)
+function initialize_mps(f, T, O::InfiniteTransferMatrix, χ::Int)
     return InfiniteMPS([
-        randn(calartype(O), ℂ^χ * _elementwise_dual(north_virtualspace(O, i)), ℂ^χ) for
-        i in 1:length(O)
+        f(T, ℂ^χ * _elementwise_dual(north_virtualspace(O, i)), ℂ^χ) for i in 1:length(O)
     ])
 end
-function initializeMPS(
-    O::MultilineTransferMatrix, virtualspaces::AbstractArray{S,2}
+function initialize_mps(
+    f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S,2}
 ) where {S}
     mpss = map(1:size(O, 1)) do r
-        return initializeMPS(O[r], virtualspaces[r, :])
+        return initialize_mps(f, T, O[r], virtualspaces[r, :])
     end
     return MPSKit.Multiline(mpss)
 end
-function initializeMPS(
-    O::MultilineTransferMatrix, virtualspaces::AbstractArray{S,1}
+function initialize_mps(
+    f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S,1}
 ) where {S}
-    return initializeMPS(O, repeat(virtualspaces, length(O), 1))
+    return initialize_mps(f, T, O, repeat(virtualspaces, length(O), 1))
 end
-function initializeMPS(O::MultilineTransferMatrix, V::ElementarySpace)
-    return initializeMPS(O, repeat([V], length(O), length(O[1])))
+function initialize_mps(f, T, O::MultilineTransferMatrix, V::ElementarySpace)
+    return initialize_mps(f, T, O, repeat([V], length(O), length(O[1])))
 end
-function initializeMPS(O::MultilineTransferMatrix, χ::Int)
-    return initializeMPS(O, repeat([ℂ^χ], length(O), length(O[1])))
+function initialize_mps(f, T, O::MultilineTransferMatrix, χ::Int)
+    return initialize_mps(f, T, O, repeat([ℂ^χ], length(O), length(O[1])))
 end
 
 @doc """

src/utility/util.jl

@@ -102,6 +102,22 @@ function absorb_s(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensor
     return u * sqrt_s, sqrt_s * vh
 end
 
+"""
+    twistdual(t::AbstractTensorMap, i)
+    twistdual!(t::AbstractTensorMap, i)
+
+Twist the i-th leg of a tensor `t` if it represents a dual space.
+"""
+function twistdual!(t::AbstractTensorMap, i::Int)
+    isdual(space(t, i)) || return t
+    return twist!(t, i)
+end
+function twistdual!(t::AbstractTensorMap, is)
+    is′ = filter(i -> isdual(space(t, i)), is)
+    return twist!(t, is′)
+end
+twistdual(t::AbstractTensorMap, is) = twistdual!(copy(t), is)
+
 # Check whether diagonals contain degenerate values up to absolute or relative tolerance
 function is_degenerate_spectrum(
     S; atol::Real=0, rtol::Real=atol > 0 ? 0 : sqrt(eps(scalartype(S)))

test/boundarymps/vumps.jl

@@ -7,12 +7,14 @@ using LinearAlgebra
 
 Random.seed!(29384293742893)
 
-const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitian=false))
+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))
 
     T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1)
-    mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)])
+    mps = initialize_mps(T, [ComplexSpace(20)])
 
     mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     N = abs(sum(expectation_value(mps, T)))
@@ -27,7 +29,7 @@ end
     psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2))
     T = PEPSKit.MultilineTransferPEPS(psi, 1)
 
-    mps = PEPSKit.initializeMPS(T, fill(ComplexSpace(20), 2, 2))
+    mps = initialize_mps(rand, scalartype(T), T, fill(ComplexSpace(20), 2, 2))
     mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     N = abs(prod(expectation_value(mps, T)))
 
@@ -37,6 +39,41 @@ end
     @test N ≈ N´ rtol = 1e-2
 end
 
+@testset "Fermionic PEPS" begin
+    D = Vect[fℤ₂](0 => 1, 1 => 1)
+    d = Vect[fℤ₂](0 => 1, 1 => 1)
+    χ = Vect[fℤ₂](0 => 10, 1 => 10)
+
+    psi = InfinitePEPS(D, d; unitcell=(1, 1))
+    n = InfiniteSquareNetwork(psi)
+    T = InfiniteTransferPEPS(psi, 1, 1)
+
+    # compare boundary MPS contraction to CTMRG contraction
+    mps = initialize_mps(T, [χ])
+    mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
+    N_vumps = abs(prod(expectation_value(mps, T)))
+
+    ctm, = leading_boundary(CTMRGEnv(psi, χ), psi)
+    N_ctm = abs(norm(psi, ctm))
+
+    @test N_vumps ≈ N_ctm rtol = 1e-2
+
+    # and again after blocking the local sandwiches
+    n´ = InfiniteSquareNetwork(map(PEPSKit.mpotensor, PEPSKit.unitcell(n)))
+    T´ = InfiniteMPO(map(PEPSKit.mpotensor, T.O))
+
+    mps´ = InfiniteMPS(randn, ComplexF64, [physicalspace(T´, 1)], [χ])
+    mps´, env´, ϵ = leading_boundary(mps´, T´, vumps_alg)
+    N_vumps´ = abs(prod(expectation_value(mps´, T´)))
+
+    ctm´, = leading_boundary(CTMRGEnv(n´, χ), n´)
+    N_ctm´ = abs(network_value(n´, ctm´))
+
+    @show N_vumps´
+    @test N_vumps´ ≈ N_vumps rtol = 1e-2
+    @test N_vumps´ ≈ N_ctm´ rtol = 1e-2
+end
+
 @testset "PEPO runthrough" begin
     function ising_pepo(beta; unitcell=(1, 1, 1))
         t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)]
@@ -64,7 +101,7 @@ end
     psi = PEPSKit.initializePEPS(O, ComplexSpace(2))
     T = InfiniteTransferPEPO(psi, O, 1, 1)
 
-    mps = PEPSKit.initializeMPS(T, [ComplexSpace(10)])
+    mps = initialize_mps(rand, scalartype(T), T, [ComplexSpace(10)])
     mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     f = abs(prod(expectation_value(mps, T)))
 
@@ -73,7 +110,7 @@ end
     psi2 = initializePEPS(O, ComplexSpace(2))
     T = InfiniteTransferPEPO(psi, O, 1, 1)
 
-    mps = PEPSKit.initializeMPS(T, [ComplexSpace(8)])
+    mps = initialize_mps(rand, scalartype(T), T, [ComplexSpace(8)])
     mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     f = abs(prod(expectation_value(mps, T)))
 end