SyTen
Tutorial

This tutorial assumes that you have successfully installed SyTen and are able to execute the resulting binaries. Note that all binaries offer the option --help which will list all options and (sometimes) offer further help.

While not all parts of this tutorial may be directly relevant to all users (especially the time-evolution and the local basis optimisation), it is strongly suggested to at least read all parts which are not obviously not applicable.

File Types

SyTen knows primarily two types of files: lattice and state files. Lattice files represent one particular physical system with one particular Hilbert space and a set of operators defined on that Hilbert space. Lattice files are generating using lattice generators, stored in lat/ and compiled into so-called lattice binaries.

State files represent a particular quantum state. State files are typically generated initially using syten-random and then acted upon with the various other tools in bin/.

Generating a Lattice

A lattice can be generated using any of the binaries compiled from .cpp files in lat/. To start, we will generate a \( \mathrm{U}(1)_{S_z} \) symmetric Heisenberg spin chain of 20 sites:

$ syten-sql-mps-spin -l 20 --sym u1 -o sp.lat

The command options are:

  • -l 20 for 20 sites
  • --sym u1 to only use the \( \mathrm{U}(1)_{S_z} \) symmetry rather than a \( \mathrm{SU}(2)_S \) symmetry (or none at all).
  • -o sp.lat to write the result into the file sp.lat

This generated file sp.lat represents this particular physical system. Calling syten-info on a lattice file will print some overview information:

$ syten-info sp.lat

Crucially, it will give a list of local site operators and global operators. Global operators are pre-defined MPOs (or BTTOs in the case of binary tree tensor networks) which can be used as-is. Local site operators are ‘recipes’ stored inside the lattice on how to generate a MPO representation on the fly. To get an operator out of the lattice, use the syntax <lattice file>:<operator expression>. Note that the whole thing usually has to be quoted for the shell. For details regarding the accepted <operator expression>s, see the FAQ. syten-info can also be called on a particular operator. In this case it will tell you how the operator transforms (per default) or its MPO bond dimensions (with the -b option). syten-print will print the MPO representation to the console.

As an example, check the MPO representations of the sz ( \( \hat s^z \)), sp ( \( \hat s^+ \)) and sm ( \( \hat s^- \)) operators, generated on sites 0, 1 and 10 (sites are zero-indexed) like so:

$ syten-print sp.lat:"sz 0 @"
$ syten-print sp.lat:"sz 1 @"
$ syten-print sp.lat:"sz 10 @"
$ syten-print sp.lat:"sp 0 @"
$ syten-print sp.lat:"sp 1 @"
$ syten-print sp.lat:"sp 10 @"
$ syten-print sp.lat:"sm 0 @"
$ syten-print sp.lat:"sm 1 @"
$ syten-print sp.lat:"sm 10 @"

Any of these commands will give a lot of output, which will be primarily the identity site tensors generated on sites 1 through 19, 0 and 2 through 19 and 0 through 9 and 11 through 19. However, when checking the output carefully, you should be able to find the active site tensors. For the sm 1 @, this is for example the following part:

Site: 1
========================================
Tensor of rank 4 with the following bases:
Leg 1:    Incoming Basis (1 Syms): [<1, [U(1)_z:{+0}@1]>]
Leg 2:    Outgoing Basis (1 Syms): [<1, [U(1)_z:{-1}@1]>]
Leg 3:    Incoming Basis (1 Syms): [<1, [U(1)_z:{-0.5}@1]>]
Leg 4:    Outgoing Basis (1 Syms): [<1, [U(1)_z:{+0.5}@1]>]
Composed of the following 1 blocks:
----------------------------------------
Tensor block of rank 4 and 1 symmetries. isZero = +0.
Reduced tensor is: (+1.0000e+00,+0.0000e+00)
Sym 1:        Idx 1 transforms as: U(1)_z:{+0}@1;    Idx 2 transforms as: U(1)_z:{-1}@1;    Idx 3 transforms as: U(1)_z:{-0.5}@1;    Idx 4 transforms as: U(1)_z:{+0.5}@1
========================================

First, we list the site and that this element is a tensor of rank 4. Then, the bases on each of the four legs (left, right, top, bottom) are enumerated. Each basis consists of a list of quantum number sectors <x, [sym1, ...]>. In this case, it consists of just one sector on each leg with a single element (the leading 1 resp. x) and U(1)_z quantum number 0, -1, -0.5 and +0.5. This is followed by a list of blocks, which transform irreducibly on each leg under every symmetry (i.e. can be assigned good quantum numbers). The tensor contains just a single block with a single 1 element. Compare this tensor to the tensor on the first site of this operator, which is the identity. Instead of mapping the +0.5 state to the -0.5 state as the above tensor does, it maps the -0.5 state to the -0.5 state and the +0.5 state to the +0.5 state.

The <operator expression> understands basic arithmetic operations such as addition or multiplication (both of operators and by scalars). For example,

$ syten-print sp.lat:"sz 0 @ sz 1 @ *"

will multiply together the two operators. sp 0 @ sm 0 @ + will add them together. The operator then does not transform irreducibly any more, resulting in a larger leftmost basis:

$ syten-info sp.lat:"sp 0 @ sm 0 @ +"
2017.04.25T16:09:08.823 5.73650e-02 I Loaded MP operator.
2017.04.25T16:09:08.823 5.75140e-02 I MP-Operator on lattice of 20 sites.
2017.04.25T16:09:08.823 5.75560e-02 I Vacuum basis (right edge): [<1, [U(1)_z:{+0}@1]>]
2017.04.25T16:09:08.823 5.75950e-02 I Sector basis (left edge):  [<1, [U(1)_z:{-1}@1]>, <1, [U(1)_z:{+1}@1]>]
2017.04.25T16:09:08.823 5.76200e-02 I Physical basis on site 1:  [<1, [U(1)_z:{-0.5}@1]>, <1, [U(1)_z:{+0.5}@1]>]

As you can see, the sector basis now contains two elements with \( S_z \) quantum numbers -1 and +1 respectively.

Once you have generated a Hamiltonian, it is a good idea to check that it is hermitian using syten-check-herm. The spin lattice we generated already defines two global operators Hp (transverse nearest-neighbour spin-spin interactions) and Hz (longitudinal nearest-neighbour spin-spin interactions). We can build the standard Heisenberg Hamiltonian from them and check its hermiticity. In comparison, sp 0 @ is certainly not Hermitian:

