Details on the implementation of Models

In this chapter, we provide some more detail on how models work, and how you might customize them. You should probably read Models (Introduction) first.

We distinguish three different ways in which the Hamiltonian can be given, and there is a base class for each one of them:

  1. In a NearestNeighborModel, we have the Hamiltonian as a sum of two-body terms, which are stored explicitly as a list of Arrays. This is the structure you need to do TEBD with the model.

  2. In a MPOModel, we have the Hamiltonian directly given as a MPO. This is the structure you need to do DMRG, ExpMPOEvolution or TDVP.

  3. In a CouplingModel, the Hamiltonian is given symbolically, in the form of terms (see terms). There are (currently) no algorithms in TeNPy that require this particular structure. We can view it more as a convenient way to specify models, which also allows us to initialize the other two structures easily.

A custom model (as well as the pre-defined models in TeNPy) should then inherit from all of the classes that are applicable.

If you define a CouplingModel structure for the model, that class offers convenient methods to initialize the other two structures, as shown in more detail below. There is a convenience class that achieves this directly, the CouplingMPOModel. It uses the same symbolical representation of the Hamiltonian, but in contrast to the plain CouplingModel, automates the initialization of the lattice and of the MPO. It also automatically initializes H_bond, if it detects that the custom model is also a subclass of NearestNeighborModel. This means that there is virtually no explicit code needed, e.g. when the TFIModel is specialized to the TFIChain.

In the rest of this intro, we introduce the classes and their ways of initializing models in more detail.

The CouplingModel: general structure

The CouplingModel provides a general, quite abstract way to specify a Hamiltonian of couplings on a given lattice. Once initialized, its methods add_onsite() and add_coupling() allow to add onsite and coupling terms repeated over the different unit cells of the lattice. In that way, it basically allows a straight-forward translation of the Hamiltonian given as a math formula \(H = \sum_{i} A_i B_{i+dx} + ...\) with onsite operators A, B,… into a model class.

The general structure for a new model based on the CouplingModel is then:

class MyNewModel3(CouplingModel,MPOModel,NearestNeighborModel):
    def __init__(self, ...):
        ...  # follow the basic steps explained below

In the initialization method __init__(self, ...) of this class you can then follow these basic steps:

  1. Read out the parameters.

  2. Given the parameters, determine the charges to be conserved. Initialize the LegCharge of the local sites accordingly.

  3. Define (additional) local operators needed.

  4. Initialize the needed Site.

    Note

    Using pre-defined sites like the SpinHalfSite is recommended and can replace steps 1-3.

  5. Initialize the lattice (or if you got the lattice as a parameter, set the sites in the unit cell).

  6. Initialize the CouplingModel with CouplingModel.__init__(self, lat).

  7. Use add_onsite() and add_coupling() to add all terms of the Hamiltonian. Here, the pairs of the lattice can come in handy, for example:

    self.add_onsite(-np.asarray(h), 0, 'Sz')
for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
    self.add_coupling(0.5*J, u1, 'Sp', u2, 'Sm', dx, plus_hc=True)
    self.add_coupling(    J, u1, 'Sz', u2, 'Sz', dx)

    Note

    The method add_coupling() adds the coupling only in one direction, i.e. not switching i and j in a \(\sum_{\langle i, j\rangle}\). If you have terms like \(c^\dagger_i c_j\) or \(S^{+}_i S^{-}_j\) in your Hamiltonian, you need to add it in both directions to get a Hermitian Hamiltonian! The easiest way to do that is to use the plus_hc option of add_onsite() and add_coupling(), as we did for the \(J/2 (S^{+}_i S^{-}_j + h.c.)\) terms of the Heisenberg model above. Alternatively, you can add the hermitian conjugate terms explicitly, see the examples in add_coupling() for more details.

    Note that the strength arguments of these functions can be (numpy) arrays for site-dependent couplings. If you need to add or multiply some parameters of the model for the strength of certain terms, it is recommended use np.asarray beforehand – in that way lists will also work fine.

  8. Finally, if you derived from the MPOModel, you can call calc_H_MPO() to build the MPO and use it for the initialization as MPOModel.__init__(self, lat, self.calc_H_MPO()).

  9. Similarly, if you derived from the NearestNeighborModel, you can call calc_H_bond() to initialize it as NearestNeighborModel.__init__(self, lat, self.calc_H_bond()). Calling self.calc_H_bond() will fail for models which are not nearest-neighbors (with respect to the MPS ordering), so you should only subclass the NearestNeighborModel if the lattice is a simple Chain.

