Skip to content

Building a Hamiltonian

A Hamiltonian in symtenet is an MPO. You build one from local operators, and tenet.models ships the standard sites so you do not write the same three matrices again — or, for the models below, the finished operator.

Named models

One function per model, each returning the MPO for an open chain of n sites, with the symmetry a parameter:

>>> from tenet.models import heisenberg, hubbard, sun_heisenberg
>>> from tenet.symmetry import SU2
>>> h = heisenberg(20)                  # U(1), J = 1, spin 1/2
>>> h_su2 = heisenberg(20, SU2)         # the same chain, one invariant term per bond
>>> len(hubbard(8, t=1.0, U=4.0)), len(sun_heisenberg(8, 3))
(8, 8)
call Hamiltonian grading
heisenberg(n, symmetry=U1, J=1.0, spin=0.5) \(J \sum_i \vec{S}_i \cdot \vec{S}_{i+1}\) U1 (default) or SU2
xxz(n, Delta=1.0, J=1.0) \(J \sum_i S^x S^x + S^y S^y + \Delta S^z S^z\) U1
transverse_field_ising(n, J=1.0, g=1.0) \(-J\left(\sum_i \sigma^z_i \sigma^z_{i+1} + g \sum_i \sigma^x_i\right)\) Z2
hubbard(n, t=1.0, U=4.0) \(-t \sum_{i\sigma} (c^\dagger_{i\sigma} c_{i+1\sigma} + h.c.) + U \sum_i n_{i\uparrow} n_{i\downarrow}\) fZ2
spinless_tv(n, t=1.0, V=1.0) \(-t \sum_i (c^\dagger_i c_{i+1} + h.c.) + V \sum_i n_i n_{i+1}\) fZ2
sun_heisenberg(n, N, J=1.0) \(J \sum_i P_{i,i+1}\) on SU(\(N\)) fundamentals SU(\(N\))

Each takes symbolic=, passed to the builder below. Every function's docstring carries its sign and ordering conventions — the basis a grading puts the site in is part of the model, so read it before comparing a number against another library. The rest of this page is what to do when your Hamiltonian is not one of these, which is the usual case.

Sites

A Site carries three things:

field what it is
phys the physical GradedSpace, ready for MPS.product / MPS.random
ops name to term operator, the rank-3 charge-leg form built by local_op
matrices the dense matrix behind each, for forms the term API does not build 
>>> from tenet.models import spin_half
>>> site = spin_half()
>>> site.phys.sectors            # U(1), the charge is 2 S^z
((U1Sector(charge=-1), 1), (U1Sector(charge=1), 1))
>>> sorted(site.ops)
['S+', 'S-', 'Sz']

The shipped sites:

call grading ops
spin_half() U1, charge \(2S^z\) Sz, S+, S- (S.S as a matrix)
spin_half(SU2) SU2, one \(j = 1/2\) multiplet S.S
spinless_fermion() fZ2 c, c+, n
spinful_fermion() fZ2, d=4 c_up, c+_up, c_dn, c+_dn, n_up, n_dn, n, n_up n_dn
hard_core_boson() U1, the occupation b, b+, n
hard_core_boson(Trivial) ungraded b, b+, n

Any other provider is refused by name. tenet.network never imports tenet.models: the driver layer decides nothing about what your operators mean, and this layer is a convenience above it. The library ships sites, which are a finite set; the lattice and the Hamiltonian stay yours.

Each entry of ops is rank 3: the two physical legs plus a D=1 leg carrying the charge the operator emits. That third leg is what makes S+ expressible as a term operator. Write your own with local_op, which takes any dense matrix:

>>> import numpy as np
>>> from tenet.network import local_op
>>> from tenet.symmetry import U1Sector
>>> sp = local_op(np.array([[0.0, 0.0], [1.0, 0.0]]), phys=site.phys, charge=U1Sector(-2))
>>> sp.ndim
3

With no charge, local_op takes an invariant k-site operator instead — a (d**k, d**k) or (d,)*2k array, the layout np.kron produces.

The four builders

you have use
a list of terms MPO.from_terms
many terms over a few patterns MPO.from_arrays
a W you are writing by hand MPO.from_entries
a W that is already a dense array MPO.from_w

The first three derive the MPO bond spaces from the operators' own charges. You declare no grading and no bond width.

from_terms — a list of terms

(coefficient, [(operator, site), ...]) tuples, identity implied on every untouched site:

>>> from tenet.network import MPO
>>> n = 8
>>> terms = []
>>> for i in range(n - 1):
...     terms.append((1.0, [(site.ops["Sz"], i), (site.ops["Sz"], i + 1)]))
...     terms.append((0.5, [(site.ops["S+"], i), (site.ops["S-"], i + 1)]))
...     terms.append((0.5, [(site.ops["S-"], i), (site.ops["S+"], i + 1)]))
>>> h = MPO.from_terms(n, terms)
>>> len(h)
8

Every term is a bond-1 MPO, direct_sum stacks them, and two compressing SVD sweeps collapse the stack to the operator Schmidt rank. cutoff= sets those sweeps' threshold; cutoff=None skips the compression and gives the exact finite-state-machine bond.

A term operator may be an invariant k-site operator, in which case its sites are given as a tuple: (1.0, [(ss, (i, i + 1))]). from_terms splits it with an SVD, and the graded MPO bond comes out of that split.

Refusal to know: two operators of one term on one site. Multiply them into a single on-site operator first.

from_arrays — patterns and index arrays

One (expr, indices, data) triple per operator pattern. expr names the operators, indices is a (T, L) integer array of sites, data is the length-T coefficient array:

>>> bonds = [(i, i + 1) for i in range(n - 1)]
>>> bonds[:3]
[(0, 1), (1, 2), (2, 3)]

bonds is the index array, and one row of it is one term's site list — one entry per operator named in expr, in the same order. ("Sz Sz", bonds, ...) therefore says: the first Sz sits on the row's first site, the second on its second, and there are n - 1 such terms. A three-operator pattern would need rows of three sites; a two-operator pattern acting twice on one site has rows like (i, i).

>>> h2 = MPO.from_arrays(n, site.ops, [
...     ("Sz Sz", bonds, [1.0] * (n - 1)),
...     ("S+ S-", bonds, [0.5] * (n - 1)),
...     ("S- S+", bonds, [0.5] * (n - 1)),
... ])
>>> len(h2)
8

Three patterns share the one bonds array, each with its own length-n - 1 coefficient array. The coefficients' dtype decides the MPO's, so a complex coupling makes a complex operator without anything else being said.

The pattern's work is done once per block in NumPy instead of once per term in Python, which is what a Hamiltonian with \(O(K^4)\) terms over a handful of patterns needs. Before the walk, three whole-array steps run: each row is sorted into site order paying the Koszul sign of every inversion of two sign-braiding operators; operators coinciding on a site are pre-multiplied into one on-site operator, and a term whose on-site product vanishes is dropped; terms agreeing on (operator labels, sites) are fused and their coefficients summed. screen= then drops merged terms below a magnitude.

Permutational symmetry is expanded by you and merged here: the eight images of \((ij \vert kl)\) are eight different operator strings.

Every operator here is rank 3 — a block gives one site index per name, so an invariant k-site operator has nowhere to put its extra indices and is refused, naming from_terms.

from_entries — a W by hand

The non-zero (i, j) entries of each site's W, one mapping per site:

>>> sz, sp2, sm = site.ops["Sz"], site.ops["S+"], site.ops["S-"]
>>> w = {                      # the Heisenberg W, its eight non-zero channels
...     (0, 0): None,          # I -- the term has not started
...     (0, 1): (0.5, sm), (1, -1): sp2,
...     (0, 2): (0.5, sp2), (2, -1): sm,
...     (0, 3): sz,        (3, -1): sz,
...     (-1, -1): None,        # I -- the term is finished
... }
>>> h3 = MPO.from_entries([w] * n)
>>> len(h3)
8

0 is a bond's IdL channel and -1 its IdR channel, at every bond, the way a lower-triangular W is printed; the open channels are 1, 2, .... An entry is None (the identity, a spectator ride on (i, i)), a number (that multiple of it), an operator, or the pair (coefficient, operator). -1 names the last index, so no bond width is declared; the charge is already on local_op's third leg, so each channel's GradedSpace is derived and the bond at a cut is the direct sum over its channels. The two boundary bonds are D = 1, so bond 0 keeps only IdL and the last bond only IdR, and the same bulk mapping serves the first and last site.

from_w — a dense W

MPO.from_w takes one dense bulk W array plus the physical space, the MPO bond space and the boundary space. Use it when the W arrives as a dense array out of a paper or another library: its entries are numbers, no charge can be recovered from them, so you supply the bond grading yourself. A wrong grading makes the build raise, and that refusal is the proof the grading is right.