$ syten-check-herm sp.lat:"Hp Hz +"
2017.04.25T16:12:44.774 3.66320e-02 I Attempting to parse A = 'sp.lat:Hp Hz +'…
2017.04.25T16:12:44.888 8.53640e-02 I Parsed file, calculating norm and |A·A-A*A|…
2017.04.25T16:12:44.942 1.39475e-01 N |A·A - A*A| = 0.00000e+00; |A·A-A*A|/|A| = 0.00000e+00: test passed, A possibly hermitian.
$ syten-check-herm sp.lat:"sp 0 @"
2017.04.25T16:12:51.550 1.95790e-02 I Attempting to parse A = 'sp.lat:sp 0 @'…
2017.04.25T16:12:51.577 4.68140e-02 I Parsed file, calculating norm and |A·A-A*A|…
2017.04.25T16:12:51.602 5.47290e-02 N |A·A - A*A| = 7.24077e+02; |A·A-A*A|/|A| = 1.00000e+00: test failed, A most likely not hermitian.

Generating a Random State

Generating a random state is done using syten-random:

$ syten-random -l sp.lat -s 0 -o rnd.state

where we wish to generate a state in the Hilbert space described by sp.lat with quantum number zero (if there are multiple quantum numbers, separate them with a comma, e.g. -s 5,0.5 in the Fermi-Hubbard case for 5 particles with spin 0.5) and write the result into a file rnd.state. By default, ten vacuum product states are generated and acted upon with random single-site operators until the global state transforms as desired.

If this generation takes a long time, it may be useful to try a different generator, e.g. the complete generator (-g c) or the vacuum generator (-g v). The latter will always generate a vacuum state, though!

Generating a Product State

Remarks
The state produced in this section will not be used later on in the tutorial, so you may skip this section.

If pyten is set up properly, the syten-mps-construct.py script can be used to construct a specific (nearly-) product state. This script reads in a file, each line of this file either describes a bond (odd lines) or a physical leg configuration (even lines). The first line describes the leftmost bond. Bond lines start with a -, site lines with a .. Following this first symbol and separated by a space is a description of which state should be taken on this site or bond. This description takes the form n@Q where Q denotes a quantum number sector and n the state within this sector to use. Alternatively, using * one may allow all states through.

To create a Néel state with the above lattice, first create a state description file:

$ cat state-desc
- 1@0
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5
- *
. 1@+0.5
- *
. 1@-0.5

and then call the script:

$ syten-mps-construct.py -i state-desc -l sp.lat -o neel.state

Applying an Operator to a State

It is possible to apply an operator (given as a MPO) to a state. For example, to apply the Hamiltonian to our random input state:

$ syten-apply-op -i rnd.state -l sp.lat:"Hz Hp +" -o H-rnd.state

Calculating Expectation Values

Given a state and a lattice, it becomes possible to calculate expectation values:

$ syten-expectation rnd.state sp.lat:'Hz Hp +'
2017.04.25T16:21:42.312 2.35330e-02 I Loaded first MP state.
2017.04.25T16:21:42.357 5.62580e-02 I Loaded MP lattice from 'sp.lat', calculating single result…
(-1.3,1.8794721221315e-18)
$ syten-expectation rnd.state sp.lat:'sz 0 @'
2017.04.25T16:21:51.793 2.37490e-02 I Loaded first MP state.
2017.04.25T16:21:51.842 5.69280e-02 I Loaded MP lattice from 'sp.lat', calculating single result…
(-0.1,-3.46944695195362e-18)

With -q, the additional information can be suppressed (only the last line will print). -r only prints the real part, -i only the imaginary part. -v allows the calculation of the variance:

$ syten-expectation -qvr rnd.state sp.lat:'Hz Hp +'
⟨O⟩:          -1.3
⟨O²⟩:         6.1875
|⟨O²⟩-⟨O⟩²|:  4.4975

which in this case is fairly large. Note that calculating the variance (or even \( \langle O^2 \rangle \)) in this way is relatively unstable if your Hamiltonian contains entries with many different orders of magnitude. This is in particular true for bosonic models, with \( N \leq 1000 \), even \( \hat n^2 \) has entries \( O(10^6) \). Now add in a small multiplicator (e.g. \( V_{ijkl} \approx 10^{-5} \)) and the square \( \hat H^2 \) and you suddendly have entries ranging from \( 10^{-10} \) to \( 10^{12} \). In this case, it may be better to apply \( \hat H - E \hat 1 \) to the state and then take the norm of the state.

There are two more ways to calculate expectation values: First, syten-expectation supports a template file input method. This is a file which contains expressions describing matrix-product operators. Those expressions are parsed and the expectation values are evaluated and placed in the same position. This means that with an input file like

$ cat input.template
Sz_0: { sz 0 @ }
Sz_1: { sz 1 @ }
Sz_2: { sz 2 @ } (Sz_1 + Sz_2): { sz 1 @ sz 2 @ + }
Hp: { Hp }
Hz: { Hz }
H: { Hz Hp + }

we can get output (numerical values are different from before as a different random state was used)

$ syten-expectation --template-file input.template -l sp.lat rnd.state 1>output.exp
2017.05.02T10:48:50.354 2.76850e-02 I Loaded first MP state.
2017.05.02T10:48:50.392 6.17180e-02 I Loaded MP lattice from 'sp.lat'.
2017.05.02T10:48:50.392 6.18540e-02 I Read template file, calculating and outputting result…
$ cat output.exp
Sz_0: (0.0999999999999998,8.64270111212828e-19)
Sz_1: (-0.1,-2.59899358718963e-18)
Sz_2: (6.2450045135165e-17,2.81838030913164e-18) (Sz_1 + Sz_2): (-0.1,2.19386721942014e-19)
Hp: (0,0)
Hz: (-0.799999999999999,2.13678452622974e-17)
H: (-0.799999999999999,5.66419973183401e-17)

Options -r, -i and -c work as normal, printing only the real or imaginary parts of both parts in column format:

$ syten-expectation --template-file input.template -l sp.lat rnd.state -r 1>output-real.exp
2017.05.02T10:50:31.357 2.74330e-02 I Loaded first MP state.
2017.05.02T10:50:31.393 6.15900e-02 I Loaded MP lattice from 'sp.lat'.
2017.05.02T10:50:31.393 6.17370e-02 I Read template file, calculating and outputting result…
$ cat output-real.exp
Sz_0: 0.0999999999999998
Sz_1: -0.1
Sz_2: 6.2450045135165e-17 (Sz_1 + Sz_2): -0.1
Hp: 0
Hz: -0.799999999999999
H: -0.799999999999999

