How to write a PXP model right

[adolepha]

I am a beginner in Python and TenPy, and I have read the notebook posted in the github https://tenpy.readthedocs.io/en/latest/ ... _tebd.html.
And I want to rewrite it to do something with the Rydberg-type model, the PXP model with open boundary condition.
$$H=X_{1} P_{2}+P_{N-1} X_{N}+\sum_{i=2}^{N-1} P_{i-1} X_{i} P_{i+1}$$
But It seems not work.

The simple model code I wrote is here, it has't consider the boundary right now:

Python: Select all

import numpy as np
import tenpy
from tenpy.models.model import CouplingModel
from tenpy.networks.site import SpinHalfSite
from tenpy.models.lattice import Chain
from tenpy.networks.mps import MPS
from tenpy.algorithms.tebd import TEBDEngine
import matplotlib.pyplot as plt
tenpy.tools.misc.setup_logging(to_stdout="INFO")
np.set_printoptions(precision=5, suppress=True, linewidth=100)
plt.rcParams['figure.dpi'] = 150

class PXPModel(CouplingModel):
    def __init__(self, L):
        # 定义二能级系统的格点（自旋 1/2）
        site = SpinHalfSite(conserve=None)  # 不守恒量子数
        # 添加自定义算符：P=|0><0|, X=泡利X
        P = np.array([[1, 0], [0, 0]])  # |0><0| 算符
        X = np.array([[0, 1], [1, 0]])  # 泡利X算符
        site.add_op('P', P)            # 注册为格点算符
        site.add_op('X', X)
        
        # 定义一维链（开边界）
        lat = Chain(L, site, bc="open", bc_MPS='finite')
        CouplingModel.__init__(self, lat)

        for i in range(L):  # 遍历所有可能中心位点
            # 定义相对位移（使用元组表示晶格方向）
            # displacements = [(-1,), (0,), (1,)]  # 左邻，中心，右邻
            self.add_multi_coupling(
            strength=1.0,
            # ops=[(0, 'P'), (1, 'X'), (2, 'P')],  # 三个算符的相对位置
            ops=[('P', [-1], 0), ('X', [0], 0), ('P', [1], 0)],  # 三个算符的相对位置
            op_string='Id',
            plus_hc=False,  # 不添加厄米共轭（H 已是厄米）
            category='PXP'
            )
            

# =================================================================
# 2. 参数设置
# =================================================================
L = 8                   # 原子链长度
dt = 0.05               # 时间步长
t_max = 10.0            # 总演化时间
trunc_params = {
    'chi_max': 100,     # MPS最大键维度
    'svd_min': 1e-10    # 奇异值截断阈值
}

