Remove InfiniteWeightPEPS and mirror_antidiag

examples/heisenberg_su/main.jl

@@ -28,7 +28,7 @@ md"""
 To construct the Heisenberg Hamiltonian as just discussed, we'll use `heisenberg_XYZ` and,
 in addition, make it real (`real` and `imag` works for `LocalOperator`s) since we want to
 use PEPS and environments with real entries. We can either initialize the Hamiltonian with
-no internal symmetries (`symm = Trivial`) or use the global $U(1)$ symmetry
+no internal symmetries (`symm = Trivial`) or use the global spin $U(1)$ symmetry
 (`symm = U1Irrep`):
 """
 
@@ -39,8 +39,9 @@ H = real(heisenberg_XYZ(ComplexF64, symm, InfiniteSquare(Nr, Nc); Jx=1, Jy=1, Jz
 md"""
 ## Simple updating
 
-We proceed by initializing a random weighted PEPS that will be evolved. First though, we
-need to define the appropriate (symmetric) spaces:
+We proceed by initializing a random PEPS that will be evolved. 
+The weights used for simple update are initialized as identity matrices.
+First though, we need to define the appropriate (symmetric) spaces:
 """
 
 Dbond = 4
@@ -57,7 +58,8 @@ else
     error("not implemented")
 end
 
-wpeps = InfiniteWeightPEPS(rand, Float64, physical_space, bond_space; unitcell=(Nr, Nc));
+peps = InfinitePEPS(rand, Float64, physical_space, bond_space; unitcell=(Nr, Nc));
+wts = SUWeight(peps);
 
 md"""
 Next, we can start the `SimpleUpdate` routine, successively decreasing the time intervals
@@ -73,8 +75,7 @@ trscheme_peps = truncerr(1e-10) & truncdim(Dbond)
 
 for (dt, tol) in zip(dts, tols)
     alg = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
-    result = simpleupdate(wpeps, H, alg; bipartite=true)
-    global wpeps = result[1]
+    global peps, wts, = simpleupdate(peps, wts, H, alg; bipartite=true)
 end
 
 md"""
@@ -83,8 +84,7 @@ md"""
 In order to compute observable expectation values, we need to converge a CTMRG environment
 on the evolved PEPS. Let's do so:
 """
-
-peps = InfinitePEPS(wpeps) ## absorb the weights
+normalize!.(peps.A, Inf)
 env₀ = CTMRGEnv(rand, Float64, peps, env_space)
 trscheme_env = truncerr(1e-10) & truncdim(χenv)
 env, = leading_boundary(
@@ -136,5 +136,5 @@ well as finite bond dimension effects and a lacking extrapolation:
 
 E_ref = -0.6675
 M_ref = 0.3767
-@show (E - E_ref) / E_ref
-@show (mean(M_norms) - M_ref) / E_ref;
+@show (E - E_ref) / abs(E_ref)
+@show (mean(M_norms) - M_ref) / M_ref;

examples/hubbard_su/main.jl

@@ -32,47 +32,46 @@ t = 1
 U = 6
 Nr, Nc = 2, 2
 H = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(Nr, Nc); t, U, mu=U / 2);
+physical_space = Vect[fℤ₂](0 => 2, 1 => 2);
 
 md"""
 ## Running the simple update algorithm
 
-Next, we'll specify the virtual PEPS bond dimension and define the fermionic physical and
-virtual spaces. The simple update algorithm evolves an infinite PEPS with weights on the
-virtual bonds, so we here need to intialize an [`InfiniteWeightPEPS`](@ref). By default,
-the bond weights will be identity. Unlike in the other examples, we here use tensors with
-real `Float64` entries:
+Suppose the goal is to use imaginary-time simple update to optimize a PEPS 
+with bond dimension D = 8, and $2 \times 2$ unit cells.
+For a challenging model like the Hubbard model, a naive evolution starting from a 
+random PEPS at D = 8 will almost always produce a sub-optimal state.
+In this example, we shall demonstrate some common practices to improve SU result.
+
+First, we shall use a small D for the random PEPS initialization, which is chosen as 4 here.
+For convenience, here we work with real tensors with `Float64` entries.
+The bond weights are still initialized as identity matrices. 
 """
 
-Dbond = 8
-physical_space = Vect[fℤ₂](0 => 2, 1 => 2)
-virtual_space = Vect[fℤ₂](0 => Dbond / 2, 1 => Dbond / 2)
-wpeps = InfiniteWeightPEPS(rand, Float64, physical_space, virtual_space; unitcell=(Nr, Nc));
+virtual_space = Vect[fℤ₂](0 => 2, 1 => 2)
+peps = InfinitePEPS(rand, Float64, physical_space, virtual_space; unitcell=(Nr, Nc));
+wts = SUWeight(peps)
 
 md"""
-Let's set the algorithm parameters: The plan is to successively decrease the time interval of
-the Trotter-Suzuki as well as the convergence tolerance such that we obtain a more accurate
-result at each iteration. To run the simple update, we call [`simpleupdate`](@ref) where we
-use the keyword `bipartite=false` - meaning that we use the full $2 \times 2$ unit cell
-without assuming a bipartite structure. Thus, we can start evolving:
+Starting from the random state, we first use a relatively large evolution time step 
+`dt = 1e-2`. After convergence at D = 4, to avoid stucking at some bad local minimum,
+we first increase D to 12, and drop it back to D = 8 after a while. 
+Afterwards, we keep D = 8 and gradually decrease `dt` to `1e-4` to improve convergence.
 """
 
-dts = [1e-2, 1e-3, 4e-4, 1e-4]
-tols = [1e-6, 1e-8, 1e-8, 1e-8]
+dts = [1e-2, 1e-2, 1e-3, 4e-4, 1e-4]
+tols = [1e-7, 1e-7, 1e-8, 1e-8, 1e-8]
+Ds = [4, 12, 8, 8, 8]
 maxiter = 20000
 
-for (n, (dt, tol)) in enumerate(zip(dts, tols))
+for (dt, tol, Dbond) in zip(dts, tols, Ds)
     trscheme = truncerr(1e-10) & truncdim(Dbond)
     alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-    global wpeps, = simpleupdate(wpeps, H, alg; bipartite=false)
+    global peps, wts, = simpleupdate(
+        peps, wts, H, alg; bipartite=false, check_interval=2000
+    )
 end
 
-md"""
-To obtain the evolved `InfiniteWeightPEPS` as an actual PEPS without weights on the bonds,
-we can just call the following constructor:
-"""
-
-peps = InfinitePEPS(wpeps);
-
 md"""
 ## Computing the ground-state energy
 
@@ -85,11 +84,11 @@ run:
 
 χenv₀, χenv = 6, 16
 env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2)
-
+normalize!.(peps.A, Inf)
 env = CTMRGEnv(rand, Float64, peps, env_space)
 for χ in [χenv₀, χenv]
     global env, = leading_boundary(
-        env, peps; alg=:sequential, tol=1e-5, trscheme=truncdim(χ)
+        env, peps; alg=:sequential, tol=1e-8, maxiter=50, trscheme=truncdim(χ)
     )
 end
 
@@ -109,4 +108,4 @@ agree:
 
 Es_exact = Dict(0 => -1.62, 2 => -0.176, 4 => 0.8603, 6 => -0.6567, 8 => -0.5243)
 E_exact = Es_exact[U] - U / 2
-@show (E - E_exact) / E_exact;
+@show (E - E_exact) / abs(E_exact);

examples/j1j2_su/main.jl

@@ -25,21 +25,21 @@ Random.seed!(29385293);
 md"""
 ## Simple updating a challenging phase
 
-Let's start by initializing an `InfiniteWeightPEPS` for which we set the required parameters
-as well as physical and virtual vector spaces. We use the minimal unit cell size
-($2 \times 2$) required by the simple update algorithm for Hamiltonians with
-next-nearest-neighbour interactions:
+Let's start by initializing an `InfinitePEPS` for which we set the required parameters
+as well as physical and virtual vector spaces. 
+The `SUWeight` used by simple update will be initialized to identity matrices.
+We use the minimal unit cell size ($2 \times 2$) required by the simple update algorithm 
+for Hamiltonians with next-nearest-neighbour interactions:
 """
 
-Dbond, χenv, symm = 4, 32, U1Irrep
-trscheme_env = truncerr(1e-10) & truncdim(χenv)
+Dbond, symm = 4, U1Irrep
 Nr, Nc, J1 = 2, 2, 1.0
 
-## random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
+## random initialization of 2x2 iPEPS (using real numbers) and SUWeight
 Pspace = Vect[U1Irrep](1//2 => 1, -1//2 => 1)
 Vspace = Vect[U1Irrep](0 => 2, 1//2 => 1, -1//2 => 1)
-Espace = Vect[U1Irrep](0 => χenv ÷ 2, 1//2 => χenv ÷ 4, -1//2 => χenv ÷ 4)
-wpeps = InfiniteWeightPEPS(rand, Float64, Pspace, Vspace; unitcell=(Nr, Nc));
+peps = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell=(Nr, Nc));
+wts = SUWeight(peps);
 
 md"""
 The value $J_2 / J_1 = 0.5$ corresponds to a [possible spin liquid phase](@cite liu_gapless_2022),
@@ -56,8 +56,7 @@ for J2 in 0.1:0.1:0.5
     H = real( ## convert Hamiltonian `LocalOperator` to real floats
         j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice=false),
     )
-    result = simpleupdate(wpeps, H, alg; check_interval)
-    global wpeps = result[1]
+    global peps, wts, = simpleupdate(peps, wts, H, alg; check_interval)
 end
 
 md"""
@@ -71,19 +70,20 @@ J2 = 0.5
 H = real(j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice=false))
 for (dt, tol) in zip(dts, tols)
     alg′ = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
-    result = simpleupdate(wpeps, H, alg′; check_interval)
-    global wpeps = result[1]
+    global peps, wts, = simpleupdate(peps, wts, H, alg′; check_interval)
 end
 
 md"""
 ## Computing the simple update energy estimate
 
 Finally, we measure the ground-state energy by converging a CTMRG environment and computing
-the expectation value, where we make sure to normalize by the unit cell size:
+the expectation value, where we first normalize tensors in the PEPS:
 """
 
-peps = InfinitePEPS(wpeps)
-normalize!.(peps.A, Inf) ## normalize PEPS with absorbed weights by largest element
+normalize!.(peps.A, Inf) ## normalize each PEPS tensor by largest element
+χenv = 32
+trscheme_env = truncerr(1e-10) & truncdim(χenv)
+Espace = Vect[U1Irrep](0 => χenv ÷ 2, 1//2 => χenv ÷ 4, -1//2 => χenv ÷ 4)
 env₀ = CTMRGEnv(rand, Float64, peps, Espace)
 env, = leading_boundary(env₀, peps; tol=1e-10, alg=:sequential, trscheme=trscheme_env);
 E = expectation_value(peps, H, env) / (Nr * Nc)

src/PEPSKit.jl

@@ -28,7 +28,6 @@ include("utility/util.jl")
 include("utility/diffable_threads.jl")
 include("utility/svd.jl")
 include("utility/rotations.jl")
-include("utility/mirror.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
 include("utility/retractions.jl")
@@ -38,7 +37,6 @@ include("networks/local_sandwich.jl")
 include("networks/infinitesquarenetwork.jl")
 
 include("states/infinitepeps.jl")
-include("states/infiniteweightpeps.jl")
 include("states/infinitepartitionfunction.jl")
 
 include("operators/infinitepepo.jl")
@@ -49,6 +47,7 @@ include("operators/models.jl")
 
 include("environments/ctmrg_environments.jl")
 include("environments/vumps_environments.jl")
+include("environments/suweight.jl")
 
 include("algorithms/contractions/ctmrg_contractions.jl")
 include("algorithms/contractions/localoperator.jl")
@@ -104,7 +103,7 @@ export su_iter, su3site_iter, simpleupdate, SimpleUpdate
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
 export InfinitePEPS, InfiniteTransferPEPS
-export SUWeight, InfiniteWeightPEPS
+export SUWeight
 export InfinitePEPO, InfiniteTransferPEPO
 export initialize_mps, initializePEPS
 export ReflectDepth, ReflectWidth, Rotate, RotateReflect

src/algorithms/time_evolution/evoltools.jl

@@ -77,20 +77,18 @@ to get the reduced tensors
 ```
         2                   1
         |                   |
-    5 - A ← 3   ====>   4 - X ← 2   1 ← a ← 3
+    5 - A - 3   ====>   4 - X ← 2   1 ← a - 3
         | ↘                 |            ↘
         4   1               3             2
 
         2                               1
         |                               |
-    5 ← B - 3   ====>   1 ← b → 3   4 → Y - 2
+    5 - B - 3   ====>   1 - b → 3   4 → Y - 2
         | ↘                  ↘          |
         4   1                 2         3
 ```
 """
 function _qr_bond(A::PEPSTensor, B::PEPSTensor)
-    # TODO: relax dual requirement on the bonds
-    @assert isdual(space(A, 3)) # currently only allow A ← B
     X, a = leftorth(A, ((2, 4, 5), (1, 3)))
     Y, b = leftorth(B, ((2, 3, 4), (1, 5)))
     @assert !isdual(space(a, 1))
@@ -125,7 +123,7 @@ $(SIGNATURES)
 
 Apply 2-site `gate` on the reduced matrices `a`, `b`
 ```
-    -1← a -← 3 -← b ← -4
+    -1← a --- 3 --- b ← -4
         ↓           ↓
         1           2
         ↓           ↓
@@ -140,14 +138,16 @@ function _apply_gate(
     gate::AbstractTensorMap{T,S,2,2},
     trscheme::TruncationScheme,
 ) where {T<:Number,S<:ElementarySpace}
+    V = space(b, 1)
+    need_flip = isdual(V)
     @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4]
-    trunc = if trscheme isa FixedSpaceTruncation
-        V = space(b, 1)
-        truncspace(isdual(V) ? V' : V)
-    else
-        trscheme
+    trunc = (trscheme isa FixedSpaceTruncation) ? truncspace(V) : trscheme
+    a, s, b, ϵ = tsvd!(a2b2; trunc, alg=TensorKit.SVD())
+    a, b = absorb_s(a, s, b)
+    if need_flip
+        a, s, b = flip_svd(a, s, b)
     end
-    return tsvd!(a2b2; trunc, alg=TensorKit.SVD())
+    return a, s, b, ϵ
 end
 
 """