Second, it is possible to calculate all expectation values of a local operator. This operation is very cheap, as it only requires sweeping once through the system. The associated tool is syten-localexpectation, it takes a state file, a lattice file and the name of the local operator, for example:

$ syten-localexpectation rnd.state sp.lat sz
2017.05.02T10:52:48.934 2.31760e-02 I Loading first state…
2017.05.02T10:52:48.938 2.76710e-02 I Loading lattice…
2017.05.02T10:52:48.985 6.93330e-02 N Overlap: (1,1.52152000758595e-17)
Site                           Real part            Imaginary part
0                 +9.999999999999952e-02    +1.598782521119432e-18
1                 -1.000000000000002e-01    +2.190961595535765e-18
2                 -1.804112415015878e-16    +8.572885622364144e-18
3                 +3.000000000000002e-01    -5.982029749571607e-18
4                 -9.999999999999995e-02    +8.116160650195126e-18
5                 -2.999999999999998e-01    -2.728928454611426e-18
6                 +6.245004513516502e-17    -9.745382243168825e-19
7                 -9.999999999999994e-02    -6.670744698797513e-18
8                 +1.000000000000000e-01    -8.880928559645446e-19
9                 -1.999999999999998e-01    -4.250448462197377e-18
10                +1.999999999999999e-01    +6.605115589029011e-18
11                +2.636779683484745e-16    +1.253370802476233e-18
12                +1.000000000000000e-01    -1.166925434123899e-18
13                -9.999999999999991e-02    +1.166925434123897e-18
14                +1.999999999999998e-01    +6.605115589029011e-18
15                -9.999999999999984e-02    -5.771968469783327e-18
16                -2.220446049250312e-16    +3.545945734620567e-19
17                +1.999999999999998e-01    +4.564560022757855e-18
18                +2.012279232133095e-16    +6.334271516214067e-19
19                -1.999999999999997e-01    -2.488725474085579e-18

Whenever possible, it is strongly advised to use syten-localexpectation over syten-expectation. By first applying an operator sp 10 @ to the state and then taking the expectation value between two different states, it is also possible to calculate a series of correlators (we use the ground state now to get non-zero long-range correlators):

$ syten-apply-op -i gs.state -l sp.lat:'sp 10 @' -o gs-sp10.state
2017.05.02T10:57:49.090 3.32850e-02 I Loaded MP state file 'gs.state'.
2017.05.02T10:57:49.090 3.33900e-02 I Loading and parsing MP operator 'sp.lat:sp 10 @'…
2017.05.02T10:57:49.118 6.14100e-02 I Applying operator naively…
2017.05.02T10:57:49.156 8.81680e-02 I Truncating resulting state to δs = 2.22045e-15 and m = 4294967295…
2017.05.02T10:57:49.225 1.45227e-01 N Incurred error: 0.
2017.05.02T10:57:49.225 1.45326e-01 I Saving to file 'gs-sp10.state'…
$ syten-localexpectation gs.state sp.lat sm gs-sp10.state -s
2017.05.02T10:57:54.101 2.15540e-02 I Loading first state…
2017.05.02T10:57:54.114 3.39400e-02 I Loading second state…
2017.05.02T10:57:54.124 4.41800e-02 I Loading lattice…
2017.05.02T10:57:54.193 1.02784e-01 N Overlap: assumed 1
Site                           Real part            Imaginary part
0                 +2.888546600063733e-02    +2.733378069566696e-17
1                 -1.936268434684807e-02    -1.700576763718013e-17
2                 +3.747563552410487e-02    +3.364953229091393e-17
3                 -2.989003183295122e-02    -2.243586596805221e-17
4                 +4.896956350885955e-02    +3.952205942574597e-17
5                 -4.448342185441484e-02    -2.550016867178118e-17
6                 +6.946753313686281e-02    +4.374065100325854e-17
7                 -7.630144227132522e-02    -2.655648876317488e-17
8                 +1.213845373388270e-01    +4.328530882650883e-17
9                 -2.425650442492071e-01    -3.139613086468641e-17
10                +5.000000001637165e-01    +1.246499055424441e-17
11                -3.492507360368306e-01    +7.871318565507556e-17
12                +1.206648079029102e-01    -6.295240819415821e-17
13                -1.255965620799623e-01    +1.090707530277025e-16
14                +6.742866565790358e-02    -4.002398715587754e-17
15                -7.909107936596899e-02    +1.972552262826181e-17
16                +4.478713451727234e-02    -2.799754276323828e-17
17                -5.782166151839252e-02    +3.933184768449661e-17
18                +2.899128653486479e-02    -2.145116317043427e-17
19                -4.369196656634142e-02    +3.375213988734672e-17
$ syten-expectation gs.state sp.lat:"sp 10 @ sm 11 @ *" -qr
-0.349250736036831
$ syten-expectation gs.state sp.lat:"sp 10 @ sm 12 @ *" -qr
0.120664807902911
$ syten-expectation gs.state sp.lat:"sp 10 @ sm 13 @ *" -qr
-0.125596562079963
$ syten-expectation gs.state sp.lat:"sp 10 @ sm 0 @ *" -qr
0.0288854660006371

As can be seen, the values calculated via syten-expectation and those obtained from syten-localexpectation agree perfectly. Note that it is necessary to pass -s/--no-normalise to syten-localexpectation to avoid the bare expectation value \( \langle 0 | s^-_i | s^+_{10} \rangle \) to be normalised by \( \langle 0 | s^+_{10} \rangle \), which would of course result in \( \pm \infty \).

Calculating a Norm

$ syten-norm rnd.state

this should print 1. The operator application above will have changed the norm, however:

$ syten-norm H-rnd.state
Loaded MP state, calculating norm…
2.48746859276656

Calculating an Overlap

Given two states, we can calculate the overlap:

$ syten-overlap rnd.state H-rnd.state
Loaded first MP State…
Loaded second MP state…
Calculating overlap…
(-1.3,-3.45763353487656e-16)

Obviously this is just the expectation value of H from before.

Running DMRG

Running DMRG is a bit tricky. Primarily, it requires a MPO acting as the Hamiltonian, a MPS to start with and a staging description. The staging description consists of one or more stages: (<stage 1>) (<stage 2>) (<stage 3>) each contained in round parantheses. Each stage can fix the number of sweeps, the truncation threshold, the maximal number of states etc. Parameters not set will be copied over from a previous stage. For example:

$ syten-dmrg -l sp.lat:"Hz Hp +" -i rnd.state -s "(m 20 x 5) (m 50) (m 100 x 10 t 1e-12) ()" -o dmrg