Note

The method add_coupling() works only for terms involving operators on 2 sites. If you have couplings involving more than two sites, you can use the add_multi_coupling() instead. A prototypical example is the exactly solvable ToricCode.

The code of the module tenpy.models.xxz_chain is included below as an illustrative example how to implement a Model. The implementation of the XXZChain directly follows the steps outline above. The XXZChain2 implements the very same model, but based on the CouplingMPOModel explained in the next section.

"""Prototypical example of a 1D quantum model: the spin-1/2 XXZ chain.

The XXZ chain is contained in the more general :class:`~tenpy.models.spins.SpinChain`; the idea of
this module is more to serve as a pedagogical example for a model.
"""
# Copyright (C) TeNPy Developers, Apache license

from .lattice import Site, Chain
from .model import CouplingModel, NearestNeighborModel, MPOModel, CouplingMPOModel
from ..linalg import np_conserved as npc
from ..tools.params import asConfig
from ..networks.site import SpinHalfSite  # if you want to use the predefined site

__all__ = ['XXZChain', 'XXZChain2']


class XXZChain(CouplingModel, NearestNeighborModel, MPOModel):
    r"""Spin-1/2 XXZ chain with Sz conservation.

    The Hamiltonian reads:

    .. math ::
        H = \sum_i \mathtt{Jxx}/2 (S^{+}_i S^{-}_{i+1} + S^{-}_i S^{+}_{i+1})
                 + \mathtt{Jz} S^z_i S^z_{i+1} \\
            - \sum_i \mathtt{hz} S^z_i

    All parameters are collected in a single dictionary `model_params`, which
    is turned into a :class:`~tenpy.tools.params.Config` object.

    Parameters
    ----------
    model_params : :class:`~tenpy.tools.params.Config`
        Parameters for the model. See :cfg:config:`XXZChain` below.

    Options
    -------
    .. cfg:config :: XXZChain
        :include: CouplingMPOModel

        L : int
            Length of the chain.
        Jxx, Jz, hz : float | array
            Coupling as defined for the Hamiltonian above.
            Defaults to ``Jxx=Jz=1`` without field ``hz=0``.
        bc_MPS : {'finite' | 'infinite'}
            MPS boundary conditions. Coupling boundary conditions are chosen appropriately.
        sort_charge : bool
            Whether to sort by charges of physical legs. `True` by default.

    """
    def __init__(self, model_params):
        # 0) read out/set default parameters
        model_params = asConfig(model_params, "XXZChain")
        L = model_params.get('L', 2, int)
        Jxx = model_params.get('Jxx', 1., 'real_or_array')
        Jz = model_params.get('Jz', 1., 'real_or_array')
        hz = model_params.get('hz', 0., 'real_or_array')
        bc_MPS = model_params.get('bc_MPS', 'finite', str)
        sort_charge = model_params.get('sort_charge', True, bool)
        # 1-3):
        USE_PREDEFINED_SITE = False
        if not USE_PREDEFINED_SITE:
            # 1) charges of the physical leg. The only time that we actually define charges!
            leg = npc.LegCharge.from_qflat(npc.ChargeInfo([1], ['2*Sz']), [1, -1])
            # 2) onsite operators
            Sp = [[0., 1.], [0., 0.]]
            Sm = [[0., 0.], [1., 0.]]
            Sz = [[0.5, 0.], [0., -0.5]]
            # (Can't define Sx and Sy as onsite operators: they are incompatible with Sz charges.)
            # 3) local physical site
            site = Site(leg, ['up', 'down'], sort_charge=sort_charge, Sp=Sp, Sm=Sm, Sz=Sz)
        else:
            # there is a site for spin-1/2 defined in TeNPy, so just we can just use it
            # replacing steps 1-3)
            site = SpinHalfSite(conserve='Sz', sort_charge=sort_charge)
        # 4) lattice
        bc = 'open' if bc_MPS == 'finite' else 'periodic'
        lat = Chain(L, site, bc=bc, bc_MPS=bc_MPS)
        # 5) initialize CouplingModel
        CouplingModel.__init__(self, lat)
        # 6) add terms of the Hamiltonian
        # (u is always 0 as we have only one site in the unit cell)
        self.add_onsite(-hz, 0, 'Sz')
        self.add_coupling(Jxx * 0.5, 0, 'Sp', 0, 'Sm', 1, plus_hc=True)
        # the `plus_hc=True` adds the h.c. term
        # see also the examples tenpy.models.model.CouplingModel.add_coupling
        self.add_coupling(Jz, 0, 'Sz', 0, 'Sz', 1)
        # 7) initialize H_MPO
        MPOModel.__init__(self, lat, self.calc_H_MPO())
        # 8) initialize H_bond (the order of 7/8 doesn't matter)
        NearestNeighborModel.__init__(self, lat, self.calc_H_bond())


