# Models¶

## What is a model?¶

Abstractly, a **model** stands for some physical (quantum) system to be described.
For tensor networks algorithms, the model is usually specified as a Hamiltonian written in terms of second quantization.
For example, let us consider a spin-1/2 Heisenberg model described by the Hamiltonian

Note that a few things are defined more or less implicitly.

The local Hilbert space: it consists of Spin-1/2 degrees of freedom with the usual spin-1/2 operators \(S^x, S^y, S^z\).

The geometric (lattice) strucuture: above, we spoke of a 1D “chain”.

The boundary conditions: do we have open or periodic boundary conditions? The “chain” suggests open boundaries, which are in most cases preferable for MPS-based methods.

The range of i: How many sites do we consider (for a 2D system: in each direction)?

Obviously, these things need to be specified in TeNPy in one way or another, if we want to define a model.

Ultimately, our goal is to run some algorithm. However, different algorithm requires the model and Hamiltonian to be specified in different forms.
We have one class for each such required form.
For example `dmrg`

requires an `MPOModel`

,
which contains the Hamiltonian written as an `MPO`

.
So a new model class suitable for DMRG should have this general structure:

```
class MyNewModel(MPOModel):
def __init__(self, model_params):
lattice = somehow_generate_lattice(model_params)
H_MPO = somehow_generate_MPO(lattice, model_params)
# initialize MPOModel
MPOModel.__init__(self, lattice, H_MPO)
```

On the other hand, if we want to evolve a state with `tebd`

we need a `NearestNeighborModel`

, in which the Hamiltonian is written in terms of
two-site bond-terms to allow a Suzuki-Trotter decomposition of the time-evolution operator:

```
class MyNewModel2(NearestNeighborModel):
"""General strucutre for a model suitable for TEBD."""
def __init__(self, model_params):
lattice = somehow_generate_lattice(model_params)
H_bond = somehow_generate_H_bond(lattice, model_params)
# initialize MPOModel
NearestNeighborModel.__init__(self, lattice, H_bond)
```

Of course, the difficult part in these examples is to generate the `H_MPO`

and `H_bond`

in the required form.
If you want to write it down by hand, you can of course do that.
But it can be quite tedious to write every model multiple times, just because we need different representations of the same Hamiltonian.
Luckily, there is a way out in TeNPy: the `CouplingModel`

. Before we describe this class, let’s
discuss the background of the `Site`

and `Site`

class.

## The Hilbert space¶

The **local Hilbert** space is represented by a `Site`

(read its doc-string!).
In particular, the Site contains the local `LegCharge`

and hence the meaning of each
basis state needs to be defined.
Beside that, the site contains the local operators - those give the real meaning to the local basis.
Having the local operators in the site is very convenient, because it makes them available by name for example when you want to calculate expectation values.
The most common sites (e.g. for spins, spin-less or spin-full fermions, or bosons) are predefined
in the module `tenpy.networks.site`

, but if necessary you can easily extend them
by adding further local operators or completely write your own subclasses of `Site`

.

The full Hilbert space is a tensor product of the local Hilbert space on each site.

Note

The `LegCharge`

of all involved sites need to have a common
`ChargeInfo`

in order to allow the contraction of tensors acting on the various sites.
This can be ensured with the function `set_common_charges()`

.

An example where `set_common_charges()`

is needed would be a coupling of different
types of sites, e.g., when a tight binding chain of fermions is coupled to some local spin degrees of freedom.
Another use case of this function would be a model with a $U(1)$ symmetry involving only half the sites, say \(\sum_{i=0}^{L/2} n_{2i}\).

Note

If you don’t know about the charges and np_conserved yet, but want to get started with models right away,
you can set `conserve=None`

in the existing sites or use
`leg = tenpy.linalg.np_conserved.LegCharge.from_trivial(d)`

for an implementation of your custom site,
where d is the dimension of the local Hilbert space.
Alternatively, you can find some introduction to the charges in the Charge conservation with np_conserved.

## The geometry : lattice class¶

