CTMRG for 2D Ising model

We test the Corner transfer matrix renormalization group (CTMRG) algorithm by setting up a well-known exact PEPS: Thermal state of the 2D classical Ising model. We match the output of the CTMRG procedure with the exact results coming from Onsager solution of the 2D Ising model.

We enforce the Z2 symmetry of the model, which prevents spontaneous symmetry breaking in the ordered phase. Enforcing this symmetry is important because, in finite systems or numerical simulations, artificial symmetry breaking can occur due to numerical inaccuracies or finite-size effects. To that end, we consider Hamiltonian \(H = - \sum_{\langle i, j \rangle} X_i X_j\), where \(X_i\) is the first Pauli matrix operating at the site \(i\), and the sum is over the nearest neighbor sites \(\langle i, j \rangle\) on a square lattice.

import numpy as np
import pytest
import yastn
import yastn.tn.fpeps as fpeps

@pytest.mark.parametrize("checkpoint_move", ['reentrant','nonreentrant', False])
def test_ctmrg_Ising(checkpoint_move, config_kwargs):
    r"""
    Use CTMRG to calculate some expectation values in classical 2D Ising model.
    Compare with analytical results.

    """
    #
    # We start by representing the partition function
    # of the model at temperature beta as a PEPS network.
    beta = 0.5
    #
    config = yastn.make_config(sym='Z2', **config_kwargs)
    if config.backend.BACKEND_ID != 'torch' and checkpoint_move != False:
        pytest.skip("checkpoint_move is not supported for this backend")
    #
    leg = yastn.Leg(config, s=1, t=(0, 1), D=(1, 1))
    T = yastn.ones(config, legs=[leg, leg, leg.conj(), leg.conj()], n=0)
    X = yastn.ones(config, legs=[leg, leg, leg.conj(), leg.conj()], n=1)
    B = yastn.zeros(config, legs=[leg, leg.conj()])
    B.set_block(ts=(0, 0), val=np.cosh(beta))
    B.set_block(ts=(1, 1), val=np.sinh(beta))
    #
    # The partition function of the system in the thermodynamic limit
    # follows from contracting infinite network:
    #
    #      |     |     |
    #    ──T──B──T──B──T──
    #      |     |     |
    #      B     B     B
    #      |     |     |
    #    ──T──B──T──B──T──
    #      |     |     |
    #      B     B     B
    #      |     |     |
    #    ──T──B──T──B──T──
    #      |     |     |
    #
    # The expectation values of spin correlation functions are obtained
    # replacing tensors T with X at sites of interest.
    #
    # We absorb bond tensors B into site tensors.
    #
    TB = yastn.ncon([T, B, B], [(-0, -1, 2, 3), [2, -2], [3, -3]])
    XB = yastn.ncon([X, B, B], [(-0, -1, 2, 3), [2, -2], [3, -3]])
    #
    # We use infinite translationally invariant square lattice with 1x1 unit cell.
    #
    geometry = fpeps.SquareLattice(dims=(1, 1), boundary='infinite')
    psi = fpeps.Peps(geometry=geometry, tensors={(0, 0): TB})
    #
    # Execute CTMRG algorithm to contract the network.
    # chi is the bond dimension of CTM environment tensors.
    # CTMRG updates are terminated based on the convergence of corner spectra.
    #
    chi = 24
    env = fpeps.EnvCTM(psi, init='eye')
    opts_svd = {"D_total": chi}
    info = env.ctmrg_(opts_svd=opts_svd, max_sweeps=200, corner_tol=1e-12, checkpoint_move=checkpoint_move)
    assert info.max_dsv < 1e-12
    assert info.converged == True
    #
    # Local magnetization is obtained by contracting
    # the CTM environment with tensor XB (normalized by the partition function).
    # It is zero by symmetry.
    #
    ev_X = env.measure_1site(XB)  # at all unique lattice sites: [(0, 0)]
    assert abs(ev_X[(0, 0)]) < 1e-10
    #
    # Nearest-neighbor XX correlator
    ev_XXnn = env.measure_nn(XB, XB)
    if beta == 0.5:
        assert abs(ev_XXnn[(0, 0), (0, 1)] - 0.872783) < 1e-6  # horizontal bond
        assert abs(ev_XXnn[(0, 0), (1, 0)] - 0.872783) < 1e-6  # vertical bond
    #
    # Exact magnetization squared in the ordered phase for beta > beta_c = 0.440686...
    #
    MX2 = (1 - np.sinh(2 * beta) ** -4) ** (1 / 4)
    #
    # We calculate the square of the spontaneous magnetization
    # from the large-distance limit of XX correlator
    #
    pairs = [((0, 0), (0, 100)),  # horizontal
             ((0, 0), (100, 0))]  # vertical
    for pair in pairs:
        ev_XXlong = env.measure_line(XB, XB, sites=pair)
        assert abs(MX2 - ev_XXlong) < 1e-10

config_kwargs = {"backend": "np"}
test_ctmrg_Ising(config_kwargs)