class XXZChain2(CouplingMPOModel, NearestNeighborModel):
    """Another implementation of the Spin-1/2 XXZ chain with Sz conservation.

    This implementation takes the same parameters as the :class:`XXZChain`, but is implemented
    based on the :class:`~tenpy.models.model.CouplingMPOModel`.

    Parameters
    ----------
    model_params : dict | :class:`~tenpy.tools.params.Config`
        See :cfg:config:`XXZChain`
    """
    default_lattice = "Chain"
    force_default_lattice = True

    def init_sites(self, model_params):
        sort_charge = model_params.get('sort_charge', True, bool)
        return SpinHalfSite(conserve='Sz', sort_charge=sort_charge)  # use predefined Site

    def init_terms(self, model_params):
        # read out parameters
        Jxx = model_params.get('Jxx', 1., 'real_or_array')
        Jz = model_params.get('Jz', 1., 'real_or_array')
        hz = model_params.get('hz', 0., 'real_or_array')
        # add terms
        for u in range(len(self.lat.unit_cell)):
            self.add_onsite(-hz, u, 'Sz')
        for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
            self.add_coupling(Jxx * 0.5, u1, 'Sp', u2, 'Sm', dx, plus_hc=True)
            self.add_coupling(Jz, u1, 'Sz', u2, 'Sz', dx)

The easiest way: the CouplingMPOModel

Since many of the basic steps above are always the same, we don’t need to repeat them all the time. So we have yet another class helping to structure the initialization of models: the CouplingMPOModel. The general structure of this class is like this:

class CouplingMPOModel(CouplingModel,MPOModel):
    default_lattice = "Chain"
    "

    def __init__(self, model_param):
        # ... follows the basic steps 1-8 using the methods
        lat = self.init_lattice(self, model_param)  # for step 4
        # ...
        self.init_terms(self, model_param) # for step 6
        # ...

    def init_sites(self, model_param):
        # You should overwrite this in most cases to ensure
        # getting the site(s) and charge conservation you want
        site = SpinSite(...)  # or FermionSite, BosonSite, ...
        return site  # (or tuple of sites)

    def init_lattice(self, model_param):
        sites = self.init_sites(self, model_param) # for steps 1-3
        # and then read out the class attribute `default_lattice`,
        # initialize an arbitrary pre-defined lattice
        # using model_params['lattice']
        # and ensure it's the default lattice if the class attribute
        # `force_default_lattice` is True.

    def init_terms(self, model_param):
        # does nothing.
        # You should overwrite this

The XXZChain2 included above illustrates, how it can be used. You need to implement steps 1-3) by overwriting the method init_sites() Step 4) is performed in the method init_lattice(), which initializes arbitrary 1D or 2D lattices; by default a simple 1D chain. If your model only works for specific lattices, you can overwrite this method in your own class. Step 6) should be done by overwriting the method init_terms(). Steps 5,7,8 and calls to the init_… methods for the other steps are done automatically if you just call the CouplingMPOModel.__init__(self, model_param).

The XXZChain and XXZChain2 work only with the Chain as lattice, since they are derived from the NearestNeighborModel. This allows to use them for TEBD in 1D (yeah!), but we can’t get the MPO for DMRG on (for example) a Square lattice cylinder - although it’s intuitively clear, what the Hamiltonian there should be: just put the nearest-neighbor coupling on each bond of the 2D lattice.

It’s not possible to generalize a NearestNeighborModel to an arbitrary lattice where it’s no longer nearest Neighbors in the MPS sense, but we can go the other way around: first write the model on an arbitrary 2D lattice and then restrict it to a 1D chain to make it a NearestNeighborModel.

Let me illustrate this with another standard example model: the transverse field Ising model, implemented in the module tenpy.models.tf_ising included below. The TFIModel works for arbitrary 1D or 2D lattices. The TFIChain is then taking the exact same model making a NearestNeighborModel, which only works for the 1D chain.

