Tutorial

- File Types
- A Word on Scripting Languages
- Finite One-Dimensional Calculations with Matrix Product States
- Generating a Matrix Product Lattice
- Generating a Random Matrix Product State
- Generating a Product ‘Matrix Product’ State
- Applying an Operator to a State
- Calculating Expectation Values
- Calculating a Norm
- Calculating an Overlap
- Running DMRG
- Local Basis Optimisations for Matrix Product States
- Reordering of Sites
- Time Evolution
- Flexible Fermionic Interactions with syten-mps-fermi-hubbard-qc
- Writing a Lattice File

- Infinite Two-Dimensional Calculations with Infinite Projected Entangled Pair States

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.

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/`

.

The compiled executables and Python scripts 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). Nils Linden has used Python 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.

The first part of this tutorial concerns matrix-product states, operators and lattices as used for 1+1-dimensional calculations.

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 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!

- 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

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

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.

- Remarks
`syten-expectation`

by default computes \( \frac{\langle a | \hat O | b \rangle}{\langle a | b \rangle} \) (with potentially the same state used for both \( |a\rangle \) and \( |b\rangle \)). This is the correct definition of expectation values but may be surprising if you don't want the normalisation by \( \langle a | b \rangle \) (or if those states are orthogonal to each other). To only compute \( \langle a | \hat O | b \rangle \), pass the`-s`

option.

`-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 \).

$ 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

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 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

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.

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.

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:

Additional reading for this section includes the [FAQ](Which time-evolution algorithm should I use?) on the available time evolution algorithms. Here, we will only discuss the `syten-krylov`

tool, though the `syten-tdvp`

tool may in fact be better in practice.

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).

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 "inc/dense/dense.h"

#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)

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"

#include "inc/util/bpo_helper.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

std::complex<double> j(-1, 0);

#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()) {

}

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).

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.

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.

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.

Infinite Projected Entangled Pair States (iPEPS) are states constructed on an infinite two-dimensional lattice by repeating a finite number of tensors of a finite-size unit cell. SyTen offers two implementations, the initial syten::IPEPS code and the newer syten::IPEPSv2 code. Both are based on the same algorithms and ideas, merely implementation details (and what is implemented) differ. The former uses type-ranked Tensors (syten::Tensor) and was used for two publications, Hubig, SciPost Phys. 5, 047 and Hubig and Cirac, SciPost Phys. 6, 031. The latter uses the newer smart tensor class (syten::STensor) and is much more straightforward to extend, adapt etc. It can also be used for variational ground-state searches using automatic differentiation and provides a more unified time-evolution interface. However, it was less extremely optimised for performance and never tested with non-abelian symmetries (which means that one may have to insert a few calls to correct CGC spaces).

In the following, we will use the second-generation code, though the overall layout of the two versions is the same: There are states (syten::IPEPSv2::State) and lattices (syten::IPEPSv2::Lattice). The former represent infinite quantum states, the latter a specific physical system. Due to using the STensor class with its uniquely identified bases, each state has to be constructed for a specific lattice – recreating the lattice (even with identical parameters) and using it together with a previously-created state will not work.

Lattice objects define a set of one-site operators. These operators can be retrieved using the syten::IPEPSv2::Lattice::get() function which also requires passing in a coordinate on which this operator (relative to the unit cell) will act. The lower left site of the \( W \times D \) unit cell has coordinate \( (0,0) \), the lower right end has coordinate \( (W-1,0) \) and the upper right corner has coordinate \( (W-1,D-1) \). syten::IPEPSv2::Lattice::get() returns an applied one-site operator syten::IPEPSv2::AOneOp. This operator can be multiplied with another to create a two-site operator syten::IPEPSv2::ATwoOp. It can also be summed together to create a syten::IPEPSv2::SumOp, which is a list of operators on a unit cell representing their sum and (in theory) being applied to all unit cells of the lattice. Exponentiating a syten::IPEPSv2::SumOp yields a syten::IPEPSv2::ProdOp, an operator which contains individual factors acting on a unit cell of a larger operator.

Expectation values of a SumOp or local expectation values of a one-site operator can be evaluated and the exponential of a SumOp can be applied to a state to perform a real or imaginary time evolution with pre-existing Python scripts. Some further scripting (see below) allows performing variational ground-state searches.

