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)
```