Infinite DMRG convergence in a critical region

How do I use this algorithm? What does that parameter do?
Post Reply
ahenry
Posts: 2
Joined: 21 Jul 2021, 01:45

Infinite DMRG convergence in a critical region

Post by ahenry »

Hi,

I'm attempting to do some infinite DMRG calculations with the Superintegrable Chiral Potts model (spin chain version). However I'm having difficulty getting DMRG to converge in the critical region. I've attached some example code, which should demonstrate the problem directly. It has some other tests which are commented out.

The model uses N-state "clocks", reducing to the Ising model for N=2. It has one parameter gamma, which controls the strength of the two terms, and is expected to have a critical region between gamma=0.493 and gamma=0.507. I can give you more details about the model if needed. There are some comments in the code too.

I find that DMRG converges well outside the critical region, and for finite systems anywhere. When it doesn't converge (infinite critical systems), it reaches an energy close to the ground state, but the energy and entropy oscillate randomly with additional sweeps. I've tried changing some DMRG parameters, including chi_mix, svd_min, mixer settings, unit cell size, and using the two-site algorithm. None of these seemed to help but it's possible I did not apply them correctly.

I find that turning the Hamiltonian's complex coefficients off i.e. setting use_sicp=0 in the example code resolves the convergence issues. In this case the model should have a critical point at gamma=0.5, rather than a critical region.

TEBD, while slow, appears to converge correctly. I was wondering if you can suggest anything else that is worth trying in terms of DMRG parameters or otherwise. Perhaps it is better to just use TEBD in the critical region.

Code: Select all

#!/usr/bin/env python3


# This code implements and tests the Superintegrable Chiral Potts model
# It contains simple tests with TEBD and DMRG
# Corresponding tests with the Ising model are included

import numpy as np
import os
import math
import logging

from tenpy.networks.site import Site
from tenpy.linalg import np_conserved as npc
from tenpy.models.model import CouplingMPOModel
from tenpy.networks.mps import MPS
from tenpy.models.spins import SpinModel
from tenpy.models.model import NearestNeighborModel

from tenpy.algorithms import dmrg
from tenpy.algorithms import tebd

from tenpy.models.tf_ising import TFIChain


# controls whether the SICP model uses the complex SICP coefficient
# if use_sicp = 0, the coefficients are just 1 (a slightly different model)
use_sicp = 1

# controls whether to tell tenpy to add hermitian conjugates
# this should not be necessary as the hamiltonian is already hermitian
use_hcs = 0

# L = 200
# bc_MPS = 'finite'
L = 2
bc_MPS = 'infinite'
TFI_model_params = {
    'L' : L,
    'J' : 0.5,
    'g' : 0.5,
    'bc_MPS' : bc_MPS,
    'conserve' : None,
    }
SICP_model_params = {
    'L' : L,
    'gamma' : 0.5,
    'bc_MPS' : bc_MPS,
    'N' : 3,
    'conserve' : 'Q',
    }
dmrg_params = {
    'mixer': True,
    'trunc_params': {
        'chi_max': 30,
        'svd_min': 1.e-10
    },
    'max_E_err': 1.e-6,
    'active_sites': 1,
}
tebd_params = {
    'order': 1,
    'delta_tau_list': [0.1, 0.01, 0.001, 1.e-4, 1.e-5, 1.e-6],
    'N_steps': 10,
    'max_error_E': 1.e-8,
    'trunc_params': {
        'chi_max': 30,
        'svd_min': 1.e-10
    },
}

class SICPSite(Site):
    def __init__(self, N=2, conserve="Q"):
        if conserve not in ["Q", "parity", None]:
            raise ValueError("invalid `conserve`: " + repr(conserve))
        self.N = N

        # The main operators we use are Z and X, which generalise the Pauli matrices

        # Z measures the clock at one site
        Z = np.zeros([N, N], dtype=np.complex64)
        for i in range(N):
            Z[i][i] = np.exp(-2.0j * i * np.pi / N)
        ZP = np.zeros([N, N], dtype=np.complex64)
        for i in range(N):
            ZP[i][i] = np.exp(2.0j * i * np.pi / N)


        # X increments a clock
        # unlike spins, clocks cycle back to the first state
        X = np.zeros([N, N])
        XP = np.zeros([N, N])
        for i in range(N):
            iP = (i + 1) % N
            XP[i][iP] = 1
            X[iP][i] = 1

        # QN operator
        # the model has a conserved quantity: the total of all the clocks mod N
        Q = np.arange(N)

        # add powers of Z and X for convenience
        ops = dict(Z=Z, ZP=ZP, X=X, XP=XP)
        for i in range(1, N):
            Zpow = np.linalg.matrix_power(Z, i)
            ops[f"Z{i}"] = Zpow
            ZPpow = np.linalg.matrix_power(ZP, i)
            ops[f"ZP{i}"] = ZPpow
            Xpow = np.linalg.matrix_power(X, i)
            ops[f"X{i}"] = Xpow
            XPpow = np.linalg.matrix_power(XP, i)
            ops[f"XP{i}"] = XPpow

        if conserve == "Q":
            chinfo = npc.ChargeInfo([N], ["Q"])
            leg = npc.LegCharge.from_qflat(chinfo, Q)
        else:
            if conserve == "parity":
                chinfo = npc.ChargeInfo([N], ["parity_Q"])
                leg = npc.LegCharge.from_qflat(chinfo, Q)
            else:
                leg = npc.LegCharge.from_trivial(N)
        self.conserve = conserve
        names = [str(i) for i in np.arange(N)]
        Site.__init__(self, leg, names, **ops)

    def __repr__(self):
        return f"SICPSite(N={self.N})"