"""Prototypical example of a quantum model: the transverse field Ising model.

Like the :class:`~tenpy.models.xxz_chain.XXZChain`, the transverse field ising chain
:class:`TFIChain` is contained in the more general :class:`~tenpy.models.spins.SpinChain`;
the idea is more to serve as a pedagogical example for a 'model'.

We choose the field along z to allow to conserve the parity, if desired.
"""
# Copyright (C) TeNPy Developers, Apache license

import numpy as np

from .model import CouplingMPOModel, NearestNeighborModel
from .lattice import Chain
from ..networks.site import SpinHalfSite

__all__ = ['TFIModel', 'TFIChain']


class TFIModel(CouplingMPOModel):
    r"""Transverse field Ising model on a general lattice.

    The Hamiltonian reads:

    .. math ::
        H = - \sum_{\langle i,j\rangle, i < j} \mathtt{J} \sigma^x_i \sigma^x_{j}
            - \sum_{i} \mathtt{g} \sigma^z_i

    Here, :math:`\langle i,j \rangle, i< j` denotes nearest neighbor pairs, each pair appearing
    exactly once.
    All parameters are collected in a single dictionary `model_params`, which
    is turned into a :class:`~tenpy.tools.params.Config` object.

    Parameters
    ----------
    model_params : :class:`~tenpy.tools.params.Config`
        Parameters for the model. See :cfg:config:`TFIModel` below.

    Options
    -------
    .. cfg:config :: TFIModel
        :include: CouplingMPOModel

        conserve : None | 'parity'
            What should be conserved. See :class:`~tenpy.networks.Site.SpinHalfSite`.
        sort_charge : bool
            Whether to sort by charges of physical legs. `True` by default.
        J, g : float | array
            Coupling as defined for the Hamiltonian above.
            Defaults to ``J=g=1``

    """
    def init_sites(self, model_params):
        conserve = model_params.get('conserve', 'parity', str)
        assert conserve != 'Sz'
        if conserve == 'best':
            conserve = 'parity'
            self.logger.info("%s: set conserve to %s", self.name, conserve)
        sort_charge = model_params.get('sort_charge', True, bool)
        site = SpinHalfSite(conserve=conserve, sort_charge=sort_charge)
        return site

    def init_terms(self, model_params):
        J = np.asarray(model_params.get('J', 1., 'real_or_array'))
        g = np.asarray(model_params.get('g', 1., 'real_or_array'))
        for u in range(len(self.lat.unit_cell)):
            self.add_onsite(-g, u, 'Sigmaz')
        for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
            self.add_coupling(-J, u1, 'Sigmax', u2, 'Sigmax', dx)
        # done


class TFIChain(TFIModel, NearestNeighborModel):
    """The :class:`TFIModel` on a Chain, suitable for TEBD.

    See the :class:`TFIModel` for the documentation of parameters.
    """
    default_lattice = Chain
    force_default_lattice = True

Automation of Hermitian conjugation

As most physical Hamiltonians are Hermitian, these Hamiltonians are fully determined when only half of the mutually conjugate terms is defined. For example, a simple Hamiltonian:

\[H = \sum_{\langle i,j\rangle, i<j} - \mathtt{J} (c^{\dagger}_i c_j + c^{\dagger}_j c_i)\]

is fully determined by the term \(c^{\dagger}_i c_j\) if we demand that Hermitian conjugates are included automatically. In TeNPy, whenever you add a coupling using add_onsite(), add_coupling(), or add_multi_coupling(), you can use the optional argument plus_hc to automatically create and add the Hermitian conjugate of that coupling term - as shown above.

Additionally, in an MPO, explicitly adding both a non-Hermitian term and its conjugate increases the bond dimension of the MPO, which increases the memory requirements of the MPOEnvironment. Instead of adding the conjugate terms explicitly, you can set a flag explicit_plus_hc in the MPOCouplingModel parameters, which will ensure two things:

  1. The model and the MPO will only store half the terms of each Hermitian conjugate pair added, but the flag explicit_plus_hc indicates that they represent self + h.c.. In the example above, only the term \(c^{\dagger}_i c_j\) would be saved.

  2. At runtime during DMRG, the Hermitian conjugate of the (now non-Hermitian) MPO will be computed and applied along with the MPO, so that the effective Hamiltonian is still Hermitian.

Note