The geometry is usually given by some kind of **lattice** structure how the sites are arranged,
e.g. implicitly with the sum over nearest neighbours \(\sum_{<i, j>}\).
In TeNPy, this is specified by a `Lattice`

class, which contains a unit cell of
a few `Site`

which are shifted periodically by its basis vectors to form a regular lattice.
Again, we have pre-defined some basic lattices like a `Chain`

,
two chains coupled as a `Ladder`

or 2D lattices like the
`Square`

, `Honeycomb`

and
`Kagome`

lattices; but you are also free to define your own generalizations.

MPS based algorithms like DMRG always work on purely 1D systems. Even if our model “lives” on a 2D lattice,
these algorithms require to map it onto a 1D chain (probably at the cost of longer-range interactions).
This mapping is also done by the lattice by defining the **order** (`order`

) of the sites.

Note

Further details on the lattice geometry can be found in Details on the lattice geometry.

## 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 forumla
\(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:

Read out the parameters.

Given the parameters, determine the charges to be conserved. Initialize the

`LegCharge`

of the local sites accordingly.Define (additional) local operators needed.

Initialize the needed

`Site`

.Note

Using pre-defined sites like the

`SpinHalfSite`

is recommended and can replace steps 1-3.Initialize the lattice (or if you got the lattice as a parameter, set the sites in the unit cell).

Initialize the

`CouplingModel`

with`CouplingModel.__init__(self, lat)`

.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 multipliy 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.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())`

.Similarly, if you derived from the

`NearestNeighborModel`

, you can call`calc_H_bond()`

to initialze 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 2018-2021 TeNPy Developers, GNU GPLv3
import numpy as np
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.
bc_MPS : {'finite' | 'infinte'}
MPS boundary conditions. Coupling boundary conditions are chosen appropriately.
"""
def __init__(self, model_params):
# 0) read out/set default parameters
model_params = asConfig(model_params, "XXZChain")
L = model_params.get('L', 2)
Jxx = model_params.get('Jxx', 1.)
Jz = model_params.get('Jz', 1.)
hz = model_params.get('hz', 0.)
bc_MPS = model_params.get('bc_MPS', 'finite')
# 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'], 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')
# 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):
return SpinHalfSite(conserve='Sz') # use predefined Site
def init_terms(self, model_params):
# read out parameters
Jxx = model_params.get('Jxx', 1.)
Jz = model_params.get('Jz', 1.)
hz = model_params.get('hz', 0.)
# 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 enure 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 Neigbors 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 2018-2021 TeNPy Developers, GNU GPLv3
import numpy as np
from .model import CouplingMPOModel, NearestNeighborModel
from .lattice import Chain
from ..tools.params import asConfig
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`.
J, g : float | array
Coupling as defined for the Hamiltonian above.
"""
def init_sites(self, model_params):
conserve = model_params.get('conserve', 'parity')
assert conserve != 'Sz'
if conserve == 'best':
conserve = 'parity'
self.logger.info("%s: set conserve to %s", self.name, conserve)
site = SpinHalfSite(conserve=conserve)
return site
def init_terms(self, model_params):
J = np.asarray(model_params.get('J', 1.))
g = np.asarray(model_params.get('g', 1.))
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:

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:

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.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 “couplig” 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]\),

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 generete 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`CouplingMPOModel`

. For example an exponentially decaying long-range interactions are not supported by the coupling model but straight-forward to include to an MPO, as demonstrated in the example`examples/mpo_exponentially_decaying.py`

. The`AKLTChain`

is another example which is directly constructed from the H_bond terms.If you want to debug or double check that the you added the correct terms to a

`CouplingMPOModel`

, you can print the coupling terms with`print(M.all_coupling_terms().to_TermList())`

, and the onsite terms with`print(M.all_onsite_terms().to_TermList())`

. More specifically, you can take only the terms from some categories, e.g. for the`TFIChain`

, you could`print(M.coupling_terms['Sigmax_i Sigmax_j'].to_TermList())`

.If the model of your interest contains Fermions, you should read the Fermions and the Jordan-Wigner transformation.

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`

.