 |
SyTen
|
|
|
SyTen
|
|
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.latThis 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.mps
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.mps. 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!
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.mps
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.mps -l sp.lat:"Hz Hp +" -o H-rnd.mps
Given a state and a lattice, it becomes possible to calculate expectation values:
$ syten-expectation rnd.mps 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.mps 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.
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.mps 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.mps 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.mps -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.mps 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.mps -l sp.lat:'sp 10 @' -o gs-sp10.mps 2017.05.02T10:57:49.090 3.32850e-02 I Loaded MP state file 'gs.mps. 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.mps'… $ syten-localexpectation gs.mps sp.lat sm gs-sp10.mps -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.mps sp.lat:"sp 10 @ sm 11 @ *" -qr -0.349250736036831 $ syten-expectation gs.mps sp.lat:"sp 10 @ sm 12 @ *" -qr 0.120664807902911 $ syten-expectation gs.mps sp.lat:"sp 10 @ sm 13 @ *" -qr -0.125596562079963 $ syten-expectation gs.mps 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.mps
this should print 1. The operator application above will have changed the norm, however:
$ syten-norm H-rnd.mps Loaded MP state, calculating norm… 2.48746859276656
Given two states, we can calculate the overlap:
$ syten-overlap rnd.mps H-rnd.mps 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.mps -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
←) or right (→)m)DMRG will result in a series of *.mps files:
$ ls -l dmrg*state -rw-rw-rw- 1 C.Hubig ls-schollwoeck 57101 2017-04-25 16:49 dmrg_1_5.mps -rw-rw-rw- 1 C.Hubig ls-schollwoeck 247530 2017-04-25 16:49 dmrg_2_4.mps -rw-rw-rw- 1 C.Hubig ls-schollwoeck 744820 2017-04-25 16:49 dmrg_3_4.mps -rw-rw-rw- 1 C.Hubig ls-schollwoeck 745628 2017-04-25 16:49 dmrg_4_2.mps
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 *.mps); do printf "${i}:\n"; syten-expectation -qrv ${i} sp.lat:"Hz Hp +"; done
rnd.mps:
⟨O⟩: -1.3
⟨O²⟩: 6.1875
|⟨O²⟩-⟨O⟩²|: 4.4975
H-rnd.mps:
⟨O⟩: -3.3070707070707
⟨O²⟩: 17.0516414141414
|⟨O²⟩-⟨O⟩²|: 6.11492475257628
dmrg_1_5.mps:
⟨O⟩: -8.68247180616517
⟨O²⟩: 75.385320972876
|⟨O²⟩-⟨O⟩²|: 4.30802299433708e-06
dmrg_2_4.mps:
⟨O⟩: -8.68247333353095
⟨O²⟩: 75.3853431910421
|⟨O²⟩-⟨O⟩²|: 3.56602924968302e-09
dmrg_3_4.mps:
⟨O⟩: -8.68247333439879
⟨O²⟩: 75.3853432025469
|⟨O²⟩-⟨O⟩²|: 8.24277846980541e-13
dmrg_4_2.mps:
⟨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 rnd $ 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}.mps 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:
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.
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.
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.
Within the namespace, we then define the main() function with the usual signature.
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.
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.
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 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.
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).
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).
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.
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.
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.
This line saves the lattice object in the file output.
We return from the syten::main() function and close the syten namespace.
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
and/or copy it into the external binary directory
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:
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.
1.9.1