The model flag explicit_plus_hc should be used in conjunction with the flag plus_hc in add_coupling() or add_multi_coupling(). If plus_hc is False while explicit_plus_hc is True the MPO bond dimension will not be reduced, but you will still pay the additional computational cost of computing the Hermitian conjugate at runtime.

Thus, we end up with several use cases, depending on your preferences. Consider the FermionModel. If you do not care about the MPO bond dimension, and want to add Hermitian conjugate terms manually, you would set model_par[‘explicit_plus_hc’] = False and write:

self.add_coupling(-J, u1, 'Cd', u2, 'C', dx)
self.add_coupling(np.conj(-J), u2, 'Cd', u1, 'C', -dx)

If you wanted to save the trouble of the extra line of code (but still did not care about MPO bond dimension), you would keep the model_par, but instead write:

self.add_coupling(-J, u1, 'Cd', u2, 'C', dx, plus_hc=True)

Finally, if you wanted a reduction in MPO bond dimension, you would need to set model_par[‘explicit_plus_hc’] = True, and write:

self.add_coupling(-J, u1, 'Cd', u2, 'C', dx, plus_hc=True)

Non-uniform terms and couplings

The CouplingModel-methods add_onsite(), add_coupling(), and add_multi_coupling() add a sum over a “coupling” term shifted by lattice vectors. However, some models are not that “uniform” over the whole lattice.

First of all, you might have some local term that gets added only at one specific location in the lattice. You can add such a term for example with add_local_term().

Second, if you have irregular lattices, take a look at the corresponding section in Details on the lattice geometry.

Finally, note that the argument strength for the add_onsite, add_coupling, and add_multi_coupling methods can not only be a numpy scalar, but also a (numpy) array. In general, the sum performed by the methods runs over the given term shifted by lattice vectors as far as possible to still fit the term into the lattice.

For the add_onsite() case this criterion is simple: there is exactly one site in each lattice unit cell with the u specified as separate argument, so the correct shape for the strength array is simply given by Ls. For example, if you want the defacto standard model studied for many-body localization, a Heisenberg chain with random , uniform onsite field \(h^z_i \in [-W, W]\),

\[H = J \sum_{i=0}^{L-1} \vec{S}_i \cdot \vec{S}_{i+1} - \sum_{i=0}^{L} h^z_i S^z_i\]

you can use the SpinChain with the following model parameters:

L = 30 # or whatever you like...
W = 5.  # MBL transition at W_c ~= 3.5 J
model_params = {
    'L': L,
    'Jx': 1., 'Jy': 1., 'Jz': 1.,
    'hz': 2.*W*(np.random.random(L) - 0.5),  # random values in [-W, W], shape (L,)
    'conserve': 'best',
}
M = tenpy.models.spins.SpinChain(model_params)

For add_coupling() and add_multi_coupling(), things become a little bit more complicated, and the correct shape of the strength array depends not only on the Ls but also on the boundary conditions of the lattice. Given a term, you can call coupling_shape() and multi_coupling_shape() to find out the correct shape for strength. To avoid any ambiguity, the shape of the strength always has to fit, at least after a tiling performed by to_array().

For example, consider the Su-Schrieffer-Heeger model, a spin-less FermionChain with hopping strength alternating between two values, say t1 and t2. You can generate this model for example like this:

L = 30 # or whatever you like...
t1, t2 = 0.5, 1.5
t_array = np.array([(t1 if i % 2 == 0 else t2) for i in range(L-1)])
model_params = {
    'L': L,
    't': t_array,
    'V': 0., 'mu': 0.,  # just free fermions, but you can generalize...
    'conserve': 'best'
}
M = tenpy.models.fermions.FermionChain(model_params)

Some random remarks on models

  • Needless to say that we have also various predefined models under tenpy.models.

  • If you want to use random parameters, you should use model.rng as a random number generator and change model_params['random_seed'] for different configurations.

  • Of course, an MPO is all you need to initialize a MPOModel to be used for DMRG; you don’t have to use the CouplingModel or CouplingMPOModel. For example, we build the model directly from an MPO in examples/mpo_exponentially_decaying.py. The AKLTChain is another example which is directly constructed from the H_bond terms.

  • We suggest writing the model to take a single parameter dictionary for the initialization, as the CouplingMPOModel does. The CouplingMPOModel converts the dictionary to a dict-like Config with some additional features before passing it on to the init_lattice, init_site, … methods. It is recommended to read out providing default values with model_params.get("key", default_value), see get().

  • When you write a model and want to include a test that it can be at least constructed, take a look at tests/test_model.py.