1D Heisenberg model with periodic BC, how to

[Albert65]

Hi All,

I am new to TeNPy and currently trying to implement a 1D Heisenberg model with periodic boundary conditions (PBC).

I have read the documentation and searched online, but I have not found a working example that supports PBC with TEBD.
My current code (see below) works fine for open boundary conditions (OBC), but fails to behave correctly when I switch to periodic.

The goal is to simulate time evolution using the TEBD engine as far as possible in time, for a simple real-coefficient 1D Hamiltonian
(Heisenberg spin-1/2 chain). I am primarily interested in understanding dynamics and entanglement spreading.

My (naive?) understanding:
* TEBD applies two-site gates sequentially to nearest neighbours, so I expected that PBC could be handled similarly by wrapping around the chain.
* I realize that issues might arise during contractions for measurements, or perhaps in how the gates are applied across the periodic boundary.
* One intuition I have is that PBC might result in more uniform bond dimension growth across the chain, rather than peaking at the center like in OBC — potentially allowing longer evolution times.

The question:
Has anyone implemented TEBD with periodic boundaries in TeNPy for the Heisenberg model (or similar)?
Is there a known workaround or example that would help?
I would be very grateful for any insights, code snippets, or clarification on whether PBC + TEBD is actually feasible in practice within TeNPy.
Thanks a lot in advance!

Python: Select all

from tenpy.networks.site import SpinHalfSite
from tenpy.models.model import MPOModel, CouplingModel, NearestNeighborModel
from tenpy.models.spins import SpinModel
from tenpy.models.lattice import Chain
from tenpy.tools.params import asConfig


class HeisenbergChain(CouplingModel, NearestNeighborModel, MPOModel):
    """
    Hamiltonian model:
        H = (sum_{i=1}^n interaction_xx[i] Sx_i kron Sx_{i+1} +
                         interaction_yy[i] Sy_i kron Sy_{i+1} +
                         interaction_zz[i] Sz_i kron Sz_{i+1}) +
            (sum_{i=1}^n h_x[i] Sx_i + h_z[i] Sz_i)

    where Sx, Sy, Sz are spin-half matrices.
    """

    def __init__(self, model_params):
        model_params = asConfig(model_params, "XXZChain")

        # All parameters must be specified.
        L = model_params.get("L", None, int)
        bc = model_params.get("bc", None, str)
        Jx = model_params.get("Jx", None, "real_or_array")
        Jy = model_params.get("Jy", None, "real_or_array")
        Jz = model_params.get("Jz", None, "real_or_array")
        hx = model_params.get("hx", None, "real_or_array")
        hz = model_params.get("hz", None, "real_or_array")
        assert L >= 2
        assert bc in ["open", "periodic"]

        site = SpinHalfSite(conserve=None)
        lat = Chain(L, site, bc=bc, bc_MPS="finite")
        CouplingModel.__init__(self, lat)

        self.add_onsite(hz, 0, "Sz")
        self.add_onsite(hx, 0, "Sx")

        # for i in range(L - 1 if bc == "open" else L):
        for i in range(L - 1):
            self.add_coupling_term(Jx, i, i + 1, "Sx", "Sx")
            self.add_coupling_term(Jy, i, i + 1, "Sy", "Sy")
            self.add_coupling_term(Jz, i, i + 1, "Sz", "Sz")

        if bc == "periodic":
            self.add_coupling_term(Jx, 0, L - 1, "Sx", "Sx")
            self.add_coupling_term(Jy, 0, L - 1, "Sy", "Sy")
            self.add_coupling_term(Jz, 0, L - 1, "Sz", "Sz")

        print(self.lat.pairs)

        MPOModel.__init__(self, lat, self.calc_H_MPO())
        NearestNeighborModel.__init__(self, lat, self.calc_H_bond())

[Johannes]

Using bc='periodic' for the Lattice does not imply periodic boundary conditions on the MPS, but just makes the add_coupling method add the additional \(S_{L-1} \cdot S_0\) coupling , similar as you did explicitly "by hand" in this case, such that the Hamiltonian implements the same couplings as on a ring.
The bc_MPS is still "finite", though, meaning that the MPS itself has open boundary conditions, ranging from 0 to L-1.
Thus, the \(S_{L-1} \cdot S_0\) is a single long-range coupling over the whole system, and the tenpy.models.model.NearestNeighborModel complains that the model you write has a term which is not "nearest neighbor" for the MPS (although it being nearest neighbor in the hamiltonian on a ring.)

See this thread with more detailed explanations and the possible workaround of using order='folded'.

Let me repeat a disclaimer/warning here again: forcing the bc=periodic Hamiltonian terms despite using open/finite bc_MPS (which we need for the efficient, local optimization with canonical form) implies twice as much entanglement from the point of view of the MPS, thus you need *significantly* larger MPS bond dimension (at least double, but for the same error order of squared bond dimension (with a fraction depending on how quickly the singular values decay), i.e. you might need to go from something like chi=100 for open chains to something like chi=5000 for bc=periodic with equally good results!!!)


Note that the TeNPy simulation setup allows to group_sites_for_algorithm.
Adjusting the minimal_TEBD.yml simulation parameters a bit, you should be able to use something like:

YAML parameters: Select all

#!/usr/bin/env -S python -m tenpy 

simulation_class : RealTimeEvolution

directory: results
output_filename : a_minimal_tebd.h5

model_class :  SpinModel   # not SpinChain!
model_params :
    L : 32
    bc_MPS : finite
    Jz : 1.

    bc_x: periodic   # explicitly force periodic coupling in hamiltonian, MPS is still open!
    order: folded   # folded order such that instead of one long-range, we have almost all next-nearest neighbour couplings

group_sites: 2  # since folded order gives at most next-nearest neighbour couplings, we're back to nearest neighbour by grouping each 2 sites.
group_to_NearestNeighborModel: True

initial_state_params:
    method : lat_product_state
    product_state : [[up], [down]]

algorithm_class: TEBDEngine
algorithm_params:
    trunc_params:
        chi_max: 120
    dt : 0.05
    N_steps : 2  # measurements every dt*N_steps

final_time :  1.

# additional measurements during the evolution
connect_measurements: 
  - - tenpy.simulations.measurement
    - m_onsite_expectation_value
    - opname: Sz

Note that the measurement functions like "m_onsite_expectation_value" automatically reverse the permutation between the "natural" order on the ring and the "folded" order of the MPS, while MPS methods like psi.expectation_value or psi.correlation_function give you their results in the folded order.

[Albert65]

Thank you so much, Johannes. I’ll make sure to follow your recommendations.