symbolic=True

from_terms, from_arrays and from_entries each take symbolic=. It decides which engine Env.heff2 runs.

At the default, the builder hands back the site tensors, and heff2 contracts them. That is what a finite-range lattice model wants: its MPO bond is five or eight wide, and the site-tensor contraction is the cheapest thing that can happen to it.

symbolic=True keeps the finite-state-machine description the terms were assembled into, and heff2 runs the prepared matvec: complementary operators assembled per bond, the sum dispatched term family by term family. That is what a Hamiltonian with a bond in the thousands needs — in practice, quantum chemistry, where it is the route that fits in memory at all.

h = MPO.from_terms(n, terms, symbolic=True)

Nothing is decided at run time and no threshold is probed. The representation the operator is in when you hand it to dmrg_ is the choice, made at build time exactly as cutoff is. The two are independent: cutoff=None with no symbolic gives exact, uncompressed site tensors.

MPO.materialize() takes an operator built symbolic=True to the site tensors — the same operator, the description dropped:

h.materialize()

Fermions

The spinful site is the \(d = 4\) basis \((\lvert 0\rangle, \lvert ud\rangle, \lvert u\rangle, \lvert d\rangle)\), the even sector first, because a dense array over a GradedSpace is laid out sector by sector. There is no Jordan-Wigner operator to place: the string is the fZ2 braiding an odd MPO bond pays when it crosses a physical line. A Hubbard chain is its terms:

>>> from tenet.models import spinful_fermion
>>> fsite, m, u = spinful_fermion(), 6, 4.0
>>> fwd = [(i, i + 1) for i in range(m - 1)]
>>> bwd = [(i + 1, i) for i in range(m - 1)]
>>> blocks = []
>>> for flavour in ("up", "dn"):
...     expr = "c+_" + flavour + " c_" + flavour
...     blocks += [(expr, fwd, [-1.0] * (m - 1)), (expr, bwd, [-1.0] * (m - 1))]
>>> blocks.append(("n_up n_dn", [(i, i) for i in range(m)], [u] * m))

Four hopping blocks — two flavours, each in both directions — and one interaction block. Writing the backward hop out explicitly is what makes the operator Hermitian; the Koszul sign of the braid is paid by the builder's own argsort, not by you.

The last block names two operators on two coincident site indices, and from_arrays multiplies them into one on-site operator before the walk. The site also ships that product under the same key, fsite.ops["n_up n_dn"], for from_terms, which places one operator per site and so needs it pre-multiplied.

>>> hub = MPO.from_arrays(m, fsite.ops, blocks)
>>> len(hub)
6

SU(2)

>>> from tenet.symmetry import SU2
>>> su2_site = spin_half(SU2)
>>> sorted(su2_site.ops)
['S.S']
>>> su2_site.ops["S.S"].ndim
4

S+ is not in the table because there is no such SU(2) operator. The rank-3 charge-leg form puts the emitted sector on a D=1 leg, and the only leg a spin-1 tensor operator could emit onto is the \(j = 1\) multiplet, whose dense dimension is 3. What SU(2) has is the invariant two-site operator, and S.S is one whole Heisenberg bond term:

>>> su2_terms = [(1.0, [(su2_site.ops["S.S"], (i, i + 1))]) for i in range(n - 1)]
>>> hsu2 = MPO.from_terms(n, su2_terms)
>>> len(hsu2)
8

from_terms splits it with an SVD and the graded MPO bond comes out of that split — the coupling lives inside the operator's own blocks, so no coupling tree is named. from_arrays cannot express this term, which is why the SU(2) site's ops is a from_terms table.

Applying and checking an operator

>>> from tenet.network import MPS
>>> psi = MPS.product(site.phys, [U1Sector(1 if i % 2 else -1) for i in range(n)])
>>> phi = h.apply(psi)           # H|psi>, exact and untruncated
>>> len(phi)
8

apply takes no chi and no cutoff: the product is exact, and truncation is MPS.compress_, by name. The cost of the exact product is the operator's bond dimension times the state's, so compress promptly on a wide operator.

variance is \(\langle\psi\vert H^2 \vert\psi\rangle/\langle\psi\vert\psi\rangle - E^2\), zero for an exact eigenstate — see DMRG for what to do with it. MPO.identity builds the D=1 identity operator, and to_dense expands the whole \(d^N \times d^N\) operator for a small chain.

Where next