# =================================================================
# 3. 初始化模型和初始态 |01010101>
# =================================================================
model = PXPModel(L)     # 创建PXP模型
# 初始态：0对应"down"（即|0>），1对应"up"（即|1>）
initial_state = ["down", "up"] * (L // 2)  # ["down", "up", "down", ...]
psi = MPS.from_product_state(model.lat.mps_sites(), initial_state, bc="finite")

It seems no bugs when you only try to compile this PXPModel

But when i use the code below to use the TEBD method almost set in TenPy, there are some bugs and I don't know how to solve it.
The Python code:

Python: Select all

tebd_params = {
    'dt': 0.05,          # 时间步长
    'N_steps': 200,      # 总步数
    'order': 2,          # 二阶Suzuki-Trotter
    'trunc_params': {
        'chi_max': 200,
        'svd_min': 1e-12
    },
    'verbose': 1
}

# 执行时间演化
eng = TEBDEngine(psi, model, tebd_params)
times = []
mag_z = []

def measurement(eng, data):
    keys = ['t', 'entropy', 'Sx', 'Sz', 'corr_XX', 'corr_ZZ', 'trunc_err']
    if data is None:
        data = dict([(k, []) for k in keys])
    data['t'].append(eng.evolved_time)
    data['entropy'].append(eng.psi.entanglement_entropy())
    data['Sx'].append(eng.psi.expectation_value('Sigmax'))
    data['Sz'].append(eng.psi.expectation_value('Sigmaz'))
    data['corr_XX'].append(eng.psi.correlation_function('Sigmax', 'Sigmax'))
    data['trunc_err'].append(eng.trunc_err.eps)
    return data

data = measurement(eng, None)
while eng.evolved_time < 1.:
    eng.run()
    measurement(eng, data)

The out put:

Python: Select all

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[1], line 187
    185 data = measurement(eng, None)
    186 while eng.evolved_time < 1.:
--> 187     eng.run()
    188     measurement(eng, data)
    189 # for _ in range(tebd_params['N_steps']):
    190 #     eng.run()
    191 #     times.append(eng.evolved_time)
   (...)
    208 # plt.tight_layout()
    209 # plt.show()

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\algorithm.py:391, in TimeEvolutionAlgorithm.run(self)
    388 start_time = time.time()
    389 Sold = np.mean(self.psi.entanglement_entropy())
--> 391 self.run_evolution(N_steps, dt)
    393 S = self.psi.entanglement_entropy()
    394 logger.info(
    395     "--> time=%(tr)3.3f %(sign)s %(tc)3.3fj, max(chi)=%(chi)d, max(S)=%(S).5f, "
    396     "avg DeltaS=%(dS).4e, since last update: %(wall_time).1fs", {
   (...)
    403         'wall_time': time.time() - start_time,
    404     })

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\algorithm.py:419, in TimeEvolutionAlgorithm.run_evolution(self, N_steps, dt)
    416 if preserve_norm:
    417     old_norm = self.psi.norm
--> 419 self.prepare_evolve(dt)
    420 trunc_err = self.evolve(N_steps, dt)
    422 self.trunc_err = self.trunc_err + trunc_err  # not += : make a copy!

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\tebd.py:277, in TEBDEngine.prepare_evolve(self, dt)
    275 order = self.options.get('order', 2, int)
    276 E_offset = self.options.get('E_offset', None, 'real')
--> 277 self.calc_U(order, dt, type_evo='real', E_offset=E_offset)

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\tebd.py:318, in TEBDEngine.calc_U(self, order, delta_t, type_evo, E_offset)
    316 self._U = []
    317 for dt in self.suzuki_trotter_time_steps(order):
--> 318     U_bond = [
    319         self._calc_U_bond(i_bond, dt * delta_t, type_evo, E_offset) for i_bond in range(L)
    320     ]
    321     self._U.append(U_bond)
    322 self.force_prepare_evolve = False

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\tebd.py:319, in <listcomp>(.0)
    316 self._U = []
    317 for dt in self.suzuki_trotter_time_steps(order):
    318     U_bond = [
--> 319         self._calc_U_bond(i_bond, dt * delta_t, type_evo, E_offset) for i_bond in range(L)
    320     ]
    321     self._U.append(U_bond)
    322 self.force_prepare_evolve = False

File E:\Pycharm\code\venv\lib\site-packages\tenpy\algorithms\tebd.py:569, in TEBDEngine._calc_U_bond(self, i_bond, dt, type_evo, E_offset)
    563 def _calc_U_bond(self, i_bond, dt, type_evo, E_offset):
    564     """Calculate exponential of a bond Hamiltonian.
    565 
    566     * ``U_bond = exp(-i dt (H_bond-E_offset_bond))`` for ``type_evo='real'``, or
    567     * ``U_bond = exp(- dt H_bond)`` for ``type_evo='imag'``.
    568     """
--> 569     h = self.model.H_bond[i_bond]
    570     if h is None:
    571         return None  # don't calculate exp(i H t), if `H` is None

AttributeError: 'PXPModel' object has no attribute 'H_bond'

My question is what's wrong with my silly code, and how can I get a right code? what is the usable code.
Or any notebook that can teach myself to construct a model step-by-step. Because I'm learning TenPy and Python by myself and no one around me can teach me those.

[Johannes]

There's some issue with your models, try reading reading the model guide for more details.
Some hints:

1. You're trying to run TEBD, which requires H_bond defined in a NearestNeighbourModel. However, the PXP model is not nearest-Neighbor, so you either need to use a different time evolution method like TDVP or ExpMPOEvolution, or use the "group" parameter of the simultations, which makes the model effectively nearest-neighbor at the expense of a larger local hilbert space (during the time evolution).
2. I would recommend to use the structure of the CouplingMPOModel rather than plain CouplingModel, since this will then also generate H_MPO which you need for DMRG or other time evolution methods.
   Further, it takes care of initializing the lattice for you, if you simply set the default_lattice on the model
3. You have a loop over i which doesn't use i. Don't do this: the add_multi_coupling already iterates over all possible 3-site couplings which fit into the lattice (given the boundary conditions). If you keep the loop, you will effectively add each term multiple times and hence get a larger prefactor which doesn't match with the boundary terms you add.
4. Use add_multi_coupling_term or add_coupling_term for the boundary terms on the very left and right, which doesn't iterate over the possible nearest-neighbor couplings.

Taking into account these hints, you should arrive at something like

Python: Select all

class PXPModel(CouplingMPOModel):
    default_lattice = 'Chain'
    force_default_lattice = True  # for the boundary terms, we implicitly assume 1D chain
    
    def init_sites(self, model_params):
        s = SpinHalfSite(conserve=None)
        s.add_op('P', np.array([[1., 0.], [0., 0.]]))
        s.add_op('X', np.array([[0., 1.], [1., 0.]]))
        return s

    def init_terms(self, model_params):
        L = model_params['L']
        J = model_params.get('J', 2.)
        self.add_multi_coupling(J, [('P', [-1], 0), ('X', [0],0), ('P', [1], 0)])
        
        if model_params['bc_x'] == 'open':
            self.add_coupling_term(J, 0, 1, 'X', 'P')
            self.add_coupling_term(J, L-2, L-1, 'P', 'X')


M = PXPModel({'L': 4})
print(M.all_coupling_terms().to_TermList())