will first sweep up to five times (x 5) with up to twenty states (m 20), then with up to fifty states (m 50) again at most five times and then with 100 states up to ten times with truncation error now reduced to \( 10^{-12} \) (m 100 x 10 t 1e-12). Finally, it will repeat the last stage if necessary (the empty description at the end). This will result in a lot of output, primarily from each iteration:

2017.04.25T16:43:01.998 8.79050e-02    1  D Stg. Swp.  #Pos            eEnergy eIt     eRes             Energy  TruncErr MixFactor    #N    #S    #R      #T    #B    eReason
2017.04.25T16:43:01.998 8.79490e-02 N  M Starting stage 1…
2017.04.25T16:43:01.998 8.79800e-02 N  M Unconverged: Minimum number of sweeps not reached.
2017.04.25T16:43:01.998 8.80080e-02 N  1 Starting  right-to-left (←) sweep 1 of stage 1 at position 19 using DMRG3S.
2017.04.25T16:43:02.000 8.92800e-02    1  ←    1    1    19  -1.7500000000e+00   2  1.7e+00  -1.7500000000e+00  0.00e+00  1.00e+00     2     2     2       2     1    invSubsp[n=0.0e+00]
2017.04.25T16:43:02.001 9.10180e-02    1  ←    1    1    18  -2.5590169944e+00   3  2.0e+00  -2.5590169944e+00  0.00e+00  1.00e+00     4     3     4       4     2    invSubsp[n=2.3e-13]
2017.04.25T16:43:02.003 9.26480e-02    1  ←    1    1    17  -2.6896926208e+00   3  2.0e-03  -2.6896926208e+00  0.00e+00  1.00e+00     4     3     4       4     2    invSubsp[n=1.1e-13]