def SICP_factor(r, N):
    u = np.exp(np.pi * 1.0j * (2.0 * r - N) / 2.0 / N) / np.sin(np.pi * r / N)
    return u

class SICP_model(CouplingMPOModel):
    def init_sites(self, model_params):
        N = model_params.get("N", 2)
        conserve = model_params.get("conserve", None)
        site = SICPSite(N, conserve)
        return site

    def init_terms(self, model_params):
        # controls the relative strength of the two terms
        g = model_params.get("gamma", 0.5)

        # number of states in the clock
        N = model_params.get("N", 2)

        # there are two terms, which directly generalise the Ising model
        # sometimes it's written with Z and X interchanged, which is equivalent
        # this way makes the QN easier to implement
        for u in range(len(self.lat.unit_cell)):
            for i in np.arange(1, N):
                a = 1
                if use_sicp: a = SICP_factor(i, N)
                self.add_onsite(-0.5 * a * g, u, f"Z{i}", plus_hc=use_hcs)

        for u1, u2, dx in self.lat.pairs["nearest_neighbors"]:
            for i in np.arange(1, N):
                a = 1
                if use_sicp: a = SICP_factor(i, N)
                self.add_coupling(
                    -0.5 * a * (1 - g), u1, f"XP{i}", u2, f"X{i}", dx, plus_hc=use_hcs
                )


# test functions
def TFI_dmrg():
    M = TFIChain(TFI_model_params)
    product_state = ["up"] * M.lat.N_sites
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    R = dmrg.run(psi, M, dmrg_params)
    E = R['E']
    if not bc_MPS == 'infinite':
        E = R['E'] / L
    S = psi.entanglement_entropy()[0]
    return(E, S)

def TFI_tebd():
    M = TFIChain(TFI_model_params)
    product_state = ["up"] * M.lat.N_sites
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    eng = tebd.TEBDEngine(psi, M, tebd_params)
    eng.run_GS()
    E = np.sum(M.bond_energies(psi)) / psi.L
    if not bc_MPS == 'infinite':
        E = np.sum(M.bond_energies(psi)) / psi.L
    # MPS.canonical_form_infinite(psi)
    S = psi.entanglement_entropy()[0]
    return(E, S)


def SICP_dmrg():
    M = SICP_model(SICP_model_params)
    product_state = ["0"] * M.lat.N_sites
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    R = dmrg.run(psi, M, dmrg_params)
    E = R['E']
    if not bc_MPS == 'infinite':
        E = R['E'] / L
    S = psi.entanglement_entropy()[0]
    return(E, S)

def SICP_tebd():
    M = SICP_model(SICP_model_params)
    M = NearestNeighborModel.from_MPOModel(M)
    product_state = ["0"] * M.lat.N_sites
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    eng = tebd.TEBDEngine(psi, M, tebd_params)
    eng.run_GS()
    E = np.sum(M.bond_energies(psi)) / psi.L
    if not bc_MPS == 'infinite':
        E = np.sum(M.bond_energies(psi)) / psi.L
    # MPS.canonical_form_infinite(psi)
    S = psi.entanglement_entropy()[0]
    return(E, S)

if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    # TFI_res_dmrg = TFI_dmrg()
    # print(TFI_res_dmrg)
    # TFI_res_tebd = TFI_tebd()
    # print(TFI_res_tebd)
    SICP_res_dmrg = SICP_dmrg()
    print("DMRG results:")
    print(SICP_res_dmrg)
    SICP_res_tebd = SICP_tebd()
    print("TEBD results:")
    print(SICP_res_tebd)

Post Reply