diff --git a/src/algorithms/ctmrg/c4v.jl b/src/algorithms/ctmrg/c4v.jl index bfefcae34..093506ac4 100644 --- a/src/algorithms/ctmrg/c4v.jl +++ b/src/algorithms/ctmrg/c4v.jl @@ -71,7 +71,16 @@ end # as expensive as `svd_vals` on CPU but much cheaper on GPU, and gives effectively # the same information. For the QR projector, off-diagonal elements are still # present until CTMRG converges, so we need a fallback for the non-diagonal case. -corner_spectrum(C::AbstractTensorMap, ::C4vCTMRG) = eigh_vals(C) +# gauge fix sign of spectrum such that the eigenvalue of largest magnitude is positive +function corner_spectrum(C::AbstractTensorMap, ::C4vCTMRG) + S = eigh_vals(C) + if -minimum(S) > maximum(S) + for (_, b) in blocks(S) + b .= .-reverse(b) # spectrum of -C: negated and, to stay ascending, reversed + end + end + return S +end """ convergence_tensors(env::CTMRGEnv, alg::C4vCTMRG) -> (corners, edges) diff --git a/test/ctmrg/c4v_sign_branches.jl b/test/ctmrg/c4v_sign_branches.jl new file mode 100644 index 000000000..10a5dca95 --- /dev/null +++ b/test/ctmrg/c4v_sign_branches.jl @@ -0,0 +1,10 @@ +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.ctmrg_c4v_sign_branches(Vector) +end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..db2c7b000 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -121,6 +121,12 @@ module CTMRGFlavors end using .CTMRGFlavors +module CTMRGC4vSignBranches + include("ctmrg/c4v_sign_branches.jl") + export ctmrg_c4v_sign_branches +end +using .CTMRGC4vSignBranches + module CTMRGGaugeFix include("ctmrg/gaugefix.jl") export ctmrg_gaugefix_asymmetric, ctmrg_gaugefix_c4v diff --git a/test/testsuite/ctmrg/c4v_sign_branches.jl b/test/testsuite/ctmrg/c4v_sign_branches.jl new file mode 100644 index 000000000..3f8311d50 --- /dev/null +++ b/test/testsuite/ctmrg/c4v_sign_branches.jl @@ -0,0 +1,61 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using PEPSKit +using Adapt +using PEPSKit: ctmrg_iteration, compute_gauge_fix_gauge, fix_relative_phases, ScramblingEnvGaugeC4v + +# For half-integer U(1) charges, (-1)^{2Q} on the ket bonds negates the sandwich, so fixed +# points come in pairs, one of which flips (C, E) to (-C, -E) every iteration. Both should converge. + +sd = 20260924 +H = heisenberg_XXZ(ComplexF64, U1Irrep, InfiniteSquare(); spin = 1 // 2) +V = U1Space(0 => 1, 1 // 2 => 1, -1 // 2 => 1) +Venv = U1Space(0 => 2, 1 // 2 => 2, -1 // 2 => 2, 1 => 1, -1 => 1) +tol = 1.0e-10 +maxiter = 400 + +# (-1)^{2Q} on a U(1) space +function charge_parity(W) + X = id(W) + for (c, b) in blocks(X) + b .*= (-1)^Int(2 * c.charge) + end + return X +end + +# sign by which one gauge-fixed iteration multiplies the edge +function iteration_sign(network, env, alg) + env′, = ctmrg_iteration(network, env, alg) + signs, = compute_gauge_fix_gauge(env′, env, ScramblingEnvGaugeC4v()) + env′ = fix_relative_phases(env′, signs) + return sign(real(dot(env.edges[1], env′.edges[1]))) +end + +function sign_twin(env) + E = env.edges[1] + g = charge_parity(space(E, 2)) + @tensor Ê[χ_out D_ket D_bra; χ_in] := E[χ_out D_ket′ D_bra; χ_in] * g[D_ket; D_ket′] + return CTMRGEnv(env.corners[1], Ê) +end + +function ctmrg_c4v_sign_branches(AT) + return @testset "C4v CTMRG converges on both sign branches of a half-integer U(1) state ($AT)" begin + Random.seed!(sd) + peps = InfinitePEPS(randn, Float64, physicalspace(H)[1, 1], V) + symmetrize!(peps, RotateReflect()) + peps = adapt(AT, InfinitePEPS(map(a -> a / norm(a), peps.A))) + network = InfiniteSquareNetwork(peps) + alg = C4vCTMRG(; tol, maxiter, projector_alg = :C4vEighProjector) + + env, info = leading_boundary(initialize_random_c4v_env(peps, Venv), peps, alg) + @test info.converged + + env_twin = sign_twin(env) + @test iteration_sign(network, env_twin, alg) == -iteration_sign(network, env, alg) + env_twin, info_twin = leading_boundary(env_twin, peps, alg) + @test info_twin.converged + @test expectation_value(peps, H, env_twin) ≈ expectation_value(peps, H, env) atol = 1.0e-8 + end +end