The first column is the current time, the second the CPU time used up to now. The third column (1) is the ID of the worker (always one here as we don't do PDMRG at the moment). The next columns are

  • the current direction, either left () or right ()
  • the current stage
  • the current sweep within that stage
  • the current site which is being updated
  • the energy resulting from the eigensolver solution
  • the number of iterations used by the eigensolver
  • one minus the overlap between the last two eigensolver iterations
  • the energy after the subspace expansion
  • the truncation error encountered
  • the subspace expansion mixing factor
  • the number of blocks in the tensor
  • the number of quantum number sectors on the MPS bond
  • the size of all reduced blocks added together (colloquially the number of states, m)
  • the total size incorporating possible non-abelian symmetries (what an equivalent abelian calculation would need)
  • the cubically-average block size
  • the exit reason from the eigensolver

DMRG will result in a series of *.state files:

$ ls -l dmrg*state
-rw-rw-rw- 1 C.Hubig ls-schollwoeck  57101 2017-04-25 16:49 dmrg_1_5.state
-rw-rw-rw- 1 C.Hubig ls-schollwoeck 247530 2017-04-25 16:49 dmrg_2_4.state
-rw-rw-rw- 1 C.Hubig ls-schollwoeck 744820 2017-04-25 16:49 dmrg_3_4.state
-rw-rw-rw- 1 C.Hubig ls-schollwoeck 745628 2017-04-25 16:49 dmrg_4_2.state

growing larger as their bond dimensions increased from stages 1 through 4. The energy and variance of these states now should result in better values:

$ for i in $(ls -1tr *.state); do printf "${i}:\n"; syten-expectation -qrv ${i} sp.lat:"Hz Hp +"; done
rnd.state:
⟨O⟩:          -1.3
⟨O²⟩:         6.1875
|⟨O²⟩-⟨O⟩²|:  4.4975
H-rnd.state:
⟨O⟩:          -3.3070707070707
⟨O²⟩:         17.0516414141414
|⟨O²⟩-⟨O⟩²|:  6.11492475257628
dmrg_1_5.state:
⟨O⟩:          -8.68247180616517
⟨O²⟩:         75.385320972876
|⟨O²⟩-⟨O⟩²|:  4.30802299433708e-06
dmrg_2_4.state:
⟨O⟩:          -8.68247333353095
⟨O²⟩:         75.3853431910421
|⟨O²⟩-⟨O⟩²|:  3.56602924968302e-09
dmrg_3_4.state:
⟨O⟩:          -8.68247333439879
⟨O²⟩:         75.3853432025469
|⟨O²⟩-⟨O⟩²|:  8.24277846980541e-13
dmrg_4_2.state:
⟨O⟩:          -8.68247333439878
⟨O²⟩:         75.3853432025467
|⟨O²⟩-⟨O⟩²|:  8.10019007028338e-13

Parallelisation in DMRG

SyTen supports four different levels of parallelisation, with three currently implemented. Of those, the --threads-tensor and --threads-super options are the most relevant for DMRG: The former allows you to parallelise over the tensor blocks in each tensor. This is helpful in all cases where you have some sort of symmetry in your system. The latter automatically spawns multiple worker threads (cf. Stoudenmire and White, 2013). For this real-space DMRG parallelisation to work correctly, it is often necessary to first do some preliminary sweeps with serial DMRG to obtain an overall relatively good ground-state approximation, e.g. with m = 200 states and five sweeps before increasing m to 5000...10000 and sweeping many more times.

Local Basis Optimisation

In bosonic systems, the local degrees of freedom are often very large. If the boson number is not conserved, the effective local degrees of freedom may often be described as a superposition of different boson numbers. If the boson number is conserved, it may still be possible to only consider a subspace (e.g. with \( n \in [ 490, 510 ] \) if there are \( n \approx 500 \) bosons per site on average). Such an optimisation is called a local basis optimisation and is probably best described in Florian Dorfner's PhD thesis. The idea is to have an extra matrix \( V_i \) on every site which translates the full basis (e.g. \( n \in [0, 1000 ]\)) into an effective local basis. That effective local basis is then treated as normal for MPS.

The tool syten-mps2lbo can translate a standard MPS into the LBO-MPS structure. When running DMRG with a LBO-MPS state as the input, DMRG will automatically switch to the LBO-MPS-DMRG mode and on every site first update the rank-3 MPS tensor, then the rank-2 matrix \( V_i \) and then the rank-3 MPS tensor again. Furthermore, it uses the secondary truncation settings st/sectruncTol, sb/sectruncMaxBlocksize, sm/sectruncMaxStates and sw/sectruncWeight to truncate the effective local basis. On some systems, this can lead to a large speed-up. Compare for example:

$ syten-mps-ssh-fermi-bose -l 4 -n 50 -o lat
$ syten-random -l lat -s 4,0 -o rnd -m 2 -g c
$ syten-dmrg -l lat:H_ges -i rnd -s "(m 100)" -o std

where the DMRG calculation itself takes about 5.4 seconds with final variance \( 4.7·10^{-11} \) to

$ syten-mps-ssh-fermi-bose -l 4 -n 50 -o lat
$ syten-random -l lat -s 4,0 -o rnd -m 2 -g c
$ syten-mps2lbo rnd
$ syten-dmrg -l lat:H_ges -i rnd.lbs -s "(m 100 sm 6)" -o lbo

where the DMRG calculation only took 3.1 seconds with variance \( 1.2·10^{-10} \).

As one can see, while the precision is slightly lower (but could in principle be increased using a larger bond dimension), the computational effort is also lower.

Note that at the moment, MPS-LBO-DMRG does not support the real-space parallelisation.

Reordering of Sites

Sometimes, the optimal position of each site is not known beforehand. For this reason, various lattice binaries have support for on-the-fly reordering of sites during lattice generation. The idea is to calculate an initial ground state using some arbitrary order of sites, then calculate the optimal order of those sites based on minimising the entanglement between sites far away from each other and then re-generate the lattice and re-run DMRG with this reordering.

This optimisation is mostly applicable to specific lattice files implementing complicated models (e.g. DMFT, quantum chemistry etc.), but we can model the process very easily with the standard Heisenberg Hamiltonian on ten sites if we add an extra-strong long-range interaction between the zeroth and eigth site. For this, first generate the standard lattice, random state and run DMRG with the additional link \( 10 \left(\hat s_0 \cdot \hat s_8\right) \).

$ syten-sql-mps-spin -l 10 -o sp
$ syten-random -l sp -o rnd
$ syten-dmrg -l sp:'H s 0 @ s 8 @ _dot 10 * +' -i rnd -o dmrg -s "(m 100 t 1e-12)" -f gs

Now, calculate the mutual information and, based on this, the optimal order of sites:

$ syten-mutualInformation -o gs -f reordering
$ cat reordering
[10]: {0  8  1  2  9  3  4  5  7  6}

This says that sites 0 and 8 should be close together, which is not really a surprise. We can also calculate some expectation values:

$ syten-expectation gs sp:"s 0 @ s 1 @ _dot" -qr
-0.0551028783463896
$ syten-expectation gs sp:"s 0 @ s 2 @ _dot" -qr
0.0343474606218299
$ syten-expectation gs sp:"s 0 @ s 8 @ _dot" -qr
-0.743944743140135
$ syten-expectation gs sp:"s 1 @ s 8 @ _dot" -qr
0.0505435789135365

As you can see, there is nearly a full singlet link between sites 0 and 8! Now, to generate the re-ordered lattice, pass the 'reordering' file to syten-sql-mps-spin:

$ syten-sql-mps-spin -l 10 -o sp-2 -f reordering
[…]
 Ordering is defined by new[old]: [0, 2, 3, 5, 6, 7, 9, 8, 1, 4].
[…]

So the old site 0 is still at position zero, while the old site 1 is at position two and the old site 8 is at position one. We can re-use the same initial state as before, since the lattice is homogeneous (if the lattice is not homogeneous, you have to regenerate the state) and run DMRG. Note that now we have to link sites 0 and 1 extra-strongly:

$ syten-dmrg -l sp-2:'H s 0 @ s 1 @ _dot 10 * +' -i rnd -o dmrg -s "(m 100 t 1e-12)" -f gs-2

On such a small system, we can represent energy expectation values essentially perfectly regardless of ordering of sites. However, we do expect that this reordering will have changed some things. For one, we have to use the mapping from before to calculate expectation values:

$ syten-expectation gs-2 sp-2:"s 0 @ s 2 @ _dot" -qr
-0.0551028765294503
$ syten-expectation gs-2 sp-2:"s 0 @ s 3 @ _dot" -qr
0.0343473933134452
$ syten-expectation gs-2 sp-2:"s 0 @ s 1 @ _dot" -qr
-0.743944740594402
$ syten-expectation gs-2 sp-2:"s 2 @ s 1 @ _dot" -qr
0.0505435698153109
$ syten-expectation gs-2 sp-2:"s 0 @ s 3 @ _dot" -qr
0.0343473933134452

Furthermore, when comparing Shannon entropies, we see a marked drop between the original system and the re-ordered system from 1.66 between sites 4 and 5 to 0.85 between sites 4 and 5.

$ syten-info -e 1 gs
2017.05.10T11:28:32.148 2.64140e-02 I Loaded MP state
2017.05.10T11:28:32.148 2.65180e-02 I MP State on lattice of 10 sites.
2017.05.10T11:28:32.148 2.65840e-02 I Vacuum basis (right edge): [<1, [SU(2)_S:{+0}@1]>]
2017.05.10T11:28:32.148 2.66160e-02 I Sector basis (left edge):  [<1, [SU(2)_S:{+0}@1]>]
2017.05.10T11:28:32.148 2.66370e-02 I Physical basis on site 1:  [<1, [SU(2)_S:{+0.5}@2]>]
2017.05.10T11:28:32.148 2.66640e-02 I Printing the Rényi entropy at every bond with α = 1
#BondL Entropy=-\sum_i s_i^2 log[\sum_i s_i^(2*1)]
0     6.93147180559945175204e-01
1     1.37854867773739608872e+00
2     1.32588531283508093139e+00
3     1.51066630975269133863e+00
4     1.66872211623284938398e+00
5     1.60379251529139876808e+00
6     1.92585050154392178889e+00
7     1.38005482573962723336e+00
8     6.93147180559945175204e-01
$ syten-info -e 1 gs-2
2017.05.10T11:28:34.002 2.57460e-02 I Loaded MP state
2017.05.10T11:28:34.002 2.58830e-02 I MP State on lattice of 10 sites.
2017.05.10T11:28:34.002 2.59530e-02 I Vacuum basis (right edge): [<1, [SU(2)_S:{+0}@1]>]
2017.05.10T11:28:34.002 2.60060e-02 I Sector basis (left edge):  [<1, [SU(2)_S:{+0}@1]>]
2017.05.10T11:28:34.002 2.60280e-02 I Physical basis on site 1:  [<1, [SU(2)_S:{+0.5}@2]>]
2017.05.10T11:28:34.002 2.60580e-02 I Printing the Rényi entropy at every bond with α = 1
#BondL Entropy=-\sum_i s_i^2 log[\sum_i s_i^(2*1)]
0     6.93147180559945175204e-01
1     4.36124402513411585369e-02
2     7.11730080105066442719e-01
3     6.54959408882573423938e-01
4     8.56123983551658862545e-01
5     7.84680267225669636488e-01
6     7.61451193413838156587e-01
7     4.76636073573340646181e-01
8     6.93147180559945175204e-01

Similarly, the decay of singular values at this bond is now much faster: Plotting

$ syten-info -s gs-2 | grep '^ *4 ' | sort -rg -k 3 > 2-4

against

$ syten-info -s gs | grep '^ *4 ' | sort -rg -k 3 > 1-4

using gnuplot:

$ plot  "1-4" u 0:3:2 w lp lt 7 lc 1 ps variable t "Singular values on bond (4,5) without reordering", "2-4" u 0:3:2 w lp lt 6 lc 2 ps variable t "Singular values on bond (4,5) with reordering"

where we use the symmetry-related multiplicity of each singular value in the second column to scale the point size. The resulting graph:

TUTORIAL-reordering-svd.svg
Singular value distribution on the fourth bond of the spin chain with and without reordering.

Time Evolution

Additional reading for this section includes the FAQ on the available time evolution algorithms. Here, we will only discuss the syten-krylov tool.

To nicely visualise the results of a time evolution, we will stick to a \( \mathrm{U}(1)-S_z \)-symmetric spin chain of twenty sites:

$ syten-sql-mps-spin -l 20 --sym u1 -o sp

our goal is then to evaluate time-dependent correlators \( \langle 0 | \hat s^z_i(t) \hat s^z_j | 0 \rangle \). For this, we first calculate the ground-state as usual:

$ syten-random -l sp -o rn
$ syten-dmrg -l sp:'Hp Hz +' -i rnd -f gs -s "(m 200)"

we then apply the \( \hat s^z_10 \) operator in the middle of the chain:

$ syten-apply-op -l sp:'sz 10 @' -i gs -o exc -t 1e-6

while discarding singular values smaller than \( 10^{-6} \). We now wish to do a real-time evolution with steps \( \delta t = 0.1 \) up to time \( T = 10 \). Since this is a pure MPS calculation, we can use one of the more advanced eigensolver modes, which are not available for BTT- or LBO-MPS calculations. Specifically, we call:

$ syten-krylov -l sp:'Hp Hz +' -i exc -d 0.1 -T 1 -t 1e-5 -e '{ mode V }'

where -i is the initial state, -l specifies the Hamiltonian with which we want to evolve the system, -d the time step, -T the maximal time, -t the truncation threshold of singular values and finally -e allows us to configure the eigensolver a bit further. See the output of syten-krylov --help for more details.

syten-krylov will then generate a Krylov subspace based on the initial vector exc and attempt to calculate the exponential \( e^{-i t \hat H} | \mathrm{exc} \rangle \) in this subspace. If this succeeds for \( t = \delta t \), it will attempt to use the same space for the next step. If this fails, it will either build an entirely new subspace or attempt to add another Krylov vector. In this particular case, three Krylov vectors will typically be sufficient to evolve to \( T = 10 \), largely because the initial state primarily consists of three eigenstates of the Hamiltonian. In larger systems, it will be necessary to re-build the subspace more frequently and take longer to time-evolve.

Now, considering the correlator above, if we fix \( i = 2 \) and \( j = 10 \), we are already halfway there. Specifically, of the expression \( \langle 0 | e^{i \hat H t} \hat s^z_2 e^{-i \hat H t} \hat s^z_{10} | 0 \rangle \), we have evaluated the right part \( e^{-i \hat H t} \hat s^z_{10} | 0 \rangle \). The left \( e^{i \hat H t} \) will simply give a phase which can straightforwardly be calculated once we know the energy:

$ syten-expectation gs sp:'Hz Hp +' --precision 17 -qr

what remains is hence the expectation value, but this is also easy to get:

$ for i in $(seq -f %g 0.1 0.1 10); do printf "${i} "; syten-expectation -qc krylov_T-${i}.state sp:'sz 2 @' gs -s; done > exp

Finally, we can use Python to read in the lines from the file exp, calculate the phase and multiply everything together:

$ (python <<EOF
from cmath import *
f = open('exp')
line = f.readline()
while line:
t, re, im = line.split()
res = (float(re) + float(im)*1j) * exp(-8.6824733343987752 * 1j * float(t))
print t, res.real, res.imag
line = f.readline()
EOF
) > res

Depending on the system, to get useful values it may be necessary to subtract the constant part \( \langle 0 | \hat s^z_i | 0 \rangle \langle 0 | \hat s^z_j | 0 \rangle \), especially if this is non-zero. From here, in principle it takes just two Fourier transformations to get to the dynamical structure factor.

A similar calculation can yield the Green's function useful for DMFT calculations. Additionally, by replacing 0.1 above by (0,-0.1), it becomes possible to do imaginary time evolution, which is often much cheaper than real-time evolution (it may then also be possible to restrict oneself to real-valued numbers only, speeding up calculations by another factor of two).

A Word on Scripting Languages

The compiled executables are intended to be the primary interface of SyTen and as such, nearly all functionality can be used from the command-line alone. It hence is possible with just a few Bash scripts to drive most ordinary calculations. Nevertheless, it may be helpful to consider using another language, such as Python, as the primary ‘scripting glue’ (e.g. because in Python, one actually can use complex numbers, as seen above). Nils Linden has used Python extremely successfully to drive difficult DMFT calculations with various levels of dynamical parallelisation over parameter regimes far beyond the scope of SyTen.

In any case, knowledge of some scripting language is extremely helpful and essentially a prerequisite, not just to call the SyTen binaries but also to manipulate the text output produced.

Additionally, it should be pointed out that most of the functionality is actually contained in the SyTen ‘library’, i.e. if desired, it is entirely possible to build C++ drivers instead of Bash scripts which directly call the relevant functions. This can even be faster, since states and lattices have to be serialised and deserialised in order for different binaries to work on them, while such serialisation and deserialisation is not necessary if all calculations take place inside a single C++ binary. However, in most applications, this overhead is very low.

In any case, the user will have the most pleasant time and the fastest results when using a language they already know, even if this language is JavaScript! As such, the only suggestion here can be to use whatever scripting glue layer with which one is most comfortable.

Flexible Fermionic Interactions with syten-mps-fermi-hubbard-qc

There are quite a few different lattice files available in lat/ and many more can be created using the existing models defined in lat/inc. However, most of them only define relatively rigid interaction and hopping terms without a lot of flexibility. There is however one lattice, syten-mps-fermi-hubbard-qc, which allows the user to supply their own interaction tensors. Specifically, the lattice reads in two text files into tensors \( t_{ij} \) and \( V_{ijkl} \) and creates two Hamiltonian defined on a Fermi-Hubbard chain of length \( L \) as

\[ \hat H_1 = \sum_{i=1}^L \sum_{j=1}^L t_{ij} \left( c^\dagger_{\uparrow, i} c_{\uparrow, j} + c^\dagger_{\downarrow, i} c_{\downarrow, j} \right) \]

and

\[ \hat H_2 = \sum_{\sigma = \uparrow, \downarrow} \sum_{\tau = \uparrow, \downarrow} \sum_{i=1}^L \sum_{j=1}^L \sum_{k=1}^L \sum_{l=1}^L \left( V_{ijkl} c^\dagger_{\sigma, i} c^\dagger_{\tau, j} c_{\tau, l} c_{\sigma, k} \right) \]

As can be seen, arbitrary one- and two-particle terms can be encoded in these tensors \( t_{ij} \) and \( V_{ijkl} \). As a first example, we will encode a nearest-neighbour interaction on a lattice of four sites and compare with the existing implementation in syten-sql-mps-fermi-hubbard. First, create a file tij which contains the textual description of the tensor \( t_{ij} \). At the moment, only the dense format is available here but there are plans to also make a ‘tabulated’ format available. tij should have contents:

$ cat tij
Start: {4, 4}/16
(0,0)  (-1,0) (0,0)  (0,0)
(-1,0) (0,0)  (-1,0) (0,0)
(0,0)  (-1,0) (0,0)  (-1,0)
(0,0)  (0,0)  (-1,0) (0,0)
End: {4, 4}/16

Note that you can replace any linebreaks by simple white space. The {4, 4}/16 are the dimensions of the tensor and its total length. This format is generated by write(), so you can generate the above output with for example the following tiny program:

#include <iostream>
int main(int, char**) {
using namespace syten;
DenseTensor<2> tij({4,4});
for(Index i(0); i != 3; ++i) {
tij[{i, (i+1)%4}] = -1;
tij[{(i+1)%4, i}] = -1;
}
return 0;
}

Afterwards, call syten-mps-fermi-hubbard-qc as

$ syten-mps-fermi-hubbard-qc -1 tij -o fh

and generate the comparison file via

$ syten-sql-mps-fermi-hubbard -l 4 -w 1 -o fh-comp --sym u1

We can already check that the operator H1 in the file fh is equal to the operator Ht in the file fh-comp:

$ syten-overlap -qn fh:H1 fh-comp:Ht
(1,0)

where -n normalises \( \mathrm{tr}\left( H_1 H_t \right) \) by \( \sqrt{ \mathrm{tr}\left( H_1 H_1 \right) \mathrm{tr}\left( H_t H_t \right)} \). These operators also have the same ground states:

$ syten-random -s 3,0.5 -l fh -o rnd
$ syten-dmrg -i rnd -l fh:H1 -f dmrg-fh -s "(m 16)"
$ syten-dmrg -i rnd -l fh-comp:Ht -f dmrg-fh-comp -s "(m 16)"
$ syten-overlap -qsn dmrg-fh dmrg-fh-comp
(0,0)

\( V_{ijkl} \) would proceed in the same fashion, though of course at much higher dimension. Both \( t_{ij} \) and \( V_{ijkl} \) are stored in row-major fashion, i.e. the second stored element is \( t_{12} \) and \( V_{1112} \) respectively. For example, to get a Hubbard- \( U \) of \( 6 \) on all sites, use the tensor

$ cat vijkl
Start: {4, 4, 4, 4}/256
(3,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (3,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (3,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (3,0)
End: {4, 4, 4, 4}/256

and generate the lattice file as

$ syten-mps-fermi-hubbard-qc -1 tij -o fh -2 vijkl

The difference between the (3,0) in the file and the effective interaction strength six comes from the double-sum over \( \sigma \) and \( \tau \) in the definition of \( \hat H_2 \). Afterwards, the operator should have equal norm with Hd 6 * from the comparison lattice and overlap close to 1:

$ syten-norm -q fh-comp:'Hd 6 *'
126.9960629311
$ syten-norm -q fh:H2
126.9960629311
$ syten-overlap -qn fh-comp:'Hd 6 *' fh:H2
(0.9999999999999999,0)

Writing a Lattice File

For many problems, appropriate lattice files already exist: square-lattice spin systems with nearest-neighbour interactions are implemented in syten-sql-mps-spin, for the Fermi-Hubbard system in syten-sql-mps-fermi-hubbard[-kns and syten-mps-fermi-hubbard-*. Spinless bosons exist in two variants with and without symmetries in syten-mps-boson-spinless-*. Additional, more specialised systems for DMFT and various embedding models also exist and may serve as an inspiration.

However, when studying a new model, it is likely necessary to write such a lattice file. In the following, we will consider the case where the underlying model (spinless bosons, fermions, spins etc.) is already implemented and we only have to take care of new interactions. In this case, it is sufficient to create a new file in lat/ with the .cpp ending. This file will be a programme which, when executed, creates the desired lattice file.

As an example, consider the U(1) model for spinless bosons, which you can also find in lat/syten-mps-boson-spinless-u1.cpp. We will go through each part of the file, explaining its purpose.

#include "lat/inc/u1.h"

The u1.h header file defines the basic bosonic model, i.e. sites with certain single-site operators. bpo_helper.h contains utilities for creating an executable and handles mostly argument parsing as well as setup of the threading environment.

namespace syten {

By placing everything in our own namespace, we can make sure that global functions, such as abs() don't take precedence over our own functions.

int main(int argc, char** argv) {

Within the namespace, we then define the main() function with the usual signature.

Index length(10);
#if SYTEN_SCALAR_COMPLEX
#else
double j{-1.};
#endif
double u(1);
Index n = 4;
std::string output;

Some variables are defined to store the options passed into the program. If complex scalars are enabled (the default), SYTEN_SCALAR_COMPLEX evaluates to true. In this case, we define a potentially complex-valued coupling j. Otherwise, the coupling is purely real, as the toolkit then does not know how to deal with complex numbers.

opt.add_options()
("length,l",boost::program_options::value(&length)->default_value(length), "length of the chain")
("nearest-neighbour,j",boost::program_options::value(&j)->default_value(j),"nearest neighbour hopping coefficient")
("on-site,u",boost::program_options::value(&u)->default_value(u),"on-site interaction coefficient")
("out,o",boost::program_options::value(&output)->default_value(""))
("max-site,n",boost::program_options::value(&n)->default_value(n),"max. number of bosons per site")
;

Here, we define the options understood by our executable. Each option requires a short and/or long name, separated by a ,. The second argument is the place where the passed-in value will be stored. Here, we store it in the variables length, j, u and output, passing the addresses of the variables (by using the & operator). With ->default_value(x) we can also set a default value for this variable which will be shown in --help. The last argument inside each () is the description of the option. The macros SYTEN_BPO_INIT and SYTEN_BPO_EXEC initialise the argument parsing variables and also inject some variables into the surrounding block.

if (bpo_help) {
std::cout << "U(1)-symmetric spinless bosonic model\n"
<< "Defaults to creating a L=10 chain in the file bo-u1-10.lat.\n"
<< "Usage: " << argv[0] << " [options]\n"
<< opt;
return 0;
}

One of these injected variables is bpo_help, it is set to true if the user passed --help or-h`. In this case, we print a short help message and exit.

if (output.empty()) {
output = "bo-u1-"+std::to_string(length)+SYTEN_EXT_MP_LATTICE;
}
if (!bpo_quiet) {
std::cerr << "Creating a lattice of length " << length
<< " to be saved into file " << output << ".\n";
}

If the user did not set an output file name, we construct one. Unless they specified -q/--quiet (in which case bpo_quiet will be true), we also tell them where we will put the lattice.

MPS::Lattice lat = MPS::Lat::U1::genBosonLattice(length, "U(1)-symmetric spinless bosonic model", n, 'n');

This is arguably the most important line. genBosonLattice() creates a lattice object of the specified length (length) with the specified number of bosons per site (n) and some description. This is already a perfectly useable lattice, it is just not particularly useful – only single-site operators are defined (and the identity).

MPS::Operator hj1(length);
MPS::Operator hj2(length);
MPS::Operator hu(length);
MPS::Operator particle_number(length);
for(Index i(0); i != length; ++i) {
if ((i+1)%50 == 0) {
hj1.truncate();
hj2.truncate();
hu.truncate();
particle_number.truncate();
}
progress(i, bpo_quiet);
if (i != length-1) {
hj1 += lat.get("c", i) * lat.get("ch", i+1);
hj2 += lat.get("c", i+1) * lat.get("ch", i);
}
hu += lat.get("n",i) * lat.get("n",i);
particle_number += lat.get("n", i);
}

The main part of the work is done here. We create four matrix-product operators hj1, hj2, hu and particle_number, each of length length and initialised to all-zero. We then loop over all sites of the lattice. If the current site is divisible by fifty, we compress the operators generated so far (to avoid linear growth of their bond dimensions). progress(i, bpo_quiet) gives nice feedback on the current status. Unless we are on the last site, we then add terms \( c^\dagger_{i+1} c_i \) and \( c^\dagger_{i} c_{i+1} \) to hj1 and hj2 respectively. The single-site operators are retrieved from the lattice object using the .get() method and multiplied together using operator*. Note the reversed order of operations – the left operand is applied first! The result of the multiplication is the added to the Hamiltonian terms hj1/hj2. Particle number and on-site interaction are generated in much the same way (but on every site of the system).

if (!bpo_quiet) { std::cerr << "\n"; }
if (!bpo_quiet) { std::cerr << "Truncating Hj1,2 and Ht…\n"; }
hu.truncate();
particle_number.truncate();
hj1.truncate();
hj2.truncate();

When adding two MPOs together, their bond dimension is the sum of the two addend bond dimensions. However, it is often possible to truncate this down to something much smaller at no loss of accuracy. The .truncate() method, called on MPOs, does this. Before saving the operators in the lattice, we hence call this method.

if (!bpo_quiet) { std::cerr << "Building and truncating full Hamiltonian and Hj…\n"; }
MPS::Operator hj = hj1 + hj2;
MPS::Operator h = j * hj1 + std::conj(j) * hj2 + u * hu;
hj.truncate();
h.truncate();

We also define two convenience-MPOs hj and h representing only hopping and hopping plus on-site interactions. u is the user-supplied interaction strength.

lat.add("Hj1","nearest neighbour hopping part 1: Σ_(i=1)^(L-1) { ch(i) c(i+1) }",hj1);
lat.add("Hj2","nearest neighbour hopping part 2: Σ_(i=1)^(L-1) { ch(i+1) c(i) }",hj2);
lat.add("Hj" ,"symmetric nearest neighbour hopping: Hj1 + Hj2",hj);
lat.add("Hu" ,"on-site squared interaction: Σ_(i=1)^L { n(i) n(i) }",hu);
lat.add("N" ,"total particle number:Σ_(i=1)^L { n(i) }",particle_number);
lat.add("H", "Hamiltonian: "+to_string_rd(j)+" Hj1 + "+to_string_rd(std::conj(j))+" Hj2 + "+to_string_rd(u)+" Hu",h);

Finally, the operators are added to the lat lattice object. The add() method takes the name of the operator (which you can use to get it in e.g. DMRG), a description shown in e.g. syten-info and finally the operator itself.

if (!bpo_quiet) { std::cerr << lat << "\n"; }

Printing the lattice object causes all the information stored in the lattice to be displayed.

save(lat, output);

This line saves the lattice object in the file output.

return 0;
}
}

We return from the syten::main() function and close the syten namespace.

int main(int argc, char** argv) { return syten::main(argc, argv); }

Finally, we define a proper main() function and tell it to directly pass control to the syten::main() function we defined above. Once the file is in place, we can compile it using

$ make lat/syten-my-new-lattice-file

and/or copy it into the external binary directory

$ make lat/syten-my-new-lattice-file--copy

From this, it should be fairly obvious to see how one could e.g. insert long-range interactions into the model: adapt the for loop to run over all pairs of (i,j) and insert the desired interaction term between sites i and j there. For an example of how to do this with variable coupling constants, see e.g. lat/syten-mps-fermi-hubbard-qc.cpp.

The primary problem is then the definition of an entirely new model, i.e. a new type of local site bases and single-site operators. These are defines in the files in lat/inc and are generally not too complicated – largely one has to define a basis with the appropriate quantum number sectors and the appropriate single-site operator components (i.e. the blocks of the tensors) by hand. For easy examples, see MPS::Lat::U1::genBosonLattice() in lat/inc/u1.cpp or MPS::Lat::SU2U1::genKondoModel() in lat/inc/su2u1.cpp which has two different types of sites: fermionic Hubbard sites and spin sites.