The unique basis IDs of each tensor unfortunately make it hard to create a proper \( 1 \times 1 \) unit cell, as this would imply that the unit cell tensor has two identical legs (one incoming and one outgoing) to the left/right and top/bottom. As a workaround, it is straightforward to create a state with one effective tensor and three identities which only re-route entanglement (cf. lat/syten-ipeps2-1x1-spin.py).

An iPEPSv2 spin lattice is straightforward to create:

$ syten-ipeps2-spin.py -w 2 -d 2 -o sp

where `-w`

and `-d`

give the width and depth of the unit cell and `-o`

the output filename. Note that `syten-info`

cannot read the lattice file, but we can load it using pyten:

$ ipython3 <snip> > from pyten import * > lat = ipeps2.Lattice("sp") > print(lat)

will print a lot of information, most of it the terms of the Hamiltonian. To get a specific one-site operator, use `.get`

:

> lat.get("sz", [0,0])

which will return

IPEPSv2::AOneOp @ {0, 0} phys_inp: [2ab79276:11f89cbf:44791511#0] phys_out: [2ab79276:11f89cbf:44791511#1] auxiliary: [5d08f2b4:5d08f2b4:f1d2891f#0] Tensor: ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--== == STensor[3;0000000000]@7ffcb91b9b50: == == Leg 1: SBasis (inc) 'aux'@[5d08f2b4:5d08f2b4:f1d2891f#0]: <{< <[U(1)_z:{+0}@1+], 1> >}> == Leg 2: SBasis (out) 'p:0,0'@[2ab79276:11f89cbf:44791511#1]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Leg 3: SBasis (inc) 'p:0,0'@[2ab79276:11f89cbf:44791511#0]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Fermionic order: [1, 2, 3] == Composed of 2 blocks. -------------------------------------------------------------------------------------- == Block 1: Tensor block of rank 3 and 1 symmetries. isZero = 0. Reduced tensor is: (-5.00000000000000000e-01,0.00000000000000000e+00) Sym 1 (no CGC): Idx 1 transforms as: U(1)_z:{+0}@1+; Idx 2 transforms as: U(1)_z:{-0.5}@1+; Idx 3 transforms as: U(1)_z:{-0.5}@1+ -------------------------------------------------------------------------------------- == Block 2: Tensor block of rank 3 and 1 symmetries. isZero = 0. Reduced tensor is: (5.00000000000000000e-01,0.00000000000000000e+00) Sym 1 (no CGC): Idx 1 transforms as: U(1)_z:{+0}@1+; Idx 2 transforms as: U(1)_z:{+0.5}@1+; Idx 3 transforms as: U(1)_z:{+0.5}@1+ -------------------------------------------------------------------------------------- ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==

The AOneOp stores the site of the unit cell on which it acts ({0,0}) as well as the IDs of the physical input and output and auxiliary tensor bases (to refer to those bases programmatically). The auxiliary basis will carry a charge, if necessary – it is zero in the example above, but in e.g. `sp`

, it is not:

> lat.get("sp", [0,0]) IPEPSv2::AOneOp @ {0, 0} phys_inp: [2ab79276:11f89cbf:44791511#0] phys_out: [2ab79276:11f89cbf:44791511#1] auxiliary: [5d08f2b4:5d08f2b4: 1f9fbe3#0] Tensor: ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--== == STensor[3;0000000000]@7ffcb91b9b50: == == Leg 1: SBasis (inc) 'aux'@[5d08f2b4:5d08f2b4: 1f9fbe3#0]: <{< <[U(1)_z:{+1}@1+], 1> >}> == Leg 2: SBasis (out) 'p:0,0'@[2ab79276:11f89cbf:44791511#1]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Leg 3: SBasis (inc) 'p:0,0'@[2ab79276:11f89cbf:44791511#0]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Fermionic order: [1, 2, 3] == Composed of 1 blocks. -------------------------------------------------------------------------------------- == Block 1: Tensor block of rank 3 and 1 symmetries. isZero = 0. Reduced tensor is: (1.00000000000000000e+00,0.00000000000000000e+00) Sym 1 (no CGC): Idx 1 transforms as: U(1)_z:{+1}@1+; Idx 2 transforms as: U(1)_z:{+0.5}@1+; Idx 3 transforms as: U(1)_z:{-0.5}@1+ -------------------------------------------------------------------------------------- ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==

Multiplying two operators on different sites together yields a ATwoOp and can be done using `*`

:

> spsm = lat.get("sp", [0,0]) * lat.get("sm", [1,0]) > spsm IPEPSv2::ATwoOp @ {0, 0} ; {1, 0} phys_inp_l: [2ab79276:11f89cbf:44791511#0] phys_inp_r: [5d08f2b4:5d08f2b4:18b56c2b#2] phys_out_l: [2ab79276:11f89cbf:44791511#1] phys_out_r: [5d08f2b4:5d08f2b4:18b56c2b#3] Tensor: ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--== == STensor[4;0000000000]@000002a780f0: == == Leg 1: SBasis (out) 'p:1,0'@[5d08f2b4:5d08f2b4:18b56c2b#3]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Leg 2: SBasis (inc) 'p:1,0'@[5d08f2b4:5d08f2b4:18b56c2b#2]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Leg 3: SBasis (inc) 'p:0,0'@[2ab79276:11f89cbf:44791511#0]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Leg 4: SBasis (out) 'p:0,0'@[2ab79276:11f89cbf:44791511#1]: <{< <[U(1)_z:{-0.5}@1+], 1>, <[U(1)_z:{+0.5}@1+], 1> >}> == Fermionic order: [1, 2, 3, 4] == Composed of 1 blocks. -------------------------------------------------------------------------------------- == Block 1: Tensor block of rank 4 and 1 symmetries. isZero = 0. Reduced tensor is: (1.00000000000000000e+00,0.00000000000000000e+00) Sym 1 (no CGC): Idx 1 transforms as: U(1)_z:{-0.5}@1+; Idx 2 transforms as: U(1)_z:{+0.5}@1+; Idx 3 transforms as: U(1)_z:{-0.5}@1+; Idx 4 transforms as: U(1)_z:{+0.5}@1+ -------------------------------------------------------------------------------------- ==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==--==

The ATwoOp has no auxiliary basis legs. For more details on how to work with those operators and create a Hamiltonian, cf. lat/syten-ipeps2-spin.py.

A random iPEPSv2 state compatible with a given lattice can be created as follows:

$ syten-ipeps2-random.py -l sp -o rnd

Because already during the lattice creation it is necessary to fix the effective quantum number per unit cell (defaulted to zero above) to set the physical basis ID on the (0,0) site properly, it is not necessary to pass this information into `syten-ipeps2-random.py`

.

IPEPS constructed from random tensors are very unwieldy, as there is no particular guarantee that the corner transfer matrix (CTM) used for the evaluation of expectation values converges. As such, it is typically a good idea to project the random state down to a bond dimension-1 state with selected states on each site:

$ syten-ipeps2-project.py -i rnd -o rnd-projected.ips -s '1@-0.5;1@0.5;1@0.5;1@-0.5' -D 1

which creates a new file `rnd-projected.ips`

of which we can evaluate expectation values:

$ syten-ipeps2-expectation.py -l sp -O H -i rnd-projected.ips 2019.06.18T16:29:04.027 2.21477e-01 N Loading operator 'H' from file 'sp'… 2019.06.18T16:29:04.029 2.23986e-01 N Loading state from file 'rnd-projected.ips'… 2019.06.18T16:29:04.030 2.24470e-01 N Building double layer… 2019.06.18T16:29:04.052 2.27937e-01 N Initialising CTM… 2019.06.18T16:29:04.073 2.48891e-01 N Working at χ = 1 2019.06.18T16:29:04.094 2.69499e-01 I Estimate 0: (-2.0000000000000000E+00,-7.8504622934188734E-17) 2019.06.18T16:29:04.115 2.91180e-01 I Estimate 1: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+1.2326e-32 2019.06.18T16:29:04.136 3.11302e-01 I Estimate 2: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+1.2326e-32 2019.06.18T16:29:04.156 3.31586e-01 I Estimate 3: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+0.0000e+00 2019.06.18T16:29:04.176 3.51840e-01 I Estimate 4: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+0.0000e+00 2019.06.18T16:29:04.176 3.51909e-01 N Converged: Δ(-1) = +0.0000e+00 < +1.0000e-08 and Δ(-2) = +0.0000e+00 < +1.0000e-08 2019.06.18T16:29:04.182 3.57968e-01 N Working at χ = 2 2019.06.18T16:29:04.194 3.69651e-01 I Estimate 0: (-2.0000000000000000E+00,-7.8504622934188734E-17) 2019.06.18T16:29:04.214 3.90020e-01 I Estimate 1: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+1.2326e-32 2019.06.18T16:29:04.234 4.10096e-01 I Estimate 2: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+1.2326e-32 2019.06.18T16:29:04.255 4.30526e-01 I Estimate 3: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+0.0000e+00 2019.06.18T16:29:04.275 4.51154e-01 I Estimate 4: (-2.0000000000000000E+00,-7.8504622934188721E-17) Δ(-1)=+0.0000e+00; Δ(-2)=+0.0000e+00 2019.06.18T16:29:04.275 4.51223e-01 N Converged: Δ(-1) = +0.0000e+00 < +1.0000e-08 and Δ(-2) = +0.0000e+00 < +1.0000e-08 2019.06.18T16:29:04.275 4.51250e-01 N Converged: δ(-1) = +0.0000e+00 < +1.0000e-08 2019.06.18T16:29:04.276 4.51275e-01 I 1 (-2.0000000000000000E+00,-7.8504622934188721E-17) +0.0000e+00 +0.0000e+00 2019.06.18T16:29:04.276 4.51293e-01 I 2 (-2.0000000000000000E+00,-7.8504622934188721E-17) +0.0000e+00 +0.0000e+00 -2.0000000000000000E+00 -7.8504622934188721E-17 0.0000E+00 2

The result of an expectation value calculation in iPEPS with CTM depends on the number of CTM steps done and the bond dimension of the CTM. To check convergence in both, `syten-ipeps2-expectation.py`

increases the bond dimension repeatedly and does a number of steps at each of them until convergence at this bond dimension (or the maximal number of steps) has been reached. In the present case, the state is unentangled and as such, the expectation value converges quickly. \( -2 \) is the correct value of the Heisenberg Hamiltonian on a \( 2 \times 2 \) unit cell with the Néel state. The last line gives the expectation value in real and imaginary part, the estimated error (from changes in the value during iterations) and the bond dimension used for this value. Note that for IPEPSv2, at the moment it is not possible to calculate expecation values of operators which are not defined as a SumOp in the lattice. You can however add arbitrary operators to an existing lattice from Python, save this lattice and use it in the calculation.

Local expectation values can likewise be evaluated:

$ syten-ipeps2-localexpectation.py -l sp -O sz -i rnd-projected.ips 2019.06.18T16:32:59.683 2.27784e-01 N Loading lattice from file 'sp'… 2019.06.18T16:32:59.685 2.29586e-01 N Loading state from file 'rnd-projected.ips'… 2019.06.18T16:32:59.685 2.30017e-01 N Building double layer… 2019.06.18T16:32:59.709 2.34243e-01 N Initialising CTM… 2019.06.18T16:32:59.729 2.54481e-01 N Working at χ = 1 2019.06.18T16:32:59.736 2.60886e-01 I Estimate 0: 2019.06.18T16:32:59.736 2.61052e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.736 2.61097e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.759 2.84730e-01 I Estimate 1 (change to previous iteration: 0.0000E+00): 2019.06.18T16:32:59.760 2.84868e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.760 2.84916e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.779 3.04158e-01 I Estimate 2 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.779 3.04293e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.779 3.04341e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.794 3.19647e-01 I Estimate 3 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.795 3.19755e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.795 3.19787e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.809 3.34255e-01 I Estimate 4 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.809 3.34420e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.809 3.34456e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.809 3.34523e-01 N Converged in CTM steps: Δ(-1) = 0.0000e+00 < 1.0000e-08 and Δ(-2) = 0.0000e+00 < 1.0000e-08 2019.06.18T16:32:59.816 3.40955e-01 N Working at χ = 2 2019.06.18T16:32:59.819 3.44384e-01 I Estimate 0: 2019.06.18T16:32:59.819 3.44510e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.819 3.44541e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.832 3.57470e-01 I Estimate 1 (change to previous iteration: 0.0000E+00): 2019.06.18T16:32:59.832 3.57570e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.832 3.57605e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.845 3.70141e-01 I Estimate 2 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.845 3.70256e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.845 3.70283e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.857 3.82610e-01 I Estimate 3 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.857 3.82708e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.857 3.82736e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.870 3.95077e-01 I Estimate 4 (change to previous iterations: 0.0000E+00, 0.0000E+00): 2019.06.18T16:32:59.870 3.95183e-01 I w = 1: (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.870 3.95207e-01 I w = 0: (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) 2019.06.18T16:32:59.870 3.95258e-01 N Converged in CTM steps: Δ(-1) = 0.0000e+00 < 1.0000e-08 and Δ(-2) = 0.0000e+00 < 1.0000e-08 2019.06.18T16:32:59.870 3.95309e-01 N Converged in χ: δ(-1) = 0.0000e+00 < 1.0000e-08 Best expectation value estimates at χ=2 (+5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) (-5.0000000000000000E-01,+0.0000000000000000E+00) (+5.0000000000000000E-01,+0.0000000000000000E+00) Associated error estimates: ±0.0000E+00 ±0.0000E+00 ±0.0000E+00 ±0.0000E+00

The tabular output is intended to resemble the unit cell layout, as we can see, we indeed have a Néel state.

The currently most stable tool in SyTen to find ground states of iPEPS is imaginary time evolution using the fast full update. Generic time evolution is implemented in bin/syten-ipeps2-update.py. This tool can proceed in multiple stages, at each stage defining a specific step size, iPEPS and environment bond dimension, minimal cost to incur during the truncation (to avoid small singular values) etc. It can also switch between different modes on-the-fly, so you can e.g. start with a simple update and then switch to a fast full update if you so wish.

A simple invocation is:

$ syten-ipeps2-update.py -l sp -H H -i rnd-projected.ips -o gs_search -s "(D 2 dt 0.1 l 5) (D 3) (D 4) (D 5)"

Here, `-l sp`

denotes the lattice file to use, `-H H`

the SumOp within that file, `-i rnd-projected.ips`

the initial state, `-o gs_search`

the output prefix and `-s (...)`

is the staging configuration. The output confirms that we have defined four stages:

2019.06.19T10:04:44.873 1.77824e-01 I Read staging description. 2019.06.19T10:04:44.873 1.77905e-01 I Stage 1: ( D 2 X 0 dt (0.1,0) tn 2 ft 1e-05 fx 40 fc 1e-14 fse 0 t 1e-06 l 5 x 1000000 c 1 ctc 0 s 1 Cs 10 Cv 1e-08 Cm WithQR Ct 1e-05 fg true fs true m Fast dS false ) 2019.06.19T10:04:44.873 1.77965e-01 I Stage 2: ( D 3 X 0 dt (0.1,0) tn 2 ft 1e-05 fx 40 fc 1e-14 fse 0 t 1e-06 l 5 x 1000000 c 1 ctc 0 s 1 Cs 10 Cv 1e-08 Cm WithQR Ct 1e-05 fg true fs true m Fast dS false ) 2019.06.19T10:04:44.873 1.78025e-01 I Stage 3: ( D 4 X 0 dt (0.1,0) tn 2 ft 1e-05 fx 40 fc 1e-14 fse 0 t 1e-06 l 5 x 1000000 c 1 ctc 0 s 1 Cs 10 Cv 1e-08 Cm WithQR Ct 1e-05 fg true fs true m Fast dS false ) 2019.06.19T10:04:44.873 1.78077e-01 I Stage 4: ( D 5 X 0 dt (0.1,0) tn 2 ft 1e-05 fx 40 fc 1e-14 fse 0 t 1e-06 l 5 x 1000000 c 1 ctc 0 s 1 Cs 10 Cv 1e-08 Cm WithQR Ct 1e-05 fg true fs true m Fast dS false ) 2019.06.19T10:04:44.873 1.78154e-01 I tau dtau E E/S ΔE/(N·δt) costfunc redD totD blkD redχ totχ blkχ #CT-U #CT-E file 2019.06.19T10:04:44.873 1.78199e-01 N Starting stage 1. 2019.06.19T10:04:44.874 1.78248e-01 N Configuration: ( D 2 X 0 dt (0.1,0) tn 2 ft 1e-05 fx 40 fc 1e-14 fse 0 t 1e-06 l 5 x 1000000 c 1 ctc 0 s 1 Cs 10 Cv 1e-08 Cm WithQR Ct 1e-05 fg true fs true m Fast dS false ) 2019.06.19T10:04:44.874 1.78282e-01 I tau dtau E E/S ΔE/(N·δt) costfunc redD totD blkD redχ totχ blkχ #CT-U #CT-E file 2019.06.19T10:04:44.949 2.44554e-01 I +0.000000e+00 +1.000000e-01 -2.000000000000e+00 -5.000000000000e-01 +nan +nan 1 1 1 1 1 1 -1 16 gs_search_T=+0.000000e+00+0.000000e+00j.ips 2019.06.19T10:04:45.373 6.68145e-01 I +1.000000e-01 +1.000000e-01 -2.297467717602e+00 -5.743669294004e-01 -2.974677e+00 +2.092276e-03 2 2 1 3 3 1 30 16 gs_search_T=+1.000000e-01+0.000000e+00j.ips 2019.06.19T10:04:45.853 1.14855e+00 I +2.000000e-01 +1.000000e-01 -2.452544327350e+00 -6.131360818374e-01 -7.753830e-01 +1.048937e-02 2 2 1 4 4 2 30 24 gs_search_T=+2.000000e-01+0.000000

The parameters of each stage are the iPEPS bond dimension `D`

, the environment bond dimension `X`

(defaults to D*D if 0), the step size `δ`

in `exp(-δH)`

, the full update convergence threshold `ft`

, the full update maximal iteration number (locally) `fx`

, the full update minimal cost `fc`

, whether the full update uses subspace expansion `fse`

, the convergence threshold in energy `t`

, the minimal number of steps `l`

, the maximal number of steps `x`

, the number of steps after which convergence is checked `c`

, whether convergence is checked for the gain in energy or the cost of the full update `ctc`

, the number of steps after which the state is changed `s`

, the number of CTM growth steps `Cs`

, the convergence threshold for the first singular value of the CTM `Cv`

, the CTM truncation method `Cm`

(either `WithQR`

, `SVDOnly`

or `SVDLiao`

), the CTM truncation threshold `Ct`

, wether the full update gauge fixes the norm tensor `fg`

, whether the gauge fixing uses the SVD `fs`

, the update mode itself `m`

(either `Fast`

, `Full`

or `Simple`

) and whether the reduced on-link tensors are created using a SVD (true) or QR (false) in `dS`

. As you can see, there are quite a few parameters. Most of these can be left at the defaults most of the time. Understanding them properly does require reading the source code and changing them may be necessary to obtain optimal results.

For now, we will simply look at the output provided:

2019.06.19T10:04:44.874 1.78282e-01 I tau dtau E E/S ΔE/(N·δt) costfunc redD totD blkD redχ totχ blkχ #CT-U #CT-E file ... 2019.06.19T10:04:59.419 1.47102e+01 I +3.100000e+00 +1.000000e-01 -2.638391745691e+00 -6.595979364229e-01 -1.204264e-06 +6.703025e-02 2 2 1 4 4 2 30 32 gs_search_T=+3.100000e+00+0.000000e+00j.ips 2019.06.19T10:04:59.872 1.51623e+01 I +3.200000e+00 +1.000000e-01 -2.638394593608e+00 -6.595986484020e-01 -8.899739e-07 +6.703623e-02 2 2 1 4 4 2 30 32 gs_search_T=+3.200000e+00+0.000000e+00j.ips 2019.06.19T10:04:59.872 1.51624e+01 N Energy decreased 8.8997e-07 < 1e-06 during last unit time, continuing to next stage… 2019.06.19T10:04:59.872 1.51630e+01 N Leaving stage, saving current state to 'gs_search_stg_1_final.ips' ... 2019.06.19T10:05:07.770 2.30583e+01 I +4.100000e+00 +1.000000e-01 -2.657822671360e+00 -6.644556678400e-01 -4.540450e-05 +4.250759e-02 3 3 2 9 9 5 30 48 gs_search_T=+4.100000e+00+0.000000e+00j.ips 2019.06.19T10:05:08.682 2.39697e+01 I +4.200000e+00 +1.000000e-01 -2.657800507935e+00 -6.644501269837e-01 +2.216343e-05 +4.324381e-02 3 3 2 9 9 5 30 48 gs_search_T=+4.200000e+00+0.000000e+00j.ips 2019.06.19T10:05:08.682 2.39698e+01 N Energy increased from -2.65782267e+00 to -2.65780051e+00, reverting steps and continuing to next stage… 2019.06.19T10:05:08.682 2.39705e+01 N Leaving stage, saving current state to 'gs_search_stg_2_final.ips' ... 2019.06.19T10:08:01.684 1.96945e+02 I +7.400000e+00 +1.000000e-01 -2.674950905381e+00 -6.687377263454e-01 -1.470723e-06 +1.202275e-02 4 4 2 16 16 6 30 72 gs_search_T=+7.400000e+00+0.000000e+00j.ips 2019.06.19T10:08:07.026 2.02286e+02 I +7.500000e+00 +1.000000e-01 -2.674954470020e+00 -6.687386175049e-01 -1.048423e-06 +1.203180e-02 4 4 2 16 16 6 30 72 gs_search_T=+7.500000e+00+0.000000e+00j.ips 2019.06.19T10:08:12.433 2.07693e+02 I +7.600000e+00 +1.000000e-01 -2.674956909664e+00 -6.687392274160e-01 -6.970412e-07 +1.203364e-02 4 4 2 16 16 6 30 72 gs_search_T=+7.600000e+00+0.000000e+00j.ips 2019.06.19T10:08:12.434 2.07693e+02 N Energy decreased 6.9704e-07 < 1e-06 during last unit time, continuing to next stage… 2019.06.19T10:08:12.436 2.07695e+02 N Leaving stage, saving current state to 'gs_search_stg_3_final.ips' .... 2019.06.19T10:16:42.222 7.17504e+02 I +1.000000e+01 +1.000000e-01 -2.676546531505e+00 -6.691366328763e-01 -1.127228e-06 +5.465222e-03 5 5 2 25 25 8 30 80 gs_search_T=+1.000000e+01+0.000000e+00j.ips 2019.06.19T10:17:01.712 7.36993e+02 I +1.010000e+01 +1.000000e-01 -2.676548798987e+00 -6.691371997468e-01 -9.069927e-07 +5.472547e-03 5 5 2 25 25 8 30 80 gs_search_T=+1.010000e+01+0.000000e+00j.ips 2019.06.19T10:17:01.712 7.36993e+02 N Energy decreased 9.0699e-07 < 1e-06 during last unit time, continuing to next stage… 2019.06.19T10:17:01.715 7.36996e+02 N Leaving stage, saving current state to 'gs_search_stg_4_final.ips'

The columns are mostly self-explanatory, `tau`

is the current value of the prefactor of the exponential, `dtau`

the current change, `E`

the expectation value of the operator, `E/S`

that value divided by the number of sites in the unit cell, `ΔE/(N·δt)`

the change in that value between the last and current check divided by the step size, `costfunc`

the cost function (either truncated singular values or difference in the full update), `redD`

, `totD`

, `blkD`

and likewise for `χ`

the reduced, total and maximal blocksize dimensions of the iPEPS state and its CTM, `CT-U`

and `-E`

the number of CTM steps done during the update and the expectation value evaluation respectively and finally the filename.

The true ground state energy of this system is -0.669437 according to Sandvik, Phys. Rev. B 56, 1997. Calculating the relative error as \( (E - E_{ex}) / E_{ex} \) gives for the above values 1.46e-2, 7.45e-3, 1.04e-3 and 4.47e-4, well in line with current best results from full update iPEPS calculations.

Expectation values can again be evaluated on the produced states, e.g.:

$ syten-ipeps2-localexpectation.py -l sp -O sz -i gs_search_stg_4_final.ips --threads-tensor 10

Note that now both more iterations at each CTM bond dimension and larger CTM bond dimensions are necessary. Also note that all of those tools understand `--threads-tensor 10`

just as well as bin/syten-dmrg.cpp does. `-x 64`

limits the maximal CTM bond dimension to 64. The above finally gives the output:

2019.06.19T10:33:42.280 5.93325e+02 I w = 0: (-3.4336800036475551E-01,-6.1590574170410547E-16) (+3.4318080287685498E-01,+6.2404373938749309E-16) 2019.06.19T10:33:42.280 5.93326e+02 N Unconverged in CTM steps: Δ(-1) = 1.5647e-07 > 1.0000e-08 or Δ(-2) = 9.7479e-07 > 1.0000e-08 2019.06.19T10:33:42.280 5.93326e+02 N Unconverged in χ: δ(-1) = 4.7985e-06 > 1.0000e-08 Best expectation value estimates at χ=64 (+3.4322932980660747E-01,+7.5972615808493723E-16) (-3.4304212935164830E-01,-7.5969176919106456E-16) (-3.4336800036475551E-01,-6.1590574170410547E-16) (+3.4318080287685498E-01,+6.2404373938749309E-16) Associated error estimates: ±2.1256E-04 ±2.0906E-04 ±2.1169E-04 ±2.0912E-04

While it complains that it is unconverged in both the CTM steps (with a change of 1e-7 in the last step) and the CTM bond dimension (with a change of 4e-6 going from 32 to 64), it does print the expectation values in tabular form with real and imaginary part and their associated estimated errors from the last three CTM bond dimension settings. As you can see, the magnetisation in this system is still much larger than expected when comparing to Monte Carlo results, which put this value at approximately 0.307.

The syten::STensor class underlying the iPEPSv2 code supports autodifferentiation. Autodifferentiation is a helpful tool to automatically obtain the gradient of the result of a specific function evaluation with respect to its input. Specifically, reverse-mode autodifferentiation is implemented for syten::STensor objects and most of the functions working on them. Only one bit is missing at the moment: the adjoint of the `U`

and `V`

tensors of the SVD for complex numbers. For real numbers, everything implemented works correctly.

In our application, we will consider the variational update of an iPEPS towards its ground state. We again use the Heisenberg square lattice model, but now with a rotated Hamiltonian such that the unit cell is one-by-one. As this is not implementable in iPEPSv2, we also add three dummy tensors which essentially just connect the tensors at position `(0,0)`

of the unit cell. This is implemented in lat/syten-ipeps2-1x1-spin.py.

For the autodifferentiation, we need a real-valued function `f: x → f(x)`

which we want to optimise with respect to `x`

. Both `x`

and `f(x)`

need to be syten::STensor objects, `f(x)`

needs to have rank 0. The pyten function `gradopt.gradient_descent_linesearch()`

implements a conjugate gradient optimisation with a line-search to find the optimal step size. An exemplary script implementing this function `f`

and the gradient descent may be:

#!/usr/bin/python3

from pyten import *

lat = ipeps2.Lattice("sp.ipl")

psi = ipeps2.State("rnd.ips")

ctm = ipeps2.CornerTransferMatrix(psi, ipeps2.CornerTransferMatrixGrowthMode.SVDLiao)

def expval(t):

trunc=1e-8

ctm.disable_autodiff()

psi[[0,0]] = t

psi.build_dl()

for i in range(0, 10):

ctm.grow(psi, Truncation(trunc, 50))

r = real(ipeps2.expectation_value(psi, ctm, lat.getSum("H")))

return r

init_e = expval(psi[[0,0]])

last_e = init_e

print("initial energy: ", init_e)

for iter in range(0, 100):

[opt_x, opt_v] = gradopt.gradient_descent_linesearch(psi[[0,0]], expval, 1, 50)

logGI("initial energy: ", init_e)

logGI("last energy: ", last_e)

logGI("After optimising: ", opt_v)

psi[[0,0]] = opt_x

psi[[0,0]].disable_autodiff()

last_e = opt_v

psi.save("trial-{}".format(iter+1))

First executing lat/syten-ipeps2-1x1-spin.py and then the Python script above will result in a series of optimisation steps being run. These are most likely not resulting in the optimal state at this bond dimension ( \( D = 3 \) by default) yet, though this is most likely a problem of the gradient descent instead of the autodifferentiation itself.

Note furthermore that the generated states sometimes behave less well than those generated by full update steps, meaning that the CTM may yield vastly wrong expectation values.